/*SIMILAR TO SECTION 3.5 OF EXERCISES*/
/*APC in Triangular Dataset*/

***CHANGE WORKING DIRECTORY**
cd "G:\APC COURSE"

clear all


***Load the dataset**
infile A5 P5 C5 D Y up Ax Px Cx in 2/l using "lung5-Mc.txt"


**Look at the first few rows of the data**
list * in 1/10


***tabulate triangular values***
ta Ax
ta Px
ta Cx

***Quick loop to draw a Lexis Diagram with the triangular subsets in Stata***
***Could make into a general program if useful?***

**MAKE SURE TO RUN THE LOOP AND THE TWOWAY COMMAND AT THE SAME TIME AS ***
***ONLY USED LOCAL MACROS***

local fn 
forvalues i=1858(5)1953 {
	local j=`i'+40
	if `j'<1943 {
		local j=1943 
	}
	local k=`i'+85
	if `k'>1993 {
		local k=1993 
	}
	local fnt (function y=x-`i', range(`j' `k') lc(bluishgray) lw(medium))
	local fn `fn' `fnt'
}

local aspect=(85-45)/(1993-1943)

twoway  `fn' (scatter Ax Px  if Ax>65 & Ax<70 & Px<1973 & Px>1968, msize(small)), ///
legend(off) aspect(`aspect') xlabel(1943(10)1993,grid glc(bluishgray) ///
 glw(thin)) ylabel(40(5)85,grid glc(bluishgray) glw(medium)) xmtick(##2,grid  ///
 glc(bluishgray) glw(medium)) xtitle(Period) ytitle(Age) title(Lexis Diagram) ///
 name(Lexis, replace)
 

***Make the spline terms using correct midpoints of triangular subsets***
capture drop *spl*

rcsgen Ax, gen(Aspl) df(6) orthog
rcsgen Px, gen(Pspl) df(6) orthog
rcsgen Cx, gen(Cspl) df(6) orthog


**Make a constant term for the age values***
gen Aspl0=1

**Fit model with splines - no set constraints by us***
glm D Aspl* Pspl* Cspl*, family(poisson) lnoffset(Y) nocons
estimates store trianglesplines


partpred Aeff, for(Aspl*) ci(Aeff_lci Aeff_uci) eform
partpred Peff, for(Pspl*) ci(Peff_lci Peff_uci) eform 
partpred Ceff, for(Cspl*) ci(Ceff_lci Ceff_uci) eform 


**Make rates per 100000py**
replace Aeff=Aeff*100000
replace Aeff_lci=Aeff_lci*100000
replace Aeff_uci=Aeff_uci*100000


***Draw Graphs***
line Aeff* Ax, sort yscale(log) name(age, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Age) ytitle(Incidence Rate (per 100,000 person-years))

line Peff* Px, sort yscale(log) name(period, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Period) ytitle(RR)

line Ceff* Cx, sort yscale(log) name(cohort, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Cohort) ytitle(RR)


***Could Choose to Centre on first Cohort and first Period**
forvalues i=1/6 {
quietly su Px
quietly gen Psplref`i'=Pspl`i' if Px==r(min)
quietly egen Psplreff`i'=max(Psplref`i')

quietly su Cx
quietly gen Csplref`i'=Cspl`i' if Cx==r(min)
quietly egen Csplreff`i'=max(Csplref`i')

gen PsplFirst`i'=Pspl`i'-Psplreff`i'
gen CsplFirst`i'=Cspl`i'-Csplreff`i'
}

glm D Aspl* PsplFirst* CsplFirst*, family(poisson) lnoffset(Y) nocons
estimates store trianglesplinesFirst

***Same Model - Different Parameterisation***
estimates stats *

***Make the appropriate predictions***
partpred AeffFirst, for(Aspl*) eform ci(AeffFirst_lci AeffFirst_uci)
partpred PeffFirst, for(PsplF*) eform ci(PeffFirst_lci PeffFirst_uci)
partpred CeffFirst, for(CsplF*) eform ci(CeffFirst_lci CeffFirst_uci)

replace AeffFirst=AeffFirst*100000
replace AeffFirst_lci=AeffFirst_lci*100000
replace AeffFirst_uci=AeffFirst_uci*100000

***Draw Graphs***
line AeffFirst* Ax, sort yscale(log) name(ageFirst, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Age) ytitle(Incidence Rate (per 100,000 person-years))


line PeffFirst* Px, sort yscale(log) name(periodFirst, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Period) ytitle(RR)

line CeffFirst* Cx, sort yscale(log) name(cohortFirst, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Cohort) ytitle(RR)




**Fit Cohort terms first to change the constraint chosen by GLM***
glm D Aspl? Cspl? Pspl?, family(poisson) lnoffset(Y) nocons
estimates store trianglesplinesCP

*Make the predictions**
partpred AeffPcons, for(Aspl?) eform ci(AeffPcons_lci AeffPcons_uci)
partpred PeffPcons, for(Pspl?) eform ci(PeffPcons_lci PeffPcons_uci)
partpred CeffPcons, for(Cspl?) eform ci(CeffPcons_lci CeffPcons_uci)

**Make rates per 100000py**
replace AeffPcons=AeffPcons*100000
replace AeffPcons_lci=AeffPcons_lci*100000
replace AeffPcons_uci=AeffPcons_uci*100000

**Plot the results***
line AeffPcons* Ax, sort yscale(log) name(agePcons, replace) ///
lc(black black black) lp(solid dash dash) legend(off)

line PeffPcons* Px, sort yscale(log) name(periodPcons, replace) ///
lc(black black black) lp(solid dash dash) legend(off)

line CeffPcons* Cx, sort yscale(log) name(cohortPcons, replace) ///
lc(black black black) lp(solid dash dash) legend(off)


*******************
***Using apcfit****
*******************
capture drop agef* perf* cohf*
apcfit, age(Ax) period(Px) cases(D) poprisk(Y) nper(100000) ///
param(ACP) refc(1900)
return list

line agefitted* Ax, sort yscale(log) name(ageapcfit, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Age) ytitle(Incidence Rate (per 100,000 person-years))


line perfitted* Px, sort yscale(log) name(perapcfit, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Period) ytitle(RR)


line cohfitted* Cx, sort yscale(log) name(cohapcfit, replace) ///
lc(black black black) lp(solid dash dash) legend(off) ///
xtitle(Cohort) ytitle(RR)


***LOOK AT THE HELP FILE FOR SYNTAX AND OTHER OPTIONS***
help apcfit