@@ -72,9 +72,9 @@ EC_IPW_OPT <- function(data,
7272 if ((! optimal_weight_flag ) && wt == 0 ) {
7373 # estimate ATE
7474 # # TODO: why do use the true propensity score?
75- temp <- df % > %
76- filter(S == 1 ) % > %
77- mutate(piA = sum(A ) / n ) % > %
75+ temp <- df | >
76+ filter(S == 1 ) | >
77+ mutate(piA = sum(A ) / n ) | >
7878 mutate(w11 = piA , w10 = 1 - piA )
7979
8080 # ## create outcomes: obs by T
@@ -112,14 +112,14 @@ EC_IPW_OPT <- function(data,
112112 # print(piS.model)
113113
114114 # estimate ATE
115- temp <- df % > %
115+ temp <- df | >
116116 mutate(
117117 piA = sum(A [S == 1 ]) / n ,
118118 piS = sum(S ) / (n + m ),
119119 piSX = predict(piS.model , newdata = df , type = " response" ),
120120 # piSX = exp(log(2) + df$X - 3 * df$U)/(1+exp(log(2) + df$X - 3 * df$U)),
121121 rx = (piSX / (1 - piSX )) * ((1 - piS ) / piS )
122- ) % > %
122+ ) | >
123123 mutate(w11 = piA , w10 = 1 - piA , w00 = rx )
124124
125125 # temp$w00[temp$w00 > 0.9] = 0.9
@@ -130,7 +130,10 @@ EC_IPW_OPT <- function(data,
130130 # ## create outcomes: obs * T
131131 Ys <- as.matrix(Y )
132132
133- potential <- data.frame ((temp $ S * temp $ A / temp $ w11 + temp $ S * (1 - temp $ A ) / temp $ w10 + (1 - temp $ S ) * temp $ w00 ) * Ys )
133+ potential <- data.frame (
134+ (temp $ S * temp $ A / temp $ w11 + temp $ S * (1 - temp $ A ) / temp $ w10 +
135+ (1 - temp $ S ) * temp $ w00 ) * Ys
136+ )
134137 mu1 <- colSums(potential [temp $ S == 1 & temp $ A == 1 , , drop = FALSE ]) / sum(temp $ S * temp $ A / temp $ w11 )
135138 mu10 <- colSums(potential [temp $ S == 1 & temp $ A == 0 , , drop = FALSE ]) / sum(temp $ S * (1 - temp $ A ) / temp $ w10 )
136139 mu00 <- colSums(potential [temp $ S == 0 , , drop = FALSE ]) / (sum((1 - temp $ S ) * temp $ w00 ))
@@ -143,8 +146,12 @@ EC_IPW_OPT <- function(data,
143146 A11 <- diag(rep(- 1 , T_follow ), nrow = T_follow )
144147 A22 <- diag(rep(- 1 , T_follow ), nrow = T_follow )
145148 A33 <- diag(rep(- mean((1 - temp $ S ) * temp $ w00 / (1 - temp $ piS )), T_follow ), nrow = T_follow )
146- A34 <- t((as.vector((1 - temp $ S ) * temp $ piSX / (temp $ piS * (1 - temp $ piSX ))) * sweep(Ys , 2 , mu00 ))) %*% model.matrix(piS.model ) / (n + m )
147- piS.beta <- t(model.matrix(piS.model )) %*% diag(- temp $ `piSX` * (1 - temp $ `piSX` )) %*% model.matrix(piS.model ) / (n + m )
149+ A34 <- t((as.vector((1 - temp $ S ) * temp $ piSX /
150+ (temp $ piS * (1 - temp $ piSX ))) * sweep(Ys , 2 , mu00 ))) %*%
151+ model.matrix(piS.model ) / (n + m )
152+ piS.beta <- t(model.matrix(piS.model )) %*%
153+ diag(- temp $ `piSX` * (1 - temp $ `piSX` )) %*%
154+ model.matrix(piS.model ) / (n + m )
148155 A44 <- piS.beta
149156
150157 n1 <- dim(A11 )[1 ]
@@ -179,7 +186,8 @@ EC_IPW_OPT <- function(data,
179186 # fit1 = lm(as.formula(paste("Y2","~",form_x)), data = filter(temp,S==1&A==0) )
180187 # sigma10 = (summary(fit1)$sigma)**2
181188 num <- sum(temp $ S * (1 - temp $ A ) / temp $ w10 ^ 2 / (sum(temp $ S * (1 - temp $ A ) / temp $ w10 ))^ 2 )
182- # + sum(temp$S*(1-temp$A))*var(temp$S*(1-temp$A)*predict(fit1,newdata =temp)/temp$w10/sum(temp$S*(1-temp$A)/temp$w10))
189+ # + sum(temp$S*(1-temp$A))*var(temp$S*(1-temp$A)* # nolint: line_length_linter.
190+ # + predict(fit1,newdata =temp)/temp$w10/sum(temp$S*(1-temp$A)/temp$w10))
183191
184192 # fit0 = lm(as.formula(paste("Y2","~",form_x)), data = filter(temp,S==0) )
185193 # sigma00 = (summary(fit0)$sigma)**2
@@ -216,8 +224,8 @@ EC_IPW_OPT <- function(data,
216224 cutoff <- qnorm(1 - alpha / 2 , lower.tail = TRUE )
217225
218226 if (Bootstrap == TRUE ) {
219- Group_ID <- df % > %
220- group_by(S , A ) % > %
227+ Group_ID <- df | >
228+ group_by(S , A ) | >
221229 mutate(group_id = cur_group_id())
222230 Group_ID <- Group_ID $ group_id
223231
@@ -267,14 +275,14 @@ EC_IPW_OPT <- function(data,
267275 lower_CI_boot = lower_CI_boot ,
268276 upper_CI_boot = upper_CI_boot
269277 )
270- return ( list (results = results , borrow_weight = wt ) )
278+ list (results = results , borrow_weight = wt )
271279 } else {
272280 results <- data.frame (
273281 point_estimates = tau ,
274282 standard_deviation = sd.tau ,
275283 lower_CI_normal = tau - sd.tau * cutoff ,
276284 upper_CI_normal = tau + sd.tau * cutoff
277285 )
278- return ( list (results = results , borrow_weight = wt ) )
286+ list (results = results , borrow_weight = wt )
279287 }
280288}
0 commit comments