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 .

Loading