Skip to content

Commit

Permalink
Merge pull request #14 from cyk0315/master
Browse files Browse the repository at this point in the history
Fix `svyglm`
  • Loading branch information
jinseob2kim authored Jan 5, 2024
2 parents dc852db + e0bff68 commit 77ca86d
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 29 deletions.
4 changes: 4 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
# jstable 1.1.4

* Fix: confidence interval calculation in `svyglm` ( thanks for `cyk0315`)

# jstable 1.1.3

* Add `addOverall` options to `CreateTableOneJS` and `svyCreateTableOneJS` to add overall column.
Expand Down
61 changes: 32 additions & 29 deletions R/svyglm.R
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
#' @rdname svyglm.display
#' @export
#' @importFrom survey svyglm
#' @importFrom stats update
#' @importFrom stats update confint


svyregress.display <- function(svyglm.obj, decimal = 2) {
Expand All @@ -25,59 +25,62 @@ svyregress.display <- function(svyglm.obj, decimal = 2) {
xs <- attr(model$terms, "term.labels")
y <- names(model$model)[1]
gaussianT <- ifelse(length(grep("gaussian", model$family)) == 1, T, F)

if (length(xs) == 0) {
stop("No independent variable")
} else if (length(xs) == 1) {
uni <- data.frame(coefNA(model))[-1, ]
uni <- cbind(data.frame(coefNA(model))[-1, ], data.frame(stats::confint(model))[-1, ])
rn.uni <- lapply(list(uni), rownames)
# uni <- data.frame(summary(survey::svyglm(as.formula(paste(y, " ~ ", xs)), design = design.model, family = model$family))$coefficients)[-1, ]
if (gaussianT) {
summ <- paste(round(uni[, 1], decimal), " (", round(uni[, 1] - 1.96 * uni[, 2], decimal), ",", round(uni[, 1] + 1.96 * uni[, 2], decimal), ")", sep = "")
summ <- paste(round(uni[, 1], decimal), " (", round(uni[,5], decimal), ",", round(uni[,6], decimal), ")", sep = "")
uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <= 0.001, "< 0.001", as.character(round(uni[, 4], decimal + 1)))), nrow = nrow(uni))
colnames(uni.res) <- c(paste("Coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value")
} else {
summ <- paste(round(exp(uni[, 1]), decimal), " (", round(exp(uni[, 1] - 1.96 * uni[, 2]), decimal), ",", round(exp(uni[, 1] + 1.96 * uni[, 2]), decimal), ")", sep = "")
summ <- paste(round(exp(uni[, 1]), decimal), " (", round(exp(uni[,5]), decimal), ",", round(exp(uni[,6]), decimal), ")", sep = "")
uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <= 0.001, "< 0.001", as.character(round(uni[, 4], decimal + 1)))), nrow = nrow(uni))
colnames(uni.res) <- c(paste("OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value")
}
rownames(uni.res) <- rownames(uni)
res <- uni.res
} else {
uni <- lapply(xs, function(v) {
data.frame(coefNA(stats::update(model, formula(paste(paste(c(". ~ .", xs), collapse = " - "), " + ", v)), design = design.model)))[-1, ]
cbind(
data.frame(coefNA(stats::update(model, formula(paste(paste(c(". ~ .", xs), collapse = " - "), " + ", v)), design = design.model)))[-1, ], data.frame(stats::confint(stats::update(model, formula(paste(paste(c(". ~ .", xs), collapse = " - "), " + ", v)), design = design.model)))[-1, ]
)
})

# uni <- lapply(xs, function(v){
# summary(survey::svyglm(as.formula(paste(y, " ~ ", v)), design = design.model))$coefficients[-1, ]
# })
rn.uni <- lapply(uni, rownames)
uni <- Reduce(rbind, uni)

if (gaussianT) {
summ <- paste(round(uni[, 1], decimal), " (", round(uni[, 1] - 1.96 * uni[, 2], decimal), ",", round(uni[, 1] + 1.96 * uni[, 2], decimal), ")", sep = "")
summ <- paste(round(uni[, 1], decimal), " (", round(uni[,5], decimal), ",", round(uni[,6], decimal), ")", sep = "")
uni.res <- t(rbind(summ, ifelse(uni[, 4] <= 0.001, "< 0.001", as.character(round(uni[, 4], decimal + 1)))))
colnames(uni.res) <- c(paste("crude coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value")
rownames(uni.res) <- rownames(uni)
# mul <- summary(model)$coefficients[-1, ]
mul <- coefNA(model)[-1, ]
mul.summ <- paste(round(mul[, 1], decimal), " (", round(mul[, 1] - 1.96 * mul[, 2], decimal), ",", round(mul[, 1] + 1.96 * mul[, 2], decimal), ")", sep = "")
mul <- cbind(coefNA(model)[-1, ], stats::confint(model)[-1, ])
mul.summ <- paste(round(mul[, 1], decimal), " (", round(mul[,5], decimal), ",", round(mul[,6], decimal), ")", sep = "")
mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <= 0.001, "< 0.001", as.character(round(mul[, 4], decimal + 1)))))
colnames(mul.res) <- c(paste("adj. coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value")
} else {
summ <- paste(round(exp(uni[, 1]), decimal), " (", round(exp(uni[, 1] - 1.96 * uni[, 2]), decimal), ",", round(exp(uni[, 1] + 1.96 * uni[, 2]), decimal), ")", sep = "")
summ <- paste(round(exp(uni[, 1]), decimal), " (", round(exp(uni[,5]), decimal), ",", round(exp(uni[,6]), decimal), ")", sep = "")
uni.res <- t(rbind(summ, ifelse(uni[, 4] <= 0.001, "< 0.001", as.character(round(uni[, 4], decimal + 1)))))
colnames(uni.res) <- c(paste("crude OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value")
rownames(uni.res) <- rownames(uni)
# mul <- summary(model)$coefficients[-1, ]
mul <- coefNA(model)[-1, ]
mul.summ <- paste(round(exp(mul[, 1]), decimal), " (", round(exp(mul[, 1] - 1.96 * mul[, 2]), decimal), ",", round(exp(mul[, 1] + 1.96 * mul[, 2]), decimal), ")", sep = "")
mul <- cbind(coefNA(model)[-1, ], stats::confint(model)[-1, ])
mul.summ <- paste(round(exp(mul[, 1]), decimal), " (", round(exp(mul[,5]), decimal), ",", round(exp(mul[,6]), decimal), ")", sep = "")
mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <= 0.001, "< 0.001", as.character(round(mul[, 4], decimal + 1)))))
colnames(mul.res) <- c(paste("adj. OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value")
}
res <- cbind(uni.res[rownames(uni.res) %in% rownames(mul.res), ], mul.res)
rownames(res) <- rownames(mul)
}

