infile logs days ch o3 year tree using "C:\My Documents\655LDA\DataLDA\sitka.raw "
rename tree id
save "C:\My Documents\655LDA\DataLDA\sitka.dta", replace

****The observed mean responses in each of the 4 chambers; for 1988 and 1989
egen mean1=mean(logs) if ch==1&days<300, by(days) 
egen mean2=mean(logs) if ch==2&days<300, by(days)
egen mean3=mean(logs) if ch==3&days<300, by(days)
egen mean4=mean(logs) if ch==4&days<300, by(days)
egen mean5=mean(logs) if ch==1&days>300, by(days) 
egen mean6=mean(logs) if ch==2&days>300, by(days)
egen mean7=mean(logs) if ch==3&days>300, by(days)
egen mean8=mean(logs) if ch==4&days>300, by(days)

label var mean1 "Ch1, 98"
label var mean2 "Ch2, 98"
label var mean3 "Ch3, 98"
label var mean4 "Ch4, 98"
label var mean5 "Ch1, 99"
label var mean6 "Ch2, 99"
label var mean7 "Ch3, 99"
label var mean8 "Ch4, 99"
 
sort days
graph7 mean* days, c(llllllll) xlab ylab s(iiiiiiii)

****The observed mean responses in each of the two treatment groups
egen meanc1=mean(logs) if o3==0&days<300, by(days) 
egen meanc2=mean(logs) if o3==0&days>300, by(days)
egen meant1=mean(logs) if o3==1&days<300, by(days)
egen meant2=mean(logs) if o3==1&days>300, by(days)
label var meanc1 "control-1988"
label var meanc2 "control-1989"
 label var meant1 "O3-treated-1988"
label var meant2 "O3-treated-1989"
sort days
graph7 meanc* meant* days, c(llll) xlab ylab s(iiii)

*** Compute the REML estimates of the covariance matrix
tsset id days

**  Plotting logsize over time, across chambers ;

** two scatterplots;
graph7 logs days if o3==0, ti("scatterplot of logsize vs. days, in normal environ") 
graph7 logs days if o3==1, ti("scatterplot of logsize vs. days, in ozone environ") 

** two  box plots;
sort days
graph7 logs if o3==0, box by(days) ti("Boxplot of logsize vs time, in normal environ") 
graph7 logs if o3==1, box by(days) ti("Boxplot of logsize vs time, in normal environ") 

**   Mean trend plot***;
xtgraph logs if (days<300) , group(o3) ti("Mean logsize vs days, 1988")  bar(se)  offset(1.5)  
xtgraph logs if (days>300) , group(o3) ti("Mean logsize vs days, 1989")  bar(se)  offset(1.5)  

****Saturated model for the mean response, ignore chamber effect
keep if days<300
anova logs days o3
predict logsres, resid

***The corresponding fitted means
predict yhat 

keep id logsres days
reshape group days 152 174 201 227 258
reshape var logsres
reshape  cons id
reshape wide
***Get the matrix in (4.6.6) pg 72 book
matrix accum Cov=logsres152 logsres174 logsres201 logsres227 logsres258, deviations noconstant
scalar adj=1/(_result(1)-4)
**get teh REML estimator of V_0 (4.6.4) pg. 71
matrix Cov=adj*Cov

 *REML estimate of covariance for 1988, P72
 matrix list Cov

**NOw get the GLS beta(hatW)
use "Z:\LDA\DataLDA\sitka.dta", clear 
keep if days<300
tab days, gen(t)
gen eta=1-o3
*covariate of time
gen gamma=days/100*eta

*We fit GEE marginal model with robust estimates of standard error
 xtgee logs t1 t2 t3 t4 t5 eta gamma, noconst i(id) corr(exc) robust
*Incorrect standard error using OLS
 reg logs t1 t2 t3 t4 t5 eta gamma, noconst
*Now create design matrix for 1988 data X
mkmat t1 t2 t3 t4 t5 eta gamma, matrix(X)
*Now create response matrix of Y for 1988 data
mkmat logs, matrix(Y)
*OLS estimate of beta  (4.3.4) pg. 60
 matrix beta=syminv(X'*X)*X'*Y
 matrix list beta
*Calculate the robust estimate of standard error of beta (4.6.1)
 ***This is W^(-1)
matrix diag=I(79) 
 matrix V=diag#Cov
***Get the estimated variance matrix (4.6.2) pg. 70
 matrix RW=syminv(X'*X)*X'*V*X*syminv(X'*X)
*Robust estimate of covariance matrix for the coefficients
 matrix list RW

* Calculate the robust standard error
matrix var=vecdiag(RW)
matrix var=var'
svmat var, name(varbeta)
keep varbeta1
drop in 8/395
gen se=sqrt(varbeta1)
 svmat beta, name(beta)
 mkmat beta1 se, matrix(Coef)
 matrix Coef=Coef'
 matrix list Coef
 matrix colnames Coef= beta1 beta2 beta3 beta4 beta5 tao gamma
*Table4.3 for 1988 P75
 matrix list Coef


***********For 1989 data, we perform the similar analysis as above
 clear
 set matsize 800
 set memory 20000 m
 use  "C:\My Documents\655LDA\DataLDA\sitka.dta", clear 
 keep if days>300
anova logs days o3
predict logsres, resid
***The corresponding fitted means
predict yhat 
keep id logsres days
reshape group days 469 496 528 556 579 613 639 674
reshape var logsres
reshape  cons id
reshape wide
matrix accum Cov=logsres469 logsres496 logsres528 logsres556 logsres579 logsres613 logsres639 logsres674, deviations noconstant
scalar adj=1/(_result(1)-4)
matrix Cov=adj*Cov

 *REML estimate of covariance for 1989, P72
 matrix list Cov
**NOw get the GLS beta(hatW)
use "C:\My Documents\655LDA\DataLDA\sitka.dta", clear 
keep if days>300
tab days, gen(t)

gen eta=1-ozone
***gen gamma=days/100*eta

*gen tao=1-o3

*We fit GEE marginal model with robust estimate of standard error
 xtgee logs t1 t2 t3 t4 t5 t6 t7 t8 eta, noconst i(id) corr(exc) robust
*Incorrect standard error using OLS
 reg logs t1 t2 t3 t4 t5 t6 t7 t8 eta, noconst
*Now create design matrix of X for 1989 data X
mkmat t1 t2 t3 t4 t5 t6 t7 t8 eta, matrix(X)
*Now create response matrix of Y for 1989 data
mkmat logs, matrix(Y)
*OLS estimate of beta
 matrix beta=syminv(X'*X)*X'*Y
 matrix list beta
*Calculate the robust estimate of standard error of beta (4.6.1)
 matrix diag=I(79) 
 matrix V=diag#Cov
 matrix RW=syminv(X'*X)*X'*V*X*syminv(X'*X)
*Robust estimate of covariance matrix for the coefficients
 matrix list RW

* Calculate the robust standard error
matrix var=vecdiag(RW)
matrix var=var'
svmat var, name(varbeta)
keep varbeta1
drop in 10/632
gen se=sqrt(varbeta1)
 svmat beta, name(beta)
 mkmat beta1 se, matrix(Coef)
 matrix Coef=Coef'
 matrix list Coef
 matrix colnames Coef= beta1 beta2 beta3 beta4 beta5 t6 t7 t8 tao 
*Table4.3 for 1989 P75
 matrix list Coef