Manli Yan
2009-Apr-11 05:25 UTC
[R] who happenly read these two paper Mohsen Pourahmadi (biometrika1999, 2000)
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 11,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~~
[[alternative HTML version deleted]]
Manli Yan
2009-Apr-11 22:16 UTC
[R] who happenly read these two paper Mohsen Pourahmadi (biometrika1999, 2000)
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~~
[[alternative HTML version deleted]]
Manli Yan
2009-Apr-12 22:14 UTC
[R] who happenly read these two paper Mohsen Pourahmadi (biometrika1999, 2000)
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>