r - 我在 lme4::lmer() 中的语法有什么问题，用于具有不平衡重复测量的拆分图设计？

Question

我正在尝试使用 R 中的 lme4 包和函数 lmer() 来为我的 split - split 图设计拟合模型。如果我没有丢失少量观察值，我会使用重复测量方差分析，但丢失的数据对于线性混合效应模型应该没有问题。

我的数据框 (data) 有一个简单的结构，包含四个因素和一个名为 all_vai 的数字结果变量。请注意，在此示例数据框中，并非所有因素的所有级别都被交叉，即使它们在我的真实数据中(缺失的观察除外)。这对我的问题来说应该无关紧要，这是为了修复有问题的语法。

collected_vai <- rnorm(125, mean = 6, sd = 1)
missing <- rep(NA, times = 3)
all_vai <- c(collected_vai, missing)
year1 <- rep(2018, times = 32)
year2 <- rep(2019, times = 32)
year3 <- rep(2020, times = 32)
year4 <- rep(2021, times = 32)
year <- c(year1, year2, year3, year4)
disturbance_severity <- rep(c(0,45,65,85), each = 32)
treatment <- rep(c("B" , "T"), each = 64)
replicate <- rep(c("A", "B", "C", "D"), each = 32)
data = data.frame(all_vai, year, disturbance_severity, treatment, replicate)
data$year <- as.factor(data$year)
data$disturbance_severity <- as.factor(data$disturbance_severity)
data$treatment <- as.factor(data$treatment)
data$replicate <- as.factor(data$replicate)

这是我为具有不同(正态分布)数值结果且没有缺失观察值的相同数据集运行的模型——也就是说，如果我现在没有不平衡的重复测量，这就是我将要运行的模型丢失数据:

VAImodel1 <- aov(all_vai ~ disturbance_severity*treatment*year + Error(replicate/disturbance_severity/treatment/year), data = data)
summary(VAImodel1)

当我运行它时，我收到错误消息:“警告消息:在 aov(mean_vai ~ disturbance_severity * treatment * Year + Error(Replicate/disturbance_severity/treatment/Year) 中:错误()模型是单一的”

我有来自不同年份的观察嵌套在不同的处理中，这些观察嵌套在不同的干扰严重性中，并且所有这些都嵌套在重复(它们是实验 block )中。所以我尝试在 lme4 中使用这个结构:

library(lme4)
library(lmerTest)

VAImodel2 <- lmer(all_vai ~ (year|replicate:disturbance_severity:treatment) + disturbance_severity*treatment*year, data = data)
summary(VAImodel2)

这是我收到的错误消息:“错误:观察次数 (=125) <= 随机效应数 (=128) for term (Year | Replicate:disturbance_severity:treatment)；随机效应参数和残差(或尺度参数)可能无法识别”

接下来，我尝试通过删除处理变量和交互项来简化我的模型，这样我就不会用尽自由度，如下所示:

VAImodel3 <- lmer(all_vai ~ (year|replicate:disturbance_severity) + disturbance_severity*year, data = data)
summary(VAImodel3)

这次我得到了一个不同的错误:“边界(奇异)拟合:见？isSingular警告信息:模型未能收敛于 1 个负特征值:-1.2e-01 "

提前感谢您的帮助。

Answer 1

您的问题是错误的数据准备!!

让我们从定义变量 year、disturbance_severity、treatment、replicate 的值开始。

library(tidyverse)

set.seed(123)
yars = 2018:2021
disturbances = c(0,45,65,85)
treatments = c("B" , "T")
replicates = c("A", "B", "C", "D")
n = length(yars)*length(disturbances)*length(treatments)*length(replicates)*1
nNA=3

请注意，我首先创建了变量 yars、disturbances、treatments 和 replicates值(value)观。

然后我计算了 n 中的数据量(您可以将乘法中的最后一个值从 1 增加到例如 10)并且确定变量 nNA 中将丢失多少个值。

关键方面是函数 expand.grid(yars, disturbances, treatments, replicates) 的使用，它将返回具有正确值分布的适当表格.

查看 expand.grid 返回的前几行。

    Var1 Var2 Var3 Var4
1   2018    0    B    A
2   2019    0    B    A
3   2020    0    B    A
4   2021    0    B    A
5   2018   45    B    A
6   2019   45    B    A
7   2020   45    B    A
8   2021   45    B    A
9   2018   65    B    A
10  2019   65    B    A
11  2020   65    B    A
12  2021   65    B    A
13  2018   85    B    A
14  2019   85    B    A
15  2020   85    B    A
16  2021   85    B    A
17  2018    0    T    A
18  2019    0    T    A

这很重要。下一步就在眼前。我们创建一个 tibble 序列并将其放入 aov 函数中。

data = tibble(sample(c(rnorm(n-nNA, mean = 6, sd = 1), rep(NA, nNA)), n)) %>% 
  mutate(expand.grid(yars, disturbances, treatments, replicates)) %>% 
  rename_with(~c("all_vai", "year", "disturbance_severity", "treatment", "replicate"))

VAImodel1 <- aov(all_vai ~ disturbance_severity*treatment*year + 
                   Error(replicate/disturbance_severity/treatment/year), data = data)
summary(VAImodel1)

输出

Error: replicate
                     Df Sum Sq Mean Sq F value Pr(>F)
disturbance_severity  1 0.1341  0.1341   0.093  0.811
treatment             1 0.0384  0.0384   0.027  0.897
Residuals             1 1.4410  1.4410               

Error: replicate:disturbance_severity
                     Df Sum Sq Mean Sq F value Pr(>F)
disturbance_severity  1 0.1391  0.1391   0.152  0.763
treatment             1 0.1819  0.1819   0.199  0.733
year                  1 1.4106  1.4106   1.545  0.431
Residuals             1 0.9129  0.9129               

Error: replicate:disturbance_severity:treatment
          Df Sum Sq Mean Sq F value Pr(>F)
treatment  1 0.4647  0.4647   0.698  0.491
year       1 0.8127  0.8127   1.221  0.384
Residuals  2 1.3311  0.6655               

Error: replicate:disturbance_severity:treatment:year
               Df Sum Sq Mean Sq F value Pr(>F)
treatment       1  2.885  2.8846   3.001  0.144
year            1  0.373  0.3734   0.388  0.560
treatment:year  1  0.002  0.0015   0.002  0.970
Residuals       5  4.806  0.9612               

Error: Within
                Df Sum Sq Mean Sq F value Pr(>F)  
treatment        1   0.03   0.031   0.039 0.8430  
year             1   1.29   1.292   1.662 0.2002  
treatment:year   1   4.30   4.299   5.532 0.0206 *
Residuals      102  79.26   0.777                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

现在没有model is singular错误!!