?pcor
parcor(X)
parcor(Cor.mat[irep,,])
nreps<-1000
p<-10
n<-100
Sigma<-rWishart(1,p,diag(rep(1,p))/p)[,,1]
Cor.mat<-PCF.mat<-array(0,c(nreps,p,p))
library(MASS)
library(ggm)
for(irep in 1:nreps){
X<-mvrnorm(n,rep(0,p),Sigma=Sigma)
Cor.mat[irep,,]<-cor(X)
PCF.mat[irep,,]<-parcor(Cor.mat[irep,,])
}
hist(PCF.mat[,1,1])
nreps
parcor(Cor.mat[irep,,])
hist(PCF.mat[,1,2])
diag(PCF.mat)<-NA
hist(PCF.mat[,1,3])
hist(PCF.mat[,1,2:3])
plot(PCF.mat[,1,2:3])
nreps<-10000
p<-10
n<-100
Sigma<-rWishart(1,p,diag(rep(1,p))/p)[,,1]
Cor.mat<-PCF.mat<-array(0,c(nreps,p,p))
library(MASS)
library(ggm)
for(irep in 1:nreps){
X<-mvrnorm(n,rep(0,p),Sigma=Sigma)
Cor.mat[irep,,]<-cor(X)
PCF.mat[irep,,]<-parcor(Cor.mat[irep,,])
}
apply(PCF.mat,2:3,range)
apply(PCF.mat,2:3,min)
apply(PCF.mat,2:3,max)
hist(PCF.mat[,8,5])
q()
1150.26+226.70-43.38 
182.32- 119.47 
21+31+3
21+31+30+8
(55/90)*86.52 + 0.21 
22+30+31+8
38*(62.85/91) 
62.85 + 95.43 + 53.08 + 26.25 
62.85/91
1333.58 + 237.61 
1571.19 - 1350
q()
install.packages('twang')
library(twang)
help(twang)
??twang
install.packages('gdata')
install.packages('ipw')
library(gdata)
nhefs <- read.xls("http://www.hsph.harvard.edu/miguel-hernan/files/2013/01/nhefs_book.xls")
nhefs <- read.csv("nhefs_book.csv",header=TRUE)
nhefs
## Prepare dataset
nhefs <- within(nhefs, {
    require(car)
    education <- car::recode(school,
                recodes = '0:8   = "1. 8th grade or less";
                           9:11  = "2. HS dropout";
                           12    = "3. HS";
                           13:15 = "4. College dropout";
                           16:30 = "5. College or more"')    
    education <- factor(education)
    university <- (school >= 16)
    age50 <- as.numeric(age > 50)
    ## Censoring indicator
    censored <- as.numeric(is.na(wt82_71))
    ## Reverse code to create non-censoring indicator
    uncensored <- as.numeric(censored == 0)
})
## Convert to factor
nhefs[,c("hbp","exercise","active","education")] <-
    lapply(nhefs[,c("hbp","exercise","active","education")], factor)
## Drop people without education
nhefs <- subset(nhefs, !is.na(education))
nhefs
## Load ipw package for IPW method
library(ipw)
## Create stabilized weights
res.st.iptw <-
    ipwpoint(exposure    = qsmk,
             family      = "binomial",
             link        = "logit",
             ## Model for numerator (unconditional probability of exposure status)
             numerator   = ~ 1,
             ## Model for denominator (conditional probability of exposure status given covariates)
             denominator = ~ sex + race + age + I(age^2) + education + smokeintensity + I(smokeintensity^2) + smokeyrs + I(smokeyrs^2) + exercise + active + hbp,
             data        = nhefs
             )
## Check stabilized weights (mean is approx 1)
summary(res.st.iptw$ipw.weights)
summary(res.st.iptw$num.mod)$coef
summary(res.st.iptw$den.mod)$coef
nhefs$st.ipcw <- res.st.ipcw$ipw.weights
nhefs$st.iptw <- res.st.iptw$ipw.weights
## Create stabilized weights
res.st.ipcw <-
    ipwpoint(exposure    = uncensored, # predict uncensored status
             family      = "binomial",
             link        = "logit",
             ## Model for numerator (conditional probability of uncensored given exposure)
             numerator   = ~ qsmk,
             ## Model for denominator (conditional probability of uncensored given exposure and covariates)
             denominator = ~ qsmk + sex + race + age + I(age^2) + education + smokeintensity + I(smokeintensity^2) + smokeyrs + I(smokeyrs^2) + exercise + active + hbp,
             data        = nhefs
             )
## Numerator model
summary(res.st.ipcw$num.mod)$coef
summary(res.st.ipcw$den.mod)$coef
nhefs$st.ipcw <- res.st.ipcw$ipw.weights
## Multiply for final weight
nhefs$st.iptcw <- with(nhefs, st.iptw * st.ipcw)
## Summarize
sapply(nhefs[!is.na(nhefs$wt82_71),c("st.iptw","st.ipcw","st.iptcw")], summary)
## Perform GEE clustered by nhefs, weighted by st.ipw, and using independet correlation structure
res.geeglm.st <- geeglm(wt82_71 ~ qsmk,
                        id      = seqn,
                        ## Explicitly drop censored patients (dropped anyway)
                        data    = subset(nhefs, censored == 0),
                        weights = st.iptcw,
                        corstr  = "independence")
sum.res.geeglm.st <- summary(res.geeglm.st)
sum.res.geeglm.st
## Create confidence intervals
confint.geeglm.st <- apply(sum.res.geeglm.st$coefficients, 1, function(VEC) {
    VEC[1] + c(-1,1) * qnorm(0.975) * VEC[2]
})
options(digits = 7)
t(confint.geeglm.st)
library(survey)
## Create a survey dataset
nhefs.svy <- svydesign(ids = ~ 1, data = subset(nhefs, censored == 0), weights = ~ st.iptcw)
## Fit regression
res.svyglm <- svyglm(formula = wt82_71 ~ qsmk, design = nhefs.svy)
## Show results
summary(res.svyglm)
options(digits = 7)
confint(res.svyglm)
?ipwpoint
haartdat
data(haartdat)
haartdat
summary(haartdat)
head(haartdat)
haartdat[1:20,]
nrow(haartdat)
hist(nhefs$wt82_71)
boxplot(wt82_71~qsmk,data=nhefs)
boxplot(wt82_71~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs)
boxplot(wt82_71~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Weight change (lbs) 1971-1982')
head(nhefs)
lm(wt82_71~qsmk,data=nhefs)
summary(lm(wt82_71~qsmk,data=nhefs))
boxplot(age~qsmk,data=nhefs)
boxplot(age~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Age')
summary(lm(age~qsmk,data=nhefs))
glm(qsmk~age,family=binomial,data=nhefs)
fq.age<-glm(qsmk~age,family=binomial,data=nhefs)
summary(fq.age)
ps.age<-fitted(fq.age)
boxplot(ps.age~qsmk,data=nhefs)
q()
install.packages('ipw')
library(gdata)
install.packages(c('ipw','gdata'))
library(ipw,gdata)
## Prepare dataset
nhefs <- within(nhefs, {
    require(car)
    education <- car::recode(school,
                recodes = '0:8   = "1. 8th grade or less";
                           9:11  = "2. HS dropout";
                           12    = "3. HS";
                           13:15 = "4. College dropout";
                           16:30 = "5. College or more"')    
    education <- factor(education)
    university <- (school >= 16)
    age50 <- as.numeric(age > 50)
    ## Censoring indicator
    censored <- as.numeric(is.na(wt82_71))
    ## Reverse code to create non-censoring indicator
    uncensored <- as.numeric(censored == 0)
})
install.packages(c('ipw','gdata','car'))
library(ipw,gdata)
nhefs <- read.csv("nhefs_book.csv",header=TRUE)
## Prepare dataset
nhefs <- within(nhefs, {
    require(car)
    education <- car::recode(school,
                recodes = '0:8   = "1. 8th grade or less";
                           9:11  = "2. HS dropout";
                           12    = "3. HS";
                           13:15 = "4. College dropout";
                           16:30 = "5. College or more"')    
    education <- factor(education)
    university <- (school >= 16)
    age50 <- as.numeric(age > 50)
    ## Censoring indicator
    censored <- as.numeric(is.na(wt82_71))
    ## Reverse code to create non-censoring indicator
    uncensored <- as.numeric(censored == 0)
})
## Convert to factor
nhefs[,c("hbp","exercise","active","education")] <-
    lapply(nhefs[,c("hbp","exercise","active","education")], factor)
## Drop people without education
nhefs <- subset(nhefs, !is.na(education))
require(ipw)
boxplot(wt82_71~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Weight change (lbs) 1971-1982')
summary(lm(wt82_71~qsmk,data=nhefs))
mean(wt82_71[qsmk==1])-mean(wt82_71[qsmk==0])
attach(nhefs)
mean(wt82_71[qsmk==1])-mean(wt82_71[qsmk==0])
q()
library(ipw,gdata)
## Prepare dataset
nhefs <- within(nhefs, {
    require(car)
    education <- car::recode(school,
                recodes = '0:8   = "1. 8th grade or less";
                           9:11  = "2. HS dropout";
                           12    = "3. HS";
                           13:15 = "4. College dropout";
                           16:30 = "5. College or more"')    
    education <- factor(education)
    university <- (school >= 16)
    age50 <- as.numeric(age > 50)
    ## Censoring indicator
    censored <- as.numeric(is.na(wt82_71))
    ## Reverse code to create non-censoring indicator
    uncensored <- as.numeric(censored == 0)
})
nhefs[,c("hbp","exercise","active","education")] <-
    lapply(nhefs[,c("hbp","exercise","active","education")], factor)
nhefs <- subset(nhefs, !is.na(education))
require(ipw)
boxplot(wt82_71~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Weight change (lbs) 1971-1982')
summary(lm(wt82_71~qsmk,data=nhefs))
mean(wt82_71[qsmk==1],na.rm=T)-mean(wt82_71[qsmk==0])
attach(nhefs)
mean(wt82_71[qsmk==1],na.rm=T)-mean(wt82_71[qsmk==0])
mean(wt82_71[qsmk==1],na.rm=T)
mean(wt82_71[qsmk==1],na.rm=T)-mean(wt82_71[qsmk==0],na.rm=T)
boxplot(age~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Age')
summary(lm(age~qsmk,data=nhefs))
summary(nhefs)
boxplot(ps.age~qsmk,data=nhefs)
summary(lm(age~qsmk,data=nhefs))
boxplot(age~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Age')
boxplot(age~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Age')
summary(lm(age~qsmk,data=nhefs))
boxplot(ps.age~qsmk,data=nhefs)
hist(ps.age)
plot(ps.age[nhefs$qsmk==1], col=rgb(0,0,1,1/4), xlim=c(0,1))
nhefs$qsmk==1
plot(ps.age[nhefs$qsmk==1], col=rgb(0,0,1,1/4))
quit.hist<-hist(ps.age[nhefs$qsmk==1])
plot(quit.hist, col=rgb(0,0,1,1/4))
ps.age[nhefs$qsmk==1]
quit.hist<-hist(ps.age[nhefs$qsmk==1])
plot(quit.hist, col=rgb(0,0,1,1/4))
quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4))
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
table(nhefs$qsmk)
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,100))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
names(nhefs)
fq.age<-glm(qsmk~age*sex,family=binomial,data=nhefs)
summary(fq.age)
ps.age<-fitted(fq.age)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
boxplot(ps.age~qsmk,data=nhefs)
fq.age<-glm(qsmk~age+exercise,family=binomial,data=nhefs)
summary(fq.age)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
boxplot(ps.age~qsmk,data=nhefs)
fq.age<-glm(qsmk~age+income,family=binomial,data=nhefs)
summary(fq.age)
ps.age<-fitted(fq.age)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
boxplot(ps.age~qsmk,data=nhefs)
boxplot(ps.age~qsmk,data=nhefs)
boxplot(ps.age~qsmk,subset=(!is.na(income),data=nhefs)
boxplot(ps.age~qsmk,subset=(!is.na(income)),data=nhefs)
boxplot(ps.age~qsmk[!is.na(income)],data=nhefs)
names(nhefs)
names(nhefs)
fq.age<-glm(qsmk~age+income+race,family=binomial,data=nhefs)
summary(fq.age)
nhefs$race
fq.age<-glm(qsmk~age+income+race,family=binomial,data=nhefs)
summary(fq.age)
ps.age<-fitted(fq.age)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
smoke.age.hist<-hist(nhefs.age[nhefs$qsmk==1])
smoke.age.hist<-hist(nhefs$age[nhefs$qsmk==1])
quit.age.hist<-hist(nhefs$age[nhefs$qsmk==1])
no.quit.age.hist<-hist(nhefs$age[nhefs$qsmk==0])
plot(quit.age.hist, col=rgb(0,0,1,1/4),xlim=range(0,100),ylim=range(0,200))
plot(no.quit.age.hist, col=rgb(1,0,1,1/4),add=T)
fq.age<-glm(qsmk~age+income+race,family=binomial,data=nhefs)
summary(fq.age)
ps.age<-fitted(fq.age)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
boxplot(ps.age~qsmk[!is.na(income)],data=nhefs)
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,0,1,1/4),add=T)
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,1,1,1/4),add=T)
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(0,0,1/2,1/4),add=T)
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200))
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200),
xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
box()
legend(0,200,c('Quit=Y','Quit=N'),col=c(rgb(0,0,1,1/4),rgb(1,0,1/2,1/4)))
legend(0,200,c('Quit=Y','Quit=N'),fill=c(rgb(0,0,1,1/4),rgb(1,0,1/2,1/4)))
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200),
xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
legend(0,200,c('Quit=Y','Quit=N'),fill=c(rgb(0,0,1,1/4),rgb(1,0,1/2,1/4)))
box()
psq<-quantile(ps.age,prob=c(1:5)/5)
psq
cut
?cut
help(cut)
ps.strata<-cut(ps.age,psq,include.lowest=TRUE)
ps.strata
ps.age
range(ps.age)
psq<-quantile(ps.age,prob=c(0:5)/5)
ps.strata<-cut(ps.age,psq,include.lowest=TRUE)
ps.strata
table(ps.strata)
sum(is.na(ps.strata))
sum(is.na(nhefs$income))
hist(nhefs$income)
boxplot(income~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Income')
nhefs$income[qsmk==0]
nhefs$income[qsmk==1]
args(cut)
cut
?cut
help.search("cut")
help("cut")
iq<-1
f0<-lm(wt82_71~qsmk,subset=(ps.strata==iq),data=nhefs)
ps.strata
q()
help(cut)
options(browser='firefox')
help(cut)
options(browser='C:\\Program Files (x86)\\Mozilla Firefox\\firefox.exe')
help(cut)
boxplot(ps.age~qsmk[!is.na(income)],data=nhefs)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200),
xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
legend(0,200,c('Quit=Y','Quit=N'),fill=c(rgb(0,0,1,1/4),rgb(1,0,1/2,1/4)))
box()
ps.strata<-cut(ps.age,psq,include.lowest=TRUE,dig.lab=3)
ps.strata
ps.strata<-cut(ps.age,psq,include.lowest=TRUE,labels=FALSE)
ps.strata
f1<-lm(wt82_71~qsmk,subset=(ps.strata==1),data=nhefs)
f2<-lm(wt82_71~qsmk,subset=(ps.strata==2),data=nhefs)
f3<-lm(wt82_71~qsmk,subset=(ps.strata==3),data=nhefs)
f4<-lm(wt82_71~qsmk,subset=(ps.strata==4),data=nhefs)
f5<-lm(wt82_71~qsmk,subset=(ps.strata==5),data=nhefs)
coef(f1)
summary(f1)
summary(lm(wt82_71~qsmk,data=nhefs))
f1<-lm(wt82_71~qsmk,subset=(ps.strata==1),data=nhefs)
f2<-lm(wt82_71~qsmk,subset=(ps.strata==2),data=nhefs)
f3<-lm(wt82_71~qsmk,subset=(ps.strata==3),data=nhefs)
f4<-lm(wt82_71~qsmk,subset=(ps.strata==4),data=nhefs)
f5<-lm(wt82_71~qsmk,subset=(ps.strata==5),data=nhefs)
summary(f2)
summary(f3)
summary(f4)
summary(f5)
coef(f1)
coef(f1)[2]
(ps.adj.est<-(coef(f1)[2]+coef(f2)[2]+coef(f3)[2]+coef(f4)[2]+coef(f5)[2])/5)
summary(lm(age~qsmk,data=nhefs))
mean(wt82_71[qsmk==1],na.rm=T)-mean(wt82_71[qsmk==0],na.rm=T)
summary(lm(wt82_71~qsmk,data=nhefs))
f1<-lm(wt82_71~qsmk+income,subset=(ps.strata==1),data=nhefs);summary(f1)
f2<-lm(wt82_71~qsmk+income,subset=(ps.strata==2),data=nhefs);summary(f2)
f3<-lm(wt82_71~qsmk+income,subset=(ps.strata==3),data=nhefs);summary(f3)
f4<-lm(wt82_71~qsmk+income,subset=(ps.strata==4),data=nhefs);summary(f4)
f5<-lm(wt82_71~qsmk+income,subset=(ps.strata==5),data=nhefs);summary(f5)
(ps.adj.est<-(coef(f1)[2]+coef(f2)[2]+coef(f3)[2]+coef(f4)[2]+coef(f5)[2])/5)
summary(lm(wt82_71~qsmk,data=nhefs))
summary(lm(wt82_71~qsmk+age,data=nhefs))
summary(lm(wt82_71~qsmk+income,data=nhefs))
f.all<-lm(wt82_71~qsmk*as.factor(ps.strata),data=nhefs);summary(f.all)
f.all<-lm(wt82_71~qsmk*as.factor(ps.strata),subset=(!is.na(income)),data=nhefs);summary(f.all)
ps.strata
length(ps.strata)
nrow(nhefs)
ps.age
length(ps.age)
ps.age<-rep(NA,nrow(fq.age))
ps.age<-rep(NA,nrow(nhefs))
ps.age[!is.na(nhefs$income)]<-fitted(fq.age)
ps.age
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200),
xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
fq.age<-glm(qsmk~age+income+race,family=binomial,data=nhefs)
summary(fq.age)
ps.age<-rep(NA,nrow(nhefs))
ps.age[!is.na(nhefs$income)]<-fitted(fq.age)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,200),
xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
legend(0,200,c('Quit=Y','Quit=N'),fill=c(rgb(0,0,1,1/4),rgb(1,0,1/2,1/4)))
box()
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,500),
xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
boxplot(ps.age~qsmk[!is.na(income)],data=nhefs)
boxplot(ps.age~qsmk,data=nhefs)
psq<-quantile(ps.age,prob=c(0:5)/5)
ps.strata<-cut(ps.age,psq,include.lowest=TRUE,labels=FALSE)
f1<-lm(wt82_71~qsmk,subset=(ps.strata==1),data=nhefs);summary(f1)
f2<-lm(wt82_71~qsmk,subset=(ps.strata==2),data=nhefs);summary(f2)
f3<-lm(wt82_71~qsmk,subset=(ps.strata==3),data=nhefs);summary(f3)
f4<-lm(wt82_71~qsmk,subset=(ps.strata==4),data=nhefs);summary(f4)
f5<-lm(wt82_71~qsmk,subset=(ps.strata==5),data=nhefs);summary(f5)
(ps.adj.est<-(coef(f1)[2]+coef(f2)[2]+coef(f3)[2]+coef(f4)[2]+coef(f5)[2])/5)
f.all<-lm(wt82_71~qsmk*as.factor(ps.strata),subset=(!is.na(income)),data=nhefs);summary(f.all)
!is.na(nhefs$income)
fq.age<-glm(qsmk~age+income+race,family=binomial,data=nhefs)
summary(fq.age)
ps.age<-rep(NA,nrow(nhefs))
ps.age[!is.na(nhefs$income)]<-fitted(fq.age)
quit.hist<-hist(ps.age[nhefs$qsmk==1])
no.quit.hist<-hist(ps.age[nhefs$qsmk==0])
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,500),
xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
legend(0,200,c('Quit=Y','Quit=N'),fill=c(rgb(0,0,1,1/4),rgb(1,0,1/2,1/4)))
box()
2
plot(quit.hist, col=rgb(0,0,1,1/4),xlim=range(0,0.5),ylim=range(0,500),
xlab='Propensity Score value',main='PS values')
plot(no.quit.hist, col=rgb(1,0,1/2,1/4),add=T)
legend(0,500,c('Quit=Y','Quit=N'),fill=c(rgb(0,0,1,1/4),rgb(1,0,1/2,1/4)))
box()
boxplot(ps.age~qsmk,data=nhefs)
psq<-quantile(ps.age,prob=c(0:5)/5)
ps.strata<-cut(ps.age,psq,include.lowest=TRUE,labels=FALSE)
f1<-lm(wt82_71~qsmk,subset=(ps.strata==1),data=nhefs);summary(f1)
f2<-lm(wt82_71~qsmk,subset=(ps.strata==2),data=nhefs);summary(f2)
f3<-lm(wt82_71~qsmk,subset=(ps.strata==3),data=nhefs);summary(f3)
f4<-lm(wt82_71~qsmk,subset=(ps.strata==4),data=nhefs);summary(f4)
f5<-lm(wt82_71~qsmk,subset=(ps.strata==5),data=nhefs);summary(f5)
(ps.adj.est<-(coef(f1)[2]+coef(f2)[2]+coef(f3)[2]+coef(f4)[2]+coef(f5)[2])/5)
summary(fq.age)
boxplot(wt82_71~qsmk,names=c('Quit=N','Quit=Y'),data=nhefs,main='Weight change (lbs) 1971-1982')
summary(lm(wt82_71~qsmk+income,data=nhefs))
summary(lm(wt82_71~qsmk,data=nhefs))
q()
