|
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
|
2
|
model <- svyglm.obj
|
|
24
|
2
|
design.model <- model$survey.design
|
|
25
|
2
|
xs <- attr(model$terms, "term.labels")
|
|
26
|
2
|
y <- names(model$model)[1]
|
|
27
|
2
|
gaussianT <- ifelse(length(grep("gaussian", model$family)) == 1, T, F)
|
|
28
|
|
|
|
29
|
2
|
if (length(xs) == 0){
|
|
30
|
0
|
stop("No independent variable")
|
|
31
|
2
|
} else if (length(xs) ==1){
|
|
32
|
2
|
uni <- data.frame(coefNA(survey::svyglm(as.formula(paste(y, " ~ ", xs)), design = design.model, family = model$family)))[-1, ]
|
|
33
|
2
|
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
|
2
|
if (gaussianT){
|
|
36
|
2
|
summ <- paste(round(uni[,1], decimal), " (", round(uni[, 1] - 1.96*uni[, 2], decimal), "," ,round(uni[, 1] + 1.96*uni[, 2], decimal), ")", sep ="")
|
|
37
|
2
|
uni.res <- matrix(cbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))), nrow = nrow(uni))
|
|
38
|
2
|
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
|
2
|
rownames(uni.res) <-rownames(uni)
|
|
45
|
2
|
res <- uni.res
|
|
46
|
|
|
|
47
|
|
} else{
|
|
48
|
2
|
uni <- lapply(xs, function(v){
|
|
49
|
2
|
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
|
2
|
rn.uni <- lapply(uni, rownames)
|
|
55
|
2
|
uni <- Reduce(rbind, uni)
|
|
56
|
|
|
|
57
|
2
|
if (gaussianT){
|
|
58
|
2
|
summ <- paste(round(uni[,1], decimal), " (", round(uni[, 1] - 1.96*uni[, 2], decimal), "," ,round(uni[, 1] + 1.96*uni[, 2], decimal), ")", sep ="")
|
|
59
|
2
|
uni.res <- t(rbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))))
|
|
60
|
2
|
colnames(uni.res) <- c(paste("crude coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value")
|
|
61
|
2
|
rownames(uni.res) <-rownames(uni)
|
|
62
|
|
#mul <- summary(model)$coefficients[-1, ]
|
|
63
|
2
|
mul <- coefNA(model)[-1, ]
|
|
64
|
2
|
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
|
2
|
mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <=0.001, "< 0.001", as.character(round(mul[, 4], decimal +1)))))
|
|
66
|
2
|
colnames(mul.res) <- c(paste("adj. coeff.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value")
|
|
67
|
|
} else{
|
|
68
|
2
|
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
|
2
|
uni.res <- t(rbind(summ, ifelse(uni[, 4] <=0.001, "< 0.001", as.character(round(uni[, 4], decimal +1)))))
|
|
70
|
2
|
colnames(uni.res) <- c(paste("crude OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "crude P value")
|
|
71
|
2
|
rownames(uni.res) <-rownames(uni)
|
|
72
|
|
#mul <- summary(model)$coefficients[-1, ]
|
|
73
|
2
|
mul <- coefNA(model)[-1, ]
|
|
74
|
2
|
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
|
2
|
mul.res <- t(rbind(mul.summ, ifelse(mul[, 4] <=0.001, "< 0.001", as.character(round(mul[, 4], decimal +1)))))
|
|
76
|
2
|
colnames(mul.res) <- c(paste("adj. OR.(", 100 - 100 * 0.05, "%CI)", sep = ""), "adj. P value")
|
|
77
|
|
}
|
|
78
|
2
|
res <-cbind(uni.res[rownames(uni.res) %in% rownames(mul.res), ], mul.res)
|
|
79
|
2
|
rownames(res) <- rownames(mul)
|
|
80
|
|
}
|
|
81
|
|
|
|
82
|
2
|
fix.all <- res
|
|
83
|
2
|
mdata <- model$model
|
|
84
|
|
## rownames
|
|
85
|
|
#fix.all.list = lapply(xs, function(x){fix.all[grepl(x, rownames(fix.all)),]})
|
|
86
|
2
|
fix.all.list <- lapply(1:length(xs), function(x){fix.all[rownames(fix.all) %in% rn.uni[[x]], ]})
|
|
87
|
2
|
varnum.mfac <- which(lapply(fix.all.list, length) > ncol(fix.all))
|
|
88
|
2
|
lapply(varnum.mfac, function(x){fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])})
|
|
89
|
2
|
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
|
2
|
rn.list <- lapply(1:length(xs), function(x){rownames(fix.all)[rownames(fix.all) %in% rn.uni[[x]]]})
|
|
93
|
2
|
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
|
2
|
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
|
2
|
if (class(fix.all.unlist)[1] == "character"){
|
|
97
|
2
|
fix.all.unlist <- t(data.frame(fix.all.unlist))
|
|
98
|
|
}
|
|
99
|
2
|
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
|
2
|
outcome.name <- names(model$model)[1]
|
|
108
|
|
|
|
109
|
|
|
|
110
|
2
|
if (gaussianT) {
|
|
111
|
2
|
first.line <- paste("Linear regression predicting ", outcome.name, sep = "", "- weighted data\n")
|
|
112
|
2
|
last.lines <- paste("No. of observations = ",
|
|
113
|
2
|
length(model$y), "\n", "AIC value = ", round(model$aic,
|
|
114
|
2
|
decimal + 2), "\n", "\n", sep = "")
|
|
115
|
|
}
|
|
116
|
|
else {
|
|
117
|
2
|
first.line <- paste("Logistic regression predicting ", outcome.name, sep = "", "- weighted data\n")
|
|
118
|
2
|
last.lines <- paste("No. of observations = ", nrow(model$model),
|
|
119
|
2
|
"\n", "\n", sep = "")
|
|
120
|
|
}
|
|
121
|
|
|
|
122
|
2
|
results <- list(first.line = first.line, table = fix.all.unlist,
|
|
123
|
2
|
last.lines = last.lines)
|
|
124
|
2
|
class(results) <- c("display", "list")
|
|
125
|
2
|
return(results)
|
|
126
|
|
|
|
127
|
|
}
|
jinseob2kim