#r #modeling #lme4 #mixed-models
#r #modeling #lme4 #mixed-models
Вопрос:
Я пытаюсь использовать пакет lme4 в R и функцию lmer (), чтобы подогнать модель для моего дизайна разделенного участка. Я бы использовал повторные измерения ANOVA, если бы у меня не отсутствовало небольшое количество наблюдений, но отсутствующие данные не должны быть проблемой с линейной моделью смешанных эффектов.
Мой фрейм данных ( data ) имеет простую структуру с четырьмя факторами и числовой переменной результата, называемой all_vai . Обратите внимание, что в этом примере фрейма данных не все уровни всех факторов пересекаются, хотя они были бы в моих реальных данных (за исключением отсутствующих наблюдений). Это не должно иметь значения для моего вопроса, который является попыткой исправить проблемный синтаксис.
collected_vai lt;- rnorm(125, mean = 6, sd = 1) missing lt;- rep(NA, times = 3) all_vai lt;- c(collected_vai, missing) year1 lt;- rep(2018, times = 32) year2 lt;- rep(2019, times = 32) year3 lt;- rep(2020, times = 32) year4 lt;- rep(2021, times = 32) year lt;- c(year1, year2, year3, year4) disturbance_severity lt;- rep(c(0,45,65,85), each = 32) treatment lt;- rep(c("B" , "T"), each = 64) replicate lt;- rep(c("A", "B", "C", "D"), each = 32) data = data.frame(all_vai, year, disturbance_severity, treatment, replicate) data$year lt;- as.factor(data$year) data$disturbance_severity lt;- as.factor(data$disturbance_severity) data$treatment lt;- as.factor(data$treatment) data$replicate lt;- as.factor(data$replicate) Вот модель, которую я запустил для идентичного набора данных с другим (нормально распределенным) числовым результатом и без пропущенных наблюдений-т. Е. Это модель, которую я бы запустил, если бы у меня сейчас не было несбалансированных повторных измерений из-за отсутствия данных:
VAImodel1 lt;- aov(all_vai ~ disturbance_severity*treatment*year Error(replicate/disturbance_severity/treatment/year), data = data) summary(VAImodel1) Когда я запускаю это, я получаю сообщение об ошибке: «Предупреждающее сообщение: В aov(mean_vai ~ нарушение постоянства * лечение * Год Ошибка(Репликация/нарушение постоянства/лечение/Год), : Модель ошибки() является сингулярной»
У меня есть наблюдения за разные годы, вложенные в разные методы лечения, которые вложены в разные степени тяжести нарушений, и все это вложено в копии (которые являются экспериментальными блоками). Поэтому я попытался использовать эту структуру в lme4:
library(lme4) library(lmerTest) VAImodel2 lt;- lmer(all_vai ~ (year|replicate:disturbance_severity:treatment) disturbance_severity*treatment*year, data = data) summary(VAImodel2) And this is the error message I get: «Error: number of observations (=125) lt;= number of random effects (=128) for term (Year | Replicate:disturbance_severity:treatment); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable»
Next I tried simplifying my model so that I was not running out of degrees of freedom, by removing the treatment variable and interaction term, like so:
VAImodel3 lt;- lmer(all_vai ~ (year|replicate:disturbance_severity) disturbance_severity*year, data = data) summary(VAImodel3) This time I get a different error: «boundary (singular) fit: see ?isSingular Warning message: Model failed to converge with 1 negative eigenvalue: -1.2e-01 »
Thank you in advance for any help.
Ответ №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)) %gt;% mutate(expand.grid(yars, disturbances, treatments, replicates)) %gt;% rename_with(~c("all_vai", "year", "disturbance_severity", "treatment", "replicate")) VAImodel1 lt;- 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(gt;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(gt;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(gt;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(gt;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(gt;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 Теперь в модели нет сингулярных ошибок!!
Комментарии:
1. Спасибо вам за вашу помощь!