#### Opening the data ####
# You need to put your address
aids<-read.table("C:/Users/laupaz/Desktop/Teaching & visiting/Causal inference (Milano Bicocca) Dicembre 2018/Lab/single.txt",header=TRUE)
# Visualizing the first 10 rows
aids[1:10,]
####### Unadjasted analysis in R of the relation among A and Y ########
chisq.test(x=aids$A,y=aids$Y)
####### Adjasted analysis in R of the relation among A and Y########
# Fitting the model in R without interaction:
summary(glm(formula=Y~A+L,family=binomial,data=aids))
# Adding an interaction term between A and L:
summary(glm(formula=Y~A+L+A*L,family=binomial,data=aids))
# STEP 1
# Fit a regression model for the outcome
fit0<-glm(formula=Y~A+L+A*L,family=binomial,data=aids)
# STEP 2
# Replace the factual exposure level with a, for each individual
# (for a=0)
aids0=aids
aids0$A=0
# STEP 3
# Estimate Pr(Y = 1|A = a; L) for each individual (for a=0)
# Installing package
library("stats", lib.loc="~/R/R-3.3.0/library")
pred0<-predict(object=fit0,newdata=aids0, type="respons")
# STEP 4
# Average these estimates over all individuals to obtain an estimate
# of Pr(Ya = 1) (for a=0)
p0<-mean(pred0)
p0
# STEP 1
# Fit a regression model for the outcome
fit1<-glm(formula=Y~A+L+A*L,family=binomial,data=aids)
# STEP 2
# Replace the factual exposure level with a, for each individual
# (for a=1)
aids1=aids
aids1$A=1
# STEP 3
# Estimate Pr(Y = 1|A = a; L) for each individual (for a=1)
pred1<-predict(object=fit1,newdata=aids1,type="respons")
# STEP 4
# Average these estimates over all individuals to obtain an estimate
# of Pr(Ya = 1) (for a=1)
p1<-mean(pred1)
p1
Marg_clor<-log((0.9625101/(1-0.9625101))/(0.3916425/(1-0.3916425)))
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))
Marg_clor
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))
Marg_clor
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))
fit<-glm(formula=A~L,family=binomial,data=aids)
# STEP 2 For each subject, use the fitted exposure model to
#estimate a subject-specific weight
pred<-predict(object=fit,type="respons")
w<-1/(aids$A*pred+(1-aids$A)*(1-pred))
# STEP 3 Use Pr(Y = 1|A = a) in the weighted sample as an
# estimate of Pr(Ya = 1)for a=0
p0<-weighted.mean(x=aids$Y[aids$A==0],w=w[aids$A==0])
p0
# STEP 3 Use Pr(Y = 1|A = a) in the weighted sample as an
# estimate of Pr(Ya = 1
p1<-weighted.mean(x=aids$Y[aids$A==1],w=w[aids$A==1])
p1
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))
Marg_clor
I(L^2)
I(aids$L^2)
aids<-read.table("C:/Users/laupaz/Desktop/Teaching & visiting/Causal inference (Milano Bicocca) Dicembre 2018/Lab/multiple.txt",header=TRUE)
# Visualizing the first 10 rows
aids[1:10,]
summary(Standard<-glm(formula=Y~A0+A1+L0,family=binomial(link = "logit"),data=aids))
summary(Sequentialt0<-glm(formula=Y~A0+L0,family=binomial,data=aids))
summary(Sequentialt1<-glm(formula=Y~A1+L0+A0+L1,family=binomial,data=aids))
# STEP 1
# Fit a regression model for the exposure at each time point, given the observed past up to that time point
fit0=glm(formula=A0~L0,family=binomial,data=aids)
fit1=glm(formula=A1~L0+A0+L1,family=binomial,data=aids)
pred0=predict(object=fit0,type="respons")
pred1=predict(object=fit1,type="respons")
w0=1/(aids$A0*pred0+(1-aids$A0)*(1-pred0))
w1=1/(aids$A1*pred1+(1-aids$A1)*(1-pred1))
w=w0*w1
summary(glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w))
#  Use this more elaborate model to estimate the joint effect
# of A0 and A1 (Beta0+Beta1) as los odds ratio
Joint_effect<-0.9510+1.4030 # 2.345
MSM<-summary(glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w))
Joint_effect<-coef(MSM)["A0"]+coef(MSM)["A1"] # 2.345
coef(MSM)["A0"]
summary(MSM)
# STEP 1
# Fit a regression model for the exposure at each time point, given the observed past up to that time point
fit0=glm(formula=A0~L0,family=binomial,data=aids)
fit1=glm(formula=A1~L0+A0+L1,family=binomial,data=aids)
# STEP 2
# For each subject, use the fitted exposure model to estimate a subject-specific weight W = 1/ Pr(A0|L0)*Pr(A1|L0;A0;L1)
pred0=predict(object=fit0,type="respons")
pred1=predict(object=fit1,type="respons")
w0=1/(aids$A0*pred0+(1-aids$A0)*(1-pred0))
w1=1/(aids$A1*pred1+(1-aids$A1)*(1-pred1))
w=w0*w1
MSM<-summary(glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w))
summary(MSM)
summary(glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w))
msm<-summary(glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w))
summary(msm)
msm<-glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w)
Joint_effect<-coef(msm)["A0"]+coef(msm)["A1"] # 2.345
msm1<-glm(formula=Y~A0+A1+A0*A1,family=quasibinomial,data=aids,weights=w)
summary(msm1)
Joint_effect<-coef(msm)["A0"]+coef(msm)["A1"]-coef(msm)["A0*A1"] # 1.5322
Joint_effect<-coef(msm1)["A0"]+coef(msm1)["A1"]-coef(msm1)["A0*A1"] # 1.5322
coef(msm1)["A0*A1"]
msm1<-glm(formula=Y~A0+A1+A0*A1,family=quasibinomial,data=aids,weights=w)
coef(msm1)["A0"]
coef(msm1)["A0*A1"]
summary(msm1)
Joint_effect<-coef(msm1)["A0"]+coef(msm1)["A1"]-coef(msm1)["A0:A1"] # 1.5322
coef(msm1)["A0:A1"]
Joint_effect<-coef(msm1)["A0"]+coef(msm1)["A1"]+coef(msm1)["A0:A1"] # 1.5322
library("sandwich", lib.loc="~/R/R-3.3.0/library")
library("lmtest", lib.loc="~/R/R-3.3.0/library")
model<-glm(Y~A0+A1,family=quasibinomial,data=aids,weights=w)
summary(model)
coeftest(model, vcov = sandwich)
aids<-read.table("C:/Users/laupaz/Desktop/Teaching & visiting/Causal inference (Milano Bicocca) Dicembre 2018/Lab/single.txt",header=TRUE)
chisq.test(x=aids$A,y=aids$Y)
summary(glm(formula=Y~A,family=binomial,data=aids))
summary(glm(formula=Y~A+L,family=binomial,data=aids))
summary(glm(formula=Y~A+L+A*L,family=binomial,data=aids))
library("AdhereR", lib.loc="~/R/R-3.3.0/library")
fit0<-glm(formula=Y~A+L+A*L,family=binomial,data=aids)
aids0=aids
aids0$A=0
# STEP 3
# Estimate Pr(Y = 1|A = a; L) for each individual (for a=0)
# Installing package
library("stats", lib.loc="~/R/R-3.3.0/library")
pred0<-predict(object=fit0,newdata=aids0, type="respons")
# of Pr(Ya = 1) (for a=0)
p0<-mean(pred0)
p0
# STEP 1
# Fit a regression model for the outcome
fit1<-glm(formula=Y~A+L+A*L,family=binomial,data=aids)
# STEP 2
# Replace the factual exposure level with a, for each individual
# (for a=1)
aids1=aids
aids1$A=1
# STEP 3
# Estimate Pr(Y = 1|A = a; L) for each individual (for a=1)
pred1<-predict(object=fit1,newdata=aids1,type="respons")
# STEP 4
# Average these estimates over all individuals to obtain an estimate
# of Pr(Ya = 1) (for a=1)
p1<-mean(pred1)
p1
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))
fit<-glm(formula=A~L,family=binomial,data=aids)
# STEP 2 For each subject, use the fitted exposure model to
#estimate a subject-specific weight
pred<-predict(object=fit,type="respons")
w<-1/(aids$A*pred+(1-aids$A)*(1-pred))
# STEP 2 For each subject, use the fitted exposure model to
#estimate a subject-specific weight
pred<-predict(object=fit,type="respons")
w<-1/(aids$A*pred+(1-aids$A)*(1-pred))
# STEP 3 Use Pr(Y = 1|A = a) in the weighted sample as an
# estimate of Pr(Ya = 1) for a=0
p0<-weighted.mean(x=aids$Y[aids$A==0],w=w[aids$A==0])
p0
# STEP 3 Use Pr(Y = 1|A = a) in the weighted sample as an
# estimate of Pr(Ya = 1) for a=1
p1<-weighted.mean(x=aids$Y[aids$A==1],w=w[aids$A==1])
p1
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))
aids<-read.table("C:/Users/laupaz/Desktop/Teaching & visiting/Causal inference (Milano Bicocca) Dicembre 2018/Lab/multiple.txt",header=TRUE)
# Visualizing the first 10 rows
aids[1:10,]
#### Ignoring L1 ####
# Standard adjustment in R fitting generalized linear models with glm
summary(Standard<-glm(formula=Y~A0+A1+L0,family=binomial(link = "logit"),data=aids))
summary(Sequentialt0<-glm(formula=Y~A0+L0,family=binomial,data=aids))
summary(Sequentialt1<-glm(formula=Y~A1+L0+A0+L1,family=binomial,data=aids))
# STEP 1
# Fit a regression model for the exposure at each time point, given the observed past up to that time point
fit0=glm(formula=A0~L0,family=binomial,data=aids)
fit1=glm(formula=A1~L0+A0+L1,family=binomial,data=aids)
pred0=predict(object=fit0,type="respons")
pred1=predict(object=fit1,type="respons")
w0=1/(aids$A0*pred0+(1-aids$A0)*(1-pred0))
w1=1/(aids$A1*pred1+(1-aids$A1)*(1-pred1))
w=w0*w1
msm<-glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w)
summary(msm)
Joint_effect<-coef(msm)["A0"]+coef(msm)["A1"] # 2.345
#  (b) Add an interaction term between a0 and a1 in the MSM.
msm1<-glm(formula=Y~A0+A1+A0*A1,family=quasibinomial,data=aids,weights=w)
summary(msm1)
#
Joint_effect<-coef(msm1)["A0"]+coef(msm1)["A1"]+coef(msm1)["A0:A1"] # 1.5322
# Sandwich estimator to have robust standard errors considering the weigthing system
# Installing packages
library("sandwich", lib.loc="~/R/R-3.3.0/library")
library("lmtest", lib.loc="~/R/R-3.3.0/library")
model<-glm(Y~A0+A1,family=quasibinomial,data=aids,weights=w)
summary(model)
coeftest(model, vcov = sandwich)
summary(model)
coeftest(model, vcov = sandwich)
msm<-glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w)
summary(msm)
