|
1
|
|
#' @title lmerExp: transform the unit of coefficients (internal function)
|
|
2
|
|
#' @description Transform the unit of coefficients to "Coeff", "OR" or "RR"
|
|
3
|
|
#' @param lmer.coef lmer coefficients.
|
|
4
|
|
#' @param family Family: "gaussian", "binomial", "poisson", "quasipoisson", etc..., Default: 'binomial'
|
|
5
|
|
#' @param dec Decimal point
|
|
6
|
|
#' @return The transforemed coefficients(95% CI), p-value
|
|
7
|
|
#' @details DETAILS
|
|
8
|
|
#' @examples
|
|
9
|
|
#' #EXAMPLE1
|
|
10
|
|
#' @rdname lmerExp
|
|
11
|
|
#' @importFrom stats pnorm
|
|
12
|
|
|
|
13
|
|
|
|
14
|
|
lmerExp = function(lmer.coef, family ="binomial", dec){
|
|
15
|
1
|
if (class(lmer.coef)[1] == "numeric"){
|
|
16
|
1
|
lmer.coef = t(data.frame(lmer.coef))
|
|
17
|
|
}
|
|
18
|
1
|
pv = 2*(1-pnorm(abs(lmer.coef[,3])))
|
|
19
|
1
|
if (family == "binomial"){
|
|
20
|
1
|
OR = paste(round(exp(lmer.coef[,1]), dec), " (", round(exp(lmer.coef[,1] - 1.96*lmer.coef[,2]), dec), ",", round(exp(lmer.coef[,1] + 1.96*lmer.coef[,2]), dec),")", sep="")
|
|
21
|
1
|
return(cbind(OR, pv))
|
|
22
|
1
|
} else if (family %in% c(NULL, "gaussian")){
|
|
23
|
1
|
coeff = paste(round(lmer.coef[,1], dec), " (", round(lmer.coef[,1] - 1.96*lmer.coef[,2], dec), ",", round(lmer.coef[,1] + 1.96*lmer.coef[,2], dec), ")", sep="")
|
|
24
|
1
|
return(cbind(coeff, pv))
|
|
25
|
0
|
} else if (family %in% c("poisson", "quasipoisson")){
|
|
26
|
0
|
RR = paste(round(exp(lmer.coef[,1]), dec), " (", round(exp(lmer.coef[,1] - 1.96*lmer.coef[,2]), dec), ",", round(exp(lmer.coef[,1] + 1.96*lmer.coef[,2]), dec),")", sep="")
|
|
27
|
0
|
return(cbind(RR, pv))
|
|
28
|
|
}
|
|
29
|
|
}
|
|
30
|
|
|
|
31
|
|
|
|
32
|
|
|
|
33
|
|
|
|
34
|
|
|
|
35
|
|
|
|
36
|
|
|
|
37
|
|
#' @title lmer.display: table for "lmerMod" or "glmerMod" object (lme4 package)
|
|
38
|
|
#' @description Make mixed effect model results from "lmerMod" or "glmerMod" object (lme4 package)
|
|
39
|
|
#' @param lmerMod.obj "lmerMod" or "glmerMod" object
|
|
40
|
|
#' @param dec Decimal, Default: 2
|
|
41
|
|
#' @param ci.ranef Show confidence interval of random effects?, Default: F
|
|
42
|
|
#' @return Table: fixed & random effect
|
|
43
|
|
#' @details DETAILS
|
|
44
|
|
#' @examples
|
|
45
|
|
#' library(geepack)
|
|
46
|
|
#' data(dietox)
|
|
47
|
|
#' dietox$Cu <- as.factor(dietox$Cu)
|
|
48
|
|
#' l1 <- lme4::lmer(Weight ~ Time + Cu + (1|Pig) + (1|Evit), data = dietox)
|
|
49
|
|
#' lmer.display(l1)
|
|
50
|
|
#' @rdname lmer.display
|
|
51
|
|
#' @export
|
|
52
|
|
#' @importFrom lme4 lmer glmer confint.merMod
|
|
53
|
|
|
|
54
|
|
lmer.display = function(lmerMod.obj, dec = 2, ci.ranef = F){
|
|
55
|
1
|
sl = summary(lmerMod.obj)
|
|
56
|
1
|
fixef = sl$coefficients[-1,]
|
|
57
|
|
|
|
58
|
1
|
mdata = lmerMod.obj@frame
|
|
59
|
1
|
forms = as.character(sl$call[[2]])
|
|
60
|
1
|
y = forms[2]
|
|
61
|
1
|
xfr = strsplit(forms[3]," \\+ ")[[1]]
|
|
62
|
1
|
xr = xfr[grepl("\\|", xfr)]
|
|
63
|
1
|
xf = xfr[!grepl("\\|", xfr)]
|
|
64
|
1
|
family.lmer = ifelse(is.null(sl$family), "gaussian", sl$family)
|
|
65
|
1
|
uni.res = ""
|
|
66
|
1
|
if (family.lmer == "gaussian"){
|
|
67
|
1
|
unis = lapply(xf, function(x){summary(lmer(as.formula(paste(y, "~", x," + ", paste(xr, collapse =" + "), sep="")), data = mdata))$coefficients})
|
|
68
|
1
|
unis2 = Reduce(rbind, unis)
|
|
69
|
1
|
uni.res <- unis2[rownames(unis2) != "(Intercept)",]
|
|
70
|
|
} else{
|
|
71
|
1
|
unis = lapply(xf, function(x){summary(glmer(as.formula(paste(y, "~", x," + ", paste(xr, collapse =" + "), sep="")), data = mdata, family = family.lmer))$coefficients})
|
|
72
|
1
|
unis2 = Reduce(rbind, unis)
|
|
73
|
1
|
uni.res <- unis2[rownames(unis2) != "(Intercept)",]
|
|
74
|
|
}
|
|
75
|
|
|
|
76
|
1
|
if(length(xf) == 1){
|
|
77
|
1
|
fix.all = lmerExp(uni.res, family=family.lmer, dec = dec)
|
|
78
|
1
|
family.label = colnames(fix.all)[1]
|
|
79
|
1
|
colnames(fix.all) = c(paste(family.label, "(95%CI)",sep = ""), "P value")
|
|
80
|
1
|
rownames(fix.all) = rownames(sl$coefficients)[-1]
|
|
81
|
|
} else{
|
|
82
|
1
|
fix.all = cbind(lmerExp(uni.res, family=family.lmer, dec = dec), lmerExp(fixef, family=family.lmer, dec = dec))
|
|
83
|
1
|
family.label = colnames(fix.all)[1]
|
|
84
|
1
|
colnames(fix.all) = c(paste("crude ", family.label, "(95%CI)",sep = ""), "crude P value", paste("adj. ", family.label, "(95%CI)",sep = ""), "adj. P value")
|
|
85
|
1
|
rownames(fix.all) = rownames(sl$coefficients)[-1]
|
|
86
|
|
}
|
|
87
|
|
|
|
88
|
1
|
ranef = data.frame(sl$varcor)[,c(1,4)]
|
|
89
|
1
|
ranef.out = cbind(as.character(round(ranef[,2], dec)), matrix(NA, nrow(ranef), ncol(fix.all) - 1 ))
|
|
90
|
1
|
rownames(ranef.out) = ranef[,1]
|
|
91
|
1
|
if (ci.ranef){
|
|
92
|
0
|
ranef.ci = confint.merMod(lmerMod.obj, parm = 1:nrow(ranef), oldNames = F)^2
|
|
93
|
0
|
ranef.paste = paste(round(ranef$vcov, dec)," (", round(ranef.ci[,1], dec), ",", round(ranef.ci[,2], dec),")", sep="")
|
|
94
|
0
|
ranef.out = cbind(ranef.paste, matrix(NA, nrow(ranef), 3))
|
|
95
|
0
|
rownames(ranef.out) = ranef[,1]
|
|
96
|
|
}
|
|
97
|
|
|
|
98
|
1
|
ranef.na = rbind(rep(NA, ncol(fix.all)), ranef.out)
|
|
99
|
1
|
rownames(ranef.na)[1] = "Random effects"
|
|
100
|
1
|
colnames(ranef.na) = colnames(fix.all)
|
|
101
|
|
|
|
102
|
|
|
|
103
|
|
## rownames
|
|
104
|
1
|
fix.all.list = lapply(xf, function(x){fix.all[grepl(x, rownames(fix.all)),]})
|
|
105
|
1
|
varnum.mfac = which(lapply(fix.all.list, length) > ncol(fix.all))
|
|
106
|
1
|
lapply(varnum.mfac, function(x){fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])})
|
|
107
|
1
|
fix.all.unlist = Reduce(rbind, fix.all.list)
|
|
108
|
|
|
|
109
|
1
|
rn.list = lapply(xf, function(x){rownames(fix.all)[grepl(x, rownames(fix.all))]})
|
|
110
|
1
|
varnum.2fac = which(lapply(xf, function(x){length(sapply(mdata, levels)[[x]])}) == 2)
|
|
111
|
0
|
lapply(varnum.2fac, function(x){rn.list[[x]] <<- paste(xf[x], ": ", levels(mdata[, xf[x]])[2], " vs ", levels(mdata[, xf[x]])[1], sep="")})
|
|
112
|
1
|
lapply(varnum.mfac, function(x){rn.list[[x]] <<- c(paste(xf[x],": ref.=", levels(mdata[, xf[x]])[1], sep=""), gsub(xf[x]," ", rn.list[[x]]))})
|
|
113
|
1
|
if (class(fix.all.unlist)[1] == "character"){
|
|
114
|
1
|
fix.all.unlist = t(data.frame(fix.all.unlist))
|
|
115
|
|
}
|
|
116
|
1
|
rownames(fix.all.unlist) = unlist(rn.list)
|
|
117
|
1
|
tb.df = as.data.frame(rbind(fix.all.unlist, ranef.na))
|
|
118
|
|
|
|
119
|
|
|
|
120
|
|
## metric
|
|
121
|
1
|
no.grp = sl$ngrps
|
|
122
|
1
|
no.obs = nrow(mdata)
|
|
123
|
1
|
ll = round(sl$logLik[[1]], dec)
|
|
124
|
1
|
aic = round(sl$AICtab, dec)[1]
|
|
125
|
1
|
met = c(NA, no.grp, no.obs, ll, aic)
|
|
126
|
1
|
met.mat = cbind(met, matrix(NA, length(met), ncol(fix.all) - 1 ))
|
|
127
|
1
|
rownames(met.mat) = c("Metrics", paste("No. of groups (", rownames(met.mat)[2:(length(no.grp) + 1)], ")", sep=""), "No. of observations", "Log-likelihood", "AIC value")
|
|
128
|
1
|
met.mat[, 1] = as.character(met.mat[, 1])
|
|
129
|
1
|
colnames(met.mat) = colnames(tb.df)
|
|
130
|
1
|
tb.lmerMod = rbind(tb.df, met.mat)
|
|
131
|
1
|
lapply(seq(2, ncol(fix.all), by =2), function(x){tb.lmerMod[, x] <<- as.numeric(as.vector(tb.lmerMod[, x]))})
|
|
132
|
|
|
|
133
|
|
## caption
|
|
134
|
1
|
cap.lmerMod = paste(sl$methTitle, " : ", y," ~ ", forms[3], sep="")
|
|
135
|
1
|
return(list(table = tb.lmerMod, caption = cap.lmerMod))
|
|
136
|
|
}
|
jinseob2kim