R: Как улучшить соответствие общей линейной смешанной модели в LME4?

Question

Вопрос переполнения стека - подбор модели и выбор для glmer() в lme4

Я все еще относительно новичок в этом, и перепробовал все, что мог придумать и погуглить, но, похоже, врезался в стену, поэтому заранее извините, если что-то из этого является базовым / неправильным. Я не знаю, следует ли разбить это на несколько вопросов, но, возможно, более полезно представить его как один процесс. Что я действительно хочу знать, так это как вы можете сказать, что созданная вами модель достаточно хороша для проверки гипотез и выработки выводов.

Я выполняю анализ большого набора данных по динамике сеянцев (20258 сеянцев с 7467 участков), который я не могу предоставить по соображениям конфиденциальности, поскольку это не мой набор данных. В этом анализе я пытаюсь определить, как несколько объясняющих переменных (длина ствола, функциональная группа, конспецифическая и гетероспецифическая плотность, а также взаимодействия между функциональными группами и плотностью) влияют на вероятность выживания для данного сеянца - используя функцию glmer() в Ime4. Перед использованием модели я масштабировал все непрерывные переменные, используя функцию scale(), и убедился, что все непостоянные переменные были перечислены как факторы, и поиграл с настройками оптимизатора и номером итерации, чтобы избавиться от некоторых стандартных предупреждений оптимизатора.

Модель структурирована так (ЗНАЧЕНИЕ AIC НИЖЕ МОДЕЛИ):

themod <- glmer(survival ~ stemlength + group + conspecific density + 
heterospecific density + group:conspecifics + group:heterospecifics 
+group:stemlength + (1|plot) + (1 + conspecifics + heterospecifics|species), 
family=binomial(link = "logit"), data = dataset, 
control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))

AIC: 14040 (<2e-16) R в квадрате установлено: 0,25441404 R в квадрате случайное число: 0,748984

Я попробовал несколько версий этих моделей на одном наборе данных и сравнил значения AIC, используя anova (по сравнению со значением базовой модели 14128,7, перечисленным ниже) и дроптермный анализ (который показал, что все переменные либо улучшили соответствие модели, либо не оказали никакого влияния либо кстати), и этот лучше всего подходил к данным. Я также изучил использование значений псевдо R в квадрате (как описано в Nakawaga et al., 2013), что, по-видимому, указывает на то, что другие модели лучше подходят, но я придерживался значений AIC. См. Другие модели ниже - каждая из них указана с квадратами AIC и R.

basemod <- glmer(survival ~ stemlength + conspecifics + heterospecifics + 
(1|plot) + (1|species),          family=binomial(link = "logit"), data = 
dataset, control=glmerControl(optimizer="bobyqa",     
optCtrl=list(maxfun=2e5)))

AIC: 14128,7 R в квадрате установлены: 0,09298017 R в квадрате случайные: 0,489361

mod1 <- glmer(survival ~ stemlength + conspecifics + heterospecifics + group 
+ group:conspecifics + group:heterospecifics + (1|plot) + (1 + conspecifics 
+ heterospecifics|species), family=binomial, data = dataset, 
control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))

AIC: 14086 (<2e-16) R в квадрате установлено: 0,27126705 R в квадрате случайно: 0,730059

mod2 <- glmer(survival ~ stemlength + conspecifics + heterospecifics + group 
+ (1|plot) + (1 + conspecifics + heterospecifics|species), family=binomial, 
data = dataset, control=glmerControl(optimizer="bobyqa", 
optCtrl=list(maxfun=2e5)))

AIC: 14128,7(<2e-16) R в квадрате: 0,35435996 R в квадрате: 0,798788

Выбрав формулу для моей модели (здесь она указана как "Теодод"), я начал проверять соответствие. Следуя совету Bolker et al., 2009, я хотел посмотреть, были ли данные однородными по категориям, но, учитывая, что у меня 7467 участков и 77 видов, разделенных на две группы для 20258 сеянцев, это кажется невозможным? Возможно, было бы хорошо исключить эти очень редкие виды, т.е. те, у которых очень мало особей? Я мог бы также сгруппировать графики, чтобы уменьшить общее число, но поскольку я использую общую линейную смешанную модель с графиком в качестве случайного эффекта, я думаю, что группировка графиков не будет иметь значения?

