-
Notifications
You must be signed in to change notification settings - Fork 2
/
skinny_effects_plot.R
66 lines (51 loc) · 1.98 KB
/
skinny_effects_plot.R
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
library(emmeans)
library(ggplot2); theme_set(theme_bw())
##
source("makestuff/makeRfuns.R") ## Will eventually be a package
startGraphics()
## two models, with and without scaled/centered focal predictor
## (centering is the important part here ...)
basemod <- lm(mpg~disp+wt,data=mtcars)
scaledmod <- lm(mpg~scale(disp)+scale(wt),data=mtcars)
## make an emmeans object and do the (weird) default plot
## Why do we need both a formula and a "keep" argument – JD
e_base <- emmeans(basemod, ~disp, cov.keep="disp", CIs=TRUE)
## Plot the emmeans object by hand
## JD: Do the tails in the middle look suspiciously flat?
## Can we compare to the base "predict" function?
## BC: can you try this?
## JD: Is it easy to supply a different (denser) range of points for x?
d_base <- as.data.frame(e_base)
pp <- (ggplot(d_base, aes(disp, emmean))
+ geom_line()
+ geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL),
colour=NA, alpha=0.1
)
)
print(pp)
## FIXME: Make this work based on variable names instead of parameters names
## zero out all but focal effects
zero_vcov <- function(m, focal_params) {
v <- vcov(m)
focal_var <- v[focal_params,focal_params]
v[] <- 0 ## set all to zero, preserving dims/names
v[focal_params, focal_params] <- focal_var
return(v)
}
## zero_vcov(basemod,"disp")
## zero_vcov(basemod,c("disp","wt"))
e_disp_base <- emmeans(basemod, ~disp, cov.keep="disp", CIs=TRUE,
vcov. = zero_vcov(basemod, "disp")
)
## works, but "anchors" to an inappropriate point (x=0)
pp %+% as.data.frame(e_disp_base)
## If we use the scaled model it works as desired
e_disp_scaled <- emmeans(scaledmod, ~disp, cov.keep="disp", CIs=TRUE,
vcov. = zero_vcov(scaledmod, "scale(disp)")
)
plot(e_disp_scaled)
print(pp %+% as.data.frame(e_disp_scaled))
## identical results here because wt is zero at the model center/reference point
e4 <- emmeans(scaledmod, ~disp, cov.keep="disp", CIs=TRUE,
vcov. = zero_vcov(scaledmod, c("scale(disp)","scale(wt)")))
pp %+% as.data.frame(e4)