# Réf. : Re: [R] Choice modelling (was:(no subject))

julien.ruiz@airfrance.fr julien.ruiz at airfrance.fr
Fri Dec 3 14:33:10 CET 2004

```

>>I could incorporate indicators of choice availability as explanotary
>>variables, but it does not seem a very good way to do it.
>>Instead, for a logit model, I have coded a likelihood computation of the
>>underlying model with varying choice set and I use optim function to get
>>the "maximum".
>>
>>
>Could you post the code?

I can but will be ashamed of my coding...

Something like the following :

calcLikelihood<-function(beta,x,y, S) {
#beta  : parameter to be estimated
#S: 0/1 matrix of choices available for each individual
#x matrix of choice caracteristic (only compute intercept)
#beta[k] : intercept of choice k

#y : choice made for each indiv
logLikObs <- 0
for (t in 1:length(y)) {
#utility of choice at t
u_ch_t= beta %*%  x[y[t],]
# overall utility of possible choices at t
u_tot_t = sum( exp( beta %*% x )) + 1
logLikObs=logLikObs + (u_ch_t - log (u_tot_t) )
}

return(logLikObs)
}

#Dummy exemple :
calcProbfromUtil <- function (utilVector) {
v<-exp(utilVector)
return(v/sum(v))
}

drawFromProb <- function (probVector) {
x<-(1:length(probVector))
u<-runif(1)
if (u > max(probVector)) {
return(min(x[probVector==max(probVector)]))
} else {
return(min(x[probVector==min(probVector[probVector>u])]))
}
}

nb_class<-4 #nb of choices

maxBeta <- 5
minBeta <-  -maxBeta
betaFree <-round(runif(nb_class-1, min=minBeta, max=maxBeta),2)
beta <- c(0,betaFree)  #parameter to be estimated

x <- diag(nb_class)

T <- (nb_class*5000)  #nb of individuals

S <- array(ifelse(runif(nb_class*T)>0.5,1,0),dim=c(T,nb_class)) #0/1 matrix of choices available for each individual

util <- t(apply(S,1,function(x) {x*beta})) #utility for each choices

prob <- t(apply(util,1,calcProbfromUtil)) #corresponding probability

y <- as.vector(apply(prob,1,drawFromProb))

#betaFree : paramÃ¨tre libre de beta

calcLikelihoodFromBetaFree <- function(betaFree) {
return(calcLikelihood(c(0,betaFree),x,0,y, S))
}

#AND FINALLY
resOptim <- optim(par=rep(0,nb_class-1),
fn=calcLikelihoodFromBetaFree,
method="L-BFGS-B",
lower=rep(minBeta*1.5,nb_class-1),
upper=rep(maxBeta*1.5,nb_class-1),
control=list(trace=0,fnscale=-1))

beta_est<-c(0,resOptim\$par)

#of course , in order to do it right beta_est should be "normalized"...

```