[R] who happenly read these two paper Mohsen Pourahmadi (biometrika1999, 2000)

Manli Yan manliyanrhelp at gmail.com
Mon Apr 13 00:14:18 CEST 2009


http://biomet.oxfordjournals.org/cgi/reprint/86/3/677   biometrika1999
http://biomet.oxfordjournals.org/cgi/reprint/94/4/1006  biometrika2000

Hi All:
  I just want to try some luck.
  I am currenly working on my project,one part of my project is to
reanalysis the kenward cattle data by using the method in Mohsen's paper,but
I found I really can get the same or close output as he did,so,any one who
have happenly read his paper,got any idea how he got the LS estimates of
gamma.
   he assumed the autoregressive phi_tj satisfy the cubic function on lag of
time,where t is from 1 to 11,j from 1 to t-1,his model is:
  phi_tj=gamma1+gamma2*(t-j)+gamma3*(t-j)^2+gamma4*(t-j)^3, #biometrika
1999,page685
 at first moment ,I thougt my design matrix should be
  1  1   1  1
  1  2  4  8
  1  3  9  27.....
but I found this is wrong,actually I should(I think) use a polynomial design
matirx with level of 10,degree=3
here is my rocode~~
y1=c(1,0.9,0.98,1.06,0.83,1.00,0.41,0.93,1.01,0.86)
y2=c(0.05,0.16,0,0.26,0.15,0.61,0.33,0.31,0.33)
y3=c(-0.23,0,0.16,-0.03,0.22,-0.03,-0.17,-0.05)
y4=c(0.04,-0.21,-0.04,-0.26,-0.03,-0.04,-0.05)
y5=c(-0.02,-0.34,0.06,-0.22,-0.11,-0.31)
y6=c(0.20,0.01,0.01,-0.26,0.01)
y7=c(-0.06,-0.14,0.39,0.23)
y8=c(0.21,0.10,0.09)
y9=c(-0.24,-0.23)
y10=c(0.13)
y=c(y1,y2,y3,y4,y5,y6,y7,y8,y9,y10) ## autoregressive paramters ,table 1
biometrika 1999,page 685
om2=matrix(0,nrow=55,ncol=3)
om1=poly(c(1:10),degree=3)             ##polynomial  design matirx  with
level =11,cubic
 om2[1:10,]=t(matrix(rep(om1[1,],10),ncol=10))
 om2[11:19,]=t(matrix(rep(om1[2,],9),ncol=9))
 om2[20:27,]=t(matrix(rep(om1[3,],8),ncol=8))
 om2[28:34,]=t(matrix(rep(om1[4,],7),ncol=7))
 om2[35:40,]=t(matrix(rep(om1[5,],6),ncol=6))
 om2[41:45,]=t(matrix(rep(om1[6,],5),ncol=5))
 om2[46:49,]=t(matrix(rep(om1[7,],4),ncol=4))
 om2[50:52,]=t(matrix(rep(om1[8,],3),ncol=3))
 om2[53:54,]=t(matrix(rep(om1[9,],2),ncol=2))
 om2[55,]=t(matrix(rep(om1[10,],1),ncol=1))
## so om2 is my orthogonal design matirx
 ## for example ,for y1=c(1,0.9,0.98,1.06,0.83,1.00,0.41,0.93,1.01,0.86),
since t-j=1,so the design matrix should be the first row of om1,and repeat
for 10 times.

t1=om2[,1]
t2=om2[,2]
t3=om2[,3]
fit<-lm(y~(t1+t2+t3))
summary(fit)


### the estimate I got is
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.09186    0.02980   3.083   0.0033 **
t1          -0.54143    0.10970  -4.936 8.94e-06 ***
t2           0.48955    0.10417   4.700 2.01e-05 ***
t3           0.56025    0.08771   6.387 5.05e-08 **

which is total different from the author's(0.18,-1.71,1.64,-1.11)

I have been tried contact the author,but he only point out that the I might
need to pay close attention to design matrix,but did not say how,I really
run out of the idea how to construct another suitable design matrix ,since
basing the paper,this should be the correct design matrix,so ,any one has
any idea???

  GREAT THANKS~~
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Pourahmadi.pdf
Type: application/pdf
Size: 159889 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20090412/82be6b70/attachment-0004.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Pourahmadi2.pdf
Type: application/pdf
Size: 142766 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20090412/82be6b70/attachment-0005.pdf>


More information about the R-help mailing list