Dear friends,
As you are probably all aware, ASREML uses the Schall method for GLMMs
In his paper, Schall supplies some cell irradiation data. I enclose
that data and an ASREML job which runs several models on it.
The 27 units represent 9 occasions, 3 dishes per occasion and 400
cells per dish.
Schall fits 3 models which in ASREML notation are
surv ~ mu
surv ~ mu !r Occ
surv ~ mu !r Occ units
ASREML gives the same variance components as Schall gives for Occ and units
in these models.
Schall also reports a residual variance which I cannot reproduce,
even in the simplest case.
The residuals returned by ASREML are (O-E)/d
where O is the observation, E is the fitted expected value and
d is the derivative (if E=np, d=np(1-p) )
and the value beside the residual in the .sln file is (in this case)
logit(p) So it is quite simple to calculate residuals on the
observed scale.
If anyone is interested and has worked out exactly what Schall is doing,
please let me know.
Arthur
<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
Arthur Gilmour PhD mailto:Arthur.Gilmour@agric.nsw.gov.au
Senior Research Scientist (Biometrics) fax: <61> 2 6391 3899
NSW Agriculture <61> 2 6391 3922
Orange Agricultural Institute telephone work: <61> 2 6391 3815
Forest Rd, ORANGE, 2800, AUSTRALIA home: <61> 2 6362 0046
ASREML is still free by anonymous ftp from pub/aar on ftp.res.bbsrc.ac.uk
or point your web browser at ftp://ftp.res.bbsrc.ac.uk/pub/aar/
To join the asreml discussion list, send the message
subscribe
mailto:asreml-request@chiswick.anprod.CSIRO.au
To send messages to the list, mailto:asreml@chiswick.anprod.CSIRO.au
Asreml list archive: http://www.chiswick.anprod.csiro.au/lists/asreml
<> <> <> <> <> <> <>
"Christ Jesus came into the world to save sinners" I Timothy 1:15.
<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
Schall cell irradiation data
occasion 9
survived !/V3
total
schall.dat !DOPART $1
!PART 1
! surv !bin !total total ~ mu !r occ
! disp cal from residuals is 32.84/26 = 1.26 [Schall gives 1.810]
surv !bin !total total ~ mu !r occ
!PART 2
! surv !bin !total total !disp ~ mu !r occ
! dispersion calc as 1.819
surv !bin !total total !disp ~ mu !r occ
!PART 3
! surv !bin !total total ~ mu !r occ units
! dispersion (see below) is 9.3882 / 26 = .36 [Schall gives 0.937]
surv !bin !total total ~ mu !r occ units
!PART 4
! surv !bin !total total !disp ~ mu !r occ units
! fitted this way is overparameterised and disp does to zero
surv !bin !total total !disp ~ mu !r occ units
!PART 5
! surv !bin !total total ~ mu
! gives dispersion calculated from residuals of 18.96 (see below)
! Schall gives 18.09
surv !bin !total total ~ mu
!PART 6
! surv !bin !total total !disp ~ mu
! gives dispersion of 18.96 = 493/26
surv !bin !total total !disp ~ mu
!PART
0
pchi_function(file='schall.sln'){
# ASREML residuals are (O-E)/npq; predicted value is logit(p)
# Want (O-E)^2/npq
resdf_read.table(file)
res_resdf$V3[resdf$V1=='Residual']
upv_resdf$V4[resdf$V1=='Residual']
pv _ exp(upv)/(1+exp(upv))
pchi_400*sum(pv*res^2*(1-pv))
}
S> pchi('schall5.sln') #Schall gets 470.34 [22.6 less]
[1] 492.97
S> pchi('schall3.sln') #Schall gets 24.36 [15 more]
[1] 9.3882
S> pchi('schall2.sln')
[1] 33.174
S> pchi('schall1.sln') #Scall gets 47.06 [14.2 more]
[1] 32.836
Models 5 1 and 3 are the one Schall fits and ASREML agrees
with respect to the extra variance components.