[R] Seeking help for using optim for MLE in R

Bert Gunter bgunter@4567 @ending from gm@il@com
Thu Jan 10 16:44:42 CET 2019


Probably: don't do this.

Use the nnet package (and there may well be others) to fit multinomial
regression. See here for a tutorial:

https://rpubs.com/rslbliss/r_logistic_ws

Cheers,
Bert


On Thu, Jan 10, 2019 at 6:18 AM Naznin Sultana <zeny.stat using gmail.com> wrote:

> Hi, I am writing a program for MLE of parameters of multinomial
> distribution
> using optim.But when I run the program outside of a function it gives me
> the
> likelihood value, but when using it for optim function it gives the error
> message "Error in X %*% beta : non-conformable arguments".
> If X, and beta are non-conformable how it gives values.
> My data has first three columns of three dependent variables and rest of
> the
> colums indicating X (indep vars).
> Please help me out. Here goes my program for k1 categories of multinomial
> distribution:
>
> #data is the data which consists of three dependent varaible in first three
> columns and rest of the columns represent covariates.
>
>
> k1<- length(unique(data[,1]))
> p<- ncol(data)-3
> beta0 <-matrix(-.00001,nrow=k1-1,ncol=(p+1)) # starting value
> beta <-as.matrix(beta0)
> beta <-as.matrix(t(beta))
>
>
>
>
> ## likelihood for y1
>
> multin.lik<- function(beta,data) ##beta is a matrix of beta's of order
> ((p+1)*(k-1))
>                 {
>                         nr<- nrow(data)
>                         nc<- ncol(data)
>
>                         y1<- data[,1]
>                         y1<- as.matrix(y1,ncol=1)
>
>                         X<-as.matrix(cbind(1,data[,4:nc])) #matrix of order
> ((n*(p+1)))
> covariates; 1 is added for intercept
>
>                         LL<- exp(X%*%beta) #LL is of order (n*(k-1))
>                         L<- as.matrix(cbind(1,LL))  #L is of order (n*k); 1
> is added for ref
> category, L0, L1, L2
>                         pi<- t(apply(L,1, function(i) i/sum(i)))
>
>
>                         lgl<- 0
>                         for (i in 1:nr)
>                                 {
>                                         if (y1[i]==0) {lgl[i]<-
> log(pi[i,1])}
>                                         else if (y1[i]==1) {lgl[i]<-
> log(pi[i,2])}
>                                         else lgl[i]<- log(pi[i,3])
>                                 lgl
>                                 }
>                         lgL<- sum(lgl)
>                         return(-lgL)
>                 }
>
>
> ## parameter estimates
> abc <-optim(beta, multin.lik,data=data,method="SANN",hessian=T)
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

	[[alternative HTML version deleted]]



More information about the R-help mailing list