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
154
155
156
157
158
159
160
|
## ----include=FALSE--------------------------------------------------------
library(knitr)
opts_chunk$set(
tidy=FALSE,fig.width=5,fig.height=5,cache=FALSE
)
## ----echo=FALSE, results='hide', include=FALSE----------------------------
#options(continue="+ ", prompt="R> ", width=76)
options(show.signif.stars=FALSE)
options(scipen=3)
## -------------------------------------------------------------------------
mvrunif <- function(n, R, min = 0, max = 1){
# method (but not code) from E. Schumann,
# "Generating Correlated Uniform Variates"
# URL:
# <http://comisef.wikidot.com/tutorial:correlateduniformvariates>
# downloaded 2015-05-21
if (!is.matrix(R) || nrow(R) != ncol(R) ||
max(abs(R - t(R))) > sqrt(.Machine$double.eps))
stop("R must be a square symmetric matrix")
if (any(eigen(R, only.values = TRUE)$values <= 0))
stop("R must be positive-definite")
if (any(abs(R) - 1 > sqrt(.Machine$double.eps)))
stop("R must be a correlation matrix")
m <- nrow(R)
R <- 2 * sin(pi * R / 6)
X <- matrix(rnorm(n * m), n, m)
X <- X %*% chol(R)
X <- pnorm(X)
min + X * (max - min)
}
gendata <- function(n = 5000, R, min = -2, max = 2, s = 1.5,
model = expression(x1 + x2 + x3)){
data <- mvrunif(n = n, min = min, max = max, R = R)
colnames(data) <- c("x1", "x2", "x3")
data <- as.data.frame(data)
data$error <- s * rnorm(n)
data$y <- with(data, eval(model) + error)
data
}
R <- function(offdiag = 0, m = 3){
R <- diag(1, m)
R[lower.tri(R)] <- R[upper.tri(R)] <- offdiag
R
}
## -------------------------------------------------------------------------
set.seed(682626)
Data.1 <- gendata(R = R(0), model = expression(x1 + x2 * x3))
round(cor(Data.1), 2)
summary(mod.1 <- lm(y ~ x1 + x2 + x3, data = Data.1))
## ----fig-contrived-1a,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
library(effects)
plot(predictorEffects(mod.1, partial.residuals=TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
rows=1, cols=3)
## ----fig-contrived-1b,include=TRUE, fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x2", "x3"), mod.1, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
## ----fig-contrived-1c,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x1", "x2"), mod.1, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
## -------------------------------------------------------------------------
set.seed(682626)
Data.2 <- gendata(R = R(0.5), model = expression(x1 + x2 * x3))
mod.2 <- lm(y ~ x1 + x2 + x3, data = Data.2)
## ----fig-contrived-2a,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(predictorEffects(mod.2, partial.residuals=TRUE),
partial.residual=list(pch=".", col="#FF00FF80",fig.show='hide'),
axes=list(x=list(rotate=45)),
rows=1, cols=3)
## ----fig-contrived-2b,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x2", "x3"), mod.2, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
## ----fig-contrived-2c,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x1", "x2"), mod.2, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80",fig.show='hide'),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
## -------------------------------------------------------------------------
set.seed(682626)
Data.3 <- gendata(R = R(0.5), model = expression(x1^2 + x2 + x3))
mod.3 <- lm(y ~ x1 + x2 + x3, data = Data.3)
## ----fig-contrived-3a,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(predictorEffects(mod.3, partial.residuals=TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
rows=1, cols=3)
## ----fig-contrived-3b,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x2", "x3"), mod.3, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
## ----fig-contrived-3c,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x1", "x2"), mod.3, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
## -------------------------------------------------------------------------
set.seed(682626)
Data.4 <- gendata(R = R(0.5), model = expression(x1^2 + x2 * x3))
mod.4 <- lm(y ~ x1 + x2 + x3, data = Data.4)
## ----fig-contrived-4a,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(predictorEffects(mod.4, partial.residuals=TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
rows=1, cols=3)
## ----fig-contrived-4b,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x2", "x3"), mod.4, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
## ----fig-contrived-4c,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x1", "x2"), mod.4, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
## ----fig-contrived-5a,include=TRUE,fig.width=5,fig.height=4,fig.show='hide'----
mod.5 <- lm(y ~ poly(x1, 2) + x2*x3, data=Data.4)
plot(Effect("x1", mod.5, partial.residuals=TRUE),
partial.residual=list(pch=".", col="#FF00FF80", span=0.2))
## ----fig-contrived-5b,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x2", "x3"), mod.5, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80"),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)), span=0.5)
## ----fig-contrived-5c,include=TRUE,fig.width=12,fig.height=4,fig.show='hide'----
plot(Effect(c("x1", "x2"), mod.5, partial.residuals = TRUE),
partial.residual=list(pch=".", col="#FF00FF80", span=0.35),
axes=list(x=list(rotate=45)),
lattice=list(layout=c(4, 1)))
|
