# See also the examples for binge in the documentation for the dataset binge.
# \donttest{
data(binge)
# matchNever creates the B-N match for the binge data.
# matchPast creates the B-P match for the binge data.
# Each match has its own propensity score, and these
# mean different things in the B-N match and the B-P
# match. The propensity score is denoted p.
# The Mahalanobis distance is the rank-based method
# described in Rosenbaum (2020, section 9.3).
# A directional caliper from Yu and Rosenbaum (2019)
# favors controls with high propensity scores, p.
#
# The two matches share most features, including:
# 1. A heavy left penalty for mismatching for female and bpRX.
# 2. Left penalties for mismatching for ageC, vigor, smokenow.
# 3. A left Mahalanobis distance for all covariates.
# 4. A heavy right penalty for mismatching smokenow.
# 5. Right penalties for education and smokeQuit.
# 6. A right Mahalanobis distance for just female, age and p.
# 7. An asymmetric, directional caliper on p.
# Because the left distance determines pairing, it emphasizes
# covariates that are out-of-balance and thought to be related
# to blood pressure. Because the right distance affects balance
# but not pairing, it worries about education and smoking in the
# distant past, as well as the propensity score which again is
# focused on covariate balance. The asymmetric caliper is
# tolerant of mismatches in the desired direction, so it too
# is placed on the right. Current smoking, smokenow, is placed
# on both the left and the right: even if we cannot always pair
# for it, perhaps we can balance it.
# In common practice, before examining outcomes, one compares
# several matched designs, fixing imperfections by adjusting
# penalties and other considerations.
#
# The B-P match is more difficult, because the P group is
# smaller than the N group. The B-P match uses the
# controlcosts feature while the B-N match does not.
# In the B-P match, there are too few young controls
# and too few controls with high propensity scores.
# Therefore, controlcosts penalize the use of controls
# with both low propensity scores (p<.4) and higher
# ages (age>42), thereby minimizing the use of these
# controls, even though the use of some of these
# controls is unavoidable in a pair match.
matchNever<-function(z,ncontrols=1){
dt<-binge[!is.na(z),]
z<-z[!is.na(z)]
dt<-dt[order(1-z,dt$SEQN),]
z<-z[order(1-z,dt$SEQN)]
rownames(dt)<-dt$SEQN
names(z)<-dt$SEQN
left<-startcost(z)
right<-startcost(z)
attach(dt)
propmod<-glm(z~age+female+education+smokenow+smokeQuit+bpRX+
bmi+vigor+waisthip,family=binomial)
p<-propmod$fitted.values
left<-addinteger(left,z,as.integer(ageC),penalty=100)
left<-addNearExact(left,z,female,penalty = 10000)
left<-addNearExact(left,z,bpRX,penalty = 10000)
left<-addNearExact(left,z,vigor,penalty=10)
left<-addinteger(left,z,smokenow,penalty=10)
left<-addMahal(left,z,cbind(age,bpRX,female,education,smokenow,
smokeQuit,bmi,vigor,waisthip))
right<-addMahal(right,z,cbind(female,age,p))
right<-addinteger(right,z,education,penalty=10)
right<-addinteger(right,z,smokenow,penalty=1000)
right<-addinteger(right,z,smokeQuit,penalty=10)
right<-addcaliper(right,z,p,caliper=c(-1,.03),penalty=10)
detach(dt)
dt<-cbind(dt,z,p)
m<-makematch(dt,left,right,ncontrols=ncontrols)
m$mset<-as.integer(m$mset)
treated<-m$SEQN[m$z==1]
treated<-as.vector(t(matrix(treated,length(treated),ncontrols+1)))
m<-cbind(m,treated)
list(m=m,dt=dt)
}
matchPast<-function(z,ncontrols=1){
dt<-binge[!is.na(z),]
z<-z[!is.na(z)]
dt<-dt[order(1-z,dt$SEQN),]
z<-z[order(1-z,dt$SEQN)]
rownames(dt)<-dt$SEQN
names(z)<-dt$SEQN
left<-startcost(z)
right<-startcost(z)
attach(dt)
propmod<-glm(z~age+female+education+smokenow+smokeQuit+bpRX+
bmi+vigor+waisthip,family=binomial)
p<-propmod$fitted.values
left<-addinteger(left,z,as.integer(ageC),penalty=100)
left<-addNearExact(left,z,female,penalty = 10000)
left<-addNearExact(left,z,bpRX,penalty = 10000)
left<-addNearExact(left,z,vigor,penalty=10)
left<-addinteger(left,z,smokenow,penalty=10)
left<-addMahal(left,z,cbind(age,bpRX,female,education,smokenow,
smokeQuit,bmi,vigor,waisthip))
right<-addMahal(right,z,cbind(female,age,p))
right<-addinteger(right,z,education,penalty=10)
right<-addinteger(right,z,smokenow,penalty=1000)
right<-addinteger(right,z,smokeQuit,penalty=10)
right<-addcaliper(right,z,p,caliper=c(-1,.03),penalty=10)
controlcosts<-((p[z==0]<.4)&(age[z==0]>42))*1000
detach(dt)
dt<-cbind(dt,z,p)
m<-makematch(dt,left,right,ncontrols=ncontrols,controlcosts=controlcosts)
m$mset<-as.integer(m$mset)
treated<-m$SEQN[m$z==1]
treated<-as.vector(t(matrix(treated,length(treated),ncontrols+1)))
m<-cbind(m,treated)
list(m=m,dt=dt)
}
z<-rep(NA,dim(binge)[1])
z[binge$AlcGroup=="B"]<-1
z[binge$AlcGroup=="P"]<-0
mPastComplete<-matchPast(z,ncontrols=1)
mPast<-mPastComplete$m
mPastComplete<-mPastComplete$dt
rm(z)
z<-rep(NA,dim(binge)[1])
z[binge$AlcGroup=="B"]<-1
z[binge$AlcGroup=="N"]<-0
mNeverComplete<-matchNever(z,ncontrols=1)
mNever<-mNeverComplete$m
mNeverComplete<-mNeverComplete$dt
rm(z)
bingeM<-rbind(mNever,mPast[mPast$z==0,])
bingeM<-bingeM[order(bingeM$treated,bingeM$AlcGroup,bingeM$SEQN),]
w<-which(colnames(bingeM)=="p")[1]
bingeM<-bingeM[,-w]
rm(binge,matchNever,matchPast,w)
old.par <- par(no.readonly = TRUE)
par(mfrow=c(1,2))
boxplot(mNever$p[mNever$z==1],mNever$p[mNever$z==0],
mNeverComplete$p[!is.element(mNeverComplete$SEQN,mNever$SEQN)],
names=c("B","mN","uN"),ylab="Propensity Score",main="Never",
ylim=c(0,.8))
boxplot(mPast$p[mPast$z==1],mPast$p[mPast$z==0],
mPastComplete$p[!is.element(mPastComplete$SEQN,mPast$SEQN)],
names=c("B","mP","uP"),ylab="Propensity Score",main="Past",
ylim=c(0,.8))
par(mfrow=c(1,2))
boxplot(mNever$age[mNever$z==1],mNever$age[mNever$z==0],
mNeverComplete$age[!is.element(mNeverComplete$SEQN,mNever$SEQN)],
names=c("B","mN","uN"),ylab="Age",main="Never",
ylim=c(20,80))
boxplot(mPast$age[mPast$z==1],mPast$age[mPast$z==0],
mPastComplete$age[!is.element(mPastComplete$SEQN,mPast$SEQN)],
names=c("B","mP","uP"),ylab="Age",main="Past",
ylim=c(20,80))
par(mfrow=c(1,2))
boxplot(mNever$education[mNever$z==1],mNever$education[mNever$z==0],
mNeverComplete$education[!is.element(mNeverComplete$SEQN,mNever$SEQN)],
names=c("B","mN","uN"),ylab="Education",main="Never",
ylim=c(1,5))
boxplot(mPast$education[mPast$z==1],mPast$education[mPast$z==0],
mPastComplete$education[!is.element(mPastComplete$SEQN,mPast$SEQN)],
names=c("B","mP","uP"),ylab="Education",main="Past",
ylim=c(1,5))
par(mfrow=c(1,2))
boxplot(mNever$bmi[mNever$z==1],mNever$bmi[mNever$z==0],
mNeverComplete$bmi[!is.element(mNeverComplete$SEQN,mNever$SEQN)],
names=c("B","mN","uN"),ylab="BMI",main="Never",
ylim=c(14,70))
boxplot(mPast$bmi[mPast$z==1],mPast$bmi[mPast$z==0],
mPastComplete$bmi[!is.element(mPastComplete$SEQN,mPast$SEQN)],
names=c("B","mP","uP"),ylab="BMI",main="Past",
ylim=c(14,70))
par(mfrow=c(1,2))
boxplot(mNever$waisthip[mNever$z==1],mNever$waisthip[mNever$z==0],
mNeverComplete$waisthip[!is.element(mNeverComplete$SEQN,mNever$SEQN)],
names=c("B","mN","uN"),ylab="Waist/Hip",main="Never",
ylim=c(.65,1.25))
boxplot(mPast$waisthip[mPast$z==1],mPast$waisthip[mPast$z==0],
mPastComplete$waisthip[!is.element(mPastComplete$SEQN,mPast$SEQN)],
names=c("B","mP","uP"),ylab="Waist/Hip",main="Past",
ylim=c(.65,1.25))
par(old.par)
# }
Run the code above in your browser using DataLab