我在lmer 4：：lmer()中的语法有什么问题？

Question

我正在尝试使用R和函数lmer()中的lme4包来为我的分裂图设计提供一个模型。如果我没有少的观测数据缺失，我会使用重复的度量方差，但缺失的数据不应该是线性混合效应模型的问题。

我的数据框架(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 *处理费*mean_vai+disturbance_severity: Error()模型是单数的)”

我有不同年份的观测结果，它们嵌套在不同的处理中，它们嵌套在不同的扰动严重程度中，所有这些都嵌套在副本中(这些是实验块)。因此，我尝试在lme4中使用这种结构：

library(lme4)
library(lmerTest)

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

这是我得到的错误信息：“误差:观测次数(=125个) <=随机效应数(=128个)用于术语(年份_

接下来，我尝试简化我的模型，通过删除处理变量和交互作用项，使我不会失去自由度，如下所示：

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

这一次，我得到了一个不同的错误：“边界(单数)拟合:见?奇异警告信息:模型未能收敛到一个负特征值：-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。

然后，我计算了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

现在没有模型是奇异的错误！！