Dear helpers,
I am using R version 2.5.1 to estimate a multinomial logit model using my
own maximum likelihood function (I work with share data and the default
function of R cannot deal with that).
However, the computer (I have an Athlon XP 3200+ with 512 GB ram) takes
quite a while to estimate the model.
With 3 categories, 5 explanatory variables and roughly 5000 observations it
takes 2-3 min. For 10 categories and 10 explanatory variables (still 5000
obs) more than 1 hour.
Is there any way I can speed up this process? (Modifying the code or
modifying some R options maybe?)
I would be really grateful if anybody could help me with this issue, I
attach my code below.
Many thanks,
Carlo
***************************************
Carlo Fezzi
Centre for Social and Economic Research
on the Global Environment (CSERGE),
School of Environmental Sciences,
University of East Anglia,
Norwich, NR4 7TJ
United Kingdom.
***************************************
# MULTILOGIT
# This function computes the estimates of a multinomial logit model
# inputs: a matrix vector of 1 and 0 (y) or of shares
# a matrix of regressors (x) - MUST HAVE COLUMN NAMES! -
# names of the variables, default = colnames(x)
# optimization methods, default = 'BFGS'
# base category, default = 1
# restrictions, default = NULL
# weights, default all equal to 1
# outputs: an object of class "multilogit.c"
# McFadden D. (1974) "Conditional logit analysis of qualitative choice
behavior", in Zarembka P. (ed.), Frontiers in Econometrics, Academic Press.
multilogit.c <- function(y, xi, xi.names = colnames(xi), c.base=1,
rest=NULL, w = rep(1,nrow(y)), method='BFGS')
{
n.obs <- sum(w)
xi<-cbind(1,xi)
colnames(xi)[1]<-"Intercept"
nx<-ncol(xi)
ny<-ncol(y)
beta<-numeric(nx*ny)
negll<- function(beta,y,xi)
{
beta[rest]<-0
beta[(((c.base-1)*nx)+1):(c.base*nx)]<-0
lli <- y * (xi%*%matrix(beta,nx,ny) - log ( apply(exp(
xi%*%matrix(beta,nx,ny)) ,1,sum ) ) )
lli<-lli*w
-sum(lli)
}
pi<- apply((y*w),2,mean)/mean(w)
ll0 <- (t(pi)%*%log(pi))*sum(w)
result<-c( optim(par = rep(0,nx*(ny)), fn = negll, y=y, xi=xi,
hessian=T, method=method),
list(varnames=xi.names, rest=rest, nx=nx, ny=ny,
npar=nx*(ny-1)-length(rest), ll0=ll0, pi=pi, xi=xi,
n.obs=n.obs,c.base=c.base,w=w))
result$par <- result$par[-(((c.base-1)*nx)+1):-(c.base*nx)]
result$hessian <-
result$hessian[-(((c.base-1)*nx)+1):-(c.base*nx),-(((c.base-1)*nx)+1):-(c.ba
se*nx)]
class(result)<-"multilogit.c"
return(result)
}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
