I am comparing growth of fish in different habitats and since length and age data is derived from repeated measurements (back-calculations) I am using a mixed model (R, nlme). I have data from ~50 fish per habitat.
My interest is a) parameters of the growth function(s) in the different habitats and b) whether those are significantly different.
From literature I adopted the following approach: Run a global model (Model 1, all data) with a single fixed effect for the parameters of the growth function and specifify the parameters as random effects on an Individual level. Also I use an AR1 autocorrelation structure. Then I run the model with fixed effects for each habitat (Model 2, could do it for each parameter but for now I test all of them at once) and compare the two models by anova.
Do you think this approach is sound?
Do I even need to specify an autocorrelation function when using mixed effects (phi is usually quite high (>0.9))?
If I select a model by information theoretic approaches, how to deal with the different levels of the model (e.g. I select a model, like Bertalanffy, Gompertz etc., based on the global model (all data) but comparing models for each of the different habitats actually reveals that different models are supported - and by quite a notable margin (tested this with nlme and gnls)? Am I overthinking this?
Thanks for your help!
The code is (fge=habitat and ID=individual):
#model 1: VBGF.base.corr <- nlme(la ~ VBGF(al, Linf, k, t0), data = growth, fixed = list(Linf + k + t0 ~ 1), random = Linf + k + t0 ~ 1, groups = ~ID, correlation = corAR1(form = ~1|ID), method = "ML", start = list(fixed=c(Linf=100,k=0.1,t0=-1.5)), control = nlmeControl(maxIter = 1024, pnlsMaxIter = 1024, msMaxIter = 1024, pnlsTol = 0.05, msTol= 0.05) ) #model 2: VBGF.fge.corr <- nlme(la ~ VBGF(al, Linf, k, t0), data = growth, fixed = list(Linf + k + t0 ~ fge-1), random = Linf + k + t0 ~ 1, groups = ~ID, correlation = corAR1(form = ~1|ID), method = "ML", start = list(fixed=c(Linf = c(101,100,105,117,114,96,107), k = c(0.067, 0.09,0.11,0.08,0.07,0.16,0.099), t0 = c(-0.99, -0.87, -0.54, -0.72, -0.95, -0.55, -0.5))), control = nlmeControl(maxIter = 1024, pnlsMaxIter = 1024, msMaxIter = 1024, pnlsTol = 0.05, msTol= 0.05) )