**************************************************************
* thallGEE.do                                                *
**************************************************************
*                                                            *  
* PURPOSE:  analysis of seizure counts                       *
*                                                            *
* AUTHOR:  P. Heagerty                                       *
*                                                            *
* DATE:  05 May 2005                                         *  
**************************************************************

infile id age tx y0 y1 y2 y3 y4  using ThallWide.raw

label variable age "Age (years)"
label variable tx "treatment"

label define tlab 0 "placebo" 1 "progabide"
label values tx tlab

***
*** some exploratory data analysis
***
summarize  age
stem age

tab tx

sort tx

*** MEAN exploratory analysis

by tx: summarize  y0 y1 y2 y3 y4

*** CORRELATION exploratory analysis

by tx: corr y0 y1 y2 y3 y4


***
*** graphical EDA
***

*** seizure by group

*graph box y0, by(tx)
*graph export c:\courses\correlated\GEE\seizure-y0box.ps, as(eps) replace

*graph box y1 y2 y3 y4, by(tx)
*graph export c:\courses\correlated\GEE\seizure-postbox.ps, as(eps) replace

***** post versus pre

gen logY4tx = ln( y4+1 )
recode logY4tx min/max=. if tx==0

gen logY0tx = ln( y0+1 )
recode logY0tx min/max=. if tx==0

gen logY4control = ln( y4+1 )
recode logY4control min/max=. if tx==1

gen logY0control = ln( y0+1 )
recode logY0control min/max=. if tx==1

*graph twoway (scatter logY4tx logY0tx, mcolor(magenta) msymbol( D )) ///
*      (scatter logY4control logY0control, mcolor(green) msymbol( Oh ))
*graph export c:\courses\correlated\GEE\seizure-PrePost.ps, as(eps) replace

*** change by tx

gen diff = y4/2 - y0/8

*histogram diff, by(tx) width(1)
*graph export c:\courses\correlated\GEE\seizure-change.ps, as(eps) replace


*** LONGITUDINAL regression models

gen logY0 = ln( y0+1 )

save ThallWide, replace

reshape long y, i(id) j(week)

gen obsTime = 2*(week>0) + 8*(week==0)
gen logObsTime = log( obsTime )

*** create some time variables

gen wk1 = (week==1)
gen wk2 = (week==2)
gen wk3 = (week==3)
gen wk4 = (week==4)

gen weekXtx = week * tx
gen wk1Xtx = wk1 * tx
gen wk2Xtx = wk2 * tx
gen wk3Xtx = wk3 * tx
gen wk4Xtx = wk4 * tx

*** GEE with all times as outcome

xtgee y week tx weekXtx, offset(logObsTime) ///
  i(id) corr(unstructured) t(week) family(poisson) link(log) robust

xtcorr

lincom tx + 4 * weekXtx
test tx weekXtx

*** DHLZ p. 165
gen post = (week>0)
gen postXtx = post * tx
xtgee y post tx postXtx, offset(logObsTime) ///
  i(id) corr(exchangeable) family(poisson) link(log) robust

xtcorr

lincom tx + postXtx
test tx postXtx


*** GEE with BASELINE as covariate, and LINEAR model for time

xtgee y week tx weekXtx logY0 if week>0, offset(logObsTime) ///
  i(id) corr(unstructured) t(week) family(poisson) link(log) robust

xtcorr

test tx weekXtx
lincom tx + 4* weekXtx

*** using jackknife for standard error calculation

jknife "xtgee y post tx postXtx, offset(logObsTime) i(id) corr(exchangeable) family(poisson) link(log) robust" _b, cluster(id)