如何在 R studio 中使用线性混合效应模型确定固定效应和随机效应？答案

【问题标题】：How to determine the fixed and random effect using linear mixed effect model in R studio?如何在 R studio 中使用线性混合效应模型确定固定效应和随机效应？
【发布时间】：2022-11-27 23:07:11
【问题描述】：

我想使用线性混合模型研究物种和类型之间眼睛大小的变化。但是，我对判断是随机效应还是固定效应有点迷惑。我尝试使用这个模型：lmer(Eyesize ~ Species + (1|Type) in R studio

这是数据：data

这是我的代码：

library(lme4)
library(lmerTest)
library(multcomp)

example_data <-read.csv("EXAMPLE 1.csv", h= TRUE)
head(example_data)
result = lmer(Eyesize ~ Species + (1|Type), data=example_data)
summary(result)

我得到了这样的分析结果

我无法得到物种A的分析结果。我认为这是因为我在这种模式下使用了错误的随机或固定效应。我还想使用 Tukey 分析来指示字母的变化。谁能帮我解决这个问题？

太感谢了。

【问题讨论】：

标签： rstudio lme4 lmertest

【解决方案1】：

Species 的第 3 层没有显示，因为它已成为参考层。

听起来您正在进行事后比较，可以使用 emmeans::emmeans 在 lme4 对象上完成。

将来，请让数据更容易访问，就像这样：

library(lme4)
library(emmeans)
# df <- read.csv("stack_q_lme4.csv")
# df$Eyesize <- as.numeric(df$Eyesize)
# df$factor <- as.factor(df$Species)
# df$Type <- as.factor(df$Type)
#dput(df)
df <- structure(list(Eyesize = c(0.428278862, 0.400995265, 0.39874335, 
                           0.393709024, 0.425950648, 0.422276073, 0.423718381, 0.436375912, 
                           0.434729533, 0.44292607, 0.474962122, 0.480422148, 0.446239189, 
                           0.456237883, 0.479825541, 0.319118225, 0.34803004, 0.392926567, 
                           0.355179776, 0.327142563, 0.292814474, 0.295765514, 0.268828342, 
                           0.269960081, 0.2834768, 0.356604986, 0.348001735, 0.360895441, 
                           0.342982001, 0.366047801, 0.317237392, 0.293341926, 0.308761142, 
                           0.30039708, 0.308900879, 0.326311003, 0.33954796, 0.316756444, 
                           0.317607287, 0.33562927), Species = c("Species A", "Species A", 
                                                                 "Species A", "Species A", "Species A", "Species A", "Species A", 
                                                                 "Species A", "Species A", "Species A", "Species A", "Species A", 
                                                                 "Species A", "Species A", "Species A", "Species B", "Species B", 
                                                                 "Species B", "Species B", "Species B", "Species B", "Species B", 
                                                                 "Species B", "Species B", "Species B", "Species B", "Species B", 
                                                                 "Species B", "Species B", "Species B", "Species B", "Species B", 
                                                                 "Species B", "Species B", "Species B", "Species C", "Species C", 
                                                                 "Species C", "Species C", "Species C"), Type = structure(c(6L, 
                                                                                                                            6L, 6L, 6L, 6L, 8L, 8L, 8L, 8L, 8L, 3L, 3L, 3L, 3L, 3L, 5L, 5L, 
                                                                                                                            5L, 5L, 5L, 7L, 7L, 7L, 7L, 7L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 
                                                                                                                            1L, 1L, 2L, 2L, 2L, 2L, 2L), .Label = c("Type B", "Type K", "Type L", 
                                                                                                                                                                    "Type P", "Type Q", "type R", "Type T", "Type W"), class = "factor"), 
               factor = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                    1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
                                    3L), .Label = c("Species A", "Species B", "Species C"), class = "factor")), row.names = c(NA, 
                                                                                                                              -40L), class = "data.frame")

这是什么我认为你追求的是边际均值和物种的成对比较。


m1 <- lmer(Eyesize ~ Species + (1|Type), data=df)

emmeans(m1, pairwise ~ Species)
# $emmeans
# Species   emmean     SE df lower.CL upper.CL
# Species A  0.436 0.0189  5    0.388    0.485
# Species B  0.323 0.0163  5    0.281    0.365
# Species C  0.327 0.0327  5    0.243    0.411
# 
# Degrees-of-freedom method: kenward-roger 
# Confidence level used: 0.95 
# 
# $contrasts
# contrast              estimate     SE df t.ratio p.value
# Species A - Species B  0.11354 0.0250  5   4.547  0.0140
# Species A - Species C  0.10919 0.0377  5   2.893  0.0742
# Species B - Species C -0.00435 0.0365  5  -0.119  0.9922
# 
# Degrees-of-freedom method: kenward-roger 
# P value adjustment: tukey method for comparing a family of 3 estimates

【讨论】：