R：在混合模型效应中绘制 RANEF
R: Plotting RANEFs in Mixed Model Effect
所以我对 R 还很陌生。我正在从事一个项目，该项目的数据集是关于鸟类年龄的。我们有来自 95 个人的超过 400,000 个观察结果。我的任务是这样的：
"This graph would be more convincing if it had some data points on it alongside the lines. You can get the datapoints from ranef within getcall in lme4. So, run the model without an age and age2 term. Then, pull out the random effect terms (the values for each individual), then plot them (y = random effect, x = age)."
这是有问题的图表（它是在 Excel 中使用我们的模型制作的）：Graph
所以我运行这个
> lmbdna<-lmer(Gs ~ (1 | Individual) + Bin + year + Mass)
> ranef(lmbdna)
>$Individual
> plot(Age, ranef(lmbdna))
我明白了
Error in xy.coords(x, y, xlabel, ylabel, log) :
'x' and 'y' lengths differ
我真的迷路了，我不知道如何绘制 运行efs 到年龄的图表。有没有办法将年龄与个人相关联以消除该错误？
这是我的一些数据：
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
感谢任何帮助。谢谢
主要问题是您提取的随机效应是针对每个人的，而您尝试绘制的年龄数据是针对每个观察结果的。您需要将其汇总到个人级别（例如，对每个人的所有观察结果取最大值以获得您正在寻找的结果。下面的示例使用您在上面提供的文本。
不确定这是否是您想要的，但主要问题是您的向量长度不同。希望下面的代码能给你一些想法。获得更大的数据样本会很棒。
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")
我不认为我们可以根据您共享的数据得到一个有意义的模型，所以让我们使用类似于 ?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() 以获得更复杂的模型，您将了解 为什么 使用数据帧结构列表。
所以我对 R 还很陌生。我正在从事一个项目，该项目的数据集是关于鸟类年龄的。我们有来自 95 个人的超过 400,000 个观察结果。我的任务是这样的：
"This graph would be more convincing if it had some data points on it alongside the lines. You can get the datapoints from ranef within getcall in lme4. So, run the model without an age and age2 term. Then, pull out the random effect terms (the values for each individual), then plot them (y = random effect, x = age)."
这是有问题的图表（它是在 Excel 中使用我们的模型制作的）：Graph
所以我运行这个
> lmbdna<-lmer(Gs ~ (1 | Individual) + Bin + year + Mass)
> ranef(lmbdna)
>$Individual
> plot(Age, ranef(lmbdna))
我明白了
Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' lengths differ
我真的迷路了，我不知道如何绘制 运行efs 到年龄的图表。有没有办法将年龄与个人相关联以消除该错误？
这是我的一些数据：
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
感谢任何帮助。谢谢
主要问题是您提取的随机效应是针对每个人的，而您尝试绘制的年龄数据是针对每个观察结果的。您需要将其汇总到个人级别（例如，对每个人的所有观察结果取最大值以获得您正在寻找的结果。下面的示例使用您在上面提供的文本。
不确定这是否是您想要的，但主要问题是您的向量长度不同。希望下面的代码能给你一些想法。获得更大的数据样本会很棒。
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")
我不认为我们可以根据您共享的数据得到一个有意义的模型，所以让我们使用类似于 ?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() 以获得更复杂的模型，您将了解 为什么 使用数据帧结构列表。