I am confused how to analyze this dataset in R: There are 2 islands, 2 sites were selected (nonrandom) per each island, 6 transects randomly were located on each site, then the number of corals was counted in three categories: bleached, partial-bleached and health. The data are collected for 2 years.

enter image description here

The transect (random effect) is nested within sites, and sites (fixed effect) within island (fixed effect), but I have counted on categorical data. I am confused which analyses should I use to compare the spatial and temporal variation in each category of health condition. I think generalized linear mixed model is perfect but multivariate abundance model may be the best. Which test (analysis) in which package is the best fit to my data? How should I deal with the categories, is it fixed effect factor?

Any suggestion would be appreciated.

I ran the model as:

model <- glmer(No.coral ~ island + health_status + (1|island:site:transect), data = data, family=poisson)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
     Family: poisson  ( log )
    Formula: No.coral ~ island + health_status + (1 | island:site:transect)
       Data: data
         AIC      BIC   logLik deviance df.resid 
     33065.5  33100.5 -16527.8  33055.5     8059 
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
     -1.733  -0.961  -0.567  -0.284 117.509 

Random effects:
 Groups               Name        Variance Std.Dev.
 island:site:transect (Intercept) 0.1081   0.3289  
Number of obs: 8064, groups:  island:site:transect, 32

Fixed effects:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -2.99554    0.10110  -29.63   <2e-16 ***
islandL         1.38643    0.12018   11.54   <2e-16 ***
health_statusM  2.29915    0.05599   41.06   <2e-16 ***
health_statusS  2.31649    0.05595   41.40   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) islndL hlth_M
islandL     -0.618              
helth_sttsM -0.503  0.000       
helth_sttsS -0.503  0.000  0.909

My questions are: Is my model written correctly? why the model output is dismissed 'health_status_H'? Also, I don't know how to show the variation (increasing/decreasing) of each health status during the 4 years of observation (not only show the significance of variations as above)? (before, I told that the samplings were conducted during 2 years, by mistake). I also need comparisons ONLY among years at each island and each site. Many thanks.


I don't know what a multivariate abundance model is, but as you thought, you should be able to use a generalised linear mixed model here. 6 is a rather low number for random intercepts, so I would also fit a non-mixed model, and fit Transect as a fixed effect, and compare the inferences from both. Hopefully they will be consistent.

Since you have count data, a good place to start is with a poisson model. If it turns out that you have excess zeros, then you can use a zero inflated model poisson model, and if you over or under dispersion, you can use a negative binomial instead of a poisson.

If you are using R then I would suggest starting with the mixed_model function in the GLMMadaptive package or the glmer function in the lme4 package.

  • $\begingroup$ Dear Robert Long Thanks for your answer. So, in the GLMM how should I introduce the column of 'health_status' to the model? Is it fixed or random factor? Also, the year, here, is covariate, right? $\endgroup$
    – fgh
    Feb 19 at 19:18
  • 1
    $\begingroup$ You're welcome :) health_status and year will both be fixed effects (covariates) $\endgroup$ Feb 19 at 19:38
  • $\begingroup$ Due to character limitation here, I add my comment as an answer. I would appreciate any help you can provide. $\endgroup$
    – fgh
    Feb 21 at 13:53
  • $\begingroup$ Please don't post an answer to add more info. You can edit your question by using the Edit facility to add more info. Anyway, (1|island:site:transect) does not make sense. It should be just (1 | transect). health_status_H is included in the intercept. For the time effect you need to include year as a fixed effect, as I said in my last comment. $\endgroup$ Feb 21 at 16:01
  • $\begingroup$ Finally, I fit the model as: 'modelF8 <- glmmTMB(No.coral ~ (1|transect) + year +island + site + health_status + genus, data = data, family = nbinom1)' . Anova shows the significant differences between sites, however, I need to show the difference of sites in term of islands, i.e. nested design (sites within islands) has to be included. I don't know how to include nested design in glmmtmb. Also, I think there is some interaction such as 'health_status*genus'. But including such interaction leads to error of non-positive-definite Hessian matrix. $\endgroup$
    – fgh
    Mar 11 at 9:54

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.