Раздутый DF в результатах lsmeans для модели lmer

Question

Я использовал lmer из пакета lme4 для запуска линейной модели смешанных эффектов. У меня есть данные о температуре за 3 года для необработанных (5) и обработанных участков (10). Модель:

modela<-lmer(ave~yr*tr+(1|pl), REML=FALSE, data=mydata)

Модель проверена на нормальность остатков; qqnorm plot Мои данные:

'data.frame':   6966 obs. of  7 variables:

$ yr : Factor w/ 3 levels "yr1","yr2","yr3": 1 1 1 1 1 1 1 1 1 1 ...

$ pl : Factor w/ 15 levels "C02","C03","C05",..: 1 1 1 1 1 1 1 1 1 1 ...

$ tr : Factor w/ 2 levels "Cont","OTC": 1 1 1 1 1 1 1 1 1 1 ...

$ ave: num  14.8 16.1 11.6 10.3 11.6 ...

Взаимодействие является значительным, поэтому я использовал lsmeans:

lsmeans(modela, pairwise~yr*tr, adjust="tukey")

В контрасты я получаю (только выдержки)

contrast                estimate        SE      df t.ratio p.value

yr1,Cont - yr2,Cont -0.727102895 0.2731808 6947.24  -2.662  0.0832

yr1,OTC - yr2,OTC   -0.990574030 0.2015650 6449.10  -4.914  <.0001

yr1,Cont - yr1,OTC  -0.005312771 0.3889335   31.89  -0.014  1.0000

yr2,Cont - yr2,OTC  -0.268783907 0.3929332   32.97  -0.684  0.9825

Мой вопрос касается высоких значений dfs для некоторых контрастов и связанных с ними, но бессмысленных низких значений p.

Может ли это быть связано с:

-присутствие АН в моем наборе данных (некоторое улучшение после удаления)

- неравные размеры выборки (например, 5 из одного лечения, 10 из другого - однако, те (год-1, продолжение-год-1, ОТС), кажется, не проблема.

Другие вопросы?

Я искал вопросы о stakverflow и проверял их.

Спасибо за любые ответы, идеи, комментарии.

Answer 1

В этом примере обработки назначаются экспериментально на участки. Наличие небольшого количества графиков, назначенных для лечения, серьезно ограничивает доступную информацию для статистического сравнения обработок. (Если бы у вас был только один участок на обработку, было бы невозможно даже сравнить обработки, потому что вы не смогли бы отделить эффект обработок от эффекта графиков.) У вас есть 10 графиков, назначенных одному лечение и 5 к другому. С точки зрения основного эффекта от лечения, таким образом, у вас есть (10-1)+(5-1) = 13 df для основного эффекта лечения, и если вы делаете

lsmeans(modela, pairwise ~ tr)

вы увидите около 13 df (возможно, меньше из-за дисбаланса и отсутствия) для этой статистики. Когда вы сравниваете комбинации лет и лечения, вы получаете примерно в 3 раза больше df, потому что есть 3 года. Однако в некоторых из этих сравнений год является одинаковым в каждой сравниваемой комбинации, и в этих сравнениях различия в графиках в основном сводятся на нет (это сравнение внутри графика); и в этих случаях df в основном происходит от остаточной ошибки для модели, которая имеет тысячи df. Из-за дисбаланса в данных эти сравнения немного загрязняются межплановыми вариациями, делая df несколько меньшим, чем остаточный дф

Похоже, вы не особенно заинтересованы в перекрестных сравнениях, таких как лечение1, год1 против лечения2, год3. Я предлагаю использовать переменные "by", чтобы сократить количество проверенных сравнений, потому что, когда вы проверяете их все, коррекция множественности излишне консервативна. Было бы что-то вроде этого:

modela.lsm = lsmeans(modela, ~ tr * yr)
pairs(modela.lsm, by = "yr")   # compare tr for each yr
pairs(modela.lsm, by = "tr")   # compare yr for each tr

Эти вызовы будут применять поправку Тьюки отдельно для каждой группы "по". Если вы хотите исправить множественность для каждой семьи, сделайте следующее:

rbind(pairs(modela.lsm, by = "yr"))
rbind(pairs(modela.lsm, by = "tr"))

По умолчанию используется многомерная коррекция t (Тьюки здесь не правильный метод). Вы даже можете сделать

rbind(pairs(modela.lsm, by = "yr"), pairs(modela.lsm, by = "tr"))

сгруппировать все сравнения в одну семью и применить многомерную поправку t.