Skip to content

Commit 0e349a8

Browse files
mattsecrestmattsecrest
authored
Merge pull request #8 from Genentech/7-clean
7 clean
2 parents 7654913 + 21c6580 commit 0e349a8

83 files changed

Lines changed: 3557 additions & 3172 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎.Rbuildignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,3 +3,4 @@
33
^LICENSE\.md$
44
^NOTES\.md$
55
^\.github$
6+
^\.lintr$

‎.lintr‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
linters: linters_with_defaults(
2+
# Suppress: domain acronyms (OLE, AIPW, IPW, DID, SCM) and
3+
# math/stats variables (X, Y, T, S, A) are not snake_case
4+
object_name_linter = NULL,
5+
# Suppress: S4 constructors (new(), setClass()) trigger false positives
6+
object_usage_linter = NULL,
7+
# Suppress: math formulas need tight spacing (e.g. x1*1, n-1)
8+
infix_spaces_linter = NULL,
9+
# Suppress: legacy commented-out code is kept intentionally for now
10+
commented_code_linter = NULL,
11+
# Suppress: continuation-line indentation style is consistent as-is
12+
indentation_linter = NULL,
13+
# 120 chars is a reasonable limit for this codebase
14+
line_length_linter(120)
15+
)
16+
exclusions: list(
17+
"inst/" = Inf
18+
)

‎DESCRIPTION‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,7 @@ Suggests:
6666
rmarkdown,
6767
lintr,
6868
spelling,
69+
styler,
6970
testthat (>= 3.0.0)
7071
Config/testthat/edition: 3
7172
Collate:

‎R/DID_EC_AIPW.R‎

