-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathModels-PresidentialIdeology.R
2386 lines (2171 loc) · 151 KB
/
Models-PresidentialIdeology.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
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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Setup, load, initial cleaning -------------------------------------------
setwd("/N/slate/rmarcano/LatinAmericapaper")
rm(list = ls( ))
#load("LatinAmerica14.RData")
#load("LatAmBayesian.RData")
library(readxl)
LatAm2014 <- read_excel("LatAm2014.xlsx", col_types = c("text", "text", "numeric",
"numeric", "text", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "text", "text",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"text", "text", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "text", "text"))
attach(LatAm2014)
unique(LatAm2014$country)
LatAm2014$gini100<-gini_equivalized*100
library(stargazer)
library(emmeans)
library(ggeffects)
library(ggplot2)
library(ggthemes)
library(tidyverse)
library(bayesplot)
library(ggridges)
library(sjPlot)
library(insight)
library(httr)
library(brms)
library(devtools)
library(tidybayes)
library(brms)
library(cowplot)
install.packages("sjmisc")
library(lme4)
theme_set(theme_tidybayes() + panel_border())
set.seed(42)
# Initial view, summary ---------------------------------------------------
summary(LatAm2014$gini100)
labels(LatAm2014)
set_labels(LatAm2014, c("Year", "Country", "CountryCode","Equivalized Gini Inequality"))
LatAm2014
# Initial model computation -----------------------------------------------
empty.model<-(lmer(gini100~(1|country),LatAm2014)) #null model
empty.model.Bayesian<-brm(data = LatAm2014, family = gaussian,
gini100~(1|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(cauchy(0,2), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
empty.model.Bayesian
Uncond.Growth.Model<-(lmer(gini100~year+(1|country),LatAm2014)) #unconditional growth model
Uncond.Growth.Model.Bayesian<-brm(data = LatAm2014, family = gaussian,
gini100~year+(1|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
PooledRace<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
PooledRace
PooledRace.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
PooledRaceArea<-mcmc_areas(
PooledRace,
pars = c("b_Blacks","b_Mulattoes","b_Mestizo","b_Amerindians","b_Asians","b_CreolesEtGafurinas","b_socexpgdp","b_educexpgdp",
"b_PresidentPartyLeftRightIndex_mean","b_wage.arrangementRegional","b_wage.arrangementSegmented","b_highlyeducatedlabor","b_logminwage","b_employment"),
prob = 0.8, # 80% intervals
prob_outer = 0.99, # 99%
point_est = "mean")
PooledDiversity<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
PooledDiversity
PooledDiversity.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
PooledHeterogeneity<-brm(data = LatAm2014, family = gaussian,
gini100 ~ HeterogeneityIndex+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
PooledHeterogeneity
PooledFractionalization<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Fractionalization+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
PooledFractionalization
PooledFractional.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
PooledFractional<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
tab_model(PooledRace,PooledHeterogeneity,PooledDiversity,PooledDiversity.lag)
tab_model(PooledFractional,PooledFractional.lag, PooledFractionalization, PooledDiversity)
# Extraction of Models for Table 1 ----------------------------------------
summary(empty.model.Bayesian,robust = TRUE)
summary(Uncond.Growth.Model.Bayesian,robust = TRUE)
round(fixef(empty.model.Bayesian,robust = TRUE),digits = 2)
round(fixef(Uncond.Growth.Model.Bayesian,robust = TRUE),digits = 2)
round(fixef(Uncond.Growth.Model.Bayesian,robust = TRUE),digits = 2)
4.12^2
a<-4.12^2
3.03 ^2
5.98^2
3.13^2
b<-3.13^2
2.93^2
3.36^2
rho<-a/(a+b)
rho # 64% of variation in Income Inequality in LatAm is explained just by country
4.15^2
a1<-4.15^2
3.06 ^2
5.99^2
2.32^2
b1<-2.32^2
2.16^2
2.49^2
#variance of residual changed from null to uncond.growth from 9.8 to 5.38
rho<-a1/(a1+b1)
rho
sigmachange<-(b-b1)/b
sigmachange
#45.40955%, that is, we may conclude 46% of within-country changes
#are systematically associated with time, year.
16.9744/(16.9744+4.6656)
9.1809-4.6656
4.5153/9.1809
tab_model(empty.model.Bayesian,show.re.var = TRUE)
tab_model(Uncond.Growth.Model.Bayesian)
tab_model(empty.model.Bayesian,Uncond.Growth.Model.Bayesian,show.re.var = TRUE)
# Extraction of Models for Table 2 ----------------------------------------
PooledRace$model
round(fixef(PooledRace,robust = TRUE),digits = 2)
round(fixef(PooledDiversity,robust = TRUE),digits = 2)
round(fixef(PooledFractional2,robust = TRUE),digits = 2)
tab_model(PooledRace,PooledDiversity,PooledFractional,PooledFractional2, PooledFractional.lag,show.re.var = TRUE)
plot_model(PooledRace,type = "diag")
plot_model(PooledRace, sort.est = TRUE,bpe.style = "dot", show.values = TRUE,bpe.color = "purple",vline.color = "red") + theme_sjplot() + coord_flip()
?plot_model
library(shinystan)
launch_shinystan(PooledRace)
#Bayes-R2 function using modelled (approximate) residual variance
#from https://avehtari.github.io/bayes_R2/bayes_R2.html
bayes_R2 <- function(fit) {
mupred <- rstanarm::posterior_epred(fit)
var_mupred <- apply(mupred, 1, var)
if (family(fit)$family == "binomial" && NCOL(y) == 1) {
sigma2 <- apply(mupred*(1-mupred), 1, mean)
} else {
sigma2 <- as.matrix(fit, pars = c("sigma"))^2
}
var_mupred / (var_mupred + sigma2)
}
#Bayes-R2 function using residuals
bayes_R2_res <- function(fit) {
y <- rstanarm::get_y(fit)
ypred <- rstanarm::posterior_epred(fit)
if (family(fit)$family == "binomial" && NCOL(y) == 2) {
trials <- rowSums(y)
y <- y[, 1]
ypred <- ypred %*% diag(trials)
}
e <- -1 * sweep(ypred, 2, y)
var_ypred <- apply(ypred, 1, var)
var_e <- apply(e, 1, var)
var_ypred / (var_ypred + var_e)
}
summary(bayes_R2(PooledRace))
summary(bayes_R2_res(PooledRace))
summary(bayes_R2(PooledDiversity))
summary(bayes_R2_res(PooledDiversity))
summary(bayes_R2(PooledFractional))
summary(bayes_R2(PooledHeterogeneity))
summary(bayes_R2_res(PooledHeterogeneity))
summary(PooledHeterogeneity,robust = TRUE)
3.77^2
summary(bayes_R2(PooledFractional))
summary(bayes_R2_res(PooledFractional))
summary(PooledFractional2,robust = TRUE)
summary(bayes_R2(PooledFractional2))
summary(bayes_R2_res(PooledFractional2))
summary(bayes_R2(PooledFractional.lag))
summary(bayes_R2_res(PooledFractional.lag))
summary(bayes_R2(PooledFractionalization))
summary(bayes_R2_res(PooledFractionalization))
summary(PooledFractionalization,robust = TRUE)
summary(PooledFractional,robust = TRUE)
summary(PooledFractional.lag,robust = TRUE)
3.79^2
summary(RandomSlopes13,robust = TRUE)
performance::variance_decomposition(PooledRace)
print(performance::icc(PooledFractional2),digits=8)
VarianceDec1<-performance::variance_decomposition(RandomSlopes13)
VarianceDec2<-performance::variance_decomposition(RandomSlopes14)
VarianceDec3<-performance::variance_decomposition(RandomSlopes15)
print(VarianceDec1,digits=4)
print(VarianceDec2,digits=4)
print(VarianceDec3,digits=4)
print(performance::icc(RandomSlopes13),digits=8)
print(performance::icc(RandomSlopes14),digits=8)
print(performance::icc(RandomSlopes15),digits=8)
PooledRace<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+left+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 10000, warmup = 2500, chains = 5, cores = 1,seed = 5)
summary(PooledRace)
tab_model(PooledRace)
PooledRaceArea<-mcmc_areas(
PooledRace,
pars = c("b_Blacks","b_Mulattoes","b_Mestizo","b_Amerindians","b_Asians","b_CreolesEtGafurinas","b_socexpgdp","b_educexpgdp",
"b_left","b_wage.arrangementRegional","b_wage.arrangementSegmented","b_highlyeducatedlabor","b_logminwage","b_employment"),
prob = 0.8, # 80% intervals
prob_outer = 0.99, # 99%
point_est = "mean")
#forest(PooledRace)
summary(LatAm2014$Blacks*100*0.93)
summary(LatAm2014$Mulattoes*100*0.06)
summary(LatAm2014$Mestizo*100*0.04)
summary(LatAm2014$Amerindians*100*0.01)
summary(LatAm2014$CreolesEtGafurinas*100*0.23)
summary(LatAm2014$Asians*100*0.72)
summary(LatAm2014$SimpsonDiversity*0.61)
summary(LatAm2014$socexpgdp*0.25)
summary(LatAm2014$educexpgdp*0.8)
summary(LatAm2014$highlyeducatedlabor*0.07)
summary(LatAm2014$logminwage*0.86)
7.138*(0.86/100)
summary(LatAm2014$employment*0.14)
PooledDiversityArea<-mcmc_areas(
PooledDiversity,
pars = c("b_SimpsonDiversity","b_socexpgdp","b_educexpgdp",
"b_PresidentPartyLeftRightIndex_mean","b_wage.arrangementRegional","b_wage.arrangementSegmented","b_highlyeducatedlabor","b_logminwage","b_employment"),
prob = 0.8, # 80% intervals
prob_outer = 0.99, # 99%
point_est = "mean")
RandInt.Race<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
RandInt.Diversity<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
RandInt.Fractional<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
RandInt.Fractional2<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
RandInt.Fractional3<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,5), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 100000, warmup = 10000, chains = 10, cores = 10,seed = 5)
RandomIntDiversity<-brmstools::forest(RandInt.Diversity,pars = c("Intercept","logminwage", "wage.arrangementregional",
"SimpsonDiversity","educexpgdp",
"PresidentPartyLeftRightIndex_mean","highlyeducatedlabor","socexpgdp","employment"))
RandomIntRace<-brmstools::forest(RandInt.Race,pars = c("Intercept","logminwage", "wage.arrangementregional",
"Blacks","Mulattoes","Mestizo","educexpgdp",
"PresidentPartyLeftRightIndex_mean","highlyeducatedlabor","socexpgdp","employment"))
RandomSlopes1<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes2<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes3<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(51.24, 4), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes4<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes5<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes6<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(51.24, 4), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes7<-brm(data = LatAm2014, family = gaussian,
gini100 ~ DiverseFractional+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes8<-brm(data = LatAm2014, family = gaussian,
gini100 ~ DiverseFractional+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes9<-brm(data = LatAm2014, family = gaussian,
gini100 ~ DiverseFractional+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(51.24, 4), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes10<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes11<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes12<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(51.24, 4), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes13<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes14<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes15<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_mean+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_mean|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes1.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes2.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes3.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(51.24, 4), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes4.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes5.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes6.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ Blacks+Mulattoes+Mestizo+Amerindians+Asians+CreolesEtGafurinas+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(51.24, 4), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes7.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ DiverseFractional+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes8.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ DiverseFractional+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes9.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ DiverseFractional+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(51.24, 4), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes10.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes11.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes12.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractLag1+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(51.24, 4), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes13.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes14.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
RandomSlopes15.lag<-brm(data = LatAm2014, family = gaussian,
gini100 ~ EthnicFractionalizationTS+socexpgdp+educexpgdp+PresidentPartyLeftRightIndex_meanlag+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+PresidentPartyLeftRightIndex_meanlag|country),
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma),
prior(lkj(2), class = cor)),
control=list(adapt_delta=0.999, max_treedepth=30),
iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
save.image("IterationsWithHighlyEducAlternateLeftsPresidentsOnly.RData")
library(shinystan)
launch_shinystan(PooledRace)
#Bayes-R2 function using modelled (approximate) residual variance
#from https://avehtari.github.io/bayes_R2/bayes_R2.html
bayes_R2 <- function(fit) {
mupred <- rstanarm::posterior_epred(fit)
var_mupred <- apply(mupred, 1, var)
if (family(fit)$family == "binomial" && NCOL(y) == 1) {
sigma2 <- apply(mupred*(1-mupred), 1, mean)
} else {
sigma2 <- as.matrix(fit, pars = c("sigma"))^2
}
var_mupred / (var_mupred + sigma2)
}
#Bayes-R2 function using residuals
bayes_R2_res <- function(fit) {
y <- rstanarm::get_y(fit)
ypred <- rstanarm::posterior_epred(fit)
if (family(fit)$family == "binomial" && NCOL(y) == 2) {
trials <- rowSums(y)
y <- y[, 1]
ypred <- ypred %*% diag(trials)
}
e <- -1 * sweep(ypred, 2, y)
var_ypred <- apply(ypred, 1, var)
var_e <- apply(e, 1, var)
var_ypred / (var_ypred + var_e)
}
summary(bayes_R2(PooledRace))
summary(bayes_R2_res(PooledRace))
summary(bayes_R2(PooledDiversity))
summary(bayes_R2_res(PooledDiversity))
summary(bayes_R2(PooledFractional2))
summary(bayes_R2_res(PooledFractional2))
summary(RandomSlopes13,robust = TRUE)
#RandomSlopesAllWage<-brm(data = LatAm2014, family = gaussian,
# gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+left+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+wage.arrangement|country),
# prior = c(prior(normal(51.24, 4), class = Intercept),
# prior(normal(0, 2), class = b),
# prior(cauchy(0,2), class = sigma),
# prior(lkj(2), class = cor)),
# control=list(adapt_delta=0.999, max_treedepth=30),
# iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
#RandomSlopesAllDiversity<-brm(data = LatAm2014, family = gaussian,
# gini100 ~ SimpsonDiversity+socexpgdp+educexpgdp+left+wage.arrangement+highlyeducatedlabor+logminwage+employment+year+(1+SimpsonDiversity|country),
# prior = c(prior(normal(51.24, 4), class = Intercept),
# prior(normal(0, 2), class = b),
# prior(cauchy(0,2), class = sigma),
# prior(lkj(2), class = cor)),
# control=list(adapt_delta=0.999, max_treedepth=30),
# iter = 100000, warmup = 20000, chains = 10, cores = 10,seed = 5)
library(effects)
#Subsets for each year
LatAm<-split(LatAm2014,LatAm2014$year)
countryyear<-c()
for (i in 1:22) {
countryyear[i]<-brm(data = LatAm[[i]], family = gaussian,
gini100~SimpsonDiversity+socexpgdp+educexpgdp+left+wage.arrangement+highlyeducatedlabor+logminwage+employment,
prior = c(prior(normal(49.43, 5.018121), class = Intercept),
prior(normal(0, 2), class = b),
prior(cauchy(0,2), class = sigma)),
control=list(adapt_delta=0.99, max_treedepth=20),
iter = 10000, warmup = 2500, chains = 4, cores = 1,seed = 5)
print(paste0("Datos para el anno ", i+1992, " listos."))
}
library(broom)
broom.mixed::tidy(M1.all)
str(rstan::extract(M1.all))
posteriors<-c()
for (i in 1:22) {
posteriors[i]<-posterior_samples(countryyear[i])
print(paste0("Posteriores para el a?o ", i+1992, " listos."))
}
post93$year <- rep(1993,times=length(post93))
post94$year <- rep(1994,times=length(post94))
post95$year <- rep(1995,times=length(post95))
post96$year <- rep(1996,times=length(post96))
post97$year <- rep(1997,times=length(post97))
post98$year <- rep(1998,times=length(post98))
post99$year <- rep(1999,times=length(post99))
post00$year <- rep(2000,times=length(post00))
post01$year <- rep(2001,times=length(post01))
post02$year <- rep(2002,times=length(post02))
post03$year <- rep(2003,times=length(post03))
post04$year <- rep(2004,times=length(post04))
post05$year <- rep(2005,times=length(post05))
post06$year <- rep(2006,times=length(post06))
post07$year <- rep(2007,times=length(post07))
post08$year <- rep(2008,times=length(post08))
post09$year <- rep(2009,times=length(post09))
post10$year <- rep(2010,times=length(post10))
post11$year <- rep(2011,times=length(post11))
post12$year <- rep(2012,times=length(post12))
post13$year <- rep(2013,times=length(post13))
post14$year <- rep(2014,times=length(post14))
post93$SimpsonDiversity_mean <- rep(mean(post93$b_SimpsonDiversity),times=30000)
post94$SimpsonDiversity_mean <- rep(mean(post94$b_SimpsonDiversity),times=30000)
post95$SimpsonDiversity_mean <- rep(mean(post95$b_SimpsonDiversity),times=30000)
post96$SimpsonDiversity_mean <- rep(mean(post96$b_SimpsonDiversity),times=30000)
post97$SimpsonDiversity_mean <- rep(mean(post97$b_SimpsonDiversity),times=30000)
post98$SimpsonDiversity_mean <- rep(mean(post98$b_SimpsonDiversity),times=30000)
post99$SimpsonDiversity_mean <- rep(mean(post99$b_SimpsonDiversity),times=30000)
post00$SimpsonDiversity_mean <- rep(mean(post00$b_SimpsonDiversity),times=30000)
post01$SimpsonDiversity_mean <- rep(mean(post01$b_SimpsonDiversity),times=30000)
post02$SimpsonDiversity_mean <- rep(mean(post02$b_SimpsonDiversity),times=30000)
post03$SimpsonDiversity_mean <- rep(mean(post03$b_SimpsonDiversity),times=30000)
post04$SimpsonDiversity_mean <- rep(mean(post04$b_SimpsonDiversity),times=30000)
post05$SimpsonDiversity_mean <- rep(mean(post05$b_SimpsonDiversity),times=30000)
post06$SimpsonDiversity_mean <- rep(mean(post06$b_SimpsonDiversity),times=30000)
post07$SimpsonDiversity_mean <- rep(mean(post07$b_SimpsonDiversity),times=30000)
post08$SimpsonDiversity_mean <- rep(mean(post08$b_SimpsonDiversity),times=30000)
post09$SimpsonDiversity_mean <- rep(mean(post09$b_SimpsonDiversity),times=30000)
post10$SimpsonDiversity_mean <- rep(mean(post10$b_SimpsonDiversity),times=30000)
post11$SimpsonDiversity_mean <- rep(mean(post11$b_SimpsonDiversity),times=30000)
post12$SimpsonDiversity_mean <- rep(mean(post12$b_SimpsonDiversity),times=30000)
post13$SimpsonDiversity_mean <- rep(mean(post13$b_SimpsonDiversity),times=30000)
post14$SimpsonDiversity_mean <- rep(mean(post14$b_SimpsonDiversity),times=30000)
post93$socexpgdp_mean <- rep(mean(post93$b_socexpgdp),times=30000)
post94$socexpgdp_mean <- rep(mean(post94$b_socexpgdp),times=30000)
post95$socexpgdp_mean <- rep(mean(post95$b_socexpgdp),times=30000)
post96$socexpgdp_mean <- rep(mean(post96$b_socexpgdp),times=30000)
post97$socexpgdp_mean <- rep(mean(post97$b_socexpgdp),times=30000)
post98$socexpgdp_mean <- rep(mean(post98$b_socexpgdp),times=30000)
post99$socexpgdp_mean <- rep(mean(post99$b_socexpgdp),times=30000)
post00$socexpgdp_mean <- rep(mean(post00$b_socexpgdp),times=30000)
post01$socexpgdp_mean <- rep(mean(post01$b_socexpgdp),times=30000)
post02$socexpgdp_mean <- rep(mean(post02$b_socexpgdp),times=30000)
post03$socexpgdp_mean <- rep(mean(post03$b_socexpgdp),times=30000)
post04$socexpgdp_mean <- rep(mean(post04$b_socexpgdp),times=30000)
post05$socexpgdp_mean <- rep(mean(post05$b_socexpgdp),times=30000)
post06$socexpgdp_mean <- rep(mean(post06$b_socexpgdp),times=30000)
post07$socexpgdp_mean <- rep(mean(post07$b_socexpgdp),times=30000)
post08$socexpgdp_mean <- rep(mean(post08$b_socexpgdp),times=30000)
post09$socexpgdp_mean <- rep(mean(post09$b_socexpgdp),times=30000)
post10$socexpgdp_mean <- rep(mean(post10$b_socexpgdp),times=30000)
post11$socexpgdp_mean <- rep(mean(post11$b_socexpgdp),times=30000)
post12$socexpgdp_mean <- rep(mean(post12$b_socexpgdp),times=30000)
post13$socexpgdp_mean <- rep(mean(post13$b_socexpgdp),times=30000)
post14$socexpgdp_mean <- rep(mean(post14$b_socexpgdp),times=30000)
post93$educexpgdp_mean <- rep(mean(post93$b_educexpgdp),times=30000)
post94$educexpgdp_mean <- rep(mean(post94$b_educexpgdp),times=30000)
post95$educexpgdp_mean <- rep(mean(post95$b_educexpgdp),times=30000)
post96$educexpgdp_mean <- rep(mean(post96$b_educexpgdp),times=30000)
post97$educexpgdp_mean <- rep(mean(post97$b_educexpgdp),times=30000)
post98$educexpgdp_mean <- rep(mean(post98$b_educexpgdp),times=30000)
post99$educexpgdp_mean <- rep(mean(post99$b_educexpgdp),times=30000)
post00$educexpgdp_mean <- rep(mean(post00$b_educexpgdp),times=30000)
post01$educexpgdp_mean <- rep(mean(post01$b_educexpgdp),times=30000)
post02$educexpgdp_mean <- rep(mean(post02$b_educexpgdp),times=30000)
post03$educexpgdp_mean <- rep(mean(post03$b_educexpgdp),times=30000)
post04$educexpgdp_mean <- rep(mean(post04$b_educexpgdp),times=30000)
post05$educexpgdp_mean <- rep(mean(post05$b_educexpgdp),times=30000)
post06$educexpgdp_mean <- rep(mean(post06$b_educexpgdp),times=30000)
post07$educexpgdp_mean <- rep(mean(post07$b_educexpgdp),times=30000)
post08$educexpgdp_mean <- rep(mean(post08$b_educexpgdp),times=30000)
post09$educexpgdp_mean <- rep(mean(post09$b_educexpgdp),times=30000)
post10$educexpgdp_mean <- rep(mean(post10$b_educexpgdp),times=30000)
post11$educexpgdp_mean <- rep(mean(post11$b_educexpgdp),times=30000)
post12$educexpgdp_mean <- rep(mean(post12$b_educexpgdp),times=30000)
post13$educexpgdp_mean <- rep(mean(post13$b_educexpgdp),times=30000)
post14$educexpgdp_mean <- rep(mean(post14$b_educexpgdp),times=30000)
post93$left_mean <- rep(mean(post93$b_left),times=30000)
post94$left_mean <- rep(mean(post94$b_left),times=30000)
post95$left_mean <- rep(mean(post95$b_left),times=30000)
post96$left_mean <- rep(mean(post96$b_left),times=30000)
post97$left_mean <- rep(mean(post97$b_left),times=30000)
post98$left_mean <- rep(mean(post98$b_left),times=30000)
post99$left_mean <- rep(mean(post99$b_left),times=30000)
post00$left_mean <- rep(mean(post00$b_left),times=30000)
post01$left_mean <- rep(mean(post01$b_left),times=30000)
post02$left_mean <- rep(mean(post02$b_left),times=30000)
post03$left_mean <- rep(mean(post03$b_left),times=30000)
post04$left_mean <- rep(mean(post04$b_left),times=30000)
post05$left_mean <- rep(mean(post05$b_left),times=30000)
post06$left_mean <- rep(mean(post06$b_left),times=30000)
post07$left_mean <- rep(mean(post07$b_left),times=30000)
post08$left_mean <- rep(mean(post08$b_left),times=30000)
post09$left_mean <- rep(mean(post09$b_left),times=30000)
post10$left_mean <- rep(mean(post10$b_left),times=30000)
post11$left_mean <- rep(mean(post11$b_left),times=30000)
post12$left_mean <- rep(mean(post12$b_left),times=30000)
post13$left_mean <- rep(mean(post13$b_left),times=30000)
post14$left_mean <- rep(mean(post14$b_left),times=30000)
post93$wage.arrangementregional_mean <- rep(mean(post93$b_wage.arrangementRegional),times=30000)
post94$wage.arrangementregional_mean <- rep(mean(post94$b_wage.arrangementRegional),times=30000)
post95$wage.arrangementregional_mean <- rep(mean(post95$b_wage.arrangementRegional),times=30000)
post96$wage.arrangementregional_mean <- rep(mean(post96$b_wage.arrangementRegional),times=30000)
post97$wage.arrangementregional_mean <- rep(mean(post97$b_wage.arrangementRegional),times=30000)
post98$wage.arrangementregional_mean <- rep(mean(post98$b_wage.arrangementRegional),times=30000)
post99$wage.arrangementregional_mean <- rep(mean(post99$b_wage.arrangementRegional),times=30000)
post00$wage.arrangementregional_mean <- rep(mean(post00$b_wage.arrangementRegional),times=30000)
post01$wage.arrangementregional_mean <- rep(mean(post01$b_wage.arrangementRegional),times=30000)
post02$wage.arrangementregional_mean <- rep(mean(post02$b_wage.arrangementRegional),times=30000)
post03$wage.arrangementregional_mean <- rep(mean(post03$b_wage.arrangementRegional),times=30000)
post04$wage.arrangementregional_mean <- rep(mean(post04$b_wage.arrangementRegional),times=30000)
post05$wage.arrangementregional_mean <- rep(mean(post05$b_wage.arrangementRegional),times=30000)
post06$wage.arrangementregional_mean <- rep(mean(post06$b_wage.arrangementRegional),times=30000)
post07$wage.arrangementregional_mean <- rep(mean(post07$b_wage.arrangementRegional),times=30000)
post08$wage.arrangementregional_mean <- rep(mean(post08$b_wage.arrangementRegional),times=30000)
post09$wage.arrangementregional_mean <- rep(mean(post09$b_wage.arrangementRegional),times=30000)
post10$wage.arrangementregional_mean <- rep(mean(post10$b_wage.arrangementRegional),times=30000)
post11$wage.arrangementregional_mean <- rep(mean(post11$b_wage.arrangementRegional),times=30000)
post12$wage.arrangementregional_mean <- rep(mean(post12$b_wage.arrangementRegional),times=30000)
post13$wage.arrangementregional_mean <- rep(mean(post13$b_wage.arrangementRegional),times=30000)
post14$wage.arrangementregional_mean <- rep(mean(post14$b_wage.arrangementRegional),times=30000)
post93$wage.arrangementsegmented_mean <- rep(mean(post93$b_wage.arrangementSegmented),times=30000)
post94$wage.arrangementsegmented_mean <- rep(mean(post94$b_wage.arrangementSegmented),times=30000)
post95$wage.arrangementsegmented_mean <- rep(mean(post95$b_wage.arrangementSegmented),times=30000)
post96$wage.arrangementsegmented_mean <- rep(mean(post96$b_wage.arrangementSegmented),times=30000)
post97$wage.arrangementsegmented_mean <- rep(mean(post97$b_wage.arrangementSegmented),times=30000)
post98$wage.arrangementsegmented_mean <- rep(mean(post98$b_wage.arrangementSegmented),times=30000)
post99$wage.arrangementsegmented_mean <- rep(mean(post99$b_wage.arrangementSegmented),times=30000)
post00$wage.arrangementsegmented_mean <- rep(mean(post00$b_wage.arrangementSegmented),times=30000)
post01$wage.arrangementsegmented_mean <- rep(mean(post01$b_wage.arrangementSegmented),times=30000)
post02$wage.arrangementsegmented_mean <- rep(mean(post02$b_wage.arrangementSegmented),times=30000)
post03$wage.arrangementsegmented_mean <- rep(mean(post03$b_wage.arrangementSegmented),times=30000)
post04$wage.arrangementsegmented_mean <- rep(mean(post04$b_wage.arrangementSegmented),times=30000)
post05$wage.arrangementsegmented_mean <- rep(mean(post05$b_wage.arrangementSegmented),times=30000)
post06$wage.arrangementsegmented_mean <- rep(mean(post06$b_wage.arrangementSegmented),times=30000)
post07$wage.arrangementsegmented_mean <- rep(mean(post07$b_wage.arrangementSegmented),times=30000)
post08$wage.arrangementsegmented_mean <- rep(mean(post08$b_wage.arrangementSegmented),times=30000)
post09$wage.arrangementsegmented_mean <- rep(mean(post09$b_wage.arrangementSegmented),times=30000)
post10$wage.arrangementsegmented_mean <- rep(mean(post10$b_wage.arrangementSegmented),times=30000)
post11$wage.arrangementsegmented_mean <- rep(mean(post11$b_wage.arrangementSegmented),times=30000)
post12$wage.arrangementsegmented_mean <- rep(mean(post12$b_wage.arrangementSegmented),times=30000)
post13$wage.arrangementsegmented_mean <- rep(mean(post13$b_wage.arrangementSegmented),times=30000)
post14$wage.arrangementsegmented_mean <- rep(mean(post14$b_wage.arrangementSegmented),times=30000)
post93$highlyeducatedlabor_mean <- rep(mean(post93$b_highlyeducatedlabor),times=30000)
post94$highlyeducatedlabor_mean <- rep(mean(post94$b_highlyeducatedlabor),times=30000)
post95$highlyeducatedlabor_mean <- rep(mean(post95$b_highlyeducatedlabor),times=30000)
post96$highlyeducatedlabor_mean <- rep(mean(post96$b_highlyeducatedlabor),times=30000)
post97$highlyeducatedlabor_mean <- rep(mean(post97$b_highlyeducatedlabor),times=30000)
post98$highlyeducatedlabor_mean <- rep(mean(post98$b_highlyeducatedlabor),times=30000)
post99$highlyeducatedlabor_mean <- rep(mean(post99$b_highlyeducatedlabor),times=30000)
post00$highlyeducatedlabor_mean <- rep(mean(post00$b_highlyeducatedlabor),times=30000)
post01$highlyeducatedlabor_mean <- rep(mean(post01$b_highlyeducatedlabor),times=30000)
post02$highlyeducatedlabor_mean <- rep(mean(post02$b_highlyeducatedlabor),times=30000)
post03$highlyeducatedlabor_mean <- rep(mean(post03$b_highlyeducatedlabor),times=30000)
post04$highlyeducatedlabor_mean <- rep(mean(post04$b_highlyeducatedlabor),times=30000)
post05$highlyeducatedlabor_mean <- rep(mean(post05$b_highlyeducatedlabor),times=30000)
post06$highlyeducatedlabor_mean <- rep(mean(post06$b_highlyeducatedlabor),times=30000)
post07$highlyeducatedlabor_mean <- rep(mean(post07$b_highlyeducatedlabor),times=30000)
post08$highlyeducatedlabor_mean <- rep(mean(post08$b_highlyeducatedlabor),times=30000)
post09$highlyeducatedlabor_mean <- rep(mean(post09$b_highlyeducatedlabor),times=30000)
post10$highlyeducatedlabor_mean <- rep(mean(post10$b_highlyeducatedlabor),times=30000)
post11$highlyeducatedlabor_mean <- rep(mean(post11$b_highlyeducatedlabor),times=30000)
post12$highlyeducatedlabor_mean <- rep(mean(post12$b_highlyeducatedlabor),times=30000)
post13$highlyeducatedlabor_mean <- rep(mean(post13$b_highlyeducatedlabor),times=30000)
post14$highlyeducatedlabor_mean <- rep(mean(post14$b_highlyeducatedlabor),times=30000)
post93$logminwage_mean <- rep(mean(post93$b_logminwage),times=30000)
post94$logminwage_mean <- rep(mean(post94$b_logminwage),times=30000)
post95$logminwage_mean <- rep(mean(post95$b_logminwage),times=30000)
post96$logminwage_mean <- rep(mean(post96$b_logminwage),times=30000)
post97$logminwage_mean <- rep(mean(post97$b_logminwage),times=30000)
post98$logminwage_mean <- rep(mean(post98$b_logminwage),times=30000)
post99$logminwage_mean <- rep(mean(post99$b_logminwage),times=30000)
post00$logminwage_mean <- rep(mean(post00$b_logminwage),times=30000)
post01$logminwage_mean <- rep(mean(post01$b_logminwage),times=30000)
post02$logminwage_mean <- rep(mean(post02$b_logminwage),times=30000)
post03$logminwage_mean <- rep(mean(post03$b_logminwage),times=30000)
post04$logminwage_mean <- rep(mean(post04$b_logminwage),times=30000)
post05$logminwage_mean <- rep(mean(post05$b_logminwage),times=30000)
post06$logminwage_mean <- rep(mean(post06$b_logminwage),times=30000)
post07$logminwage_mean <- rep(mean(post07$b_logminwage),times=30000)
post08$logminwage_mean <- rep(mean(post08$b_logminwage),times=30000)
post09$logminwage_mean <- rep(mean(post09$b_logminwage),times=30000)
post10$logminwage_mean <- rep(mean(post10$b_logminwage),times=30000)
post11$logminwage_mean <- rep(mean(post11$b_logminwage),times=30000)
post12$logminwage_mean <- rep(mean(post12$b_logminwage),times=30000)
post13$logminwage_mean <- rep(mean(post13$b_logminwage),times=30000)
post14$logminwage_mean <- rep(mean(post14$b_logminwage),times=30000)
post93$employment_mean <- rep(mean(post93$b_employment),times=30000)
post94$employment_mean <- rep(mean(post94$b_employment),times=30000)
post95$employment_mean <- rep(mean(post95$b_employment),times=30000)
post96$employment_mean <- rep(mean(post96$b_employment),times=30000)
post97$employment_mean <- rep(mean(post97$b_employment),times=30000)
post98$employment_mean <- rep(mean(post98$b_employment),times=30000)
post99$employment_mean <- rep(mean(post99$b_employment),times=30000)
post00$employment_mean <- rep(mean(post00$b_employment),times=30000)
post01$employment_mean <- rep(mean(post01$b_employment),times=30000)
post02$employment_mean <- rep(mean(post02$b_employment),times=30000)
post03$employment_mean <- rep(mean(post03$b_employment),times=30000)
post04$employment_mean <- rep(mean(post04$b_employment),times=30000)
post05$employment_mean <- rep(mean(post05$b_employment),times=30000)
post06$employment_mean <- rep(mean(post06$b_employment),times=30000)
post07$employment_mean <- rep(mean(post07$b_employment),times=30000)
post08$employment_mean <- rep(mean(post08$b_employment),times=30000)
post09$employment_mean <- rep(mean(post09$b_employment),times=30000)
post10$employment_mean <- rep(mean(post10$b_employment),times=30000)
post11$employment_mean <- rep(mean(post11$b_employment),times=30000)
post12$employment_mean <- rep(mean(post12$b_employment),times=30000)
post13$employment_mean <- rep(mean(post13$b_employment),times=30000)
post14$employment_mean <- rep(mean(post14$b_employment),times=30000)
colnames(fits)
library(plyr)
fits = rbind.fill(post93, post94, post95, post96,post97,post98,post99,post00,
post01,post02,post03,post04,post05,post06,post07,post08,post09,
post10,post11,post12,post13,post14)
save.image("/geode2/home/u070/rmarcano/Carbonate/Documents/Latin America paper/IterationsWithHighlyEducAlternateLefts.RData")
library(ggdist)
Heterogeneity<-ggplot(fits, aes(y = (as.year), x = SimpsonDiversity_mean)) +
stat_halfeye(aes(x=b_SimpsonDiversity),
.width = c(.90, .5),
slab_fill= "deepskyblue1",
point_color = "gold1",
point_fill = "darkslategray",
interval_color = "gray0") +
labs(x="Effect on Inequality from Racial Heterogeneity")+
geom_vline(xintercept = 0, linetype = "dashed")
Heterogeneity<-Heterogeneity+ coord_flip()+geom_hline(yintercept = 0,colour="red")+ylab(NULL)+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
Heterogeneity
SocExpGDP<-ggplot(fits, aes(y = (as.year), x = socexpgdp_mean)) +
stat_halfeye(aes(x=b_socexpgdp),.width = c(.90, .5)) +
labs(x="Effect on Inequality from Social Expenditure as a share of GDP")+
geom_vline(xintercept = 0, linetype = "dashed")
SocExpGDP<-SocExpGDP+ coord_flip()+ylab(NULL)+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
EmployIneq
SocExpGDP
EduExpGDP<-ggplot(fits, aes(y = (as.year), x = educexpgdp_mean)) +
stat_halfeye(aes(x=b_educexpgdp),.width = c(.90, .5)) +
labs(x="Effect on Inequality from Educational Expenditure as a share of GDP")+