3
$\begingroup$

I am using gls in nlme. My response variable is spatial so I am using gls with correlation structure. I am determining which structure to use based on Zuur 2009, comparing AIC scores of models with all relevent fixed effects and only differing in the correlation structures. I am doing this in REML as at this stage all models have identical fixed effects. Then I am using Likelihood ratio tests in ML to select best model (in terms of fixed effects), again following Zuur 2009. Once the best model is determined, I am re-running it in REML for my unbiased parameter estimates. My problem is that I am attempting to test several hypotheses so ultimately need to compare the best model for several different hypotheses (eg. a prey model compared to a habitat model compared to a topography model). I understand that I cannot compare these using AIC in REML because the fixed effects are different even though the random effects are the same.

Can I compare the AIC scores of the best models as they are determined using ML (ie BEFORE I re-run them in REML for unbiased parameter estimates)? I have read a fair amount of the literature (including some related topics here), but I have not been able to determine once and for all if this is suitable as most comments I have seen are concerned with finding the best model and not comparing different models.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ This is a very good question. But as part of your question you say you are using model selection that is unblinded to Y. This will result in biased parameter estimates and badly distorted standard errors, confidence intervals, and P-values. $\endgroup$
    – Frank Harrell
    May 3, 2013 at 12:45
  • $\begingroup$ Hi Frank. Sorry for my ignorance, but what do you mean by "unblinded to Y"? The model selection I employ to get best model (eg prey) starts with the full set of prey-related variables and removes non-significant in a backward stepwise fashion (anova() in ML). Then I compare with best habitat (for example) model, obtained with same procedure. In final model comparison, I want to use AIC (from ML - see above question) to compare all (ie prey alone, habitat alone and prey+habitat). I obviously dont want biased parameter estimates etc. so want to be clear on what you have pointed out. Thanks. $\endgroup$
    – A.Kittle
    May 5, 2013 at 15:51
  • 1
    $\begingroup$ What gave you the idea to remove non-significant variables? The process you described is fully unblinded to $Y$ [unlike data reduction methods such as principal components, variable clustering, and redundancy analysis]. It will bias all aspects of statistical inference and will in no way make predictions more accurate. Think of it this way: if you bias the estimate of $\sigma^2$ low, everything that comes out of the model is anti-conservative. You are effectively using $\hat{\sigma}^2 = \sum(error^2)$ [from a model dredged from several models] / $(n-p-1)$ where $p$ hides phantom d.f. $\endgroup$
    – Frank Harrell
    May 5, 2013 at 16:23
  • $\begingroup$ This doesn't sound good. I originally wanted to compare models with full suite of relevent variables (ie. model w/ all relevent prey variables to model w/all relevent habitat variables), but was advised it was better to undertake "model competition within model competition" whereby I determine "best" prey model and "best" habitat model and compare those. $\endgroup$
    – A.Kittle
    May 5, 2013 at 16:50
  • $\begingroup$ That kind of data dredging creates multiple problems, and the data are incapable of reliably telling you which model is "correct". If you had only two pre-specified models to choose between, things would be better. $\endgroup$
    – Frank Harrell
    May 5, 2013 at 17:48

1 Answer 1

Reset to default
2
$\begingroup$

Can I compare the AIC scores of the best models as they are determined using ML (ie BEFORE I re-run them in REML for unbiased parameter estimates)? Yes, but be careful that comparison of non-nested models has to be based on the information criteria AIC and BIC. Check visually at least your sample variogram estimates also.

I can't seem to locate the perfect reference point but Pinheiro & Bates in Mixed Effects Models in S and S-Plus, sub-chapter 5.4 Fitting Extended Linear Models with gls does what you describe. Importantly the anova() method can be used to compare gls() and lme() objects. They stress multiple times though that the choice between a gls and and lme should be take other factors in consideration besides the information criteria.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.