[R] logLIk(lme(...))?

[Spencer Graves]

I found my problem:


The following function gave llGp1, llGp2 and ll22:


logLik_lm <- function(object){
   res <- resid(object)
   n <- length(res)
   s2MLE <- sum(res^2)/n
   lglk <- (-n/2)*(log(2*pi*s2MLE)+1)
   lglk
}
logLik(fitGp1)
logLik(fitGp1)-logLik_lm(fitGp1)
llGp1 - logLik_lm(fitGp1)


logLik(fitGp2)
logLik(fitGp2)-logLik_lm(fitGp2)
llGp2 - logLik_lm(fitGp2)

logLik(fit22)
logLik(fit22)-logLik_lm(fit22)
ll22 -logLik_lm(fit22)


	  These differences were all 0 to within roundoff error.  That 
confirmed for me that I could safely compare logLik.lm and logLik.lme.


	  What I thought should have been a linear operation wasn't.  Please 
excuse the waste of your time.


	  Thanks,
	  Spencer Graves


On 8/29/23 11:15 AM, Spencer Graves wrote:
> Hello, all:
> 
> 
>        I have a dataset with 2 groups.  I want to estimate 2 means and 2 
> standard deviations.  I naively think I should be able to use lme to do 
> that, e.g., lme(y~gp, random=y~1|gp, method='ML').  I think I should get 
> the same answer as from lm(y~1, ...) within each level of group.  I can 
> get the same means, but I don't know how to extract the within-gp 
> standard deviations, and the sum of logLik for the latter two does not 
> equal the former.
> 
> 
> TOY EXAMPLE:
> 
> 
> library(nlme)
> 
> 
> set.seed(1)
> 
> 
> lmePblm <- data.frame(y=c(rnorm(5, 1, 2), rnorm(5,3,5)),
>                        gp=factor(rep(1:2, each=5)))
> 
> 
> fit22 <- lme(y~gp, lmePblm, random=~1|gp, method='ML')
> 
> 
> fitGp1 <- lm(y~1, lmePblm[lmePblm$gp==1, ])
> 
> 
> fitGp2 <- lm(y~1, lmePblm[lmePblm$gp==2, ])
> 
> 
> (ll22 <- logLik(fit22))
> 
> 
> (llGp1 <- logLik(fitGp1))
> 
> 
> (llGp2 <- logLik(fitGp2))
> 
> 
> # Why isn't (ll22 = llGp1+llGp2)?
> 
> 
> (ll22 - llGp1-llGp2)
> 
> 
> # And secondarily, how can I get the residual standard deviations
> # within each gp from fit22?
> 
> 
>        Thanks,
>        Spencer Graves
> 
> ______________________________________________
> 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.