CODEWITHSUNDEEP
rloopsmaxlog-likelihood
Bruce
Bruce

Reputation: 81

How to save global of maximum log-likelihood data from loop?

I am doing this:

## Example data
library(poLCA)
data(gss82)
f <- cbind(PURPOSE,ACCURACY,UNDERSTA,COOPERAT)~1

mlmat <- matrix(NA,nrow=500,ncol=4)

## Do repreat 500 times LCA for best llik
for (i in 1:500) {
  gss.lc <- poLCA(f,gss82,nclass=3,maxiter=3000,tol=1e-7)
  mlmat[i,1] <- gss.lc$llik
  o <- order(gss.lc$probs$UNDERSTA[,1],decreasing=T)
  mlmat[i,-1] <- gss.lc$P[o]
 }

It's come from poLCA paper (http://www.sscnet.ucla.edu/polisci/faculty/lewis/pdf/poLCA-JSS-final.pdf p14. table1)

This is result: enter image description here

I use loop to do 500 times just show one objective on data panel of R.

I want know how to just save the maximumlog-likelihood is -2754.545 output, because it's global of maximum log-likelihood.

Upvotes: 0

Views: 41

Answers (2)

dajah&#246;
dajah&#246;

Reputation: 11

See Linzer and Lewis (2011, p. 16):

To avoid these local maxima, a user should always either 1) call poLCA at least a couple of times; or 2) utilize the nrep argument to attempt to locate the parameter values that globally maximize the log-likelihood function.

Both your answer and Edo's answer use the first approach, while the second approach might be more swift, because (Linzer and Lewis 2011, p. 12):

Setting nrep > 1 automates the search for the global – rather than just a local – maximum of the log-likelihood function. poLCA returns only the parameter estimates corresonding to the model producing the greatest likelihood.

Retrieve maximum log-likelihood with

gss.lc$llik

Upvotes: 1

Edo
Edo

Reputation: 7858

Based on the comment you left, you don't need a for loop and you don't need mlmat.

You can run 500 times the same code and save the results into a list. Then you select the item of the list with the highest llik. At that point you can just select the item of the list that you are looking for.

library(purrr)

# run 500 times
gss.lcs <- map(1:500, \(i) poLCA(f,gss82,nclass=3,maxiter=3000,tol=1e-7))

# get the position of the highest one
pos <- gss.lcs |> map_dbl("llik") |> which.max()
pos

# get the element with the highest llik
gss.lc <- gss.lcs[[pos]]

gss.lc

Upvotes: 1

Related Questions