R：在混合模型效应中绘制 RANEF

Question

所以我对 R 很陌生。我正在做一个项目，其数据集是鸟类时代的外观。我们有来自 95 个人的 >400,000 次观察。我的任务是这样的：

“如果这个图在线条旁边有一些数据点，会更有说服力。你可以在 lme4 的 getcall 中从 ranef 获取数据点。因此，运行没有年龄和年龄 2 项的模型。然后，拉出随机效应项（每个人的值），然后绘制它们（y = 随机效应，x = 年龄）。"

这是有问题的图表（它是使用我们的模型在 Excel 中制作的）：Graph

所以我运行了这个

> lmbdna<-lmer(Gs ~ (1 | Individual) + Bin + year + Mass)
> ranef(lmbdna)
>$Individual

> plot(Age, ranef(lmbdna))

我得到了这个

xy.coords(x, y, xlabel, ylabel, log) 中的错误： “x”和“y”长度不同

我真的迷路了，我不知道如何将 ranefs 绘制到年龄。有没有办法将年龄与个人联系起来以消除该错误？

这是我的一些数据：

Indiv   Age Mass   MaxDepth Depth       Gs       PDBA  Bin year
69903   12  1015    3.806   3.025   0.1854302   92.7151 A   N
52712   20  957.5   3.806   3.025   0.204678    102.339 A   T
55969   19  1002.5  3.806   3.025   0.222338    111.169 A   T
64442   15  1040    3.806   3.025   0.1872954   93.6477 A   T
76252   11  940     3.806   3.025   0.223136    111.568 A   T
53391   21  1022.5  3.806   3.025   0.234452    117.226 A   E
53391   21  1022.5  3.806   3.025   0.299438    149.719 A   E
60117   18  937.5   3.806   3.025   0.1469442   73.4721 A   E
60151   18  970     3.806   3.025   0.1941052   97.0526 A   E
52712   20  957.5   3.855   3.025   0.1812926   90.6463 A   T
52712   20  957.5   3.855   3.025   0.25101     125.505 A   T
64442   15  1040    3.855   3.025   0.1850976   92.5488 A   T
64442   15  1040    3.855   3.025   0.1026478   51.3239 A   T
76252   11  940     3.855   3.025   0.235822    117.911 A   T
78712   10  880     3.855   3.025   0.1638106   81.9053 A   T
87819   7   1000    3.855   3.025   0.166391    83.1955 A   T
90281   6   957.5   3.855   3.025   0.1493948   74.6974 A   T
60151   18  970     3.855   3.025   0.1904232   95.2116 A   E
69944   12  915     3.904   3.025   0.256504    128.252 A   N
3260    24  960     3.904   3.025   0.168019    84.0095 A   T
52712   20  957.5   3.904   3.025   0.270704    135.352 A   T
64442   15  1040    3.904   3.025   0.1507102   75.3551 A   T
71432   12  970     3.904   3.025   0.1238154   61.9077 A   T
80538   15  917.5   3.904   3.025   0.236976    118.488 A   E
76583   14  870     3.952   3.025   0.295982    147.991 A   N
84420   7   1005    3.952   3.025   0.1861876   93.0938 A   N
87819   7   1000    3.952   3.025   0.178179    89.0895 A   T
53391   21  1022.5  3.952   3.025   0.1917954   95.8977 A   E
53391   21  1022.5  3.952   3.025   0.1482036   74.1018 A   E
53391   21  1022.5  3.952   3.025   0.1999868   99.9934 A   E
53391   21  1022.5  3.952   3.025   0.276334    138.167 A   E
60151   18  970     3.952   3.025   0.1776108   88.8054 A   E
80538   15  917.5   3.952   3.025   0.188733    94.3665 A   E
69944   12  915     4.001   3.025   0.2596      129.8   A   N
3260    24  960     4.001   3.025   0.1824546   91.2273 A   T

感谢任何帮助。谢谢

Answer 1

主要问题是您提取的随机效应是针对每个人的，而您尝试绘制的年龄数据是针对每个观察的。您需要将其汇总到个人级别（例如，在每个人的所有观察中取最大值以获得您正在寻找的结果。下面的示例使用您在上面提供的文本。

不确定这是否是您所追求的，但主要问题是您的向量长度不同。希望下面的代码能给你一些想法。获得更大的数据样本会很棒。

require(lme4)

# the data snapshot supplied above only contains one level of the factor bin, let's add another
dta[dta$Indiv < 7000, "Bin"] <- "B"

# estimate model
lmbdna <- lmer(Gs ~ (1 | Indiv) + Bin + year + Mass, dta)

# random effects are for indivuals, so need to aggregate individuals' ages
plt_dta <- aggregate(dta$Age, by = list(dta$Indiv), max)

# create one dataset by merging with random effects
plt_dta <- merge(plt_dta, ranef(lmbdna)$Indiv, by.x = 1, by.y = 0)

# plot data
plot(plt_dta[,2], plt_dta[,3], xlab = "age", ylab = "ranef")

Answer 2

我认为我们无法根据您共享的数据获得有意义的模型，所以让我们使用类似于 ?ranef 中的示例：

library("lme4")
fm1 = lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)

my_ranef = ranef(fm1)
str(my_ranef)
# List of 1
#  $ Subject:'data.frame':  18 obs. of  1 variable:
#   ..$ (Intercept): num [1:18] 40.78 -77.85 -63.11 4.41 10.22 ...
#  - attr(*, "class")= chr "ranef.mer"

str() 可用于查看对象的结构。在这里我们可以看到my_ranef是一个包含1个项目的列表，而那个项目是一个数据框（名为Subject），有18行1列（名为(Intercept)）。让我们仔细看看数据框：

head(my_ranef$Subject)
#     (Intercept)
# 308   40.783710
# 309  -77.849554
# 310  -63.108567
# 330    4.406442
# 331   10.216189
# 332    8.221238

所以数据框有对应每个个体主题的行名，然后随机效应截距在（截距）列。所以我们可以像这样绘制随机效应：

plot(row.names(my_ranef$Subject), my_ranef$Subject[, "(Intercept)"])

您遇到的问题是您将整个列表或数据框作为y 参数提供给您的绘图，而您只需要提取向量。

我们也可以提取随机效应

est_ranef = my_ranef$Subject[, "(Intercept)"]

（请注意，括号使 "(Intercept)" 成为非标准列名，这增加了提取它的难度，简单的 $ 不起作用a la my_ranef$Subject$(Intercept)，但是我们可以使用[ 和上面引用的列名，或者反引号可以与$ 一起使用，如下所示：

my_ranef$Subject$`(Intercept)` 

如果您查看ranef() 以获得更复杂的模型，您会看到为什么使用数据帧结构列表。