6
$\begingroup$

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?

$\endgroup$

1 Answer 1

5
$\begingroup$

There are 2 main problems here:

  1. 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.

  2. 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.

  3. 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))

enter image description here

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.

$\endgroup$
6
  • $\begingroup$ Hi Robert, the code worked well! Thank you. Please can you explain me what is the meaning of the variance of fixed and random effects? has site effect a higher effect in the model than season (0.003>0.00057). How can I explain the percentaje of the variance for these data. Why crop habitat in fixed effect does not appear? Thank you in advance $\endgroup$ Jul 29, 2020 at 9:49
  • $\begingroup$ No problem, I don't understand this question. The random effects have a variance. You can see that in the output. The fixed effects are fixed. They don't have a variance. They have a standard error. Crop habitat is there, in the fixed effects. What is your research question ? $\endgroup$ Jul 29, 2020 at 10:35
  • $\begingroup$ Is crop habitat the Intercept? I would like to know what percentage of the variance the different factor explain. For that I supposted that I must take the total variance and make a correspondence with the variance effects. According to your statement that the fixed effect, this fuction can not extracted the variance of the fixed effects: vc_fixed <- as.matrix(vcov(m0)) var_fixed <- diag(vc_fixed); var_fixed? These codes are explained in my question. Thank you Robert $\endgroup$ Jul 29, 2020 at 10:41
  • $\begingroup$ Ahh I see what you mean. Yes, you can do that. Another approach is to fit a variance components model: lmer(Inc.~ (1|Habitat) + (1|Season) + (1|Site) and this will directly tell you the variances, although with only 3 seasons and 4 habitats this may not be reliable. And yes, Crop is the reference level for Habitat and is included in the intercept. $\endgroup$ Jul 29, 2020 at 10:48
  • $\begingroup$ No, I mean, what is the meaningful that site have a variance of 0.0031294. Is it the 0.3% of the total variance? Sorry for my naive questions. Thank you $\endgroup$ Jul 29, 2020 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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