# Problems with bimodal variable: lmers give boundary is singular, inverse gaussian glmers?

Sorry for the fomatting. I have a non-independent “sisters and brothers” (animals from a single nest) variable (cha) following a bimodal distribution (see below). I am trying to “explain” it with: condition="treatment" (2 of them, either 0 or 1) and sex (either male or female coded as 0 or 1, respectively). This is an ecological study, so there is no control over variables. There is a problem: all of condition 0 corresponds to males except for a single female, while condition 1 corresponds to all females. I have tried linear mixed models (lmer) with the lmerTest (lme4) package. I usually diagnose graphically, with cAIC4 and performance packages (see below). All combinations of condition (Tx) and sex give: “boundary is singular” problem. This is due to the nest (grouping) variable apparently not contributing to the variance (see below). Thus, I also tried generalized linear mixed models (glmer). The one that seems to work “best” is guassian with inverse link (see below). However, I suspect this is still incorrect. Any help would be appreciated.

> c=read_excel("CHA.xlsx")

> View(c)

Tx nest sex gender cha
0 2 1 female 0.053
0 2 0 male 0.018
0 2 0 male 0.109
0 2 0 male 0.004
0 3 0 male 0.030
0 3 0 male 0.026
0 3 0 male 0.164
0 3 0 male 0.008
1 15 1 female 0.116
1 15 1 female 0.067
1 15 1 female 0.093
1 15 1 female 0.053
1 16 1 female 0.094
1 16 1 female 0.092
1 16 1 female 0.206
1 16 1 female 0.138
> cha1<-lmer(sqrt(cha)~(1|nest)+Tx+sex,data=c,REML=F)

boundary (singular) fit: see ?isSingular

> cha1

Linear mixed model fit by maximum likelihood  ['lmerModLmerTest']

Formula: sqrt(cha) ~ (1 | nest) + Tx + sex

Data: c

|   AIC   |    BIC  | logLik| deviance|df.resid|
|---------|---------|-------|---------|--------|
|-139.5950|-127.9414|74.7975|-149.5950|   71   |

Random effects:

| Groups |   Name    | Std.Dev.|
|--------|-----------|---------|
| nest   |(Intercept)| 0.00000 |
|Residual|           | 0.09044 |

Number of obs: 76, groups:  nest, 20

Fixed Effects:

|(Intercept)|   Tx     |   sex   |
|-----------|----------|---------|
|  0.205551 | 0.009486 | 0.024667|

optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

> cha2<-lmer(sqrt(cha)~(1|nest)+Tx*sex,data=c,REML=F)

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

boundary (singular) fit: see ?isSingular

> cha3<-lmer(sqrt(cha)~(1|nest)+Tx:sex,data=c,REML=F)

boundary (singular) fit: see ?isSingular

> cha4<-lmer(sqrt(cha)~(1|nest)+Tx,data=c,REML=F)

boundary (singular) fit: see ?isSingular

> cha5<-lmer(sqrt(cha)~(1|nest)+sex,data=c,REML=F)

boundary (singular) fit: see ?isSingular

> cha6<-lmer(sqrt(cha)~(1|nest),data=c,REML=F)

boundary (singular) fit: see ?isSingular

> compare_performance(cha1,cha2,cha3,cha4,cha5,cha6)

Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.

# Comparison of Model Performance Indices
|Name| Model         |   AIC  |AIC wei|   BIC  |BIC wei|R2 (cond.)|R2 (marg.)| RMSE |Sigma |
|----|---------------|--------|-------|--------|-------|----------|----------|------|------|---
|cha1|lmerModLmerTest|-139.595| 0.084 |-127.941| 0.021 |          |   0.035  |0.090 | 0.090|
|cha2|lmerModLmerTest|-139.595| 0.084 |-127.941| 0.021 |          |   0.035  |0.090 | 0.090|
|cha3|lmerModLmerTest|-141.523| 0.220 |-132.200| 0.176 |          |   0.034  |0.090 | 0.090|
|cha4|lmerModLmerTest|-141.523| 0.220 |-132.200| 0.176 |          |   0.034  |0.090 | 0.090|
|cha5|lmerModLmerTest|-141.584| 0.227 |-132.261| 0.181 |          |   0.034  |0.090 | 0.090|
|cha6|lmerModLmerTest|-140.961| 0.166 |-133.969| 0.426 |          |   0.000  |0.092 | 0.092|

