*********************************************************************************************; ****** Autoimmunity, AFCR.raw data set *****************************; ****** Two treatments (AZ or AZ + MP) effect on AFCR over 18 months *****; ****** LDA Lab -11 -- Focus on logistic regression in LDA *************; ****** Stata version - 7 ************************************; *********************************************************************************************; ** Before doing analysis, please specify where you want to save your results ; * Turn off -more- pause; clear; set more off; * Save log file on disk; capture log close; log using "F:\lab11\afcr.log", replace; * Extend linesize for log; set linesize 100; **********************************************************************************************; ** Step -1, to read in data, use infile command ** ; **************************************************************; infile id time sqafcr group ptreat age using "C:\Documents and Settings\Sorina\Desktop\LDA\Lab11\afcr.raw", clear; codebook ** change group indicator to 0 and 1; replace group=group-1 ** generate binary outcome for longitudinal data **; ** Y=1 mean AFCR become lower - so treatment is good; summ sqafcr, detail gen Y=. replace Y=1 if sqafcr<=11 & sqafcr!=. replace Y=0 if sqafcr>11 & sqafcr!=. drop sqafcr **Data is already in long format.Recall the reshape commands: *reshape wide Y, i(id) j(time) *reshape long Y, i(id) j(time) save "C:\Documents and Settings\Sorina\Desktop\LDA\Lab11\afcr.dta", replace ** Step -2, EDA and model inference, based on cross-sectional data ** ; ** Y = the outcome at time=18 month, regress on baseline covariates time=0; *****************************************************************************************; ** EDA on the relationship between Y and Z & X, at baseline time **; tabulate Y group if time==18 csi 48 39 12 20 tabulate Y ptreat if time==18 csi 48 39 24 8 sort Y graph7 age if time==18, box by(Y) by Y: summ age if time==18 ** Logistic regression between Y and Z & X, at baseline time **; glm Y group age ptreat if time==18, f(bin) l(logit) glm, eform logit Y group age ptreat if time==18 ** same thng can be obtained with: logistic Y group age ptreat if time==18, coef ** directly report odds ratio estimates, rather than coefficients; logistic Y group age ptreat if time==18 ** Step -3, Logistic regression for Longitudinal data analysis setting ** ; ** (1) Marginal Logistic regression Model, without random effect *****; ****************************************************************************************; ** convert an ordinary data into a longitudinal dataset, specifying subject index and time index; tsset id time ** describe the pattern of data, the distribution of covariates (all categorical); xtdes xttab group xttab ptreat xttab Y sort Y by Y: xtsum age ** model based inference xi: xtgee Y time group age i.ptreat , nolog fam(bin) link(logit) corr(uns) test group xtcorr xi: xtgee Y time group age i.ptreat , nolog fam(bin) link(logit) corr(ind) test group xtcorr xi: xtgee Y time group age i.ptreat , nolog fam(bin) link(logit) corr(exc) test group xtcorr xi: xtgee Y time group age i.ptreat , nolog fam(bin) link(logit) corr(ar1) test group xtcorr ** using robust option to check inference ; xi: xtgee Y time group age i.ptreat , nolog fam(bin) link(logit) corr(uns) robust test group xtcorr xi: xtgee Y time group age i.ptreat , nolog fam(bin) link(logit) corr(ind) robust test group xtcorr xi: xtgee Y time group age i.ptreat , nolog fam(bin) link(logit) corr(exc) robust test group xtcorr ** using robust option to check inference, with interaction between group*time ; xi: xtgee Y time age ptreat i.group*time, nolog fam(bin) link(logit) corr(uns) robust; test _IgroXtime_1; xtcorr; xi: xtgee Y time age ptreat i.group*time, nolog fam(bin) link(logit) corr(ind) robust; test _IgroXtime_1; xtcorr; xi: xtgee Y time age ptreat i.group*time, nolog fam(bin) link(logit) corr(exc) robust; test _IgroXtime_1; xtcorr; ** Step -4, Logistic regression for Longitudinal data analysis setting ** ; ** (2) Logistic model with random effects, random intercept *****; ****************************************************************************************; ** default: the random intercept ~ N(0, \sigma^2) ** Do not need to specify correlation structure, since we assum Y is independent to each other conditional on random intercept; xi: xtlogit Y time group age i.ptreat , nolog i(id) re test group gllamm Y time group age ptreat , nolog fam(bin) link(logit) i(id ) ** Checking: interaction between group*time xi: xtlogit Y time group age i.ptreat i.group*time, nolog i(id) re test _IgroXtime_1 log close