Post hoc для двоичного GLMM (lme4) и сюжета - PullRequest
Новая
       79

Post hoc для двоичного GLMM (lme4) и сюжета

0 голосов
/ 26 сентября 2018

Так что я новичок, пытающийся сделать GLMM и провести специальный анализ ... помогите!Я собрал двоичные данные по 9 красоткам при 6 уровнях освещенности, 1 = реакция на движение барабана оптомотора, 0 = нет реакции.Мои данные были импортированы в R с заголовками «Animal_ID, light_intensity, response».Идентификатор животного (1-9), повторяемый для каждой интенсивности света (3.36-0.61) (см. Ниже)

Используя следующий код (пакет lme4), я выполнил GLMM и нашел уровень освещенности, чтобы иметьзначительное влияние на ответ: 

d = data.frame(id = data[,1], var = data$Light_Intensity, Response = data$Response)
model <- glmer(Response~var+(1|id),family="binomial",data=d)
summary(model)

Возвращает 

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: Response ~ var + (1 | Animal_ID)
Data: d

 AIC      BIC   logLik deviance df.resid 
  66       72      -30       60       51 

Scaled residuals: 
Min      1Q  Median      3Q     Max 
-3.7704 -0.6050  0.3276  0.5195  1.2463 

Random effects:
Groups    Name        Variance Std.Dev.
Animal_ID (Intercept) 1.645    1.283   
Number of obs: 54, groups:  Animal_ID, 9

Fixed effects:
        Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -1.7406     1.0507  -1.657   0.0976 .
var           1.1114     0.4339   2.561   0.0104 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr)
var -0.846

Затем работает: 

m1 <- update(model, ~.-var)
anova(model, m1, test = 'Chisq')

Возвращает 

Data: d
Models:
m1: Response ~ (1 | Animal_ID)
model: Response ~ var + (1 | Animal_ID)
      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
m1     2 72.555 76.533 -34.278   68.555                            
model  3 66.017 71.983 -30.008   60.017 8.5388      1   0.003477 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Я установилпакеты multcomp и lsmeans в попытке выполнить Tukey post hoc, чтобы увидеть разницу, но столкнулись с трудностями в обоих случаях.

Запуск: 

summary(glht(m1,linfct=mcp("Animal_ID"="Tukey")))

Возвращает: «Ошибкав mcp2matrix (модель, linfct = linfct): переменные 'Animal_ID' были указаны в 'linfct', но не могут быть найдены в 'model'! "

Running: 

lsmeans(model,pairwise~Animal_ID,adjust="tukey")

