1 ```#' @title geeUni: The coefficient of univariate gee (internal function) ``` 2 ```#' @description Extract the coefficients of univariate gee using geeglm function (geepack package). ``` 3 ```#' @param y Dependant variable ``` 4 ```#' @param x Independent variable ``` 5 ```#' @param data Data ``` 6 ```#' @param id.vec Vector of id (should be ordered) ``` 7 ```#' @param family Family: "gaussian", "binomial", "poisson", "quasipoisson", etc... ``` 8 ```#' @param cor.type Correlation structure, Default: 'exchangeable' ``` 9 ```#' @return coefficient, standard error, p-value ``` 10 ```#' @details DETAILS ``` 11 ```#' @examples ``` 12 ```#' library(geepack) ``` 13 ```#' data(dietox) ``` 14 ```#' dietox\$Cu <- as.factor(dietox\$Cu) ``` 15 ```#' gee.uni <- geeUni("Weight", "Time", data = dietox, id.vec = dietox\$Pig, ``` 16 ```#' family = "gaussian", cor.type = "exchangeable") ``` 17 ```#' @rdname geeUni ``` 18 ```#' @importFrom geepack geeglm ``` 19 ```#' @importFrom stats as.formula ``` 20 ```#' @export ``` 21 22 23 ```geeUni = function(y, x, data, id.vec, family, cor.type = "exchangeable"){ ``` 24 1 ``` form <- as.formula(paste(y, "~", x)) ``` 25 1 ``` res <- geepack::geeglm(form, data = data, family = family, id = id.vec, corstr = cor.type) ``` 26 1 ``` coef <- summary(res)\$coefficients[-1, -3] ``` 27 1 ``` return(coef) ``` 28 ```} ``` 29 30 31 32 ```#' @title geeExp: transform the unit of coefficients (internal function) ``` 33 ```#' @description Transform the unit of coefficients to "Coeff", "OR" or "RR" ``` 34 ```#' @param gee.coef geeUni object. ``` 35 ```#' @param family Family: "gaussian", "binomial", "poisson", "quasipoisson", etc..., Default: 'binomial' ``` 36 ```#' @param dec Decimal point ``` 37 ```#' @return The transforemed coefficients(95% CI), p-value ``` 38 ```#' @details DETAILS ``` 39 ```#' @examples ``` 40 ```#' library(geepack) ``` 41 ```#' data(dietox) ``` 42 ```#' dietox\$Cu <- as.factor(dietox\$Cu) ``` 43 ```#' gee.uni <- geeUni("Weight", c("Time", "Cu"), data = dietox, id.vec = dietox\$Pig, ``` 44 ```#' family = "gaussian", cor.type = "exchangeable") ``` 45 ```#' gee.exp <- geeExp(gee.uni, "binomial", 2) ``` 46 ```#' @rdname geeExp ``` 47 ```#' @export ``` 48 49 50 ```geeExp = function(gee.coef, family ="binomial", dec){ ``` 51 1 ``` if (family == "binomial"){ ``` 52 1 ``` OR <- paste(round(exp(gee.coef[,1]), dec), " (", round(exp(gee.coef[,1] - 1.96*gee.coef[,2]), dec), ",", round(exp(gee.coef[,1] + 1.96*gee.coef[,2]), dec),")", sep="") ``` 53 1 ``` return(cbind(OR, gee.coef[,3])) ``` 54 1 ``` } else if (family == "gaussian"){ ``` 55 1 ``` coeff <- paste(round(gee.coef[,1], dec), " (", round(gee.coef[,1] - 1.96*gee.coef[,2], dec), ",", round(gee.coef[,1] + 1.96*gee.coef[,2], dec), ")", sep="") ``` 56 1 ``` return(cbind(coeff, gee.coef[,3])) ``` 57 0 ``` } else if (family %in% c("poisson", "quasipoisson")){ ``` 58 0 ``` RR <- paste(round(exp(gee.coef[,1]), dec), " (", round(exp(gee.coef[,1] - 1.96*gee.coef[,2]), dec), ",", round(exp(gee.coef[,1] + 1.96*gee.coef[,2]), dec),")", sep="") ``` 59 0 ``` return(cbind(RR, gee.coef[,3])) ``` 60 ``` } ``` 61 ```} ``` 62 63 64 65 ```#' @title geeglm.display ``` 66 ```#' @description Make gee results from "geeglm" object ``` 67 ```#' @param geeglm.obj "geeglm" object ``` 68 ```#' @param decimal Decimal, Default: 2 ``` 69 ```#' @return List: caption, main table, metrics table ``` 70 ```#' @details DETAILS ``` 71 ```#' @examples ``` 72 ```#' library(geepack) ``` 73 ```#' data(dietox) ``` 74 ```#' dietox\$Cu <- as.factor(dietox\$Cu) ``` 75 ```#' gee01 <- geeglm (Weight ~ Time + Cu , id =Pig, data = dietox, ``` 76 ```#' family=gaussian,corstr="ex") ``` 77 ```#' geeglm.display(gee01) ``` 78 ```#' @seealso ``` 79 ```#' \code{\link[data.table]{data.table-package}} ``` 80 ```#' \code{\link[stats]{complete.cases}} ``` 81 ```#' @rdname geeglm.display ``` 82 ```#' @export ``` 83 ```#' @importFrom data.table data.table ``` 84 ```#' @importFrom stats complete.cases ``` 85 86 ```geeglm.display = function(geeglm.obj, decimal = 2){ ``` 87 1 ``` family.gee <- geeglm.obj\$family[[1]] ``` 88 1 ``` corstr.gee <- geeglm.obj\$corstr ``` 89 1 ``` y <- as.character(geeglm.obj\$terms[[2]]) ``` 90 1 ``` xs <- names(geeglm.obj\$model)[-1] ``` 91 ``` ## rownames ``` 92 1 ``` geeglm.obj\$data <- data.table::data.table(geeglm.obj\$data) ``` 93 1 ``` nomiss <- stats::complete.cases(geeglm.obj\$data[, c(y, xs), with = F]) ``` 94 1 ``` gee.uni.list <- lapply(xs, function(x){geeUni(y, x, data = geeglm.obj\$data[nomiss, ], id.vec = geeglm.obj\$id, family = family.gee, cor.type = corstr.gee)}) ``` 95 1 ``` rn.uni <- lapply(gee.uni.list, rownames) ``` 96 1 ``` gee.uni <- Reduce(rbind, gee.uni.list) ``` 97 ``` ``` 98 1 ``` if (length(xs) ==1){ ``` 99 1 ``` gee.res <- geeExp(gee.uni, family = family.gee, dec = decimal) ``` 100 1 ``` gee.res.list <- lapply(1:length(xs), function(x){gee.res[rownames(gee.uni) %in% rn.uni[[x]], ]}) ``` 101 1 ``` varnum.mfac <- which(lapply(gee.res.list, length) > ncol(gee.res)) ``` 102 0 ``` lapply(varnum.mfac, function(x){gee.res.list[[x]] <<- rbind(rep(NA, ncol(gee.res)), gee.res.list[[x]])}) ``` 103 1 ``` gee.res.modi <- Reduce(rbind, gee.res.list) ``` 104 1 ``` if (nrow(gee.uni) == 1){ ``` 105 1 ``` gee.res.modi <- t(data.frame(gee.res.modi)) ``` 106 ``` } ``` 107 108 1 ``` family.label <- colnames(gee.res.modi)[1] ``` 109 1 ``` colnames(gee.res.modi) <- c(paste(family.label, "(95%CI)",sep = ""), "P value") ``` 110 ``` } else{ ``` 111 1 ``` gee.multi <- summary(geeglm.obj)\$coefficients[-1, -3] ``` 112 1 ``` gee.res <- cbind(geeExp(gee.uni, family = family.gee, dec = decimal), geeExp(gee.multi, family = family.gee, dec = decimal)) ``` 113 1 ``` gee.res.list <- lapply(1:length(xs), function(x){gee.res[rownames(gee.uni) %in% rn.uni[[x]], ]}) ``` 114 1 ``` varnum.mfac <- which(lapply(gee.res.list, length) > ncol(gee.res)) ``` 115 1 ``` lapply(varnum.mfac, function(x){gee.res.list[[x]] <<- rbind(rep(NA, ncol(gee.res)), gee.res.list[[x]])}) ``` 116 1 ``` gee.res.modi <- Reduce(rbind, gee.res.list) ``` 117 1 ``` family.label <- colnames(gee.res.modi)[1] ``` 118 1 ``` colnames(gee.res.modi) <- c(paste("crude ", family.label, "(95%CI)",sep = ""), "crude P value", paste("adj. ", family.label, "(95%CI)",sep = ""), "adj. P value") ``` 119 ``` } ``` 120 ``` ``` 121 ``` ``` 122 ``` ``` 123 1 ``` rn.list <- lapply(1:length(xs), function(x){rownames(gee.uni)[rownames(gee.uni) %in% rn.uni[[x]]]}) ``` 124 1 ``` varnum.2fac <- which(lapply(xs, function(x){length(geeglm.obj\$xlevels[[x]])}) == 2) ``` 125 0 ``` lapply(varnum.2fac, function(x){rn.list[[x]] <<- paste(xs[x], ": ", geeglm.obj\$xlevels[[xs[x]]][2], " vs ", geeglm.obj\$xlevels[[xs[x]]][1], sep="")}) ``` 126 1 ``` lapply(varnum.mfac, function(x){rn.list[[x]] <<- c(paste(xs[x],": ref.=", geeglm.obj\$xlevels[[xs[x]]][1], sep=""), gsub(xs[x]," ", rn.list[[x]]))}) ``` 127 ``` ``` 128 ``` ``` 129 1 ``` rownames(gee.res.modi) <- unlist(rn.list) ``` 130 ``` #out = as.data.frame(gee.res.modi) ``` 131 ``` #lapply(seq(2, ncol(gee.res), by =2), function(x){out[, x] <<- as.numeric(as.vector(out[, x]))}) ``` 132 1 ``` out <- gee.res.modi ``` 133 1 ``` pv.colnum <- which(colnames(out) %in% c("P value", "crude P value", "adj. P value")) ``` 134 1 ``` for (i in pv.colnum){ ``` 135 1 ``` out[, i] <- ifelse(as.numeric(out[, i]) < 0.001, "< 0.001", round(as.numeric(out[, i]), decimal + 1)) ``` 136 ``` } ``` 137 ``` ``` 138 ``` ## Metric ``` 139 1 ``` info.gee <- as.character(c(NA, round(as.numeric(summary(geeglm.obj)\$corr[1]), decimal+1), length(unique(geeglm.obj\$id)), length(geeglm.obj\$y))) ``` 140 1 ``` info.df <- cbind(info.gee, matrix(NA, 4, ncol(gee.res) -1)) ``` 141 1 ``` colnames(info.df) <- colnames(out) ``` 142 ``` #lapply(seq(2, ncol(gee.res), by =2), function(x){info.df[, x] <<- as.numeric(info.df[, x])}) ``` 143 1 ``` rownames(info.df) <- c("","Estimated correlation parameters", "No. of clusters", "No. of observations") ``` 144 ``` ``` 145 ``` ## Caption ``` 146 1 ``` cap.gee = paste("GEE(", family.gee, ") predicting ", y, " by ", paste(xs, collapse = ", "), " - Group ", as.character(geeglm.obj\$call\$id)[length(as.character(geeglm.obj\$call\$id))], sep="") ``` 147 ``` ``` 148 1 ``` return(list(caption = cap.gee, table = out, metric = info.df)) ``` 149 ```} ```

Read our documentation on viewing source code .

Loading