Seeking help for using optim for MLE in R

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Seeking help for using optim for MLE in R

Naznin Sultana
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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Seeking help for using optim for MLE in R

Bert Gunter-2
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 <[hidden email]> 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]]
>
> ______________________________________________
> [hidden email] 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]]

______________________________________________
[hidden email] 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.