3
$\begingroup$

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:

enter image description here

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

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)

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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 $\endgroup$ Commented Jul 7, 2020 at 18:23
  • $\begingroup$ 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? $\endgroup$ Commented Jul 8, 2020 at 9:37
  • $\begingroup$ 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 $\endgroup$ Commented Jul 8, 2020 at 10:56
  • $\begingroup$ 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). $\endgroup$ Commented Jul 8, 2020 at 11:12

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.