Question

I have been using the code below to successfully modify the 'Zt', 'L', and 'A' slots of models fit using lme4 versions <1.0. I just updated to lme4 1.0-4 today and found that the model objects are different. Can anyone provide insight/guidance as to how to modify these slots in the new lmer model objects?

dat<-structure(list(pop1 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 10L, 
                    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 
                    4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 7L, 
                    7L, 8L, 8L, 9L), pop2 = c(2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 1L, 
                    3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
                    5L, 6L, 7L, 8L, 9L, 10L, 6L, 7L, 8L, 9L, 10L, 7L, 8L, 9L, 10L, 
                    8L, 9L, 10L, 9L, 10L, 10L), X = c(0.49136, 0.75587, 0.93952, 
                    0.61278, 0.79934, 1.07918, 1.13354, 1.15836, 1.2014, 0.43136, 
                    0.77815, 0.716, 0.93952, 1.13672, 1.16137, 1.18184, 1.21748, 
                    0.65321, 0.86332, 1.04922, 1.19866, 1.20412, 1.22272, 1.24797, 
                    0.89763, 1.08991, 1.19033, 1.15836, 1.17319, 1.18752, 0.64345, 
                    0.93952, 0.98227, 1.01703, 1.07188, 0.78533, 0.94939, 0.99564, 
                    1.06819, 0.64345, 0.716, 0.85126, -0.04576, 0.4624, 0.30103), 
           Y = c(0.491694, 0.394703, 0.113303, 0.156597, 0.450924, 0.487845, 
                 0.21821, 0.129027, -0.131522, 0.35156, -0.116826, 0.18941, 
                 0.306608, 0.258401, 0.008552, -0.024369, -0.305258, -0.013628, 
                 0.215715, 0.13783, 0.467272, 0.088882, 0.084295, -0.172337, 
                 -0.206725, -0.084339, -0.191651, -0.001586, -0.079501, -0.195094, 
                 0.232045, 0.17102, 0.003742, -0.023688, -0.26085, 0.205326, 
                 0.172809, 0.133219, -0.159054, 0.082231, 0.011025, -0.238611, 
                 0.732679, 0.478058, 0.325698)), .Names = c("pop1", "pop2", 
                  "X", "Y"), class = "data.frame", row.names = c(NA, -45L))

library(lme4) # lme4 versions >1.0 have different model output    

# Specify the model formula 
    lmer_mod <- as.formula("Y ~ X + (1|pop1)") 

# Create the Zl and ZZ matrices 
    Zl <- lapply(c("pop1","pop2"), function(nm) Matrix:::fac2sparse(dat[[nm]], "d", drop=FALSE)) 
    ZZ <- Reduce("+", Zl[-1], Zl[[1]]) 

# Fit lmer model to the data 
    mod <- lmer(lmer_mod, data = dat, REML = TRUE) 

# Replace the following slots in the fitted model
# These slots don't exist in this form in the new lmerMod objects
    mod@Zt <- ZZ 
    mod@A <- ZZ 
    mod@L <- Cholesky(tcrossprod(ZZ), LDL=FALSE, Imult=1) 

# Refit the model to the same response data 
    Final.mod <- refit(mod, dat[,Y])

Any help or insight as to how to modify these slots will be greatly appreciated. In the meantime, I guess I will stick to using an older version of lme4 for these models.

Was it helpful?

Solution

Does this do what you want? (This follows ?modular pretty closely ...)

Create the Zl and ZZ matrices:

Zl <- lapply(c("pop1","pop2"),
         function(nm) Matrix:::fac2sparse(dat[[nm]], "d", drop=FALSE)) 
ZZ <- Reduce("+", Zl[-1], Zl[[1]]) 

Construct the random-effects data structures:

lf <- lFormula(Y ~ X + (1|pop1), data=dat)

Modify them:

lf$reTrms$Zt <- ZZ

Proceed through the remaining model-construction and fitting steps:

dfun <- do.call(mkLmerDevfun,lf)   ## create dev fun from modified lf
opt <- optimizeLmer(dfun)          ## fit the model
## make the results into a 'merMod' object
fit <- mkMerMod(environment(dfun), opt, lf$reTrms,
         fr = lf$fr)
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top