Cross Validated Asked on January 1, 2022

My question is related with my previous post Extract variance of the fixed effect in a glmm.

However, in this case I change the model that the GLMM follow.

It follows a log family and as there are many zeros in my dataset,so I used a zero inflation method.

I would like to get the variation (variance component) in incidence (inc.) within each habitat while being mindful of random factors such as season and site

This is my data set:

```
## Incidence:
Incidence <- data.frame(Inc. = c(0.4400, 0.5102, 0.2979, 0.2667, 0.0000, 0.0000,
0.0200, 0.0213, 0.0000, 0.0238, 0.0256, 0.0000,
0.0000, 0.1538, 0.0417, 0.0000, 0.0734, 0.0000,
0.0000, 0.0000, 0.1293, 0.0072, 0.0000, 0.0078,
0.0000, 0.0000, 0.0000, 0.0068, 0.0000, 0.0000,
0.0068),
Habitat = c("Crop", "Crop", "Crop", "Crop", "Edge", "Edge",
"Edge", "Edge", "Edge", "Edge", "Edge", "Edge",
"Edge", "Edge", "Edge", "Oakwood", "Oakwood",
"Oakwood", "Oakwood", "Oakwood", "Oakwood",
"Oakwood", "Oakwood", "Wasteland", "Wasteland",
"Wasteland", "Wasteland", "Wasteland", "Wasteland",
"Wasteland", "Wasteland"),
Season = c("Summer", "Summer", "Summer", "Summer", "Autumn",
"Autumn", "Autumn", "Autumn", "Spring", "Spring",
"Spring", "Spring", "Summer", "Summer", "Summer",
"Autumn", "Autumn", "Autumn", "Autumn", "Spring",
"Spring", "Spring", "Spring", "Autumn", "Autumn",
"Autumn", "Autumn", "Spring", "Spring", "Spring",
"Spring"),
Site = c("M1", "M2", "M3", "M4", "L1", "L2", "L3", "L4",
"L1", "L2", "L3", "L4", "L1", "L2", "L3", "Q1",
"Q2", "Q3", "Q4", "Q1", "Q2", "Q3", "Q4", "E1",
"E2", "E3", "E4", "E1", "E2", "E3", "E4"))
```

With the aim to get the variation I check previously with a shapiro wilk test how is the distribution of my dataset by Rstudio.

```
shapiro.test(x = Incidence$Inc.):
Shapiro-Wilk normality test
data: Incidence$Incidence
W = 0.56708, p-value = 2.092e-08
```

Moreover I got the homocedasticity with a levene test:

```
leveneTest(y = Incidence$Inc., group = Incidence$Habitat, center = "median")
Levene's Test for Homogeneity of Variance (center = "median")
Df F value Pr(>F)
group 3 6.3481 0.002129 **
27
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```

Afterward I check how is the distribution using:

```
Input_2 <- Incidence$Inc.
library(rriskDistributions)
Prueba <- fit.cont(as.vector(t(Input_2)))
```

and I got a log distribution

Then I performed a glmm of this dataset in R:

```
GlM_habitats <- glmmTMB(Inc.~ Habitat + (1|Season)+ (1|Site),
data = Incidence,
ziformula = ~1,
family = poisson(link = "log"))
#Warning message:
#In glmmTMB(Inc.~ Habitat + (1 | Season) + (1 | Site), data = Incidence, :
#non-integer counts in a poisson model
Anova(GlM_habitats)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: Incidence
Chisq Df Pr(>Chisq)
Habitat 3.0632 3 0.382
summary(GlM_habitats)
Family: poisson ( log )
Formula: Inc.~ Habitat + (1 | Season) + (1 | Site)
Zero inflation: ~1
Data: Incidence
AIC BIC logLik deviance df.resid
23.5 33.5 -4.7 9.5 24
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Season (Intercept) 5.656e-13 7.52e-07
Site (Intercept) 1.176e-13 3.43e-07
Number of obs: 31, groups: Season, 3; Site, 16
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9710 0.8125 -1.195 0.232
HabitatEdge -2.6780 2.0382 -1.314 0.189
HabitatOakwood -2.6696 2.3290 -1.146 0.252
HabitatWasteland -4.9528 6.8841 -0.720 0.472
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -24.1 43216.9 -0.001 1
```

Then as in previous post they anwered to me I tried to extract the variance of fixed effect:

