2013-09-30 63 views
1

我一直在使用下面的代碼來成功修改適合使用lme4版本< 1.0的模型的'Zt','L'和'A'插槽。我剛剛更新到lme4 1.0-4今天,發現模型對象是不同的。任何人都可以提供有關如何在新lmer模型對象中修改這些插槽的見解/指導?如何修改插槽lme4> 1.0

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 
    [email protected] <- ZZ 
    [email protected] <- ZZ 
    [email protected] <- Cholesky(tcrossprod(ZZ), LDL=FALSE, Imult=1) 

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

任何幫助或瞭解如何修改這些插槽將不勝感激。同時,我想我會堅持使用舊版本的lme4來處理這些模型。

+0

讀取新的'lme4'包'modular'頁面一開始... –

回答

1

這是做你想做的嗎? (在此之前,?modular非常密切......)

創建Z1和ZZ矩陣:

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

構建隨機效應的數據結構:

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

對其進行修改:

lf$reTrms$Zt <- ZZ 

繼續完成剩餘的模型構建和安裝步驟:

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) 
+0

這看起來像它有潛力,但偏差函數產生這個錯誤:在Lambdat%*%UT錯誤: Cholmod錯誤'A和B內部尺寸必須匹配'文件../MatrixOps/cholmod_ssmult.c,第82行 –

+0

這很有趣,它適用於我。 (例如'VarCorr(fit)= pop1(截距)= 0.10208,殘差= 0.14750')。你可以給'sessionInfo()'嗎?這與lme4 1.1-0,Matrix 1.0-14 ...你可能需要確保你最近重新安裝了矩陣和RcppEigen - 在編譯的基礎庫必須匹配... –

+0

我的歉意,本。您的解決方案奏效我沒有在我的原始代碼中正確指定我的Zl矩陣(應該是'dat',而不是我的代碼中的'data')。使用正確的Zl和ZZ矩陣,您的代碼可以正常工作,並且結果與預期輸出相匹配。謝謝! (我編輯了我的原始代碼,這樣可以複製) –