Я также хотел посмотреть, были ли "отклики преобразованных данных линейными по отношению к непрерывным предикторам" - из Bolker et al., 2009, и я решил, что непрерывные параметры линейны в пространстве логарифмов - или что результат модель (logodds), нанесенная на график против каждой из трех непрерывных переменных (длина ствола, специфическая плотность и гетероспецифическая плотность), показывает линейную зависимость. Для этого я использовал этот код для создания этих графиков (нажмите на ссылку):

modvalues <- as.data.frame(cbind(predict(themod))) 
plot(modvalues$V1, dataset$stemlength, , xlab = "Log odds", ylab = "Stem 
length")
abline(lm(dataset$stemlength~modvalues$V1), col="red") 
plot(modvalues$V1, dataset$conspecifics, xlab = "Log odds", ylab = 
"Conspecific densities")
abline(lm(dataset$conspecifics~modvalues$V1), col="red")  
plot(modvalues$V1, dataset$heterospecifics, xlab = "Log odds", ylab = 
"Heterospecific densities") 
abline(lm(dataset$heterospecifics~modvalues$V1), col="red")

Рисунки, показывающие взаимосвязь между логарифмическими коэффициентами и непрерывными параметрами длины ствола, конспецифической плотности и гетероспецифической плотности

Тем не менее, они показывают, возможно, почти линейную зависимость для длины стебля и удельной плотности, но, конечно, не для гетероспецифической плотности. В-третьих, я хотел посмотреть, есть ли остаточные выбросы, используя дополнительные диагностические графики, сгенерированные с использованием этого кода:

plot(themod, id=0.05,idLabels=~.obs) # shows us outliers

Я удалил эти выбросы, используя этот код:

dfs <- dataset
dfs$resid <- resid(themod)
dfs$fitted <- fitted(themod)
z <- resid(m3)
z <- z[abs(z-mean(z))<2*sd(z)]
newdat <- dfs[dfs$resid %in% z, ] 
themod_noout <- update(themod, data=newdat) 

Это улучшило ситуацию, но оставило некоторые дополнительные выбросы, которые я пытался удалить с помощью этого кода:

newdat2 <- newdat
newdat2 <- newdat2[order(newdat2$resid), ]
newdat2 <- tail(newdat2, -10)
newdat2 <- head(newdat2, -10)
themod_noout2 <- update(themod_noout,data=newdat2)

Но это привело к очень странной схеме при рассмотрении остатков? что улучшило линейность в пространстве логарифмов для длины ствола и, возможно, конспецифичной плотности - но не гетероспецифической плотности (см. рисунки - нажмите ссылку ниже).

Остаточные графики, показывающие выбросы исходной модели (их), модель после одной загрузки выбросов и модель после двух нагрузок.

Наконец, я хотел знать, соответствует ли распределение внутри групп предполагаемому распределению (биномиальному) - что опять-таки решает проблему наличия очень большого набора данных со многими категориями, но, учитывая, что выживание может быть только биномиальным, возможно, это не проблема?

Но я поставлю это как второй вопрос, так как больше не могу публиковать ссылки на изображения.

(Я также ознакомился с проработанными примерами из Bolker, 2015 г. (я не могу опубликовать ссылку, поскольку у меня нет репутации - но все они предназначены для гораздо меньших наборов данных, и я изо всех сил стараюсь применить весь этот метод к мой процесс).

Итак, мои вопросы: как выбрать из моей исходной модели формул установки / структуры - используя AIC, dropterm, R в квадрате или диагностические графики остаточных распределений?

Как я должен сравнивать однородность по категориям и должен ли я уменьшить свою выборочную сторону или исключить редкие виды?

Что происходит с моими выбросами?

Достаточно ли хороши какие-либо из этих моделей, и как мне сказать?

Наконец, если нет - что я могу сделать, чтобы улучшить их, и как я мог бы реструктурировать свою модель, чтобы сделать это?

Огромное спасибо всем, у кого есть время и усилия, чтобы пройти через все это, и еще раз, заранее извиняюсь за мое невежество. Ура!