1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
#' @title Plot the results of a simulation of the random effects
#' @name plotREsim
#' @description Plot the simulated random effects on a ggplot2 chart. Points that
#' are distinguishable from zero (i.e. the confidence band based on \code{level}
#' does not cross the red line) are highlighted. Currently, the plots are ordered
#' according to the grouping factor.
#' @param data a data.frame generated by \code{\link{REsim}} with simulations of
#' the random effects of a \code{\link{merMod}}
#' @param level the width of the confidence interval
#' @param stat a character value indicating the variable name in data of the
#' midpoint of the estimated interval, e.g. "mean" or "median"
#' @param sd a logical indicating whether or not to plot error bars around
#' the estimates (default is TRUE). Calculates the width of the error bars
#' based on \code{level} and the variable named "sd" in \code{data}
#' @param sigmaScale a numeric value to divide the estimate and the standard
#' deviation by in the case of doing an effect size calculation
#' @param oddsRatio logical, should the parameters be converted to odds ratios
#' before plotting
#' @param labs logical, include the labels of the groups on the x-axis
#' @param facet Accepts either logical (\code{TRUE}) or \code{list} to specify which
#' random effects to plot. If \code{TRUE}, facets by both \code{groupFctr} and \code{term}.
#' If \code{list} selects the panel specified by the named elements of the list
#' @return a ggplot2 plot of the coefficient effects
#' @examples
#' \donttest{
#' fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
#' (p1 <- plotREsim(REsim(fm1)))
#' #Plot just the random effects for the Days slope
#' (p2 <- plotREsim(REsim(fm1), facet= list(groupFctr= "Subject", term= "Days")))
#' }
#' @export
#' @import ggplot2
plotREsim <- function(data, level = 0.95, stat = "median", sd = TRUE,
sigmaScale = NULL, oddsRatio = FALSE, labs = FALSE,
facet= TRUE){
# error checking
plot_sim_error_chks(type= "RE", level= level, stat= stat, sd= sd, sigmaScale= sigmaScale,
oddsRatio= oddsRatio, labs= labs, facet= facet)
# check for faceting
facet_logical <- is.logical(facet)
if (!facet_logical) {
data <- data[data$groupFctr == facet[[1]] & data$term == facet[[2]], ]
}
if(!missing(sigmaScale)){
data[, "sd"] <- data[, "sd"] / sigmaScale
data[, stat] <- data[, stat] / sigmaScale
}
data[, "sd"] <- data[, "sd"] * qnorm(1-((1-level)/2))
data[, "ymax"] <- data[, stat] + data[, "sd"]
data[, "ymin"] <- data[, stat] - data[, "sd"]
data[, "sig"] <- data[, "ymin"] > 0 | data[, "ymax"] < 0
hlineInt <- 0
if(oddsRatio == TRUE){
data[, "ymax"] <- exp(data[, "ymax"])
data[, stat] <- exp(data[, stat])
data[, "ymin"] <- exp(data[, "ymin"])
hlineInt <- 1
}
data <- data[order(data[,"groupFctr"], data[,"term"], data[,stat]),]
rownames(data) <- 1:nrow(data)
data[,"xvar"] <- factor(paste(data$groupFctr, data$groupID, sep=""),
levels=unique(paste(data$groupFctr,data$groupID, sep="")),
ordered=TRUE)
if(labs == TRUE){
xlabs.tmp <- element_text(face = "bold", angle=90, vjust=.5)
} else {
data[,"xvar"] <- as.numeric(data[,"xvar"])
xlabs.tmp <- element_blank()
}
p <- ggplot(data, aes(x = .data[["xvar"]], y = .data[[stat]], ymax = .data[["ymax"]], ymin = .data[["ymin"]])) +
geom_hline(yintercept = hlineInt, color = I("red"), linewidth = I(1.1)) +
geom_point(color="gray75", alpha=1/(nrow(data)^.33), size=I(0.5)) +
geom_point(data=subset(data, sig==TRUE), size=I(3)) +
labs(x = "Group", y = "Effect Range", title = "Effect Ranges") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = xlabs.tmp,
axis.ticks.x = element_blank())
if (sd) {
p <- p +
geom_pointrange(alpha = 1/(nrow(data)^.33)) +
geom_pointrange(data=subset(data, sig==TRUE), alpha = 0.25)
}
# check facet
if (facet_logical) {
return(p + facet_grid(term ~ groupFctr, scales = "free_x"))
} else {
return(p)
}
}
#' @title Plot the results of a simulation of the fixed effects
#' @name plotFEsim
#' @description Plot the simulated fixed effects on a ggplot2 chart
#' @param data a data.frame generated by \code{\link{FEsim}} with simulations of
#' the fixed effects of a \code{\link{merMod}}
#' @param level the width of the confidence interval
#' @param stat a character value indicating the variable name in data of the
#' midpoint of the estimated interval, e.g. "mean" or "median"
#' @param sd logical, indicating whether or not to plot error bars around
#' the estimates (default is TRUE). Calculates the width of the error bars
#' based on \code{level} and the variable named "sd" in \code{data}
#' @param intercept logical, should the intercept be included, default is FALSE
#' @param sigmaScale a numeric value to divide the estimate and the standard
#' deviation by in the case of doing an effect size calculation
#' @param oddsRatio logical, should the parameters be converted to odds ratios
#' before plotting
#' @return a ggplot2 plot of the coefficient effects
#' @examples
#' fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
#' (p1 <- plotFEsim(FEsim(fm1)))
#' @export
#' @import ggplot2
plotFEsim <- function(data, level=0.95, stat = "median", sd = TRUE,
intercept = FALSE, sigmaScale = NULL, oddsRatio = FALSE){
# error checking
plot_sim_error_chks(type= "FE", level= level, stat= stat, sd= sd, sigmaScale= sigmaScale,
oddsRatio= oddsRatio, labs= TRUE, facet= TRUE)
if(!missing(sigmaScale)){
data[, "sd"] <- data[, "sd"] / sigmaScale
data[, stat] <- data[, stat] / sigmaScale
}
if(intercept == FALSE){
data <- data[data$term != "(Intercept)", ]
}
data[, "sd"] <- data[, "sd"] * qnorm(1-((1-level)/2))
data[, "ymax"] <- data[, stat] + data[, "sd"]
data[, "ymin"] <- data[, stat] - data[, "sd"]
hlineInt <- 0
if(oddsRatio == TRUE){
data[, "ymax"] <- exp(data[, "ymax"])
data[, stat] <- exp(data[, stat])
data[, "ymin"] <- exp(data[, "ymin"])
hlineInt <- 1
}
xvar <- "term"
data$term <- as.character(data$term)
data$term <- factor(data$term , levels = data[order(data[, stat]), 1])
p <- ggplot(aes(x = .data[[xvar]], y = .data[[stat]], ymax = .data[["ymax"]], ymin = .data[["ymin"]]), data = data) +
geom_hline(yintercept = hlineInt, color = I("red")) +
geom_point(size=I(3)) +
coord_flip() +
theme_bw()
if (sd) {
p <- p + geom_errorbar(width = 0.2)
}
p
}
|
