2

我正在尝试在论文中绘制用于视觉表示的二元混合效果模型的结果。

我使用 lme 来拟合混合模型：

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

Pass = 二进制 1/0 zone.time, length & Fat = 连续

产生：

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

4

1 回答 1

2

有点令人困惑，但我收集到您想从混合效应模型中提取拟合结果，然后绘制它们。那是对的吗？

创建相似数据

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 轴上使用数据集中的变量，例如 Fat，在 y 轴上使用拟合值。

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

plot(fitted(M2))

拟合值与变量 Fat 的关系图

如您所见，对于这些虚构数据，关系是线性的。

于 2014-09-02T15:57:13.300 回答