如何为视觉表示绘制二元混合效果模型

Question

我想绘制一个二进制混合效果模型的结果在视觉表示中的一篇论文。

我使用LME以适应混合模型：

M2 <- lme(Pass ~ zone.time + length + Fat, 
     random =~ 1 | Year) 

通行证=二进制1/0 zone.time，长度&脂肪=连续

得到：

Linear mixed-effects model fit by maximum likelihood 
Data: DF1 
    AIC  BIC logLik 
39.05604 47.25981 -13.52802 

Random effects: 
Formula: ~1 | Year 
     (Intercept) Residual 
StdDev: 5.03879e-06 0.3857927 

Fixed effects: Pass ~ zone.time + length + Fat 
       Value Std.Error DF t-value p-value 
(Intercept) 4.549716 1.2384118 24 3.673832 0.0012 
zone.time 0.299438 0.1239111 24 2.416559 0.0236 
length  -0.006718 0.0019492 24 -3.446603 0.0021 
Fat   -0.051460 0.0213211 24 -2.413563 0.0238 
Correlation: 
      (Intr) zon.tm length 
zone.time 0.045    
length -0.979 -0.168  
Fat  -0.447 -0.191 0.330 

Standardized Within-Group Residuals: 
     Min   Q1  Med   Q3  Max 
-1.9097237 -0.7802111 -0.1410353 0.5683329 2.0908188 

Number of Observations: 29 
Number of Groups: 2 

我然后去计算预测值和标准误差：

MyData <- expand.grid(zone.time = seq(1,3.6, length = 10), 
        length = seq(525, 740, length = 10), 
        Fat = seq(3.7, 17, length = 10)) 
X <- model.matrix(~zone.time + length + Fat, data = MyData) 

提取参数和参数协方差矩阵

betas <- fixef(M2) 

用于样本数据使用

betas<- structure(c(4.54971638246632, 0.299438350935228, -0.00671801197327911,-0.0514597408192487), .Names = c("(Intercept)", "zone.time", "length","Fat")) 

。

Covbetas <- vcov(M2) 

样本数据使用：

Covbetas <- structure(c(1.32212400759181, 0.0059001955657893, -0.00203725210229123, 
-0.0101822039057957, 0.0059001955657893, 0.0132361635192455, 
-3.50672281561515e-05, -0.000434188193496185, -0.00203725210229123, 
-3.50672281561515e-05, 3.27522409259271e-06, 1.18250356798504e-05, 
-0.0101822039057957, -0.000434188193496185, 1.18250356798504e-05, 
0.000391886154502855), .Dim = c(4L, 4L), .Dimnames = list(c("(Intercept)", 
"zone.time", "length", "Fat"), c("(Intercept)", "zone.time", 
"length", "Fat"))) 

计算在预测器规模

MyData$eta <- X %*% betas 
MyData$Pi <- exp(MyData$eta)/(1 + exp(MyData$eta)) 

计算上的预测器功能的规模的SE拟合值

MyData$se <- sqrt(diag(X %*% Covbetas %*% t(X))) 
MyData$SeUp <- exp(MyData$eta + 1.96 *MyData$se)/(1 + exp(MyData$eta + 1.96 *MyData$se)) 
MyData$SeLo <- exp(MyData$eta - 1.96 *MyData$se)/(1 + exp(MyData$eta - 1.96 *MyData$se)) 

head(MyData) 

这是计算预测值的正确方法？

我该如何绘制这个视觉呈现？

我应该使用像

library(effects) 
plot(allEffects(M2, default.levels=50)) 

或GGPLOT2

Answer 1

类的混乱，但我现在找你要提取您的混合效应模型拟合结果，然后绘制出来。那是对的吗？

创建类似的数据

set.seed(64) 

fooDF <- data.frame(Pass = rbinom(n = 100, size = 1, prob = 0.5), zone.time = rnorm(n = 100), length = rnorm(n = 100), 
        Fat = rnorm(n = 100), Year = seq(1913, 2012)) 

M2 <- lme(Pass ~ zone.time + length + Fat, 
       random =~ 1 | Year, data = fooDF) 

您可以通过

head(fitted(M2, level = 0)) 

    1913  1914  1915  1916  1917  1918 
    0.4948605 0.7506069 0.5317316 0.5429997 0.6584630 0.7555496 

得到整体预测的结果，您可以简单地绘制适合像这样

plot(fitted(M2, level = 0)) 

您还可以在使用变量你的数据集，比如脂肪，在x轴上的拟合值在y上。

plotDF <- data.frame(fat = fooDF$Fat, fitted = fitted(M2, level = 0)) 
plot(plotDF) 

plot(fitted(M2)) 

正如你可以看到，这个虚拟数据的关系是线性的。