Edits round 1 #12
There are no files selected for viewing
Binary file not shown.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -7,20 +7,33 @@ | |
| # (Miettinen/Schouten approach to directly estimating relative risks) | ||
| 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 | ||
| data <- data |> | ||
|
There was a problem hiding this comment. see earlier comment about base R pipe |
||
| dplyr::mutate(.clusterid = dplyr::row_number()) |> | ||
| dplyr::rename(outc = dplyr::one_of(!!yvar)) |> | ||
| tidyr::uncount(outc + 1) |> | ||
|
There was a problem hiding this comment. The see the elegance in this approach! However, in the spirit of my comment in #10, it would be nice not to rely on tidyverse functions down the road. Let us focus on testing and bug fixes for now. |
||
| dplyr::group_by(.clusterid) |> | ||
| dplyr::mutate(outc = floor(outc / dplyr::row_number())) |> | ||
| dplyr::ungroup() |> | ||
| dplyr::rename(!!yvar := "outc") | ||
|
|
||
| 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" | ||
| ) | ||
|
|
||
| return(fit) | ||
| } | ||
|
|
||
|
|
@@ -46,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), | ||
| "%" | ||
| ) | ||
| 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 | ||
| ) | ||
| ) | ||
| )) | ||
| ci[] <- cf[parm] + ses %o% fac | ||
| return(ci) | ||
| } | ||
|
|
@@ -77,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, ... | ||
| ) { | ||
| # 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")) | ||
| 1 | ||
| else if (df.r > 0) { | ||
| est.disp <- TRUE | ||
|
|
@@ -111,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) | ||
| ) | ||
| ###### end changes ######## | ||
| covmat <- dispersion * covmat.unscaled | ||
| var.cf <- diag(covmat) | ||
|
|
@@ -123,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|)") | ||
| ) | ||
| } | ||
| 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|)") | ||
| ) | ||
| } | ||
| else { | ||
| coef.table <- cbind(coef.p, NaN, NaN, NaN) | ||
| dimnames(coef.table) <- list( | ||
| names(coef.p), c(dn, "t value", "Pr(>|t|)") | ||
| ) | ||
| } | ||
| df.f <- NCOL(Qr$qr) | ||
| } | ||
| else { | ||
| coef.table <- matrix(data = NA, 0L, 4L) | ||
| dimnames(coef.table) <- list(NULL, c("Estimate", "Std. Error", "t value", "Pr(>|t|)")) | ||
| covmat.unscaled <- covmat <- matrix(data = NA, 0L, 0L) | ||
|
There was a problem hiding this comment. Are you certain that |
||
| df.f <- length(aliased) | ||
| } | ||
| keep <- match(c("call", "terms", "family", "deviance", "aic", | ||
|
|
||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
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