Skip to content

Commit a7d8f92

Browse files
clarkliminggithub-actions[bot]danielinteractive
authored
fix robin_lm issue (#101)
* fix issue * [skip style] [skip vbump] Restyle files * trigger action * trigger action * add additional tests * prevent vcovhc being used in non gaussian * Update NEWS.md Co-authored-by: Daniel Sabanes Bove <danielinteractive@users.noreply.github.com> --------- Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Daniel Sabanes Bove <danielinteractive@users.noreply.github.com>
1 parent ee428a8 commit a7d8f92

File tree

7 files changed

+134
-24
lines changed

7 files changed

+134
-24
lines changed

‎NEWS.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@
99
### Bug Fixes
1010

1111
* Fixed a bug in covariate-adjusted stratified survival function estimation in `robin_surv` which could occur when there are character covariates with values only appearing in one stratum, which could have failed or lead to incorrect results.
12+
* Fixed a issue in `robin_lm` that variance method does not apply correctly.
13+
* Fixed a issue in `robin_glm` that `vcovHC` could previously be used for non-Gaussian family.
1214

1315
* Fixed another bug in covariate-adjusted stratified survival function estimation in `robin_surv`, which resulted from design matrices separately derived per stratum. Now the design matrix is created once including the stratum indicator, and then the stratum-specific parts are extracted as needed.
1416

‎R/robin_glm.R‎

Lines changed: 20 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -24,30 +24,36 @@
2424
#' treatment = treatment ~ s1, contrast = "difference"
2525
#' )
2626
robin_glm <- function(
27-
formula,
28-
data,
29-
treatment,
30-
contrast = c("difference", "risk_ratio", "odds_ratio", "log_risk_ratio", "log_odds_ratio"),
31-
contrast_jac = NULL,
32-
vcov = "vcovG",
33-
family = gaussian(),
34-
vcov_args = list(),
35-
pair,
36-
...) {
27+
formula,
28+
data,
29+
treatment,
30+
contrast = c("difference", "risk_ratio", "odds_ratio", "log_risk_ratio", "log_odds_ratio"),
31+
contrast_jac = NULL,
32+
vcov = "vcovG",
33+
family = gaussian(),
34+
vcov_args = list(),
35+
pair,
36+
...
37+
) {
3738
attr(formula, ".Environment") <- environment()
38-
# check if using negative.binomial family with NA as theta.
39-
# If so, use MASS::glm.nb instead of glm.
39+
if (is.function(family)) {
40+
family <- family()
41+
}
4042
assert_subset(all.vars(formula), names(data))
4143
assert_subset(all.vars(treatment), names(data))
4244
has_interaction <- h_interaction(formula, treatment)
4345
use_vcovhc <- identical(vcov, "vcovHC") || identical(vcov, vcovHC)
44-
if (use_vcovhc && (has_interaction || !identical(contrast, "difference"))) {
46+
if (test_character(contrast)) {
47+
contrast <- match.arg(contrast)
48+
}
49+
if (use_vcovhc && (has_interaction || !identical(contrast, "difference") || !identical(family$family, "gaussian"))) {
4550
stop(
4651
"Huber-White variance estimator is ONLY supported when the expected outcome difference is estimated",
4752
"using a linear model without treatment-covariate interactions; see the 2023 FDA guidance."
4853
)
4954
}
50-
55+
# check if using negative.binomial family with NA as theta.
56+
# If so, use MASS::glm.nb instead of glm.
5157
if (identical(family$family, "Negative Binomial(NA)")) {
5258
fit <- MASS::glm.nb(formula, data = data, ...)
5359
} else {
@@ -59,7 +65,6 @@ robin_glm <- function(
5965
}
6066

6167
if (test_character(contrast)) {
62-
contrast <- match.arg(contrast)
6368
if (contrast %in% c("odds_ratio", "risk_ratio")) {
6469
warning(
6570
c(

‎R/robin_lm.R‎

Lines changed: 11 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -17,18 +17,17 @@
1717
#' treatment = treatment ~ s1
1818
#' )
1919
robin_lm <- function(
20-
formula,
21-
data,
22-
treatment,
23-
vcov = "vcovG",
24-
vcov_args = list(),
25-
pair,
26-
...) {
20+
formula,
21+
data,
22+
treatment,
23+
vcov = "vcovG",
24+
vcov_args = list(),
25+
pair,
26+
...
27+
) {
2728
attr(formula, ".Environment") <- environment()
2829
assert_subset(all.vars(formula), names(data))
2930
assert_subset(all.vars(treatment), names(data))
30-
fit <- lm(formula, data = data, ...)
31-
pc <- predict_counterfactual(fit, treatment, data, variance = vcov, vcov_args = vcov_args)
3231
has_interaction <- h_interaction(formula, treatment)
3332
use_vcovhc <- identical(vcov, "vcovHC") || identical(vcov, vcovHC)
3433
if (use_vcovhc && has_interaction) {
@@ -37,6 +36,9 @@ robin_lm <- function(
3736
without treatment-covariate interactions; see the 2023 FDA guidance."
3837
)
3938
}
39+
fit <- lm(formula, data = data, ...)
40+
# same as robin_glm, to allow vcov as function to be used, also allow the function name to be captured.
41+
pc <- eval(bquote(predict_counterfactual(fit, treatment, data, vcov = .(substitute(vcov)), vcov_args = vcov_args)))
4042
if (missing(pair)) {
4143
pair <- pairwise(names(pc$estimate))
4244
}

‎tests/testthat/_snaps/robin_glm.md‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,3 +42,25 @@
4242
---
4343
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
4444

45+
---
46+
47+
Code
48+
robin_glm(y_b ~ treatment + s1, data = glm_data, treatment = treatment ~ s1)
49+
Output
50+
Model : y_b ~ treatment + s1
51+
Randomization: treatment ~ s1 ( Simple )
52+
Variance Type: vcovG
53+
Marginal Mean:
54+
Estimate Std.Err 2.5 % 97.5 %
55+
pbo 0.366010 0.033864 0.299638 0.4324
56+
trt1 0.580976 0.035027 0.512324 0.6496
57+
trt2 0.610163 0.034473 0.542598 0.6777
58+
59+
Contrast : difference
60+
Estimate Std.Err Z Value Pr(>|z|)
61+
trt1 v.s. pbo 0.214966 0.048633 4.4202 9.862e-06 ***
62+
trt2 v.s. pbo 0.244153 0.048244 5.0608 4.175e-07 ***
63+
trt2 v.s. trt1 0.029187 0.049064 0.5949 0.5519
64+
---
65+
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
66+

‎tests/testthat/_snaps/robin_lm.md‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
# robin_lm works correctly
2+
3+
Code
4+
robin_lm(y_b ~ treatment + s1, data = glm_data, treatment = treatment ~ s1)
5+
Output
6+
Model : y_b ~ treatment + s1
7+
Randomization: treatment ~ s1 ( Simple )
8+
Variance Type: vcovG
9+
Marginal Mean:
10+
Estimate Std.Err 2.5 % 97.5 %
11+
pbo 0.366010 0.033864 0.299638 0.4324
12+
trt1 0.580976 0.035027 0.512324 0.6496
13+
trt2 0.610163 0.034473 0.542598 0.6777
14+
15+
Contrast : h_diff
16+
Estimate Std.Err Z Value Pr(>|z|)
17+
trt1 v.s. pbo 0.214966 0.048633 4.4202 9.862e-06 ***
18+
trt2 v.s. pbo 0.244153 0.048244 5.0608 4.175e-07 ***
19+
trt2 v.s. trt1 0.029187 0.049064 0.5949 0.5519
20+
---
21+
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
22+

‎tests/testthat/test-robin_glm.R‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,4 +67,7 @@ test_that("robin_glm works for glm.nb", {
6767
test_that("robin_glm can be printed correctly", {
6868
expect_snapshot(robin_res1)
6969
expect_snapshot(robin_res2)
70+
expect_snapshot(
71+
robin_glm(y_b ~ treatment + s1, data = glm_data, treatment = treatment ~ s1)
72+
)
7073
})

‎tests/testthat/test-robin_lm.R‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,4 +27,58 @@ test_that("robin_lm works correctly", {
2727
"Huber-White variance estimator is ONLY supported when using a linear model
2828
without treatment-covariate interactions; see the 2023 FDA guidance."
2929
)
30+
expect_snapshot(
31+
robin_lm(y_b ~ treatment + s1, data = glm_data, treatment = treatment ~ s1)
32+
)
33+
})
34+
35+
36+
test_that("robin_lm give same result as robin_glm", {
37+
expect_silent(
38+
f1 <- robin_lm(
39+
y ~ treatment * s1,
40+
data = glm_data,
41+
treatment = treatment ~ pb(s1),
42+
)
43+
)
44+
expect_silent(
45+
f2 <- robin_glm(
46+
y ~ treatment * s1,
47+
data = glm_data,
48+
treatment = treatment ~ pb(s1),
49+
)
50+
)
51+
expect_equal(
52+
f1$marginal_mean$estimate,
53+
f2$marginal_mean$estimate
54+
)
55+
expect_equal(
56+
f1$contrast$estimate,
57+
f2$contrast$estimate
58+
)
59+
60+
expect_silent(
61+
f1 <- robin_lm(
62+
y ~ treatment + s1,
63+
data = glm_data,
64+
treatment = treatment ~ pb(s1),
65+
vcov = vcovHC
66+
)
67+
)
68+
expect_silent(
69+
f2 <- robin_glm(
70+
y ~ treatment + s1,
71+
data = glm_data,
72+
treatment = treatment ~ pb(s1),
73+
vcov = vcovHC
74+
)
75+
)
76+
expect_equal(
77+
f1$marginal_mean$estimate,
78+
f2$marginal_mean$estimate
79+
)
80+
expect_equal(
81+
f1$contrast$estimate,
82+
f2$contrast$estimate
83+
)
3084
})

0 commit comments

Comments
 (0)