Dear mixed-modellers;

I have couple of questions regarding mixed-model specification in R (nlme) and (lme4).

My data have hierarchical structure.  They come from standard forest inventory plots (i.e. all trees are measured in a fixed area circular plot of size 400 m2) and some plots are re-measured more than two cycles i.e. more than one measurement periods (e.g. 1999-2004, 2000-2005, 2001-2006, 2002-2007, 2003-2008, 2004-2009).  Most of the plots were measured only one measurement period, some plots that were established in 1999 were measured in 2004 and 2009 (three cycles or two measurement periods). I have nested measurement periods (MP) within Plot and trees are nested in the Plot.

My approach:
a) First, I was debating whether should I use nonlinear or log-linearized model.  I chose to use log-linearized model as one of the predictor variable is dummy variable so that I can fit and test its (dummy variable) effects in the model. 
b) Despite its slow convergence, I chose nlme package as it allows me to specify covariance structure (AR(1)) and specify weighting in order to stabilize non-constant variance if there exists any.

I have specified mixed-effects model as follow;

m.0 <- lme(dds.log ~ dbh.log+dbh+BAL+BA+SITES,
            random=~1|PLot/MP, # measured periods (MP) are nested within Plots
            data=tree, method="ML")
m.AR1 <- update(m.0, correlation = corAR1()) 
m.0d <-  update(m.AR1, weights = varPower(form = ~dbh))
m.0w <-  update(m.AR1, weights = varPower())
m.0w.REML <-  update(m.0w, method="REML")
m.0w.NOSE <-  m.0w fitted without SITES covariate

where, response variable is log linearized individual tree increment, predictors are dbh and BAL (individual tree variables), whereas BA is plot level variable varies by MP and Plot, and SITES is categorical Plot level variable (assumed to be constant for next few decades).  MP is measured period e.g. 1999-2004 or 2004-2009 and Plot is plot ID. 

My questions are:
a) residual variance after weighting: when I used weighted model (m.0w) - residual variance modelled as a power function of mean fitted value, then residual variance increases as expected (weight is power of negative delta), but variance due to plot and MP remained more or less constant.
    (i) Since this is log-linearized model, I would like to apply bias correction factor as dds = exp(predicted_dds.log+(s^2(Plot)+S^2(MP)+S^2(res))/2).   Any effects from weighting to my log-linearized bias correction factor?  or S^2(residual variance changes as expected when applying weight), should I be using the residual variance after weighting to correct the log-linearized bias? Any suggestions?

> anova(m.0, m.AR1, m.0d, m.0w)
      Model df      AIC      BIC    logLik   Test  L.Ratio p-value
m.0       1 31 48602.89 48857.54 -24270.44                        
m.AR1     2 32 48600.60 48863.47 -24268.30 1 vs 2  4.28633  0.0384
m.0d      3 33 48569.17 48840.25 -24251.58 2 vs 3 33.43440  <.0001
m.0w      4 33 46205.77 46476.85 -23069.88              

The data fits m.0w model better than other null models.  Finally, I chose m.0w and used REML to parameterize and estimate variance components.

>print(m.0w) # final model with varPower weighting as a function of mean fitted value
Random effects:
 Formula: ~1 | Plot
        (Intercept)
StdDev:   0.3974489

 Formula: ~1 | MP %in% Plot
        (Intercept) Residual
StdDev:   0.2081098 1.283800

Correlation Structure: AR(1)
 Formula: ~1 | Plot/MP 
 Parameter estimate(s):
        Phi 
0.007990434 
Variance function:
 Structure: Power of variance covariate
 Formula: ~fitted(.) 
 Parameter estimates:
     power 
-0.8423614 
Number of Observations: 27298
Number of Groups: 
        Plot MP %in% Plot 
        1182         1210 


>print(m.AR1) # model without weighting
Random effects:
 Formula: ~1 | Plot
        (Intercept)
StdDev:   0.3928492

 Formula: ~1 | MP %in% Plot
        (Intercept)  Residual
StdDev:   0.2008393 0.5616311

Correlation Structure: AR(1)
 Formula: ~1 | Plot/MP 
 Parameter estimate(s):
       Phi 
0.01321357 
Number of Observations: 27298
Number of Groups: 
        Plot MP %in% Plot 
        1182         1210 


   (ii) I am trying to calculate RMSE to estimate what proportion of variance in basal area increment is explained by Plot and/or MP?, should I be using variance estimated after weighting?

b) I would like to calculate what proportion of variation in basal area increment is explained by SITES?  I fit the model without SITES (m.0w.NOSE) and Likelihood ratio test suggests that SITES is highly significant predictor in the model. I still would like to quantify how much variation in basal area increment is explained by the SITES as a predictor variable (for this I calculated RMSE and compare model with and without SITES).  Mostly the model fitted with SITES reduces Plot level variance (as expected since this is plot level variable).  If I calculate the percent of RMSE reduced by SITES before and after weighting, they are different because only residual variance is increased due to weighting (negative delta value). Any suggestions?

> anova(m.0w, m0w.NOSE)
            Model df      AIC      BIC    logLik   Test  L.Ratio p-value
m.0w.NOSE     1 10 46532.47 46614.61 -23256.23                        
m.0w            2 33 46205.77 46476.85 -23069.88 1 vs 2 372.7003  <.0001

> print(m.0w.NOSE)
Random effects:
 Formula: ~1 | Plot
        (Intercept)
StdDev:    0.502633

 Formula: ~1 | MP %in% Plot
        (Intercept) Residual
StdDev:   0.2096941 1.292326

Correlation Structure: AR(1)
 Formula: ~1 | Plot/MP 
 Parameter estimate(s):
        Phi 
0.007753674 
Variance function:
 Structure: Power of variance covariate
 Formula: ~fitted(.) 
 Parameter estimates:
     power 
-0.8486992 
Number of Observations: 27298
Number of Groups: 
        Plot MP %in% Plot 
        1182         1210 


Sorry for a long email. Any pointer you could provide me on these questons would be very very helpful. Thank you and look forward to hearing from you!

Regards,
Bharat
	[[alternative HTML version deleted]]


