Модель смешанных эффектов в R с использованием lme
Я анализирую данные для наблюдательного исследования роста форели. Таким образом, в датафрейме много NA. Кроме того, разные методы лечения имеют разное количество наблюдений. Так я думаю это можно назвать несбалансированным дизайном?
Year Site FishID Size Cover AGE L1
2010 LT1 10_LT1 _ 11 Large Heavy 2 5.88
2010 LT3 10_LT3 _ 14 Large Heavy 2 5.228571429
2010 SO4 10_SO4 _ 8 Small Open 0 NA
2010 SO5 10_SO5 _ 22 Small Open 0 NA
2011 LT1 11_LT1 _ 14 Large Heavy 1 6.44
2011 LT1 11_LT1 _ 15 Large Heavy 1 6.25
2011 LT1 11_LT1 _ 16 Large Heavy 1 6.421052632
2011 LT1 11_LT1 _ 18 Large Heavy 1 7.74
2011 SO5 11_SO5 _ 6 Small Open 1 7.7625
2011 SO5 11_SO5 _ 8 Small Open 1 6.914285714
2011 SO5 11_SO5 _ 13 Small Open 1 6.5
2011 SO5 11_SO5 _ 16 Small Open 1 7.2
2011 ST1 11_ST1 _ 21 Small Heavy 0 NA
2011 ST2 11_ST2 _ 10 Small Heavy 0 NA
2011 ST2 11_ST2 _ 5 Small Heavy 0 NA
2011 ST3 11_ST3 _ 20 Small Heavy 0 NA
2011 ST4 11_ST4 _ 5 Small Heavy 0 NA
2011 ST1 11_ST1 _ 9 Small Heavy 1 7.521428571
2011 ST1 11_ST1 _ 17 Small Heavy 1 8.169230769
2011 ST1 11_ST1 _ 20 Small Heavy 1 7.03125
Мои фиксированные эффекты: размер потока: 2 уровня покрытия потока: 2 уровня, год: 2 уровня и возраст: 3 уровня.
Я также случайно выбрал несколько сайтов для проверки размера потока и эффектов покрытия потока. Таким образом, я предполагаю, что сайт классифицируется как случайный эффект?
L1 является прокси для скорости роста.
На этом этапе моя максимальная модель выглядит так:
m1=lme(L1~Year*Age*Size*Cover, random=~1|Site ,data=Trout_Growth,method="ML",na.action=na.exclude)
qqnorm(residuals(m1))
qqline(residuals(m1)) # Normality OK
plot(density(residuals(m1),na.rm=T)) ***## heteroscedasticity present (lme should handle this??) ##***
summary(m1)
anova(m1,test=T,type="marginal") ***## marginal used because of unbalanced design?? ##***
Я, конечно, хорошо посмотрел этот веб-сайт, и это здорово, но я все еще не уверен в правильности обозначения для включения случайного фактора сайта. (Кажется, что я и Лмер используют разные обозначения. Это правильно?)
Полагаю, я спрашиваю, является ли это правильной отправной точкой для уточнения моей модели. Я начинаю по-настоящему наслаждаться R, но без регулярного одобрения в отношении кода R возникает соблазн вернуться к удобному программному обеспечению для статистики на базе Windows.
Любые мнения или советы будут с благодарностью.
Diarm