Как исправить «единственное соответствие» с glmer (lme4) в R? - PullRequest
0 голосов
/ 12 апреля 2019

Задача

Я пытаюсь подобрать модели glmer с переменными, варьирующимися от 0 до 1, используя lme4 в R, но я всегда получаю ошибку «единственного соответствия». Я пробовал разные вещи, но до сих пор было невозможно избавиться от этой ошибки. (У меня точно такая же ошибка для другой модели, в которой переменные варьируются от -1 до 1)

Фоны и вещи, которые я пробовал

Я хочу предсказать, как движение между птицами в день X (mvt_index) предсказать индекс ассоциации в день X + 1 (sbs_nextday). У меня есть 6 стад (вольеров) из 28 птиц из 3 линий (высокий, низкий, дикий тип) на период 28 дней. Вот 5 случайно выбранных строк моего набора данных (fulldata_mod1.2):


             date pair_id aviaries lines sbs_nextday mvt_index
57552  2017-06-15   70_83       G3     W 0.007614837 0.0000000
118075 2017-06-24 158_159       G6     H 0.001187648 0.0000000
13727  2017-06-28 125_143       G2     L 0.001333333 0.0000000
106103 2017-06-18   65_82       G3     W 0.002484913 0.1666667
24790  2017-07-06   41_52       G5     L 0.004279601 0.2500000

Индекс ассоциации (sbs_nextday) и индекс движения (mvt_index) варьируются между 0 и 1. В следующей модели (mod1.2) каждая точка данных представляет значение индекса ассоциации и индекса движения по диаде (pair_id) и по дням ( Дата). Первая модель, которую я запускал, была такой:

mod1.2 <- lmer(sbs_nextday ~ mvt_index  + (1|pair_id) + (1|date) + (1|aviaries) + (1|lines),data=fulldata_mod1.2, lmerControl(optCtrl = list(maxfun = 10000)))


Error in if (REML) p else 0L : argument is not interpretable as logical
In addition: Warning message:
In if (REML) p else 0L :
  the condition has length > 1 and only the first element will be used

Подумав немного о моей модели, я решил убрать дату как случайный эффект, потому что, поскольку у меня есть одно значение на диаду и на день, pair_id и date будут объяснять одну и ту же дисперсию, если оба будут использоваться в качестве случайных эффектов. Затем каждый вольер имеет свой набор индивидуальных идентификаторов (и pair_id). Следовательно, я использовал pair_id как вложенный случайный эффект в вольерах. Наконец, поскольку я обрезал данные (от 0 до 1), я не могу использовать lmer, потому что он будет пытаться соответствовать гауссовскому распределению, и это невозможно с моими данными. Кроме того, у меня есть набор данных с нулевым раздувом, поэтому я использовал модель glmer с биномиальным семейством и функцией logit link следующим образом:

mod1.2 <- glmer(cbind(sbs_nextday, (1-sbs_nextday)) ~ mvt_index  + (1|aviaries/pair_id), data=fulldata_mod1.2, family=binomial(link='logit'))

boundary (singular) fit: see ?isSingular
Warning message:
In eval(family$initialize, rho) : non-integer counts in a binomial glm!

summary(mod1.2)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: cbind(sbs_nextday, (1 - sbs_nextday)) ~ mvt_index + (1 | aviaries/pair_id)
   Data: fulldata_mod1.2

     AIC      BIC   logLik deviance df.resid 
   392.5    426.7   -192.2    384.5    39069 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.0406 -0.0076  0.0700  0.4228 30.1497 

Random effects:
 Groups           Name        Variance Std.Dev.
 pair_id:aviaries (Intercept) 0        0       
 aviaries         (Intercept) 0        0       
Number of obs: 39073, groups:  pair_id:aviaries, 2254; aviaries, 6

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -7.6670     0.2554 -30.024   <2e-16 ***
mvt_index     1.2567     0.5989   2.098   0.0359 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
mvt_index -0.577
convergence code: 0
boundary (singular) fit: see ?isSingular

Я посмотрел? IsSingular и проверил матрицу дисперсии-ковариации:

VarCorr(mod1.2)
 Groups           Name        Std.Dev.
 pair_id:aviaries (Intercept) 0       
 aviaries         (Intercept) 0   

Затем я проверил, сбалансирован ли мой набор данных для каждой переменной в модели:

lines:
table(fulldata_mod1.2$lines)
    H     L     W 
12171 13988 12914

aviaries:
table( fulldata_mod1.2$aviaries)
  G1   G2   G3   G4   G5   G6 
5176 6822 7193 5721 7166 6995

date:
table( fulldata_mod1.2$date)
2017-06-12 2017-06-13 2017-06-14 2017-06-15 2017-06-16 2017-06-17 2017-06-18 2017-06-19 2017-06-20 2017-06-21 
      1113       1914       1995       1839       1861       1788       1652       1703       1574       1398 
2017-06-22 2017-06-23 2017-06-24 2017-06-25 2017-06-26 2017-06-27 2017-06-28 2017-06-29 2017-06-30 2017-07-01 
      1854       1497       1578       1457       1204       1494       1179       1128        826       1220 
2017-07-02 2017-07-03 2017-07-04 2017-07-05 2017-07-06 2017-07-07 2017-07-08 2017-07-09 
       486       1300       1049       1087       1471       1054       1256       1096 

Затем я использовал glm, чтобы проверить, работает ли модель без случайных эффектов:

mod1.4 <- glm(cbind(sbs_nextday, (1-sbs_nextday)) ~ mvt_index   ,data=fulldata_mod1.2,family=binomial(link='logit'))
Warning message:
In eval(family$initialize) : non-integer counts in a binomial glm!

summary(mod1.4)

Call:
glm(formula = cbind(sbs_nextday, (1 - sbs_nextday)) ~ mvt_index, 
    family = binomial(link = "logit"), data = fulldata_mod1.2)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.25747  -0.17050  -0.13572  -0.04642   2.17994  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.22387    0.04689 -90.083  < 2e-16 ***
mvt_index    0.83360    0.12672   6.578 4.76e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2011.4  on 39072  degrees of freedom
Residual deviance: 1973.3  on 39071  degrees of freedom
AIC: 1488.1

Number of Fisher Scoring iterations: 7

Так что, когда я сравниваю резюме glm и glmer, я получаю примерно одинаковые результаты. Но я хочу контролировать случайные эффекты, и модель все еще не работает. Я буду очень признателен, если вы сможете дать мне несколько советов, и я уже благодарю вас за то время, которое вы потратите на это.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...