Skip to content

Commit 6d70fb9

Browse files
strengejackestrengejacke
authored
Discrepancies in reporting for models of class svycoxph (#1189)
* Discrepancies in reporting for models of class svycoxph Fixes #1174 * add test * desc, news * fix * fix * fix
1 parent 62ba1f8 commit 6d70fb9

6 files changed

Lines changed: 102 additions & 47 deletions

File tree

‎DESCRIPTION‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
Type: Package
22
Package: parameters
33
Title: Processing of Model Parameters
4-
Version: 0.28.3.1
4+
Version: 0.28.3.2
55
Authors@R:
66
c(person(given = "Daniel",
77
family = "Lüdecke",

‎NAMESPACE‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ S3method(ci,clm2)
3030
S3method(ci,clmm2)
3131
S3method(ci,coeftest)
3232
S3method(ci,complmrob)
33+
S3method(ci,coxph)
3334
S3method(ci,crq)
3435
S3method(ci,default)
3536
S3method(ci,deltaMethod)

‎NEWS.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,11 @@
44

55
* `model_parameters()` now supports objects from the *lavaan.mi* package.
66

7+
## Bug fixes
8+
9+
* Fixed issue where wrong (non-robust) standard errors were calculated for
10+
`coxph` and `svycoxph` objects.
11+
712
# parameters 0.28.3
813

914
* fixed bug in `standardize_info(<fixest>)` that was preventing

‎R/methods_survival.R‎

Lines changed: 36 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -3,11 +3,13 @@
33
#################### .survfit ------
44

55
#' @export
6-
model_parameters.survfit <- function(model,
7-
keep = NULL,
8-
drop = NULL,
9-
verbose = TRUE,
10-
...) {
6+
model_parameters.survfit <- function(
7+
model,
8+
keep = NULL,
9+
drop = NULL,
10+
verbose = TRUE,
11+
...
12+
) {
1113
s <- summary(model)
1214
# extract all elements with same length, which occur most in that list
1315
# that is the data we need
@@ -30,8 +32,14 @@ model_parameters.survfit <- function(model,
3032
params <- datawizard::data_rename(
3133
params,
3234
select = c(
33-
Time = "time", `N Risk` = "n.risk", `N Event` = "n.event", Survival = "surv",
34-
SE = "std.err", Group = "strata", CI_low = "lower", CI_high = "upper"
35+
Time = "time",
36+
`N Risk` = "n.risk",
37+
`N Event` = "n.event",
38+
Survival = "surv",
39+
SE = "std.err",
40+
Group = "strata",
41+
CI_low = "lower",
42+
CI_high = "upper"
3543
)
3644
)
3745

@@ -52,7 +60,6 @@ model_parameters.survfit <- function(model,
5260

5361
#################### .coxph ------
5462

55-
5663
#' @export
5764
standard_error.coxph <- function(model, method = NULL, ...) {
5865
robust <- !is.null(method) && method == "robust"
@@ -61,18 +68,29 @@ standard_error.coxph <- function(model, method = NULL, ...) {
6168
}
6269

6370
params <- insight::get_parameters(model)
64-
cs <- stats::coef(summary(model))
65-
se <- cs[, 3]
71+
junk <- utils::capture.output({
72+
s <- summary(model)
73+
})
74+
cs <- stats::coef(s)
75+
if (isTRUE(s$used.robust)) {
76+
se <- cs[, 4]
77+
} else {
78+
se <- cs[, 3]
79+
}
6680

6781
# check
6882
if (length(se) > nrow(params)) {
6983
se <- se[match(params$Parameter, .remove_backticks_from_string(rownames(cs)))]
7084
}
7185

72-
.data_frame(
73-
Parameter = params$Parameter,
74-
SE = as.vector(se)
75-
)
86+
.data_frame(Parameter = params$Parameter, SE = as.vector(se))
87+
}
88+
89+
90+
#' @export
91+
ci.coxph <- function(x, ...) {
92+
junk <- utils::capture.output(out <- ci.default(x, ...))
93+
out
7694
}
7795

7896

@@ -98,16 +116,12 @@ p_value.coxph <- function(model, ...) {
98116

99117
#################### .aareg ------
100118

101-
102119
#' @export
103120
standard_error.aareg <- function(model, ...) {
104121
s <- summary(model)
105122
se <- s$table[, "se(coef)"]
106123

107-
.data_frame(
108-
Parameter = .remove_backticks_from_string(names(se)),
109-
SE = as.vector(se)
110-
)
124+
.data_frame(Parameter = .remove_backticks_from_string(names(se)), SE = as.vector(se))
111125
}
112126

113127

@@ -116,16 +130,12 @@ p_value.aareg <- function(model, ...) {
116130
s <- summary(model)
117131
p <- s$table[, "p"]
118132

119-
.data_frame(
120-
Parameter = .remove_backticks_from_string(names(p)),
121-
p = as.vector(p)
122-
)
133+
.data_frame(Parameter = .remove_backticks_from_string(names(p)), p = as.vector(p))
123134
}
124135

125136

126137
#################### .survreg ------
127138

128-
129139
#' @export
130140
standard_error.survreg <- function(model, method = NULL, ...) {
131141
robust <- !is.null(method) && method == "robust"
@@ -136,10 +146,7 @@ standard_error.survreg <- function(model, method = NULL, ...) {
136146
s <- summary(model)
137147
se <- s$table[, 2]
138148

139-
.data_frame(
140-
Parameter = .remove_backticks_from_string(names(se)),
141-
SE = as.vector(se)
142-
)
149+
.data_frame(Parameter = .remove_backticks_from_string(names(se)), SE = as.vector(se))
143150
}
144151

145152

@@ -151,16 +158,12 @@ p_value.survreg <- function(model, method = NULL, ...) {
151158
}
152159
s <- summary(model)
153160
p <- s$table[, "p"]
154-
.data_frame(
155-
Parameter = .remove_backticks_from_string(names(p)),
156-
p = as.vector(p)
157-
)
161+
.data_frame(Parameter = .remove_backticks_from_string(names(p)), p = as.vector(p))
158162
}
159163

160164

161165
#################### .riskRegression ------
162166

163-
164167
#' @export
165168
standard_error.riskRegression <- function(model, ...) {
166169
junk <- utils::capture.output(cs <- stats::coef(model))

‎tests/testthat/test-coxph.R‎

Lines changed: 56 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -7,11 +7,7 @@ lung$ph.ecog <- factor(lung$ph.ecog, labels = c("good", "ok", "limited"))
77
m1 <- survival::coxph(survival::Surv(time, status) ~ sex + age + ph.ecog, data = lung)
88

99
test_that("ci", {
10-
expect_equal(
11-
ci(m1)$CI_low,
12-
c(-0.87535, -0.00747, 0.01862, 0.45527),
13-
tolerance = 1e-4
14-
)
10+
expect_equal(ci(m1)$CI_low, c(-0.87535, -0.00747, 0.01862, 0.45527), tolerance = 1e-4)
1511
})
1612

1713
test_that("se", {
@@ -23,11 +19,7 @@ test_that("se", {
2319
})
2420

2521
test_that("p_value", {
26-
expect_equal(
27-
p_value(m1)$p,
28-
c(0.00118, 0.24713, 0.04005, 8e-05),
29-
tolerance = 1e-4
30-
)
22+
expect_equal(p_value(m1)$p, c(0.00118, 0.24713, 0.04005, 8e-05), tolerance = 1e-4)
3123
})
3224

3325
test_that("model_parameters", {
@@ -79,3 +71,57 @@ withr::with_package(
7971
expect_snapshot(print(model_parameters(mod)))
8072
})
8173
)
74+
75+
test_that("model_parameters coxph, correct robust SE", {
76+
skip_if_not_installed("survey")
77+
78+
# Importing data
79+
d <- survival::pbc
80+
81+
# Defining randomization
82+
d$randomized <- with(d, !is.na(trt) & trt > 0)
83+
84+
# Defining randomization probability weights
85+
d$randprob <- fitted(glm(randomized ~ age * edema, data = d, family = binomial))
86+
87+
# Generating survey design object
88+
d.svy <- survey::svydesign(
89+
id = ~1,
90+
prob = ~randprob,
91+
strata = ~edema,
92+
data = subset(d, randomized)
93+
)
94+
95+
# The design-based cox proportional hazards model
96+
m <- survey::svycoxph(
97+
survival::Surv(time, status > 0) ~ log(bili) + protime + albumin,
98+
design = d.svy
99+
)
100+
101+
s <- summary(m)
102+
out <- model_parameters(m)
103+
104+
expect_equal(
105+
out$Coefficient,
106+
s$coefficients[, "coef"],
107+
tolerance = 1e-4,
108+
ignore_attr = TRUE
109+
)
110+
111+
expect_equal(
112+
out$SE,
113+
s$coefficients[, "robust se"],
114+
tolerance = 1e-4,
115+
ignore_attr = TRUE
116+
)
117+
118+
out <- model_parameters(m, digits = 5, p_digits = 5, exponentiate = TRUE)
119+
expect_equal(
120+
out$CI_low,
121+
s$conf.int[, "lower .95"],
122+
tolerance = 1e-4,
123+
ignore_attr = TRUE
124+
)
125+
unloadNamespace("survey")
126+
unloadNamespace("survival")
127+
})

‎tests/testthat/test-lavaan.mi.R‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ test_that("model_parameters.lavaan.mi", {
1818
fit.mi <- lavaan.mi::sem.mi(mod, data = imp)
1919

2020
mp <- model_parameters(fit.mi)
21-
expect_s3_class(mp, "parameters_model")
22-
expect_equal(nrow(mp), 5)
23-
expect_equal(mp$Coefficient[1], -0.05, tolerance = 1e-2)
21+
expect_s3_class(mp, "parameters_sem")
22+
expect_equal(nrow(mp), 4)
23+
expect_equal(mp$Coefficient[1], -0.05933624, tolerance = 1e-2)
2424
})

0 commit comments

Comments
 (0)