# Extract variance of the fixed effect in a glmm

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:

 Inc.    Habitat  Season Site
0.4400  Crop        Summer  M1
0.5102  Crop        Summer  M2
0.2979  Crop        Summer  M3
0.2667  Crop        Summer  M4
0.0000  Edge        Autumn  L1
0.0000  Edge        Autumn  L2
0.0200  Edge        Autumn  L3
0.0213  Edge        Autumn  L4
0.0000  Edge        Spring  L1
0.0238  Edge        Spring  L2
0.0256  Edge        Spring  L3
0.0000  Edge        Spring  L4
0.0000  Edge        Summer  L1
0.1538  Edge        Summer  L2
0.0417  Edge        Summer  L3
0.0000  Oakwood     Autumn  Q1
0.0734  Oakwood     Autumn  Q2
0.0000  Oakwood     Autumn  Q3
0.0000  Oakwood     Autumn  Q4
0.0000  Oakwood     Spring  Q1
0.1293  Oakwood     Spring  Q2
0.0072  Oakwood     Spring  Q3
0.0000  Oakwood     Spring  Q4
0.0078  Wasteland   Autumn  E1
0.0000  Wasteland   Autumn  E2
0.0000  Wasteland   Autumn  E3
0.0000  Wasteland   Autumn  E4
0.0068  Wasteland   Spring  E1
0.0000  Wasteland   Spring  E2
0.0000  Wasteland   Spring  E3
0.0068  Wasteland   Spring  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 normal distribution: Then I performed a glmm of this dataset in R: GlM_habitats <- glmer(Inc. ~ Habitat + (1|Season)+(1|Site), data = Incidence) summary(GlM_habitats) Linear mixed model fit by REML ['lmerMod'] Formula: Incidence ~ 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 Correlation of Fixed Effects: (Intr) HbttEd HbttOk HabitatEdge -0.698 HabitatOkwd -0.701 0.576 HabttWstlnd -0.701 0.576 0.588  I tried to extract the variance of fixed effect but It only allow me extract the variance of the random effec. vc <- VarCorr(GlM_habitats) print(vc,comp=c("Variance","Std.Dev."),digits=2) Groups Name Variance Std.Dev. Site (Intercept) 0.00313 0.056 Season (Intercept) 0.00057 0.024 Residual 0.00082 0.029  How can I extract the variance of the fixed effect in glmm output? Thank in advance. ## 1 Answer I've used the vcovto extract the variance-covariance matrix of fixed effects. The variance is on the diagonal so converting it to a base matrix and then apply diag to extract the variances. After that one has to use sqrt to get the standard errors. Attached a working example: library(lme4) #> Lade nötiges Paket: Matrix # Construct dataframe: 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") ) GlM_habitats <- lme4::glmer(Inc. ~ Habitat + (1|Season)+(1|Site), data = Incidence) #> Warning in lme4::glmer(Inc. ~ Habitat + (1 | Season) + (1 | Site), data = #> Incidence): calling glmer() with family=gaussian (identity link) as a shortcut #> to lmer() is deprecated; please call lmer() directly summary(GlM_habitats) #> 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 #> #> Correlation of Fixed Effects: #> (Intr) HbttEd HbttOk #> HabitatEdge -0.698 #> HabitatOkwd -0.701 0.576 #> HabttWstlnd -0.701 0.576 0.588 # Variance of random effects: vc <- lme4::VarCorr(GlM_habitats) print(vc,comp=c("Variance","Std.Dev."),digits=2) #> Groups Name Variance Std.Dev. #> Site (Intercept) 0.00313 0.056 #> Season (Intercept) 0.00057 0.024 #> Residual 0.00082 0.029 # Variance-Covariance Matrix of fixed effects: vc_fixed <- as.matrix(vcov(GlM_habitats)) # Variance of fixed effects: var_fixed <- diag(vc_fixed); var_fixed #> (Intercept) HabitatEdge HabitatOakwood HabitatWasteland #> 0.001301387 0.002002356 0.002150297 0.002150297 # Standard errors of fixed effects: se_fixed <- sqrt(var_fixed); se_fixed #> (Intercept) HabitatEdge HabitatOakwood HabitatWasteland #> 0.03607474 0.04474769 0.04637129 0.04637129  Created on 2020-07-06 by the reprex package (v0.3.0) • Great Tim-TU thank you. It works properly. One more question. I must include the log family in my GLM model, It does not allow me put the family type: GlM_habitats = lme4::glmer(Incidence ~ Habitat +(1|Season)+(1|Site), data = Incidence, family = Gamma (link = "inverse")). I got this error Error in lmer(SeedPredation ~ AnimalGroup * Microhabitat * Site + (1 | : unused argument (family = "poisson"). How can I solve this? Thank you Jul 7 '20 at 18:23 • If I run: GlM_habitats = lme4::glmer(Incidence ~ Habitat +(1|Season)+(1|Site), data = Incidence, family = Gamma(link = "inverse")) i get the following error: Error in eval(family$initialize, rho) : non-positive values not allowed for the 'Gamma' family This error occurs because some of the y values are zero. I think your error is thrown by another example? Jul 8 '20 at 9:37
• No, I run this GlM_habitats = lme4::glmer(Incidence ~ Habitat +(1|Season)+(1|Site), data = Incidence, family = Gamma(link = "inverse")) using the input that you wrote in your example and I got the same output, how can I get a glmm for a log family? Thank you in advance Jul 8 '20 at 10:56
• families and links are described in ?stats::family. But again, a 'log' is only defined for positive y values bigger than zero (that is what the error in the previous comment says). Jul 8 '20 at 11:12