1

2
#' @title svyregress.display: table for svyglm.object
3
#' @description table for svyglm.object (survey package).
4
#' @param svyglm.obj svyglm.object
5
#' @param decimal digit, Default: 2
6
#' @return table
7
#' @details DETAILS
8
#' @examples 
9
#'  library(survey);data(api)
10
#'  apistrat$tt = c(rep(1, 20), rep(0, nrow(apistrat) -20))
11
#'  dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
12
#'  ds <- svyglm(api00~ell+meals+cname+mobility, design=dstrat)
13
#'  ds2 <- svyglm(tt~ell+meals+cname+mobility, design=dstrat, family = quasibinomial())
14
#'  svyregress.display(ds)
15
#'  svyregress.display(ds2)
16
#' @rdname svyglm.display
17
#' @export 
18
#' @importFrom survey svyglm
19

20

21
svyregress.display <- function(svyglm.obj, decimal = 2){
22
  
23 1
  model <- svyglm.obj
24 1
  design.model <- model$survey.design
25 1
  xs <- attr(model$terms, "term.labels")
26 1
  y <- names(model$model)[1]
27 1
  gaussianT <- ifelse(length(grep("gaussian", model$family)) == 1, T, F)
28
  
29 1
  if (length(xs) == 0){
30 0
    stop("No independent variable")
31 1
  } else if (length(xs) ==1){
32 1
    uni <- data.frame(coefNA(survey::svyglm(as.formula(paste(y, " ~ ", xs)), design = design.model, family = model$family)))[-1, ]
33 1
    rn.uni <- lapply(list(uni), rownames)
34
    #uni <- data.frame(summary(survey::svyglm(as.formula(paste(y, " ~ ", xs)), design = design.model, family = model$family))$coefficients)[-1, ]
35 1
    if (gaussianT){
36 1
      summ <- paste(round(uni[,1], decimal), " (", round(uni[, 1] - 1.96*uni[, 2], decimal), "," ,round(uni[, 1] + 1.96*uni[, 2], decimal), ")", sep ="")
37 1
      uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))), nrow = nrow(uni))
38 1
      colnames(uni.res) <- c(paste("Coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value")
39
    } else{
40 0
      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 ="")
41 0
      uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))), nrow = nrow(uni))
42 0
      colnames(uni.res) <- c(paste("OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "P value")
43
    }
44 1
    rownames(uni.res) <-rownames(uni)
45 1
    res <- uni.res
46
    
47
  } else{
48 1
    uni <- lapply(xs, function(v){
49 1
      data.frame(coefNA(survey::svyglm(as.formula(paste(y, " ~ ", v)), design = design.model, family = model$family)))[-1, ]
50
    })
51
    #uni <- lapply(xs, function(v){
52
    #  summary(survey::svyglm(as.formula(paste(y, " ~ ", v)), design = design.model))$coefficients[-1, ]
53
    #})
54 1
    rn.uni <- lapply(uni, rownames)
55 1
    uni <- Reduce(rbind, uni)
56
    
57 1
    if (gaussianT){
58 1
      summ <- paste(round(uni[,1], decimal), " (", round(uni[, 1] - 1.96*uni[, 2], decimal), "," ,round(uni[, 1] + 1.96*uni[, 2], decimal), ")", sep ="")
59 1
      uni.res <- t(rbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))))
60 1
      colnames(uni.res) <- c(paste("crude coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value")
61 1
      rownames(uni.res) <-rownames(uni)
62
      #mul <- summary(model)$coefficients[-1, ]
63 1
      mul <- coefNA(model)[-1, ]
64 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 ="")
65 1
      mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <=0.001, "< 0.001", as.character(round(mul[, 4], decimal +1)))))
66 1
      colnames(mul.res) <- c(paste("adj. coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value")
67
    } else{
68 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 ="")
69 1
      uni.res <- t(rbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))))
70 1
      colnames(uni.res) <- c(paste("crude OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value")
71 1
      rownames(uni.res) <-rownames(uni)
72
      #mul <- summary(model)$coefficients[-1, ]
73 1
      mul <- coefNA(model)[-1, ]
74 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 ="")
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. OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value")
77
    }
78 1
    res <-cbind(uni.res[rownames(uni.res) %in% rownames(mul.res), ], mul.res)
79 1
    rownames(res) <- rownames(mul)
80
  }
81
  
82 1
  fix.all <- res
83 1
  mdata <- model$model
84
  ## rownames
85
  #fix.all.list = lapply(xs, function(x){fix.all[grepl(x, rownames(fix.all)),]})     
86 1
  fix.all.list <- lapply(1:length(xs), function(x){fix.all[rownames(fix.all) %in% rn.uni[[x]], ]})
87 1
  varnum.mfac <- which(lapply(fix.all.list, length) > ncol(fix.all))
88 1
  lapply(varnum.mfac, function(x){fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])})
89 1
  fix.all.unlist <- Reduce(rbind, fix.all.list)
90
  
91
  #rn.list = lapply(xs, function(x){rownames(fix.all)[grepl(x, rownames(fix.all))]})
92 1
  rn.list <- lapply(1:length(xs), function(x){rownames(fix.all)[rownames(fix.all) %in% rn.uni[[x]]]})
93 1
  varnum.2fac <- which(xs %in% names(model$xlevels)[lapply(model$xlevels, length) == 2])
94 0
  lapply(varnum.2fac, function(x){rn.list[[x]] <<- paste(xs[x], ": ", model$xlevels[[xs[x]]][2], " vs ", model$xlevels[[xs[x]]][1], sep="")})
95 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]]))})
96 1
  if (class(fix.all.unlist)[1] == "character"){
97 1
    fix.all.unlist <- t(data.frame(fix.all.unlist))
98
  }
99 1
  rownames(fix.all.unlist) <- unlist(rn.list)
100
  
101
  #pv.colnum = which(colnames(fix.all.unlist) %in% c("P value", "crude P value", "adj. P value"))
102
  #for (i in pv.colnum){
103
  #  fix.all.unlist[, i] = ifelse(as.numeric(fix.all.unlist[, i]) < 0.001, "< 0.001", round(as.numeric(fix.all.unlist[, i]), decimal + 1))
104
  #}
105
  
106
  
107 1
  outcome.name <- names(model$model)[1]
108
  
109
  
110 1
  if (gaussianT) {
111 1
    first.line <- paste("Linear regression predicting ", outcome.name, sep = "", "- weighted data\n")
112 1
    last.lines <- paste("No. of observations = ",
113 1
                        length(model$y), "\n", "AIC value = ", round(model$aic,
114 1
                                                                     decimal + 2), "\n", "\n", sep = "")
115
  }
116
  else {
117 1
    first.line <- paste("Logistic regression predicting ", outcome.name, sep = "", "- weighted data\n")
118 1
    last.lines <- paste("No. of observations = ", nrow(model$model),
119 1
                        "\n", "\n", sep = "")
120
  }
121

122 1
  results <- list(first.line = first.line, table = fix.all.unlist,
123 1
                last.lines = last.lines)
124 1
  class(results) <- c("display", "list")
125 1
  return(results)
126

127
}

Read our documentation on viewing source code .

Loading