fix.all <- res
mdata <- model$model
## rownames
Expand All @@ -90,7 +93,7 @@ svyregress.display <- function(svyglm.obj, decimal = 2) {
fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])
})
fix.all.unlist <- Reduce(rbind, fix.all.list)

# rn.list = lapply(xs, function(x){rownames(fix.all)[grepl(x, rownames(fix.all))]})
rn.list <- lapply(1:length(xs), function(x) {
rownames(fix.all)[rownames(fix.all) %in% rn.uni[[x]]]
Expand All @@ -106,37 +109,37 @@ svyregress.display <- function(svyglm.obj, decimal = 2) {
fix.all.unlist <- t(data.frame(fix.all.unlist))
}
rownames(fix.all.unlist) <- unlist(rn.list)

# pv.colnum = which(colnames(fix.all.unlist) %in% c("P value", "crude P value", "adj. P value"))
# for (i in pv.colnum){
# fix.all.unlist[, i] = ifelse(as.numeric(fix.all.unlist[, i]) < 0.001, "< 0.001", round(as.numeric(fix.all.unlist[, i]), decimal + 1))
# }


outcome.name <- names(model$model)[1]


if (gaussianT) {
first.line <- paste("Linear regression predicting ", outcome.name, sep = "", "- weighted data\n")
last.lines <- paste("No. of observations = ",
length(model$y), "\n", "AIC value = ", round(
model$aic,
decimal + 2
), "\n", "\n",
sep = ""
length(model$y), "\n", "AIC value = ", round(
model$aic,
decimal + 2
), "\n", "\n",
sep = ""
)
} else {
first.line <- paste("Logistic regression predicting ", outcome.name, sep = "", "- weighted data\n")
last.lines <- paste("No. of observations = ", nrow(model$model),
"\n", "\n",
sep = ""
"\n", "\n",
sep = ""
)
}

results <- list(
first.line = first.line, table = fix.all.unlist,
last.lines = last.lines
)
class(results) <- c("display", "list")
return(results)
}
}

0 comments on commit 77ca86d

Please sign in to comment.