Create your Gitee Account
Explore and code with more than 5 million developers,Free private repositories !:)
Sign up
This repository doesn't specify license. Without author's permission, this code is only for learning and cannot be used for other purposes.
Clone or download
chapter12.do 10.31 KB
Copy Edit Web IDE Raw Blame History
/***************************************************************
Stata code for Causal Inference: What If by Miguel Hernan & Jamie Robins
Date: 10/10/2019
Author: Eleanor Murray
For errors contact: ejmurray@bu.edu
****************************************************************/
/***************************************************************
PROGRAM 12.1
Descriptive statistics from NHEFS data (Table 12.1)
****************************************************************/
use "nhefs.dta", clear
*-Provisionally ignore subjects with missing values for follow-up weight
*-Sample size after exclusion: N = 1566
drop if wt82==.
*- Calculate mean weight change in those with and without smoking cessation
label define qsmk 0 "No smoking cessation" 1 "Smoking cessation"
label values qsmk qsmk
by qsmk, sort: egen years = mean(age) if age < .
label var years "Age, years"
by qsmk, sort: egen male = mean(100 * (sex==0)) if sex < .
label var male "Men, %"
by qsmk, sort: egen white = mean(100 * (race==0)) if race < .
label var white "White, %"
by qsmk, sort: egen university = mean(100 * (education == 5)) if education < .
label var university "University, %"
by qsmk, sort: egen kg = mean(wt71) if wt71 < .
label var kg "Weight, kg"
by qsmk, sort: egen cigs = mean(smokeintensity) if smokeintensity < .
label var cigs "Cigarettes/day"
by qsmk, sort: egen meansmkyrs = mean(smokeyrs) if smokeyrs < .
label var kg "Years smoking"
by qsmk, sort: egen noexer = mean(100 * (exercise == 2)) if exercise < .
label var noexer "Little/no exercise"
by qsmk, sort: egen inactive = mean(100 * (active==2)) if active < .
label var inactive "Inactive daily life"
*-Output table
foreach var of varlist years male white university kg cigs meansmkyrs noexer inactive {
tabdisp qsmk, cell(`var') format(%3.1f)
}
*------------------
* PROGRAM 12.2
* Estimating IP weights for Section 12.2
* Data from NHEFS
*------------------
*-Fit a logistic model for the IP weights
logit qsmk sex race c.age##c.age ib(last).education c.smokeintensity##c.smokeintensity ///
c.smokeyrs##c.smokeyrs ib(last).exercise ib(last).active c.wt71##c.wt71
*-Output predicted conditional probability of quitting smoking for each individual
predict p_qsmk, pr
*-Generate nonstabilized weights as P(A=1|covariates) if A = 1 and 1-P(A=1|covariates) if A = 0
gen w=.
replace w=1/p_qsmk if qsmk==1
replace w=1/(1-p_qsmk) if qsmk==0
*-Check the mean of the weights; we expect it to be close to 2.0
summarize w
*-Fit marginal structural model in the pseudopopulation
*-Weights assigned using pweight = w
*-Robust standard errors using cluster() option where 'seqn' is the ID variable
regress wt82_71 qsmk [pweight=w], cluster(seqn)
/**************************************************************
PROGRAM 12.3
Estimating stabilized IP weights for Section 12.3
Data from NHEFS
***************************************************************/
*-Fit a logistic model for the denominator of the IP weights and predict the conditional probability of smoking
logit qsmk sex race c.age##c.age ib(last).education c.smokeintensity##c.smokeintensity ///
c.smokeyrs##c.smokeyrs ib(last).exercise ib(last).active c.wt71##c.wt71
predict pd_qsmk, pr
*-Fit a logistic model for the numerator of ip weights and predict Pr(A=1)
logit qsmk
predict pn_qsmk, pr
*-Generate stabilized weights as f(A)/f(A|L)
gen sw_a=.
replace sw_a=pn_qsmk/pd_qsmk if qsmk==1
replace sw_a=(1-pn_qsmk)/(1-pd_qsmk) if qsmk==0
*-Check distribution of the stabilized weights
summarize sw_a
*-Fit marginal structural model in the pseudopopulation
regress wt82_71 qsmk [pweight=sw_a], cluster(seqn)
/**********************************************************
FINE POINT 12.2
Checking positivity
**********************************************************/
*-Check for missing values within strata of covariates, for example:
tab age qsmk if race==0 & sex==1 & wt82!=.
tab age qsmk if race==1 & sex==1 & wt82!=.
/***************************************************************
PROGRAM 12.4
Estimating the parameters of a marginal structural mean model
with a continuous treatment Data from NHEFS
Section 12.4
***************************************************************/
drop sw_a
*-Analysis restricted to subjects reporting <=25 cig/day at baseline: N = 1162
keep if smokeintensity <=25
*-Fit a linear model for the denominator of the IP weights and calculate the mean expected smoking intensity
regress smkintensity82_71 sex race c.age##c.age ib(last).education c.smokeintensity##c.smokeintensity ///
c.smokeyrs##c.smokeyrs ib(last).exercise ib(last).active c.wt71##c.wt71
quietly predict p_den
*-Generate the denisty of the denomiator expectation using the mean expected smoking intensity and the residuals, assuming a normal distribution
*-Note: The regress command in STATA saves the root mean squared error for the immediate regression as e(rmse), thus there is no need to calculate it again.
gen dens_den = normalden(smkintensity82_71, p_den, e(rmse))
*-Fit a linear model for the numerator of ip weights, calculate the mean expected value, and generate the density
quietly regress smkintensity82_71
quietly predict p_num
gen dens_num = normalden( smkintensity82_71, p_num, e(rmse))
*-Generate the final stabilized weights from the estimated numerator and denominator, and check the weights distribution
gen sw_a=dens_num/dens_den
summarize sw_a
*-Fit a marginal structural model in the pseudopopulation
regress wt82_71 c.smkintensity82_71##c.smkintensity82_71 [pweight=sw_a], cluster(seqn)
*-Output the estimated mean Y value when smoke intensity is unchanged from baseline to 1982
lincom _b[_cons]
*-Output the estimated mean Y value when smoke intensity increases by 20 from baseline to 1982
lincom _b[_cons] + 20*_b[smkintensity82_71 ] +400*_b[c.smkintensity82_71#c.smkintensity82_71]
/*************************************************************
PROGRAM 12.5
Estimating the parameters of a marginal structural logistic model
Data from NHEFS
Section 12.4
**************************************************************/
clear
use "nhefs.dta"
*-Provisionally ignore subjects with missing values for follow-up weight
*-Sample size after exclusion: N = 1566
drop if wt82==.
*-Estimate the stabilized weights for quitting smoking as in PROGRAM 12.3
*-Fit a logistic model for the denominator of the IP weights and predict the conditional probability of smoking
logit qsmk sex race c.age##c.age ib(last).education c.smokeintensity##c.smokeintensity ///
c.smokeyrs##c.smokeyrs ib(last).exercise ib(last).active c.wt71##c.wt71
predict pd_qsmk, pr
*-Fit a logistic model for the numerator of ip weights and predict Pr(A=1)
logit qsmk
predict pn_qsmk, pr
*-Generate stabilized weights as f(A)/f(A|L)
gen sw_a=.
replace sw_a=pn_qsmk/pd_qsmk if qsmk==1
replace sw_a=(1-pn_qsmk)/(1-pd_qsmk) if qsmk==0
summarize sw_a
*-Fit marginal structural model in the pseudopopulation
*-NOTE: Stata has two commands for logistic regression, logit and logistic
*-Using logistic allows us to output the odds ratios directly
*-We can also output odds ratios from the logit command using the or option (default logit output is regression coefficients
logistic death qsmk [pweight=sw_a], cluster(seqn)
***************************************************************
*PROGRAM 12.6
*Assessing effect modification by sex using a marginal structural mean model
*Data from NHEFS
*Section 12.5
**************************************************************
drop pd_qsmk pn_qsmk sw_a
*-Check distribution of sex
tab sex
*-Fit logistc model for the denominator of IP weights, as in PROGRAM 12.3
logit qsmk sex race c.age##c.age ib(last).education c.smokeintensity##c.smokeintensity ///
c.smokeyrs##c.smokeyrs ib(last).exercise ib(last).active c.wt71##c.wt71
predict pd_qsmk, pr
*-Fit logistic model for the numerator of IP weights, no including sex
logit qsmk sex
predict pn_qsmk, pr
*-Generate IP weights as before
gen sw_a=.
replace sw_a=pn_qsmk/pd_qsmk if qsmk==1
replace sw_a=(1-pn_qsmk)/(1-pd_qsmk) if qsmk==0
summarize sw_a
*-Fit marginal structural model in the pseudopopulation, including interaction term between quitting smoking and sex
regress wt82_71 qsmk##sex [pw=sw_a], cluster(seqn)
*****************************************************************
* PROGRAM 12.7
* Estimating IP weights to adjust for selection bias due to censoring
* Data from NHEFS
* Section 12.6
****************************************************************
clear
use "nhefs.dta"
*-Analysis including all individuals regardless of missing wt82 status: N=1629
*-Generate censoring indicator: C = 1 if wt82 missing
gen byte cens = (wt82 == .)
*-Check distribution of censoring by quitting smoking and baseline weight
tab cens qsmk, column
bys cens: summarize wt71
*-Fit logistic regression model for the denominator of IP weight for A
logit qsmk sex race c.age##c.age ib(last).education c.smokeintensity##c.smokeintensity ///
c.smokeyrs##c.smokeyrs ib(last).exercise ib(last).active c.wt71##c.wt71
predict pd_qsmk, pr
*-Fit logistic regression model for the numerator of IP weights for A
logit qsmk
predict pn_qsmk, pr
*-Fit logistic regression model for the denominator of IP weights for C, including quitting smoking
logit cens qsmk sex race c.age##c.age ib(last).education c.smokeintensity##c.smokeintensity ///
c.smokeyrs##c.smokeyrs ib(last).exercise ib(last).active c.wt71##c.wt71
predict pd_cens, pr
*-Fit logistic regression model for the numerator of IP weights for C, including quitting smoking
logit cens qsmk
predict pn_cens, pr
*-Generate the stabilized weights for A (sw_a)
gen sw_a=.
replace sw_a=pn_qsmk/pd_qsmk if qsmk==1
replace sw_a=(1-pn_qsmk)/(1-pd_qsmk) if qsmk==0
*-Generate the stabilized weights for C (sw_c)
*-NOTE: the conditional probability estimates generated by our logistic models for C represent the conditional probability of being censored (C=1)
*-We want weights for the conditional probability of bing uncensored, Pr(C=0|A,L)
gen sw_c=.
replace sw_c=(1-pn_cens)/(1-pd_cens) if cens==0
*-Generate the final stabilized weights and check distribution
gen sw=sw_a*sw_c
summarize sw
*-Fit marginal structural model in the pseudopopulation
regress wt82_71 qsmk [pw=sw], cluster(seqn)

Comment ( 0 )

Sign in for post a comment