```
# Variance of random effects:
vc <- lme4::VarCorr(GlM_habitats)
print(vc,comp=c("Variance","Std.Dev."),digits=2)
Conditional model:
Groups Name Variance Std.Dev.
Season (Intercept) 5.7e-13 7.5e-07
Site (Intercept) 1.2e-13 3.4e-07
# Variance-Covariance Matrix of fixed effects:
vc_fixed <- as.matrix(vcov(GlM_habitats))
# Variance of fixed effects:
var_fixed <- diag(vc_fixed); var_fixed
[[1]]
(Intercept) HabitatEdge HabitatOakwood HabitatWasteland
(Intercept) 0.660153 -0.660153 -0.660153 -0.660153
HabitatEdge -0.660153 4.154245 0.660153 0.660153
HabitatOakwood -0.660153 0.660153 5.424338 0.660153
HabitatWasteland -0.660153 0.660153 0.660153 47.390362
# Standard errors of fixed effects:
se_fixed <- sqrt(var_fixed); se_fixed
```

When I perform this analysis I got this

```
Error in sqrt(var_fixed) : non-numeric argument to mathematical function
```

I would like to know how to interpret this result and to know if they have been performed OK. I can not believe that `Season`

and `Site`

have very low variance and the ANOVA results give a p value that is not significant. Moreover, I do not know why the Standard errors of fixed effects do not work.

What Am I doing wrong?

There are 2 main problems here:

As with other linear models there is no requirement for the outcome variable to be normally distributed in a linear mixed effects model. So

`shapiro.test(x = Incidence$Inc.)`

is a waste of time and so is any procedure that tries to find the distribution of the outcome, such as`fit.cont`

that you use - such things might be of interest to theoreticians but they are of very limited value to applied research. We would, however, like the residuals to be, at least approximately, normally distributed.You have fitted a poisson model. Poisson models are for data with a count (integer) outcome. You have a numeric variable so the first model to fit is a standard linear mixed effects model.

You have only 3 levels of

`Season`

. This should probably be a fixed effect.

So, with your data we can fit:

```
> m0 <- lmer(Inc.~ Habitat + (1|Season)+ (1|Site),
+ data = Incidence)
> summary(m0)
Linear mixed model fit by REML ['lmerMod']
Formula: Inc. ~ Habitat + (1 | Season) + (1 | Site)
Data: Incidence
REML criterion at convergence: -78.9
Scaled residuals:
Min 1Q Median 3Q Max
-1.45229 -0.30319 -0.01575 0.20558 2.53994
Random effects:
Groups Name Variance Std.Dev.
Site (Intercept) 0.0031294 0.05594
Season (Intercept) 0.0005702 0.02388
Residual 0.0008246 0.02872
Number of obs: 31, groups: Site, 16; Season, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.35450 0.03607 9.827
HabitatEdge -0.32669 0.04475 -7.301
HabitatOakwood -0.31616 0.04637 -6.818
HabitatWasteland -0.33973 0.04637 -7.326
```

and then we can inspect the residuals histgram:

```
hist(residuals(m0))
```

which looks fine. There is no need to perform a statistical test for normality.

Note that you should probably model `Season`

as a fixed effect, not random.

Answered by Robert Long on January 1, 2022

0 Asked on January 28, 2021 by chua-s-yang

0 Asked on January 28, 2021 by visionenthusiast

0 Asked on January 28, 2021 by zqq

lasso machine learning mathematical statistics ridge regression self study

1 Asked on January 28, 2021 by nlapidot

0 Asked on January 28, 2021 by tricostume

1 Asked on January 27, 2021 by mathews24

0 Asked on January 26, 2021 by thecity2

2 Asked on January 26, 2021

0 Asked on January 26, 2021 by bonesones

0 Asked on January 26, 2021 by kplauritzen

2 Asked on January 26, 2021 by user369210

0 Asked on January 26, 2021

causality confounding difference in difference econometrics panel data

1 Asked on January 25, 2021 by kiril-e-proykov

0 Asked on January 24, 2021 by oja-niva

classification machine learning random forest recommender system regression

0 Asked on January 24, 2021 by mephisto73

0 Asked on January 23, 2021 by valjean

1 Asked on January 23, 2021 by tjt

5 Asked on January 23, 2021 by user239903

0 Asked on January 21, 2021 by pooza

Get help from others!

Recent Answers

- Joshua Engel on Why fry rice before boiling?
- haakon.io on Why fry rice before boiling?
- Jon Church on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?
- Peter Machado on Why fry rice before boiling?

Recent Questions

- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?
- Does Google Analytics track 404 page responses as valid page views?

© 2023 AnswerBun.com. All rights reserved. Sites we Love: PCI Database, MenuIva, UKBizDB, Menu Kuliner, Sharing RPP, SolveDir