How to test for significant differences between growth curves in repeated measures mixed-models?

Question

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.

Questions:

1. Do you think this approach is sound?

2. Do I even need to specify an autocorrelation function when using mixed effects (phi is usually quite high (>0.9))?

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