In PMassicotte/makroalgaeindicators: MacroalgaeIndicators

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)