Возвращает: «Ошибка в lsmeans.character.ref.grid (object = new (« ref.grid », model.info = list (: нет переменной с именем Animal_ID в эталонной сетке»

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

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

Заранее большое спасибо!

ОБНОВЛЕНИЕ:

Новый код - 

Light <- c("3.36","3.36","3.36","3.36","3.36","3.36","3.36","3.36","3.36","2.98","2.98","2.98","2.98","2.98","2.98","2.98","2.98","2.98","2.73","2.73","2.73","2.73","2.73","2.73","2.73","2.73","2.73","2.15","2.15","2.15","2.15","2.15","2.15","2.15","2.15","2.15","1.72","1.72","1.72","1.72","1.72","1.72","1.72","1.72","1.72","0.61","0.61","0.61","0.61","0.61","0.61","0.61","0.61","0.61")
Subject <- c("1","2","3","4","5","6","7","8","9","1","2","3","4","5","6","7","8","9","1","2","3","4","5","6","7","8","9","1","2","3","4","5","6","7","8","9","1","2","3","4","5","6","7","8","9","1","2","3","4","5","6","7","8","9")
Value <- c("1","0","1","0","1","1","1","0","1","1","0","1","1","1","1","1","1","1","0","1","1","1","1","1","1","0","1","0","0","1","1","1","1","1","1","1","0","0","0","1","0","0","1","0","1","0","0","0","1","1","0","1","0","0")

data <- data.frame(Light, Subject, Value)

library(lme4)

model <- glmer(Value~Light+(1|Subject),family="binomial",data=data)
summary(model)

Возвраты: 

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: Value ~ Light + (1 | Subject)
   Data: data

     AIC      BIC   logLik deviance df.resid 
    67.5     81.4    -26.7     53.5       47 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6564 -0.4884  0.2193  0.3836  1.2418 

Random effects:
 Groups  Name        Variance Std.Dev.
 Subject (Intercept) 2.687    1.639   
Number of obs: 54, groups:  Subject, 9

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.070e+00  1.053e+00  -1.016   0.3096  
Light1.72   -7.934e-06  1.227e+00   0.000   1.0000  
Light2.15    2.931e+00  1.438e+00   2.038   0.0416 *
Light2.73    2.931e+00  1.438e+00   2.038   0.0416 *
Light2.98    4.049e+00  1.699e+00   2.383   0.0172 *
Light3.36    2.111e+00  1.308e+00   1.613   0.1067  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr) Lg1.72 Lg2.15 Lg2.73 Lg2.98
Light1.72 -0.582                            
Light2.15 -0.595  0.426                     
Light2.73 -0.595  0.426  0.555              
Light2.98 -0.534  0.361  0.523  0.523       
Light3.36 -0.623  0.469  0.553  0.553  0.508

Затем выполняется: 

m1 <- update(model, ~.-Light)
anova(model, m1, test= 'Chisq')

Возвраты: 

Data: data
Models:
m1: Value ~ (1 | Subject)
model: Value ~ Light + (1 | Subject)
      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
m1     2 72.555 76.533 -34.278   68.555                           
model  7 67.470 81.393 -26.735   53.470 15.086      5       0.01 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Наконец, запустив: 

library(lsmeans)
lsmeans(model,list(pairwise~Light),adjust="tukey")

Возвращает (теперь это действительно работает!): 

$`lsmeans of Light`
 Light    lsmean       SE df  asymp.LCL asymp.UCL
 0.61  -1.070208 1.053277 NA -3.1345922 0.9941771
 1.72  -1.070216 1.053277 NA -3.1345997 0.9941687
 2.15   1.860339 1.172361 NA -0.4374459 4.1581244
 2.73   1.860332 1.172360 NA -0.4374511 4.1581149
 2.98   2.978658 1.443987 NA  0.1484964 5.8088196
 3.36   1.040537 1.050317 NA -1.0180467 3.0991215

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 

$`pairwise differences of contrast`
 contrast         estimate       SE df z.ratio p.value
 0.61 - 1.72  7.933829e-06 1.226607 NA   0.000  1.0000
 0.61 - 2.15 -2.930547e+00 1.438239 NA  -2.038  0.3209
 0.61 - 2.73 -2.930539e+00 1.438237 NA  -2.038  0.3209
 0.61 - 2.98 -4.048866e+00 1.699175 NA  -2.383  0.1622
 0.61 - 3.36 -2.110745e+00 1.308395 NA  -1.613  0.5897
 1.72 - 2.15 -2.930555e+00 1.438239 NA  -2.038  0.3209
 1.72 - 2.73 -2.930547e+00 1.438238 NA  -2.038  0.3209
 1.72 - 2.98 -4.048874e+00 1.699175 NA  -2.383  0.1622
 1.72 - 3.36 -2.110753e+00 1.308395 NA  -1.613  0.5897
 2.15 - 2.73  7.347728e-06 1.357365 NA   0.000  1.0000
 2.15 - 2.98 -1.118319e+00 1.548539 NA  -0.722  0.9793
 2.15 - 3.36  8.198019e-01 1.302947 NA   0.629  0.9889
 2.73 - 2.98 -1.118326e+00 1.548538 NA  -0.722  0.9793
 2.73 - 3.36  8.197945e-01 1.302947 NA   0.629  0.9889
 2.98 - 3.36  1.938121e+00 1.529202 NA   1.267  0.8029

Результаты приведены по шкале отношения шансов в логах (не ответа).Корректировка значения P: метод Тьюки для сравнения семейства из 6 оценок

1 Ответ

0 голосов
/ 29 сентября 2018

Ваша модель указывает Animal_ID как случайный эффект.Функции glht и lsmeans работают только для сравнений с фиксированным эффектом.

...