Adding to Wim's comments:
The first REML is failing primarily because the AI algorithm does not
deal well
with negative variance components. As Wim showed, the Fisher scoring
algorithm
deals much better with these cases. The AI method does get to the right
result
if you let it run to 50 iterations but, because it has explored areas
outside the
valid parameter space (gamma < 0.11), still thinks it has not achieved
convergence. Using very small initial values, eg. 0.01, gives better
behaviour
in this case. It is usually helpful in these cases to print out
monitoring.
The discrepancy with the ANOVA analysis happens because the two methods
use different covariance structures. If Z is the design matrix for
factor groups,
then REML uses variance structure:
V = s2*( gZZ' + I )
whereas ANOVA uses
V = s2*( gZinv(Z'Z)Z + I )
These two forms are equivalent for groups with equal numbers, but not in
unbalanced cases. I would guess Nick generated the data using the REML
variance structure, as this fits better. If you put the ANOVA variance
structure into the model, you get out the ANOVA residual variance, see
below.
Hope this helps,
Sue Welham
Biomathematics and Bioinformatics
Rothamsted Research
Harpenden UK AL5 2JQ
E-mail: [log in to unmask]
Tel: +44 (0)1582 763133 ext 2278
Fax: +44 (0)1582 467116
Rothamsted Research is a company limited by guarantee, registered in
England under the registration number 2393175 and a not for profit
charity number 802038.
factor [levels = 10] group
read [setnvalues = yes] group, y
1 -1.03554
1 1.07015
1 1.05308
1 -0.442382
1 -1.64084
1 1.51995
2 -0.233871
2 -1.40783
2 -0.0713538
2 2.23295
2 0.0634799
2 0.698427
2 0.739394
2 0.493626
2 -2.44862
3 -0.512397
3 0.328796
3 0.526371
3 1.25548
3 -0.177776
3 0.325434
3 1.12841
3 -0.296016
4 3.1801
5 -1.70055
5 0.694679
5 0.15136
6 0.170231
6 -1.48726
6 -1.17832
6 -1.27697
6 1.99981
6 0.303725
6 0.224378
6 -0.422661
6 0.290695
7 -0.557846
7 0.235275
7 -0.0452684
7 0.602743
8 -1.20067
8 0.464304
8 0.74811
9 1.66716
9 0.957378
9 0.569089
10 -0.480984
10 0.204959
10 -0.0565612
10 0.0678034
10 -1.16203
10 -0.316222
10 1.03135 :
treatments group
anova [p=#,stratum] y
aplot fitted,normal,halfnormal,histogram
vcomponents random = group; init=0.01
reml [print = model, components,mon; max=50] y
vcomponents random = group
reml [print = model, components,mon; meth=fish] y
vkeep group; des=Dg
symm Vz
calc Vz = Dg*+inv(t(Dg)*+Dg)*+t(Dg)
prin Vz; dec=3; fie=8
flrv [p=r] t(Dg)*+Dg; lrv=lZ
calc Dgh = lZ[1]*+sqrt(1/lZ[2])*+t(lZ[1])
calc Zg = Dg*+Dgh
" fit ANOVA variance model "
vcomponents random = Zg
reml [print = model, components,mon] y
" look at deviance across values of variance ratio "
vari gam; val=!(-0.1,-0.09...-0.01,0.01,0.02...0.1)
vari [nval=gam] dev,S2
for [index=i] g=#gam
vcomponents random = group; init=g,1; con=fix,pos
reml [print = *] y
vkeep [dev=dv; sigma2=s2]
calc dev$[i]=dv
calc S2$[i]=s2
endfor
dgraph dev; gam
dgraph S2; gam
prin gam,S2,dev
|