1 2 ```#' @title coefNA: make coefficient table with NA ``` 3 ```#' @description Make coefficient table with NA ``` 4 ```#' @param model glm object (gaussian or binomial) ``` 5 ```#' @return coefficient table with NA ``` 6 ```#' @details DETAILS ``` 7 ```#' @examples ``` 8 ```#' coefNA(glm(mpg ~ wt + qsec, data = mtcars)) ``` 9 ```#' @rdname coefNA ``` 10 ```#' @export ``` 11 ```#' @importFrom stats coef ``` 12 13 ```coefNA <- function(model){ ``` 14 1 ``` coef.rownames <- merge(coef(summary(model)), model\$coefficients, by =0, all= T) ``` 15 1 ``` coef.matrix <- as.matrix(coef.rownames[, -c(1, ncol(coef.rownames))]) ``` 16 1 ``` rownames(coef.matrix) <- coef.rownames[, "Row.names"] ``` 17 1 ``` return(coef.matrix[names(model\$coefficients), ]) ``` 18 ```} ``` 19 20 21 ```#' @title glmshow.display: Show summary table of glm object. ``` 22 ```#' @description Show summary table of glm object(regression, logistic). ``` 23 ```#' @param glm.object glm.object ``` 24 ```#' @param decimal digits, Default: 2 ``` 25 ```#' @return table ``` 26 ```#' @details DETAILS ``` 27 ```#' @examples ``` 28 ```#' glmshow.display(glm(mpg ~ wt + qsec, data = mtcars)) ``` 29 ```#' @seealso ``` 30 ```#' \code{\link[stats]{glm}} ``` 31 ```#' @rdname glmshow.display ``` 32 ```#' @export ``` 33 ```#' @importFrom stats glm cor predict ``` 34 35 ```glmshow.display <- function (glm.object, decimal = 2){ ``` 36 1 ``` model <- glm.object ``` 37 0 ``` if(!any(class(model) %in% c("lm", "glm"))){stop("Model not from GLM")} ``` 38 ``` ``` 39 1 ``` xs <- attr(model\$terms, "term.labels") ``` 40 1 ``` y <- names(model\$model)[1] ``` 41 1 ``` gaussianT <- ifelse(length(grep("gaussian", model\$family)) == 1, T, F) ``` 42 1 ``` data <- model\$data ``` 43 ``` ``` 44 ``` ## table ``` 45 1 ``` if (length(xs) == 0){ ``` 46 0 ``` stop("No independent variable") ``` 47 1 ``` } else if (length(xs) ==1){ ``` 48 1 ``` uni <- data.frame(coefNA(stats::glm(as.formula(paste(y, " ~ ", xs)), data = data, family = model\$family)))[-1, ] ``` 49 1 ``` rn.uni <- lapply(list(uni), rownames) ``` 50 1 ``` if (gaussianT){ ``` 51 1 ``` summ <- paste(round(uni[,1], decimal), " (", round(uni[, 1] - 1.96*uni[, 2], decimal), "," ,round(uni[, 1] + 1.96*uni[, 2], decimal), ")", sep ="") ``` 52 1 ``` uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))), nrow = nrow(uni)) ``` 53 1 ``` colnames(uni.res) <- c(paste("Coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value") ``` 54 ``` } else{ ``` 55 1 ``` 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 ="") ``` 56 1 ``` uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))), nrow = nrow(uni)) ``` 57 1 ``` colnames(uni.res) <- c(paste("OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value") ``` 58 ``` } ``` 59 1 ``` rownames(uni.res) <-rownames(uni) ``` 60 1 ``` res <- uni.res ``` 61 ``` ``` 62 ``` } else{ ``` 63 1 ``` uni <- lapply(xs, function(v){ ``` 64 1 ``` data.frame(coefNA(stats::glm(as.formula(paste(y, " ~ ", v)), data = data, family = model\$family)))[-1, ] ``` 65 ``` }) ``` 66 1 ``` rn.uni <- lapply(uni, rownames) ``` 67 1 ``` uni <- Reduce(rbind, uni) ``` 68 1 ``` if (gaussianT){ ``` 69 1 ``` summ <- paste(round(uni[,1], decimal), " (", round(uni[, 1] - 1.96*uni[, 2], decimal), "," ,round(uni[, 1] + 1.96*uni[, 2], decimal), ")", sep ="") ``` 70 1 ``` uni.res <- t(rbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1))))) ``` 71 1 ``` colnames(uni.res) <- c(paste("crude coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value") ``` 72 1 ``` rownames(uni.res) <-rownames(uni) ``` 73 1 ``` mul <- coefNA(model)[-1, ] ``` 74 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 ="") ``` 75 1 ``` mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <=0.001, "< 0.001", as.character(round(mul[, 4], decimal +1))))) ``` 76 1 ``` colnames(mul.res) <- c(paste("adj. coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value") ``` 77 ``` } else{ ``` 78 1 ``` 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 ="") ``` 79 1 ``` uni.res <- t(rbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1))))) ``` 80 1 ``` colnames(uni.res) <- c(paste("crude OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value") ``` 81 1 ``` rownames(uni.res) <-rownames(uni) ``` 82 1 ``` mul <- coefNA(model)[-1, ] ``` 83 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 ="") ``` 84 1 ``` mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <=0.001, "< 0.001", as.character(round(mul[, 4], decimal +1))))) ``` 85 1 ``` colnames(mul.res) <- c(paste("adj. OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value") ``` 86 ``` } ``` 87 ``` ``` 88 1 ``` res <-cbind(uni.res[rownames(uni.res) %in% rownames(mul.res), ], mul.res) ``` 89 1 ``` rownames(res) <- rownames(mul) ``` 90 ``` } ``` 91 ``` ``` 92 ``` ## label ``` 93 1 ``` fix.all <- res ``` 94 ``` ``` 95 ``` ## rownames ``` 96 1 ``` fix.all.list <- lapply(1:length(xs), function(x){fix.all[rownames(fix.all) %in% rn.uni[[x]], ]}) ``` 97 1 ``` varnum.mfac <- which(lapply(fix.all.list, length) > ncol(fix.all)) ``` 98 1 ``` lapply(varnum.mfac, function(x){fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])}) ``` 99 1 ``` fix.all.unlist <- Reduce(rbind, fix.all.list) ``` 100 ``` ``` 101 1 ``` rn.list <- lapply(1:length(xs), function(x){rownames(fix.all)[rownames(fix.all) %in% rn.uni[[x]]]}) ``` 102 1 ``` varnum.2fac <- which(xs %in% names(model\$xlevels)[lapply(model\$xlevels, length) == 2]) ``` 103 0 ``` lapply(varnum.2fac, function(x){rn.list[[x]] <<- paste(xs[x], ": ", model\$xlevels[[xs[x]]][2], " vs ", model\$xlevels[[xs[x]]][1], sep="")}) ``` 104 1 ``` lapply(varnum.mfac, function(x){rn.list[[x]] <<- c(paste(xs[x],": ref.=", model\$xlevels[[xs[x]]][1], sep=""), gsub(xs[x]," ", rn.list[[x]]))}) ``` 105 1 ``` if (class(fix.all.unlist)[1] == "character"){ ``` 106 1 ``` fix.all.unlist <- t(data.frame(fix.all.unlist)) ``` 107 ``` } ``` 108 1 ``` rownames(fix.all.unlist) <- unlist(rn.list) ``` 109 ``` ``` 110 ``` #pv.colnum = which(colnames(fix.all.unlist) %in% c("P value", "crude P value", "adj. P value")) ``` 111 ``` #for (i in pv.colnum){ ``` 112 ``` # fix.all.unlist[, i] = ifelse(as.numeric(fix.all.unlist[, i]) < 0.001, "< 0.001", round(as.numeric(fix.all.unlist[, i]), decimal + 1)) ``` 113 ``` #} ``` 114 ``` ``` 115 ``` ``` 116 1 ``` outcome.name <- y ``` 117 ``` ``` 118 ``` ``` 119 1 ``` if (gaussianT) { ``` 120 1 ``` first.line <- paste("Linear regression predicting ", outcome.name, sep = "", "\n") ``` 121 1 ``` last.lines <- paste("No. of observations = ", ``` 122 1 ``` length(model\$y), "\n", "R-squared = ", round(cor(model\$y, predict(model))^2, decimal + 2), "\n", ``` 123 1 ``` "AIC value = ", round(model\$aic, decimal + 2), "\n", "\n", sep = "") ``` 124 ``` } ``` 125 ``` else { ``` 126 1 ``` first.line <- paste("Logistic regression predicting ", outcome.name, sep = "", "\n") ``` 127 1 ``` last.lines <- paste("No. of observations = ", ``` 128 1 ``` length(model\$y), "\n", ``` 129 1 ``` "AIC value = ", round(model\$aic, decimal + 2), "\n", "\n", sep = "") ``` 130 ``` } ``` 131 ``` ``` 132 1 ``` results <- list(first.line = first.line, table = fix.all.unlist, ``` 133 1 ``` last.lines = last.lines) ``` 134 1 ``` class(results) <- c("display", "list") ``` 135 1 ``` return(results) ``` 136 ``` ``` 137 ```} ``` 138

Read our documentation on viewing source code .