Hi Ben.
Following your answer on Stackoverflow, I tried to fix parameter as you described.
library(lme4) rm(list = ls()) set.seed(101) ss <- sleepstudy[sample(nrow(sleepstudy),size=round(0.9*nrow(sleepstudy))),] m1 <- lmer(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), ss) summary(m1)
This model gives me the following estimates on fixed effects:
coef(summary(m1))
To test your answer, I tried to fix the estimated variances of m1
. Doing this, I was expecting to find exactly the same estimates on the fixed effects.
dd <- as.function(m1) # These values are those originally estimated ff <- dd(c(sqrt(703.38),sqrt(35.34)))
Now, let's compare original estimates with the new ones.
coef(summary(m1)) environment(dd)$pp$beta(1) ## new parameters
As we can see, these are slightly different. What I am missing?
Thank you for your help, Philippe
opt <- list(par = c(0, 0), fval = ff, conv = 0) lmod <- lFormula(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), ss) m1X <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr, mc = quote(hacked_lmer())) buildMM <- function(theta) { dd <- as.function(m1) ff <- dd(theta) opt <- list(par = c(0, 0), fval = ff, conv = 0) mm <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr, mc = quote(hacked_lmer())) return(mm) } objfun <- function(x,target=c(700,30)) { mm <- buildMM(sqrt(x)) return(sum((unlist(VarCorr(mm))-target)^2)) } s0 <- c(700,30)/sigma(m1)^2 opt <- optim(fn=objfun,par=s0) mm_final <- buildMM(sqrt(opt$par)) summary(mm_final)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.