categories with count data, random and fixed effect in nested model

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.

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)
summary(model)

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.
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.
• You're welcome :) health_status and year will both be fixed effects (covariates) – Robert Long Feb 19 at 19:38
• 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. – Robert Long Feb 21 at 16:01