-
Notifications
You must be signed in to change notification settings - Fork 2
/
mediate_bin_preds.R
104 lines (89 loc) · 3.19 KB
/
mediate_bin_preds.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
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
library(shellpipes)
library(vareffects); varefftheme()
library(ggpubr)
library(ggplot2)
library(dplyr)
loadEnvironments()
startGraphics()
## Not mediated
### Not corrected
pred_notmediated_none <- varpred(mod_notmediated_bin
, "x"
, bias.adjust="none"
, modelname="none"
)
pred_notmediated_none_mean <- getmeans(pred_notmediated_none, what="estimate")
### Bias corrected
pred_notmediated_pop <- varpred(mod_notmediated_bin
, "x"
, bias.adjust="population"
, modelname="bias corrected"
)
pred_notmediated_pop_mean <- getmeans(pred_notmediated_pop, what="estimate")
### Binned obs
binned_df <- binfun(mod_notmediated_bin, focal="x", bins=50, groups=NULL)
### Combine all predictions
vlist <- list(pred_notmediated_none, pred_notmediated_pop)
pred_notmediated_plots <- (comparevarpred(vlist=vlist
, lnames=NULL
, plotit=TRUE
, addmarginals=FALSE
, ci=FALSE
)
+ geom_hline(data=pred_notmediated_none_mean, aes(yintercept=fit, colour=model, lty=model))
+ geom_hline(data=pred_notmediated_pop_mean, aes(yintercept=fit, colour=model, lty=model))
+ geom_hline(data=observed_df_med, aes(yintercept=zbin, colour="observed", lty="observed"))
+ geom_vline(data=observed_df_med, aes(xintercept=x), lty=2, col="grey")
+ geom_point(data=binned_df, aes(x=x, y=zbin), colour="grey")
+ scale_colour_manual(breaks = c("observed", "none", "bias corrected")
, values=c("observed"="red", "none"="blue", "bias corrected"="black")
)
+ scale_linetype_manual(values=c("observed"=2, "none"=3, "bias corrected"=4))
+ labs(colour="model", linetype="model", title="Non-mediated")
+ theme(legend.position="bottom")
)
## Mediated
### Not corrected
pred_mediated_none <- varpred(mod_mediated_bin
, "x"
, bias.adjust="none"
, modelname="none"
)
pred_mediated_none_mean <- getmeans(pred_mediated_none, what="estimate")
### Bias corrected
pred_mediated_pop <- varpred(mod_mediated_bin
, "x"
, bias.adjust="population"
, modelname="bias corrected"
)
pred_mediated_pop_mean <- getmeans(pred_mediated_pop, what="estimate")
### Binned obs
binned_df <- binfun(mod_mediated_bin, focal="x", bins=50, groups=NULL)
### Combine all predictions
vlist <- list(pred_mediated_none, pred_mediated_pop)
pred_mediated_plots <- (comparevarpred(vlist=vlist
, lnames=NULL
, plotit=TRUE
, addmarginals=FALSE
, ci=FALSE
)
+ geom_hline(data=pred_mediated_none_mean, aes(yintercept=fit, colour=model, lty=model))
+ geom_hline(data=pred_mediated_pop_mean, aes(yintercept=fit, colour=model, lty=model))
+ geom_hline(data=observed_df_med, aes(yintercept=zbin, colour="observed", lty="observed"))
+ geom_vline(data=observed_df_med, aes(xintercept=x), lty=2, col="grey")
+ geom_point(data=binned_df, aes(x=x, y=zbin), colour="grey")
+ scale_colour_manual(breaks = c("observed", "none", "bias corrected")
, values=c("observed"="red", "none"="blue", "bias corrected"="black")
)
+ scale_linetype_manual(values=c("observed"=2, "none"=3, "bias corrected"=4))
+ labs(colour="model", linetype="model", title="Mediated")
+ theme(legend.position="bottom")
)
pred_mediate_bin_plots <- ggarrange(pred_notmediated_plots
, pred_mediated_plots + rremove("ylab")
, common.legend=TRUE
, legend="bottom"
, ncol=2
)
print(pred_mediate_bin_plots)
saveVars(pred_mediate_bin_plots)