# Re: How to compute AIC from ASReml fitted model on R?

From: <brian.cullis_at_DPI.NSW.GOV.AU>
Date: Wed, 22 Jul 2009 12:48:00 +1000

Dear saidou
You raise a number of issues which I must address.

The logic I was using for the suggestion of how to compute the AIC for
model selection in lmms was based solely on the assumption that the fixed
model is constant. It is well known that REML likelihoods cannot be
compared when the fixed effects change. Therefore if you are interested in
fixed effects then the approach I documented in the previous email no
longer applies. This is the danger of email lists in that the full story

There is doubt in the literature about the RIC (residual information
criteria) proposed by Shi and Tsai (2002, JRSSB) - see recent article by
Chenlei Leng who suggest that it is not appropriate but present a
corrected RIC

I have not investigated these at all.

Further you are correct that we do not have an update() or a step() method
for ASReml-R, but certainly the former is on the agenda, modulo keeping
the fixed model the same.

If it helps I have had a consulting job which req'd model selection, and I
used a forward selection approach with a FDR approach for selection. The R
code is given below
This is purely indicative but may help you to create some code for
yourself.
There is also the issue of small sample adjustments for Wald statistics -
which I didnt use here but for most situations I would recommend it (see
Kenward and Roger, Biometrics paper and a more recetn article I believe??)

HTH

##############################
# scan individual variables
###########################

xlist <- names(GP)[c(8,11,14,18,19,26,27,28,29,31,32,33,41,46,57,59)]
ww <- c()
for(x in xlist){
ff <-formula(paste('y~1',x,sep='+'))
GP.asr <- asreml(fixed=ff,na.method.X='include',data=GP)
ww <- c(ww,wald(GP.asr)\$Pr[2])}

first.in.df <- data.frame(Variate=xlist,Wald.p=ww)

#################################
# choose tlend
# and then scan variables again
#################################

xlist <- names(GP)[c(8,11,14,18,19,26,27,28,29,31,32,33,41,57,59)]
ww <- c()
for(x in xlist){
ff <-formula(paste('y~1+tlend',x,sep='+'))
GP.asr <- asreml(fixed=ff,na.method.X='include',data=GP)
ww <- c(ww,wald(GP.asr)\$Pr[3])}

second.in.df <- data.frame(Variate=xlist,Wald.p=ww)

######################
# etc etc
#
########################

warm regards

Brian Cullis
Senior Principal Research Scientist
NSW Department of Primary Industries
Wagga Wagga Agricultural Institute

Visiting Professorial Fellow
School of Mathematics and Applied Statistics
Faculty of Informatics
University of Wollongong

Professor,
Faculty of Agriculture, Food & Natural Resources
The University of Sydney