Lines changed: 135 additions & 126 deletions
Original file line numberDiff line numberDiff line change
@@ -5,12 +5,14 @@
55
#' @param model_form_mu0_ext character. The model formula for the outcome model in the external data (mu0_ext).
66
#' @param data data.frame. The input data containing the outcome, trial status, treatment, and covariates.
77
#' @param outcome_col_name character. The column name for the outcome variable in the data.
8-
#' @param trial_status_col_name character. The column name for the trial status variable in the data (indicating RCT vs external control).
8+
#' @param trial_status_col_name character. The column name for the trial
9+
#' status variable in the data (indicating RCT vs external control).
910
#' @param treatment_col_name character. The column name for the treatment variable in the data.
1011
#' @param covariates_col_name character vector. The column names for the covariates in the data.
1112
#' @param Bootstrap logical. Whether to use bootstrap for inference.
1213
#' @param R numeric. The number of bootstrap replications.
13-
#' @param bootstrap_CI_type character. The type of bootstrap confidence interval to compute (e.g., "bca", "norm", "perc", "basic", "stud").
14+
#' @param bootstrap_CI_type character. The type of bootstrap confidence
15+
#' interval to compute (e.g., "bca", "norm", "perc", "basic", "stud").
1416
#' @param alpha numeric. The significance level for confidence intervals.
1517
#' @param T_cross numeric. The time point that separates the placebo-control period and the follow-up period.
1618
#' @param quiet Logical. If \code{TRUE}, suppress printed output.
@@ -19,141 +21,148 @@
1921
#' @return tau and standard deviation
2022
#' @export
2123
#'
22-
DID_EC_AIPW = function(data,
23-
outcome_col_name,
24-
trial_status_col_name,
25-
treatment_col_name,
26-
covariates_col_name,
27-
T_cross,
28-
model_form_piS = "",
29-
model_form_piA = "",
30-
model_form_mu0_ext = "",
31-
Bootstrap = F,
32-
R = 5e2,
33-
bootstrap_CI_type = "bca",
34-
alpha = 0.05,
35-
quiet = TRUE
36-
){
37-
38-
Y = subset(data, select = outcome_col_name)
39-
S = subset(data, select = trial_status_col_name)
40-
A = subset(data, select = treatment_col_name)
41-
X = subset(data, select = covariates_col_name)
42-
names(S) = "S"
43-
names(A) = "A"
44-
45-
model_form_piS = unlist(strsplit(model_form_piS, "~"))
46-
model_form_piS[1] = "S"
47-
model_form_piS = paste0(model_form_piS, collapse = "~")
48-
49-
model_form_piA = unlist(strsplit(model_form_piA, "~"))
50-
model_form_piA[1] = "A"
51-
model_form_piA = paste0(model_form_piA, collapse = "~")
52-
53-
df = data.frame(Y, S = S, A = A, X)
54-
N = nrow(df)
55-
T_follow = ncol(Y)
56-
T_pc = T_cross # placebo-control period
57-
n = sum(df$S) # RCT sample size
58-
m = sum(1 - df$S) # external control sample size
59-
pi.S = n/N
60-
61-
#
62-
piS = glm(as.formula(model_form_piS), data = df, family = "binomial")
63-
piSX = predict(piS, newdata = filter(df), type = "response")
64-
if (model_form_piA == ""){
65-
piAX = sum(A)/n
66-
}else{
67-
piA = glm(as.formula(model_form_piA), data = filter(df, S == 1), family = "binomial")
68-
piAX = predict(piA, newdata = filter(df), type = "response")
24+
DID_EC_AIPW <- function(data,
25+
outcome_col_name,
26+
trial_status_col_name,
27+
treatment_col_name,
28+
covariates_col_name,
29+
T_cross,
30+
model_form_piS = "",
31+
model_form_piA = "",
32+
model_form_mu0_ext = "",
33+
Bootstrap = FALSE,
34+
R = 5e2,
35+
bootstrap_CI_type = "bca",
36+
alpha = 0.05,
37+
quiet = TRUE) {
38+
Y <- subset(data, select = outcome_col_name)
39+
S <- subset(data, select = trial_status_col_name)
40+
A <- subset(data, select = treatment_col_name)
41+
X <- subset(data, select = covariates_col_name)
42+
names(S) <- "S"
43+
names(A) <- "A"
44+
45+
model_form_piS <- unlist(strsplit(model_form_piS, "~"))
46+
model_form_piS[1] <- "S"
47+
model_form_piS <- paste0(model_form_piS, collapse = "~")
48+
49+
model_form_piA <- unlist(strsplit(model_form_piA, "~"))
50+
model_form_piA[1] <- "A"
51+
model_form_piA <- paste0(model_form_piA, collapse = "~")
52+
53+
df <- data.frame(Y, S = S, A = A, X)
54+
N <- nrow(df)
55+
T_follow <- ncol(Y)
56+
T_pc <- T_cross # placebo-control period
57+
n <- sum(df$S) # RCT sample size
58+
m <- sum(1 - df$S) # external control sample size
59+
pi.S <- n / N
60+
61+
#
62+
piS <- glm(as.formula(model_form_piS), data = df, family = "binomial")
63+
piSX <- predict(piS, newdata = filter(df), type = "response")
64+
if (model_form_piA == "") {
65+
piAX <- sum(A) / n
66+
} else {
67+
piA <- glm(as.formula(model_form_piA), data = filter(df, S == 1), family = "binomial")
68+
piAX <- predict(piA, newdata = filter(df), type = "response")
6969
}
70-
70+
7171
# predict Y0 from outcome regression models
72-
model_list_ext = lapply(1:T_follow, function(x){
73-
assign(paste0("m.ext", x), lm(as.formula(model_form_mu0_ext[x]), data = filter(df, S==0)))
72+
model_list_ext <- lapply(1:T_follow, function(x) {
73+
assign(paste0("m.ext", x), lm(as.formula(model_form_mu0_ext[x]), data = filter(df, S == 0)))
7474
})
75-
Y0 = data.frame(sapply(1:T_follow, function(x){predict(model_list_ext[[x]], newdata = filter(df))}))
76-
colnames(Y0) = paste0("y", 1:T_follow, "_0")
75+
Y0 <- data.frame(sapply(1:T_follow, function(x) {
76+
predict(model_list_ext[[x]], newdata = filter(df))
77+
}))
78+
colnames(Y0) <- paste0("y", 1:T_follow, "_0")
7779
# for residual
78-
Yr = Y - Y0
79-
colnames(Yr) = paste0("y", 1:T_follow, "_r")
80-
81-
temp = df %>% cbind(., Y0, Yr) %>%
82-
mutate(piAX = piAX, #sum(A)/sum(S)
83-
piSX = piSX,
84-
rx = piSX * (1 - pi.S)/(1 - piSX)/pi.S) %>%
85-
mutate(w11 = 1/piAX, w10 = 1/(1 - piAX), w00 = rx)
86-
80+
Yr <- Y - Y0
81+
colnames(Yr) <- paste0("y", 1:T_follow, "_r")
82+
83+
temp <- cbind(df, Y0, Yr) |>
84+
mutate(
85+
piAX = piAX,
86+
piSX = piSX,
87+
rx = piSX * (1 - pi.S) / (1 - piSX) / pi.S
88+
) |>
89+
mutate(w11 = 1 / piAX, w10 = 1 / (1 - piAX), w00 = rx)
90+
8791
# create outcomes
88-
Ys = as.matrix(Yr)
89-
potential = (temp$S*temp$A*temp$w11/sum(temp$S*temp$A*temp$w11) +
90-
temp$S*(1-temp$A)*temp$w10/sum(temp$S*(1-temp$A)*temp$w10) +
91-
(1-temp$S)*temp$w00/sum((1-temp$S)*temp$w00)) * Ys
92-
93-
94-
mu_S1A1 = colSums(data.frame(potential[temp$S == 1 & temp$A == 1, (T_pc+1):T_follow]))
95-
mu_S0A0 = colSums(data.frame(potential[temp$S == 0, (T_pc+1):T_follow]))
96-
bias = sum(rowMeans(data.frame(potential[temp$S==1 & temp$A==0, 1:T_pc]))) -
97-
sum(rowMeans(data.frame(potential[temp$S==0, 1:T_pc])))
98-
99-
100-
tau = mu_S1A1 - mu_S0A0 - bias
101-
102-
names(tau) = paste0("tau", (T_pc+1):T_follow)
103-
92+
Ys <- as.matrix(Yr)
93+
potential <- (temp$S * temp$A * temp$w11 / sum(temp$S * temp$A * temp$w11) +
94+
temp$S * (1 - temp$A) * temp$w10 / sum(temp$S * (1 - temp$A) * temp$w10) +
95+
(1 - temp$S) * temp$w00 / sum((1 - temp$S) * temp$w00)) * Ys
96+
97+
98+
mu_S1A1 <- colSums(data.frame(potential[temp$S == 1 & temp$A == 1, (T_pc + 1):T_follow]))
99+
mu_S0A0 <- colSums(data.frame(potential[temp$S == 0, (T_pc + 1):T_follow]))
100+
bias <- sum(rowMeans(data.frame(potential[temp$S == 1 & temp$A == 0, 1:T_pc]))) -
101+
sum(rowMeans(data.frame(potential[temp$S == 0, 1:T_pc])))
102+
103+
104+
tau <- mu_S1A1 - mu_S0A0 - bias
105+
106+
names(tau) <- paste0("tau", (T_pc + 1):T_follow)
107+
104108
# summarize results
105-
106-
cutoff = qnorm(1-alpha/2, lower.tail = T)
107-
108-
if(Bootstrap){
109-
Group_ID = df %>% group_by(S, A) %>% mutate(group_id = cur_group_id())
110-
Group_ID = Group_ID$group_id
111-
112-
boot.ci.type = switch (bootstrap_CI_type,
113-
norm = "normal",
114-
bca = "bca",
115-
stud = "student",
116-
perc = "percent",
117-
basic = "basic"
109+
110+
cutoff <- qnorm(1 - alpha / 2, lower.tail = TRUE)
111+
112+
if (Bootstrap) {
113+
Group_ID <- df |>
114+
group_by(S, A) |>
115+
mutate(group_id = cur_group_id())
116+
Group_ID <- Group_ID$group_id
117+
118+
boot.ci.type <- switch(bootstrap_CI_type,
119+
norm = "normal",
120+
bca = "bca",
121+
stud = "student",
122+
perc = "percent",
123+
basic = "basic"
124+
)
125+
126+
boot.out <- boot(
127+
data = df,
128+
statistic = DID_EC_AIPW_bootstrap,
129+
outcome_col_name = outcome_col_name,
130+
trial_status_col_name = trial_status_col_name,
131+
treatment_col_name = treatment_col_name,
132+
covariates_col_name = covariates_col_name,
133+
T_cross = T_cross,
134+
model_form_piS = model_form_piS,
135+
model_form_piA = model_form_piA,
136+
model_form_mu0_ext = model_form_mu0_ext,
137+
R = R,
138+
strata = Group_ID
139+
)
140+
141+
lower_CI_boot <- sapply(
142+
1:(T_follow - T_pc),
143+
function(x) {
144+
ci <- boot.ci(boot.out, conf = 1 - alpha, type = bootstrap_CI_type, index = x)
145+
ci[[boot.ci.type]][4]
146+
}
147+
)
148+
149+
upper_CI_boot <- sapply(
150+
1:(T_follow - T_pc),
151+
function(x) {
152+
ci <- boot.ci(boot.out, conf = 1 - alpha, type = bootstrap_CI_type, index = x)
153+
ci[[boot.ci.type]][5]
154+
}
118155
)
119-
120-
boot.out <- boot(data = df,
121-
statistic = DID_EC_AIPW_bootstrap,
122-
outcome_col_name = outcome_col_name,
123-
trial_status_col_name = trial_status_col_name,
124-
treatment_col_name = treatment_col_name,
125-
covariates_col_name = covariates_col_name,
126-
T_cross = T_cross,
127-
model_form_piS = model_form_piS,
128-
model_form_piA = model_form_piA,
129-
model_form_mu0_ext = model_form_mu0_ext,
130-
R = R,
131-
strata = Group_ID)
132-
133-
lower_CI_boot = sapply(1:(T_follow - T_pc),
134-
function(x) {
135-
ci = boot.ci(boot.out, conf = 1 - alpha, type = bootstrap_CI_type, index = x)
136-
ci[[boot.ci.type]][4]
137-
})
138-
139-
upper_CI_boot = sapply(1:(T_follow - T_pc),
140-
function(x) {
141-
ci = boot.ci(boot.out, conf = 1 - alpha, type = bootstrap_CI_type, index = x)
142-
ci[[boot.ci.type]][5]
143-
})
144-
145-
sd.boot = sqrt(diag(var(boot.out$t)))
146-
147-
results = data.frame(
156+
157+
sd.boot <- sqrt(diag(var(boot.out$t)))
158+
159+
results <- data.frame(
148160
point_estimates = tau,
149161
lower_CI_boot = lower_CI_boot,
150162
upper_CI_boot = upper_CI_boot
151163
)
152-
return(results)
153-
154-
}else{
164+
results
165+
} else {
155166
stop("No other inference methods defined!")
156167
}
157-
158168
}
159-

0 commit comments

Comments
 (0)