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.

$\endgroup$
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)

$\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$ Jul 7 '20 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$
    – Tim-TU
    Jul 8 '20 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$ Jul 8 '20 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$
    – Tim-TU
    Jul 8 '20 at 11:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.