boundary (singular) fit: see ?isSingular

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

boundary (singular) fit: see ?isSingular

boundary (singular) fit: see ?isSingular

boundary (singular) fit: see ?isSingular

> cAIC(cha12)

Error: \$ operator is invalid for atomic vectors
Warning message:
In deleteZeroComponents.merMod(m) :
Model has no random effects variance components larger than zero.

> anocAIC(cha7,cha8,cha9,cha10,cha11)

|           model                 |  cll |   df   | cAIC  |Refit|
|:--------------------------------|:----:|:------:|:-----:|----:|
|sqrt(cha) ~ (1 / nest) + Tx:sex  | 74.76| 1073.12|1996.71|FALSE|
|sqrt(cha) ~ (1 / nest) + Tx * sex| 74.80|  768.75|1387.90|FALSE|
|sqrt(cha) ~ (1 / nest) + Tx + sex| 74.80|  768.75|1387.90|FALSE|
|sqrt(cha) ~ (1 / nest) + Tx      | 74.76| 1073.12|1996.71|FALSE|
|sqrt(cha) ~ (1 / nest) + sex     | 74.79|  887.67|1625.76|FALSE|

> compare_performance(cha7,cha8,cha9,cha10,cha11,cha12)

Random effect variances not available. Returned R2 does not account for random effects.

# Comparison of Model Performance Indices
|Name|  Model  | AIC  |AIC wei|  BIC |BIC wei|R2 (cond.)|R2 (marg.)| RMSE|Sigma| ICC     |
|----|---------|------|-------|------|-------|----------|----------|-----|-----|---------|------
|cha7|glmerMod |8.622 | 0.156 |17.945| 0.084 |  0.934   |  0.934   |0.090|0.090|5.868e-10|
|cha8|glmerMod |10.622| 0.057 |22.275| 0.010 |  0.936   |  0.936   |0.090|0.090|5.830e-08|
|cha9|glmerMod |10.622| 0.057 |22.275| 0.010 |  0.936   |  0.936   |0.090|0.090|5.830e-08|
|cha10|glmerMod|8.622 | 0.156 |17.945| 0.084 |  0.934   |  0.934   |0.090|0.090|5.868e-10|
|cha11|glmerMod|8.622 | 0.156 |17.945| 0.084 |  0.936   |  0.936   |0.090|0.090|8.675e-10|
|cha12|glmerMod|6.643 | 0.419 |13.636| 0.728 |          |  0.000   |0.092|0.092|

> shapiro.test(residuals(cha1))

Shapiro-Wilk normality test

data:  residuals(cha1)

W = 0.95864, p-value = 0.01433

> shapiro.test(residuals(cha2))

Shapiro-Wilk normality test

data:  residuals(cha2)

W = 0.95864, p-value = 0.01433

> shapiro.test(residuals(cha3))

Shapiro-Wilk normality test

data:  residuals(cha3)

W = 0.95996, p-value = 0.01713

> shapiro.test(residuals(cha4))

Shapiro-Wilk normality test

data:  residuals(cha4)

W = 0.95996, p-value = 0.01713

> shapiro.test(residuals(cha5))

Shapiro-Wilk normality test

data:  residuals(cha5)

W = 0.95819, p-value = 0.01349

> shapiro.test(residuals(cha6))

Shapiro-Wilk normality test

data:  residuals(cha6)

W = 0.96318, p-value = 0.02658

> shapiro.test(residuals(cha7))

Shapiro-Wilk normality test

data:  residuals(cha7)

W = 0.95996, p-value = 0.01714

> shapiro.test(residuals(cha8))

Shapiro-Wilk normality test

data:  residuals(cha8)

W = 0.95864, p-value = 0.01433

> shapiro.test(residuals(cha9))

Shapiro-Wilk normality test

data:  residuals(cha9)

W = 0.95864, p-value = 0.01433

> shapiro.test(residuals(cha10))

Shapiro-Wilk normality test

data:  residuals(cha10)

W = 0.95996, p-value = 0.01714

> shapiro.test(residuals(cha11))

Shapiro-Wilk normality test

data:  residuals(cha11)

W = 0.95819, p-value = 0.01349

> shapiro.test(residuals(cha12))

Shapiro-Wilk normality test

data:  residuals(cha12)

W = 0.96318, p-value = 0.02658