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

[julien.ruiz@airfrance.fr]

>>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"...