# The two examples below are from Rosenbaum (2024).
###########
# Tighten match from 1-3 to 1-2 to adjust for BMI.
# BMI might be affected by alcohol consumption, so
# the primary comparison did not adjust for it.
# See Rosenbaum (1984).
###########
data(aHDLt)
result<-tighten(aHDLt,aHDLt$z,aHDLt$block,
x=cbind(aHDLt$age,aHDLt$education),
f=cbind(aHDLt$ibmi,(aHDLt$bmi>22.5)+(aHDLt$bmi>27.5)+(aHDLt$bmi>32.5)),
ncontrols=2)
omit<-aHDLt[!is.element(aHDLt$SEQN,result$SEQN),]
boxplot(result$bmi[result$z==1],result$bmi[result$z==0],omit$bmi,
names=c("D","C","O"),ylab="BMI")
boxplot(result$hdl[result$z==1],result$hdl[result$z==0],omit$hdl,
names=c("D","C","O"),ylab="HDL Cholesterol")
# Tightened blocks
table(aHDLt$z)
table(result$z)
table(table(aHDLt$block))
table(table(result$block))
table(table(result$mset))
##################
# Dentist in 1 Year Differential Effect
##################
#
# Example of tightening in support of a differential
# comparison to address a generic bias; see Rosenbaum (2006).
#
#
# First, create a data frame in which each block contains a
# differential comparison.
#
# Identify blocks in which the treated individual did not go
# to the dentist in the past year, but at lease one control did.
# Vector dife picks out the blocks that have this pattern.
data(aHDLt)
dif1<-tapply(((aHDLt$z==1)&(aHDLt$dentist1Y==0)),aHDLt$block,sum)==1
dif0<-tapply(((aHDLt$z==0)&(aHDLt$dentist1Y==1)),aHDLt$block,sum)>=1
dif<-dif1&dif0
dife<-as.vector(rbind(dif,dif,dif,dif))
elig<-((aHDLt$z==1)&(aHDLt$dentist1Y==0))|((aHDLt$z==0)&(aHDLt$dentist1Y==1))
# Now, form a data.frame for the people eligible for this comparison
# The tightened match will occur in this data.frame.
aHDLd<-cbind(aHDLt,elig)[dife,]
aHDLe<-aHDLd[((aHDLd$z==1)&(aHDLd$dentist1Y==0))|
((aHDLd$z==0)&(aHDLd$dentist1Y==1)),]
rm(dif,dif0,dif1,dife,elig,aHDLd)
# With data frame aHDLe, several tightened block designs
# are constructed.
x<-cbind(aHDLe$age,aHDLe$education)
# No pairs are deleted. Education is not quite balanced.
m1<-tighten(aHDLe,aHDLe$z,aHDLe$block,x=x,f=aHDLe$education)
table(m1$z,m1$education)
tapply(m1$age,m1$z,summary)
# 6 pairs are deleted. Education is alsmost balanced.
m2<-tighten(aHDLe,aHDLe$z,aHDLe$block,x=x,f=aHDLe$education,subset=150)
table(m2$z,m2$education)
tapply(m2$age,m2$z,summary)
# 13 pairs are deleted. Perhaps too many. Education is balanced.
m3<-tighten(aHDLe,aHDLe$z,aHDLe$block,x=x,f=aHDLe$education,subset=50)
table(m3$z,m3$education)
tapply(m3$age,m3$z,summary)
oldpar<-par(mfrow=c(1,3))
barplot(t(table(m1$education,1-m1$z)),beside=TRUE,ylim=c(0,50),
ylab="Count",xlab="Education Level",main="All 118 Pairs",
col=gray.colors(2))
barplot(t(table(m2$education,1-m2$z)),beside=TRUE,ylim=c(0,50),
ylab="Count",xlab="Education Level",main="Best 112 Pairs")
barplot(t(table(m3$education,1-m3$z)),beside=TRUE,ylim=c(0,50),
ylab="Count",xlab="Education Level",main="Best 105 Pairs")
par(oldpar)
Run the code above in your browser using DataLab