Skip to content
Open
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 48 additions & 28 deletions R/caseduplicate.R
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
estimate_duplicate <- function(formula, data, ...) {
yvar <- as.character(all.vars(formula)[1])

#TODO the following assumes that outc is coded as 0/1. We should throw an
# error when this is not true
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

100%! How about here:

risks/R/estimate_risk.R

Lines 216 to 217 in f48f0a0

...) {
implausible <- 0.99999

data <- data |>
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see earlier comment about base R pipe

dplyr::mutate(.clusterid = dplyr::row_number()) |>
dplyr::rename(outc = dplyr::one_of(!!yvar)) |>
Expand All @@ -17,12 +19,21 @@ estimate_duplicate <- function(formula, data, ...) {
dplyr::ungroup() |>
dplyr::rename(!!yvar := "outc")

fit <- eval(substitute(stats::glm(formula = formula,
family = binomial(link = "logit"),
data = data)))
fit <- eval(substitute(
stats::glm(
formula = formula,
family = binomial(link = "logit"),
data = data
)
))
class(fit) <- c("duplicate", class(fit))
fit <- estimate_maxprob(fit = fit, formula = formula, data = data,
link = "logit")
fit <- estimate_maxprob(
fit = fit,
formula = formula,
data = data,
link = "logit"
)

return(fit)
}

Expand All @@ -48,19 +59,24 @@ confint.duplicate <- function(object, parm = NULL, level = 0.95, ...) {
# parm <- pnames[parm]
a <- (1 - level)/2
a <- c(a, 1 - a)
pct <- paste(format(100 * a, trim = TRUE, scientific = FALSE, digits = 3),
"%")
pct <- paste(
format(100 * a, trim = TRUE, scientific = FALSE, digits = 3),
"%"
)
fac <- qnorm(a)
ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, pct))
# Robust covariance accounting for clustering
obj_sandwich <- object
class(obj_sandwich) <- "glm"
ses <- sqrt(diag(sandwich::sandwich(
x = obj_sandwich,
bread. = sandwich::bread(obj_sandwich),
meat. = sandwich::meatCL(obj_sandwich,
type = "HC0",
cluster = object$data$.clusterid))))
ses <- sqrt(diag(
sandwich::sandwich(
x = obj_sandwich,
bread. = sandwich::bread(obj_sandwich),
meat. = sandwich::meatCL(
obj_sandwich, type = "HC0", cluster = object$data$.clusterid
)
)
))
ci[] <- cf[parm] + ses %o% fac
return(ci)
}
Expand All @@ -79,13 +95,14 @@ confint.duplicate <- function(object, parm = NULL, level = 0.95, ...) {
#' @param symbolic.cor Not used
#' @param ... Other arguments, not used
#' @export
summary.duplicate <- function(object, dispersion = NULL, correlation = FALSE,
symbolic.cor = FALSE, ...) {
summary.duplicate <- function(
object, dispersion = NULL, correlation = FALSE, symbolic.cor = FALSE, ...
) {
# a modification of summary.glm()
est.disp <- FALSE
df.r <- object$df.residual
if (is.null(dispersion))
dispersion <- if (object$family$family %in% c("poisson", "binomial"))
if(is.null(dispersion))
dispersion <- if(object$family$family %in% c("poisson", "binomial"))
1
else if (df.r > 0) {
est.disp <- TRUE
Expand Down Expand Up @@ -113,9 +130,9 @@ summary.duplicate <- function(object, dispersion = NULL, correlation = FALSE,
covmat.unscaled <- sandwich::sandwich(
x = obj_sandwich,
bread. = sandwich::bread(x = obj_sandwich),
meat. = sandwich::meatCL(x = obj_sandwich,
type = "HC0",
cluster = object$data$.clusterid))
meat. = sandwich::meatCL(
x = obj_sandwich, type = "HC0", cluster = object$data$.clusterid)
)
###### end changes ########
covmat <- dispersion * covmat.unscaled
var.cf <- diag(covmat)
Expand All @@ -125,26 +142,29 @@ summary.duplicate <- function(object, dispersion = NULL, correlation = FALSE,
if (!est.disp) {
pvalue <- 2 * pnorm(-abs(tvalue))
coef.table <- cbind(coef.p, s.err, tvalue, pvalue)
dimnames(coef.table) <- list(names(coef.p), c(dn,
"z value", "Pr(>|z|)"))
dimnames(coef.table) <- list(
names(coef.p), c(dn, "z value", "Pr(>|z|)")
)
}
else if (df.r > 0) {
pvalue <- 2 * pt(-abs(tvalue), df.r)
coef.table <- cbind(coef.p, s.err, tvalue, pvalue)
dimnames(coef.table) <- list(names(coef.p), c(dn,
"t value", "Pr(>|t|)"))
dimnames(coef.table) <- list(
names(coef.p), c(dn, "t value", "Pr(>|t|)")
)
}
else {
coef.table <- cbind(coef.p, NaN, NaN, NaN)
dimnames(coef.table) <- list(names(coef.p), c(dn,
"t value", "Pr(>|t|)"))
dimnames(coef.table) <- list(
names(coef.p), c(dn, "t value", "Pr(>|t|)")
)
}
df.f <- NCOL(Qr$qr)
}
else {
coef.table <- matrix(, 0L, 4L)
coef.table <- matrix(data = NA, 0L, 4L)
dimnames(coef.table) <- list(NULL, c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
covmat.unscaled <- covmat <- matrix(, 0L, 0L)
covmat.unscaled <- covmat <- matrix(data = NA, 0L, 0L)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you certain that data = NA is exactly the same as leaving this space blank? This is verbatim code from stats::glm(), which we know works. I see the aesthetics point, but why change a winning horse?

df.f <- length(aliased)
}
keep <- match(c("call", "terms", "family", "deviance", "aic",
Expand Down