Age heaping calculation
Init
options(
digits = 2
)
library(kirkegaard)## Loading required package: tidyverse## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()## Loading required package: magrittr##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extract## Loading required package: weights## Loading required package: Hmisc## Loading required package: lattice## Loading required package: survival## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:base':
##
## format.pval, units## Loading required package: assertthat##
## Attaching package: 'assertthat'## The following object is masked from 'package:tibble':
##
## has_name## Loading required package: psych##
## Attaching package: 'psych'## The following object is masked from 'package:Hmisc':
##
## describe## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha##
## Attaching package: 'kirkegaard'## The following object is masked from 'package:psych':
##
## rescale## The following object is masked from 'package:assertthat':
##
## are_equal## The following objects are masked from 'package:purrr':
##
## is_logical, is_numeric## The following object is masked from 'package:base':
##
## +load_packages(
readxl,
betareg,
GGally,
bayesbr,
patchwork
)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2theme_set(theme_bw())
ah_y = scale_y_continuous("Age heaping (% pass)", labels = scales::percent)Data
Data comes in a weird format.
#raw
ah = read_excel("data/Numeracy(Total)_Compact.xlsx",
skip = 0,
guess_max = 9999)
#add ISO
ah$ISO = pu_translate(ah$`country name`, messages = F)
#compute missing by case, and get rid of cases with no data
years = (1500:2022) %>% as.character()
ah %>% select(!!years) %>% miss_by_case(reverse = T) -> missingness_case
missingness_case %>% table2()ah_notempty = ah[missingness_case != 0, ]
#same for years
ah_notempty %>% select(!!years) %>% miss_by_var(reverse = T) -> missingness_year
missingness_year %>% table2()#which years have some data?
missingness_year[missingness_year != 0] %>% names() -> years_with_data
years_with_data## [1] "1500" "1520" "1570" "1590" "1620" "1650" "1670" "1680" "1690" "1700"
## [11] "1710" "1720" "1730" "1740" "1750" "1760" "1770" "1780" "1790" "1800"
## [21] "1810" "1820" "1830" "1840" "1850" "1860" "1870" "1880" "1890" "1900"
## [31] "1910" "1920" "1930" "1940" "1950" "1960" "1970"ah_notempty %>% select(!!years_with_data) %>% {bind_cols(ISO = ah_notempty$ISO,
.)} -> ah_final
#transform
for (year in years_with_data) {
ah_final[[year]] = (ah_final[[year]] / 100) %>% winsorise(upper = .99, lower = .01)
}
#does any country appear in every series?
(countries_datapoints = ah_final %>% miss_by_case(reverse = T) %>% set_names(ah_final$ISO) %>% sort())## HRV TCD ETH GAB FSM NER YEM ZWE MMR DZA BGD CAF PAK SRB SVK AFG BHS BRB BLZ BEN
## 3 4 4 4 4 4 4 4 5 6 6 6 6 6 6 7 7 7 7 7
## BIH BGR BFA BDI KHM CMR CPV CHL COG CIV CUB GRC GRD GTM GNB GUY HND IRN IRQ ISR
## 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## JAM KWT LBR LBY MKD MDG MYS MLI MLT MHL MUS MNG MAR NPL NGA PRY PER SGP LCA SUR
## 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## SYR TZA GMB TON TTO TUN VUT VNM ZMB BRN GHA LUX PHL TGO UGA DOM SLV HTI MDV NIC
## 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 9 9 9 9 9
## SWZ THA BHR BWA ISL KEN PAN BOL FJI IDN ZAF KOR LKA AUS CRI NZL SVN ARM AZE BLR
## 9 9 10 10 10 10 10 11 11 11 11 11 11 12 12 12 12 13 13 13
## COL EST GEO IND LVA TJK TKM UZB NLD CZE JPN KAZ KGZ LTU MDA ROU TUR UKR BEL BRA
## 13 13 13 13 13 13 13 13 14 15 15 15 15 15 15 15 15 15 16 16
## EGY URY CAN ECU IRL NOR POL FIN AUT CHN DNK PRT ESP SWE ITA CHE USA RUS DEU HUN
## 16 16 17 17 17 17 17 18 19 19 19 19 19 19 20 21 21 22 23 23
## GBR MEX FRA ARG
## 23 24 26 27#other datasets
#becker data
becker = read_excel("data/NIQ-DATASET-V1.3.3/NIQ-DATA (V1.3.3).xlsx", sheet = 2, range = "A2:N205") %>%
df_legalize_names()
#smart fraction data for UN regions
sf = read_rds("data/smart fraction data_out.rds")
#add national IQ from Lynn and Rindermann to compare
d = full_join(
becker,
sf %>% select(ISO, UN_macroregion),
by = "ISO"
) %>%
full_join(
tibble(ISO = names(countries_datapoints), datapoints = countries_datapoints),
by = "ISO"
)
#countries with IQ dataModels
The data are multiplied by 100, so divide by 100 first. Then, winsorise values to .01 and .99.
Z-scores and index country
#no but France appears in almost all and did not change its population much in this period, unlike Argentina which appears slightly more
#so we use France as index
#use reverse normal transformaiton to get the implied z score
#then subtract France's z score to set all times on the same scale relative to France's score
#copy initial
ah_z = ah_final
#transform to France z
for (year in years_with_data) {
#to z score
ah_z[[year]] = ah_final[[year]] %>% winsorise(upper = .99, lower = .01) %>% qnorm()
#subtract France
ah_z[[year]] = ah_z[[year]] - ah_z[[year]][ah_final$ISO == "FRA"]
}
#which years have ceiling issues?
yearly_data = tibble(
year = years_with_data %>% as.numeric(),
ceiling = ah_final[, years_with_data] %>% map_dbl(~mean(.>=.99, na.rm=T)),
countries = missingness_year[years_with_data]
)
yearly_data %>%
ggplot(aes(year, ceiling)) +
geom_line() +
scale_y_continuous("% countries at ceiling", labels = scales::percent) ->
p1
yearly_data %>%
ggplot(aes(year, countries)) +
geom_line() +
scale_y_continuous("Countries with data") ->
p2
p1 / p2
GG_save("figs/year_ceiling_data.png")
#average of z scores in the 1800s
ah_z %>%
select(`1800`:`1890`) %>%
df_standardize(exclude_range_01 = F) %>%
rowMeans(na.rm = T) ->
mean_zs_1800s
mean_zs_1800s %>% describe2()#model as beta function
ah_final_long = ah_final %>%
pivot_longer(
cols = -ISO
) %>%
mutate(
year = name %>% as.numeric(),
ISO = ISO %>% fct_relevel("GBR")
)
#join with Rindermann IQ for showcasing the issue
full_join(
ah_z %>% select(ISO, `1910`) %>% df_legalize_names(),
becker %>% select(ISO, R),
) %>%
GG_scatter("x1910", "R", case_names = "ISO") +
labs(x = "Age heaping z-score in 1910", y = "Rindermann national IQ")## Joining, by = "ISO"
## `geom_smooth()` using formula 'y ~ x'
GG_save("figs/ah_1910z_R.png")## `geom_smooth()` using formula 'y ~ x'#plot real data showing growth
model_visual = expand_grid(
ISO = c("GBR", "DNK", "ARG", "FRA", "BRA", "TUR", "HTI"),
year = ah_final_long$year %>% unique()
) %>% mutate(ID = as.character(ISO) + year) %>%
left_join(
ah_final_long %>% mutate(ID = as.character(ISO) + year, obs = value) %>% select(ID, obs)
)## Joining, by = "ID"ggplot(model_visual, aes(year, obs, color = ISO)) +
geom_point() +
geom_smooth() +
ah_y## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 131 rows containing non-finite values (stat_smooth).## Warning: Removed 131 rows containing missing values (geom_point).
GG_save("figs/ah_examples.png")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 131 rows containing non-finite values (stat_smooth).
## Removed 131 rows containing missing values (geom_point).Beta regression
#fit simplest model
betareg::betareg(
value ~ year + ISO,
data = ah_final_long,
) -> beta_fit
#summary
beta_fit##
## Call:
## betareg::betareg(formula = value ~ year + ISO, data = ah_final_long)
##
## Coefficients (mean model with logit link):
## (Intercept) year ISOAFG ISOARG ISOARM ISOAUS
## -21.8198 0.0136 -4.6518 -1.8905 -2.7920 -0.3891
## ISOAUT ISOAZE ISOBDI ISOBEL ISOBEN ISOBFA
## 0.2450 -2.8988 -3.3775 0.3366 -3.8728 -3.5573
## ISOBGD ISOBGR ISOBHR ISOBHS ISOBIH ISOBLR
## -5.2322 -0.7544 -3.8377 -0.9517 -0.8064 -1.9625
## ISOBLZ ISOBOL ISOBRA ISOBRB ISOBRN ISOBWA
## -1.2291 -3.0054 -2.0301 -1.0120 -0.5852 -1.8877
## ISOCAF ISOCAN ISOCHE ISOCHL ISOCHN ISOCIV
## -3.3229 -0.4262 0.3560 -2.0612 0.3023 -2.3749
## ISOCMR ISOCOG ISOCOL ISOCPV ISOCRI ISOCUB
## -3.6627 -0.8474 -2.6181 -1.5004 -2.1954 -0.7306
## ISOCZE ISODEU ISODNK ISODOM ISODZA ISOECU
## -0.2240 0.2215 0.9001 -3.0455 -2.4858 -2.8478
## ISOEGY ISOESP ISOEST ISOETH ISOFIN ISOFJI
## -5.3069 -0.3158 -0.3868 -2.6740 0.1276 -2.2507
## ISOFRA ISOFSM ISOGAB ISOGEO ISOGHA ISOGMB
## -0.3082 -1.2169 -1.2764 -2.6418 -3.2678 -4.0025
## ISOGNB ISOGRC ISOGRD ISOGTM ISOGUY ISOHND
## -3.7084 -2.1233 -1.7641 -3.2038 -1.4121 -1.8931
## ISOHRV ISOHTI ISOHUN ISOIDN ISOIND ISOIRL
## -0.4680 -3.3765 -0.7259 -3.7604 -4.3781 -1.1622
## ISOIRN ISOIRQ ISOISL ISOISR ISOITA ISOJAM
## -4.0094 -3.5592 0.0255 -1.3732 -0.4853 -1.1348
## ISOJPN ISOKAZ ISOKEN ISOKGZ ISOKHM ISOKOR
## 0.5240 -1.5872 -3.3354 -1.9426 -1.4781 -0.4673
## ISOKWT ISOLBR ISOLBY ISOLCA ISOLKA ISOLTU
## -4.5013 -3.1935 -2.9067 -1.6996 -2.3506 -1.7287
## ISOLUX ISOLVA ISOMAR ISOMDA ISOMDG ISOMDV
## -0.0577 -1.0217 -4.8597 -1.7313 -2.9491 -2.1670
## ISOMEX ISOMHL ISOMKD ISOMLI ISOMLT ISOMMR
## -1.9317 -1.0383 -0.8026 -3.6804 -0.7445 -2.7262
## ISOMNG ISOMUS ISOMYS ISONER ISONGA ISONIC
## -0.8026 -0.8418 -2.4193 -3.5094 -4.0528 -3.1918
## ISONLD ISONOR ISONPL ISONZL ISOPAK ISOPAN
## 1.3903 -0.1186 -4.2442 -0.6115 -5.3725 -2.2100
## ISOPER ISOPHL ISOPOL ISOPRT ISOPRY ISOROU
## -2.9601 -3.0602 -1.0829 -0.7790 -1.7543 -0.9090
## ISORUS ISOSGP ISOSLV ISOSRB ISOSUR ISOSVK
## -1.5310 -0.9929 -3.1547 -3.0181 -0.6152 -1.4691
## ISOSVN ISOSWE ISOSWZ ISOSYR ISOTCD ISOTGO
## 0.1218 0.2589 -2.1500 -3.5707 -3.6199 -3.0169
## ISOTHA ISOTJK ISOTKM ISOTON ISOTTO ISOTUN
## -1.0200 -2.5877 -2.3194 -0.6663 -2.1139 -2.4604
## ISOTUR ISOTZA ISOUGA ISOUKR ISOURY ISOUSA
## -3.6690 -3.1278 -3.3949 -1.6534 -1.3442 -0.8422
## ISOUZB ISOVNM ISOVUT ISOYEM ISOZAF ISOZMB
## -2.1589 -0.8026 -1.3031 -4.2544 -2.4480 -2.0558
## ISOZWE
## -1.5007
##
## Phi coefficients (precision model with identity link):
## (phi)
## 28.6#coefs
get_ISO_coefs = function(x) {
x %>% coef() %>% {
.[str_detect(names(.), "^ISO")] -> x
names(x) = names(x) %>% str_replace("ISO", "")
x
#add GBR, which is the baseline level
c(GBR = 0, x)
}
}
beta_fit %>% get_ISO_coefs() %>% sort()## PAK EGY BGD MAR AFG KWT IND YEM NPL NGA IRN
## -5.373 -5.307 -5.232 -4.860 -4.652 -4.501 -4.378 -4.254 -4.244 -4.053 -4.009
## GMB BEN BHR IDN GNB MLI TUR CMR TCD SYR IRQ
## -4.002 -3.873 -3.838 -3.760 -3.708 -3.680 -3.669 -3.663 -3.620 -3.571 -3.559
## BFA NER UGA BDI HTI KEN CAF GHA GTM LBR NIC
## -3.557 -3.509 -3.395 -3.377 -3.377 -3.335 -3.323 -3.268 -3.204 -3.194 -3.192
## SLV TZA PHL DOM SRB TGO BOL PER MDG LBY AZE
## -3.155 -3.128 -3.060 -3.045 -3.018 -3.017 -3.005 -2.960 -2.949 -2.907 -2.899
## ECU ARM MMR ETH GEO COL TJK DZA TUN ZAF MYS
## -2.848 -2.792 -2.726 -2.674 -2.642 -2.618 -2.588 -2.486 -2.460 -2.448 -2.419
## CIV LKA TKM FJI PAN CRI MDV UZB SWZ GRC TTO
## -2.375 -2.351 -2.319 -2.251 -2.210 -2.195 -2.167 -2.159 -2.150 -2.123 -2.114
## CHL ZMB BRA BLR KGZ MEX HND ARG BWA GRD PRY
## -2.061 -2.056 -2.030 -1.962 -1.943 -1.932 -1.893 -1.891 -1.888 -1.764 -1.754
## MDA LTU LCA UKR KAZ RUS ZWE CPV KHM SVK GUY
## -1.731 -1.729 -1.700 -1.653 -1.587 -1.531 -1.501 -1.500 -1.478 -1.469 -1.412
## ISR URY VUT GAB BLZ FSM IRL JAM POL MHL LVA
## -1.373 -1.344 -1.303 -1.276 -1.229 -1.217 -1.162 -1.135 -1.083 -1.038 -1.022
## THA BRB SGP BHS ROU COG USA MUS BIH VNM MNG
## -1.020 -1.012 -0.993 -0.952 -0.909 -0.847 -0.842 -0.842 -0.806 -0.803 -0.803
## MKD PRT BGR MLT CUB HUN TON SUR NZL BRN ITA
## -0.803 -0.779 -0.754 -0.745 -0.731 -0.726 -0.666 -0.615 -0.612 -0.585 -0.485
## HRV KOR CAN AUS EST ESP FRA CZE NOR LUX GBR
## -0.468 -0.467 -0.426 -0.389 -0.387 -0.316 -0.308 -0.224 -0.119 -0.058 0.000
## ISL SVN FIN DEU AUT SWE CHN BEL CHE JPN DNK
## 0.026 0.122 0.128 0.221 0.245 0.259 0.302 0.337 0.356 0.524 0.900
## NLD
## 1.390#plot model predictions and real data
#add model predictions
model_visual$pred = predict(
beta_fit,
newdata = model_visual
)
model_visual %>%
ggplot(aes(year, color = ISO)) +
geom_point(aes(y = obs)) +
geom_line(aes(y = pred)) +
ah_y## Warning: Removed 131 rows containing missing values (geom_point).
GG_save("figs/ah model 1.png")## Warning: Removed 131 rows containing missing values (geom_point).#fit 2nd model
#50 countries with most data
most_data_ISO = ah_final %>% miss_by_case() %>% set_names(ah_final$ISO) %>% sort() %>% .[1:142] %>% names()
betareg::betareg(
value ~ year * ISO,
data = ah_final_long %>% filter(ISO %in% most_data_ISO),
# link = "cloglog",
# type = "BR",
control = betareg:::betareg.control(
maxit = 999999,
hessian = T
)
) -> beta_fit_2
#summary
beta_fit_2 %>% summary()## Warning in sqrt(diag(object$vcov)): NaNs produced##
## Call:
## betareg::betareg(formula = value ~ year * ISO, data = ah_final_long %>%
## filter(ISO %in% most_data_ISO), control = betareg:::betareg.control(maxit = 999999,
## hessian = T))
##
## Standardized weighted residuals 2:
## Min 1Q Median 3Q Max
## -5.157 -0.447 0.013 0.521 6.948
##
## Coefficients (mean model with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.38e+01 NaN NaN NaN
## year 9.25e-03 NaN NaN NaN
## ISOAFG -1.12e+01 NaN NaN NaN
## ISOARG -1.65e+01 3.40e-01 -48.60 <2e-16 ***
## ISOARM -7.69e+01 1.54e+00 -49.92 <2e-16 ***
## ISOAUS 1.83e+01 NaN NaN NaN
## ISOAUT -3.77e-01 NaN NaN NaN
## ISOAZE -5.91e+01 NaN NaN NaN
## ISOBDI -5.23e+01 NaN NaN NaN
## ISOBEL -2.45e+01 NaN NaN NaN
## ISOBEN -2.96e+01 NaN NaN NaN
## ISOBFA -1.11e+01 NaN NaN NaN
## ISOBGD -3.89e+01 NaN NaN NaN
## ISOBGR 2.84e+00 NaN NaN NaN
## ISOBHR -1.22e+02 NaN NaN NaN
## ISOBHS 4.99e+00 NaN NaN NaN
## ISOBIH 1.75e+01 NaN NaN NaN
## ISOBLR -5.61e+01 9.37e-01 -59.85 <2e-16 ***
## ISOBLZ -1.88e+01 NaN NaN NaN
## ISOBOL -5.68e+01 NaN NaN NaN
## ISOBRA -3.84e+00 NaN NaN NaN
## ISOBRB 1.19e+00 NaN NaN NaN
## ISOBRN 1.80e+01 NaN NaN NaN
## ISOBWA -4.09e+01 NaN NaN NaN
## ISOCAF -3.60e+01 NaN NaN NaN
## ISOCAN -3.12e+00 NaN NaN NaN
## ISOCHE -7.37e-02 NaN NaN NaN
## ISOCHL -2.58e+01 NaN NaN NaN
## ISOCHN 1.08e+01 NaN NaN NaN
## ISOCIV -1.93e+01 NaN NaN NaN
## ISOCMR -3.16e+01 NaN NaN NaN
## ISOCOG 3.99e+01 NaN NaN NaN
## ISOCOL -1.23e+01 NaN NaN NaN
## ISOCPV 1.11e+01 NaN NaN NaN
## ISOCRI -2.74e+01 1.14e+00 -23.90 <2e-16 ***
## ISOCUB 6.42e+00 NaN NaN NaN
## ISOCZE -1.31e+01 NaN NaN NaN
## ISODEU -4.19e+00 NaN NaN NaN
## ISODNK 3.28e+00 NaN NaN NaN
## ISODOM -3.23e+01 NaN NaN NaN
## ISODZA -2.62e+01 NaN NaN NaN
## ISOECU -1.31e+01 NaN NaN NaN
## ISOEGY -1.35e+01 NaN NaN NaN
## ISOESP -7.23e-01 NaN NaN NaN
## ISOEST -8.13e+00 NaN NaN NaN
## ISOETH 1.33e+02 NaN NaN NaN
## ISOFIN 1.61e+01 NaN NaN NaN
## ISOFJI -6.02e+01 1.11e+00 -54.04 <2e-16 ***
## ISOFRA -1.11e+01 NaN NaN NaN
## ISOFSM 3.98e+01 NaN NaN NaN
## ISOGAB 2.19e+01 NaN NaN NaN
## ISOGEO -4.88e+01 2.27e+00 -21.49 <2e-16 ***
## ISOGHA -1.17e+01 NaN NaN NaN
## ISOGMB -6.17e+00 NaN NaN NaN
## ISOGNB -1.65e+01 NaN NaN NaN
## ISOGRC -5.32e+01 NaN NaN NaN
## ISOGRD -1.99e+01 NaN NaN NaN
## ISOGTM -1.43e+01 NaN NaN NaN
## ISOGUY -1.06e+01 NaN NaN NaN
## ISOHND 2.38e+01 NaN NaN NaN
## ISOHTI -2.08e+01 NaN NaN NaN
## ISOHUN -1.52e+01 4.11e-01 -36.97 <2e-16 ***
## ISOIDN -5.02e+00 NaN NaN NaN
## ISOIND 5.06e+00 NaN NaN NaN
## ISOIRL -1.70e+01 NaN NaN NaN
## ISOIRN -2.87e+01 NaN NaN NaN
## ISOIRQ -2.92e+01 NaN NaN NaN
## ISOISL 5.35e+00 NaN NaN NaN
## ISOISR -8.01e+01 NaN NaN NaN
## ISOITA -1.17e+00 NaN NaN NaN
## ISOJAM 1.91e+01 NaN NaN NaN
## ISOJPN 1.99e+00 NaN NaN NaN
## ISOKAZ -3.80e+01 1.57e+00 -24.24 <2e-16 ***
## ISOKEN -3.05e+01 NaN NaN NaN
## ISOKGZ -3.90e+01 7.48e-01 -52.15 <2e-16 ***
## ISOKHM 1.17e+01 NaN NaN NaN
## ISOKOR 1.80e+01 NaN NaN NaN
## ISOKWT -5.95e+01 NaN NaN NaN
## ISOLBR -1.61e+00 NaN NaN NaN
## ISOLBY -3.68e+01 NaN NaN NaN
## ISOLCA -2.05e+01 NaN NaN NaN
## ISOLKA -3.17e+01 1.38e+00 -22.94 <2e-16 ***
## ISOLTU -5.40e+01 8.36e-01 -64.60 <2e-16 ***
## ISOLUX 1.93e+01 NaN NaN NaN
## ISOLVA -4.04e+01 NaN NaN NaN
## ISOMAR -3.31e+01 NaN NaN NaN
## ISOMDA -4.79e+01 8.34e-01 -57.48 <2e-16 ***
## ISOMDG -2.37e+01 NaN NaN NaN
## ISOMDV 4.31e+01 NaN NaN NaN
## ISOMEX 2.36e+00 NaN NaN NaN
## ISOMHL 2.12e+00 NaN NaN NaN
## ISOMKD 1.80e+01 NaN NaN NaN
## ISOMLI 3.25e+00 NaN NaN NaN
## ISOMLT 2.37e+01 NaN NaN NaN
## ISOMMR 3.61e+00 NaN NaN NaN
## ISOMNG 1.80e+01 NaN NaN NaN
## ISOMUS 3.75e+00 NaN NaN NaN
## ISOMYS -5.47e+00 3.19e+00 -1.72 0.086 .
## ISONER 1.07e+02 NaN NaN NaN
## ISONGA 8.87e+00 NaN NaN NaN
## ISONIC -2.21e+01 NaN NaN NaN
## ISONLD 3.77e-02 NaN NaN NaN
## ISONOR 4.50e-01 NaN NaN NaN
## ISONPL -1.18e+01 NaN NaN NaN
## ISONZL 5.60e+00 NaN NaN NaN
## ISOPAK -6.82e+01 NaN NaN NaN
## ISOPAN -3.14e+01 NaN NaN NaN
## ISOPER -5.33e+01 NaN NaN NaN
## ISOPHL -5.68e+01 NaN NaN NaN
## ISOPOL -4.67e+00 NaN NaN NaN
## ISOPRT -2.09e+00 2.52e-01 -8.30 <2e-16 ***
## ISOPRY -3.59e+01 NaN NaN NaN
## ISOROU -2.79e+01 NaN NaN NaN
## ISORUS -1.51e+01 2.27e-01 -66.40 <2e-16 ***
## ISOSGP 1.18e+01 NaN NaN NaN
## ISOSLV -2.09e+01 NaN NaN NaN
## ISOSRB -5.18e+01 NaN NaN NaN
## ISOSUR 2.74e+01 NaN NaN NaN
## ISOSVK -4.86e+00 6.91e+00 -0.70 0.482
## ISOSVN 1.71e+01 NaN NaN NaN
## ISOSWE -3.07e+00 NaN NaN NaN
## ISOSWZ -1.31e+01 NaN NaN NaN
## ISOSYR -4.94e+01 NaN NaN NaN
## ISOTCD 1.01e+02 NaN NaN NaN
## ISOTGO -1.54e+01 NaN NaN NaN
## ISOTHA -7.48e+00 NaN NaN NaN
## ISOTJK -3.67e+01 2.00e+00 -18.38 <2e-16 ***
## ISOTKM -4.09e+01 9.06e-01 -45.09 <2e-16 ***
## ISOTON 1.80e+01 NaN NaN NaN
## ISOTTO -1.16e+01 NaN NaN NaN
## ISOTUN -5.29e+01 2.26e+00 -23.44 <2e-16 ***
## ISOTUR -2.82e+01 NaN NaN NaN
## ISOTZA 6.13e+00 NaN NaN NaN
## ISOUGA -2.04e+01 NaN NaN NaN
## ISOUKR -4.03e+01 7.57e-01 -53.31 <2e-16 ***
## ISOURY -2.85e+01 7.59e-01 -37.52 <2e-16 ***
## ISOUSA -8.09e+00 NaN NaN NaN
## ISOUZB -2.61e+01 8.13e-01 -32.13 <2e-16 ***
## ISOVNM 1.80e+01 NaN NaN NaN
## ISOVUT -8.35e+00 NaN NaN NaN
## ISOYEM -6.62e+00 NaN NaN NaN
## ISOZAF -1.69e+01 1.47e+00 -11.48 <2e-16 ***
## ISOZMB -7.00e+00 NaN NaN NaN
## year:ISOAFG 3.61e-03 NaN NaN NaN
## year:ISOARG 8.12e-03 2.75e-04 29.52 <2e-16 ***
## year:ISOARM 3.96e-02 8.85e-04 44.80 <2e-16 ***
## year:ISOAUS -9.42e-03 NaN NaN NaN
## year:ISOAUT 2.37e-04 NaN NaN NaN
## year:ISOAZE 3.01e-02 NaN NaN NaN
## year:ISOBDI 2.55e-02 NaN NaN NaN
## year:ISOBEL 1.39e-02 NaN NaN NaN
## year:ISOBEN 1.36e-02 NaN NaN NaN
## year:ISOBFA 4.13e-03 NaN NaN NaN
## year:ISOBGD 1.77e-02 NaN NaN NaN
## year:ISOBGR -1.50e-03 NaN NaN NaN
## year:ISOBHR 6.16e-02 NaN NaN NaN
## year:ISOBHS -2.64e-03 NaN NaN NaN
## year:ISOBIH -9.03e-03 NaN NaN NaN
## year:ISOBLR 2.91e-02 5.47e-04 53.16 <2e-16 ***
## year:ISOBLZ 9.57e-03 NaN NaN NaN
## year:ISOBOL 2.84e-02 NaN NaN NaN
## year:ISOBRA 1.02e-03 NaN NaN NaN
## year:ISOBRB -9.42e-04 NaN NaN NaN
## year:ISOBRN -9.25e-03 NaN NaN NaN
## year:ISOBWA 2.07e-02 NaN NaN NaN
## year:ISOCAF 1.73e-02 NaN NaN NaN
## year:ISOCAN 1.62e-03 NaN NaN NaN
## year:ISOCHE 2.83e-04 NaN NaN NaN
## year:ISOCHL 1.27e-02 NaN NaN NaN
## year:ISOCHN -5.87e-03 NaN NaN NaN
## year:ISOCIV 9.05e-03 NaN NaN NaN
## year:ISOCMR 1.48e-02 NaN NaN NaN
## year:ISOCOG -2.08e-02 NaN NaN NaN
## year:ISOCOL 5.28e-03 NaN NaN NaN
## year:ISOCPV -6.16e-03 NaN NaN NaN
## year:ISOCRI 1.35e-02 6.10e-04 22.20 <2e-16 ***
## year:ISOCUB -3.29e-03 NaN NaN NaN
## year:ISOCZE 7.23e-03 NaN NaN NaN
## year:ISODEU 2.37e-03 NaN NaN NaN
## year:ISODNK -1.28e-03 NaN NaN NaN
## year:ISODOM 1.55e-02 NaN NaN NaN
## year:ISODZA 1.26e-02 NaN NaN NaN
## year:ISOECU 5.62e-03 NaN NaN NaN
## year:ISOEGY 4.43e-03 NaN NaN NaN
## year:ISOESP -4.66e-05 1.33e-04 -0.35 0.726
## year:ISOEST 4.40e-03 NaN NaN NaN
## year:ISOETH -6.93e-02 NaN NaN NaN
## year:ISOFIN -8.28e-03 NaN NaN NaN
## year:ISOFJI 3.07e-02 6.01e-04 51.11 <2e-16 ***
## year:ISOFRA 6.00e-03 NaN NaN NaN
## year:ISOFSM -2.04e-02 NaN NaN NaN
## year:ISOGAB -1.13e-02 NaN NaN NaN
## year:ISOGEO 2.47e-02 1.34e-03 18.49 <2e-16 ***
## year:ISOGHA 4.61e-03 NaN NaN NaN
## year:ISOGMB 1.33e-03 NaN NaN NaN
## year:ISOGNB 6.91e-03 NaN NaN NaN
## year:ISOGRC 2.72e-02 NaN NaN NaN
## year:ISOGRD 9.82e-03 NaN NaN NaN
## year:ISOGTM 6.03e-03 NaN NaN NaN
## year:ISOGUY 5.09e-03 NaN NaN NaN
## year:ISOHND -1.32e-02 NaN NaN NaN
## year:ISOHTI 9.32e-03 NaN NaN NaN
## year:ISOHUN 8.08e-03 2.58e-04 31.36 <2e-16 ***
## year:ISOIDN 8.41e-04 NaN NaN NaN
## year:ISOIND -4.86e-03 NaN NaN NaN
## year:ISOIRL 8.68e-03 NaN NaN NaN
## year:ISOIRN 1.31e-02 NaN NaN NaN
## year:ISOIRQ 1.36e-02 NaN NaN NaN
## year:ISOISL -2.62e-03 NaN NaN NaN
## year:ISOISR 4.18e-02 NaN NaN NaN
## year:ISOITA 3.36e-04 NaN NaN NaN
## year:ISOJAM -1.01e-02 NaN NaN NaN
## year:ISOJPN -8.46e-04 NaN NaN NaN
## year:ISOKAZ 1.96e-02 7.81e-04 25.08 <2e-16 ***
## year:ISOKEN 1.44e-02 NaN NaN NaN
## year:ISOKGZ 1.99e-02 4.14e-04 48.03 <2e-16 ***
## year:ISOKHM -6.68e-03 NaN NaN NaN
## year:ISOKOR -9.25e-03 NaN NaN NaN
## year:ISOKWT 2.90e-02 NaN NaN NaN
## year:ISOLBR -6.18e-04 NaN NaN NaN
## year:ISOLBY 1.79e-02 NaN NaN NaN
## year:ISOLCA 1.01e-02 NaN NaN NaN
## year:ISOLKA 1.56e-02 7.14e-04 21.90 <2e-16 ***
## year:ISOLTU 2.80e-02 4.70e-04 59.61 <2e-16 ***
## year:ISOLUX -1.04e-02 NaN NaN NaN
## year:ISOLVA 2.13e-02 NaN NaN NaN
## year:ISOMAR 1.50e-02 NaN NaN NaN
## year:ISOMDA 2.48e-02 4.65e-04 53.31 <2e-16 ***
## year:ISOMDG 1.11e-02 NaN NaN NaN
## year:ISOMDV -2.35e-02 NaN NaN NaN
## year:ISOMEX -2.52e-03 NaN NaN NaN
## year:ISOMHL -1.22e-03 NaN NaN NaN
## year:ISOMKD -9.25e-03 NaN NaN NaN
## year:ISOMLI -3.42e-03 NaN NaN NaN
## year:ISOMLT -1.27e-02 NaN NaN NaN
## year:ISOMMR -3.34e-03 NaN NaN NaN
## year:ISOMNG -9.25e-03 NaN NaN NaN
## year:ISOMUS -2.03e-03 NaN NaN NaN
## year:ISOMYS 1.80e-03 1.66e-03 1.08 0.279
## year:ISONER -5.62e-02 NaN NaN NaN
## year:ISONGA -6.61e-03 NaN NaN NaN
## year:ISONIC 1.01e-02 NaN NaN NaN
## year:ISONLD 5.84e-04 NaN NaN NaN
## year:ISONOR -9.56e-05 NaN NaN NaN
## year:ISONPL 4.13e-03 NaN NaN NaN
## year:ISONZL -2.91e-03 NaN NaN NaN
## year:ISOPAK 3.29e-02 NaN NaN NaN
## year:ISOPAN 1.55e-02 NaN NaN NaN
## year:ISOPER 2.69e-02 NaN NaN NaN
## year:ISOPHL 2.85e-02 NaN NaN NaN
## year:ISOPOL 2.07e-03 NaN NaN NaN
## year:ISOPRT 5.16e-04 2.33e-04 2.22 0.027 *
## year:ISOPRY 1.82e-02 NaN NaN NaN
## year:ISOROU 1.47e-02 NaN NaN NaN
## year:ISORUS 7.49e-03 1.88e-04 39.82 <2e-16 ***
## year:ISOSGP -6.13e-03 NaN NaN NaN
## year:ISOSLV 9.51e-03 NaN NaN NaN
## year:ISOSRB 2.69e-02 NaN NaN NaN
## year:ISOSUR -1.44e-02 NaN NaN NaN
## year:ISOSVK 1.87e-03 3.87e-03 0.48 0.629
## year:ISOSVN -8.82e-03 NaN NaN NaN
## year:ISOSWE 1.97e-03 NaN NaN NaN
## year:ISOSWZ 5.97e-03 NaN NaN NaN
## year:ISOSYR 2.42e-02 NaN NaN NaN
## year:ISOTCD -5.31e-02 NaN NaN NaN
## year:ISOTGO 6.70e-03 NaN NaN NaN
## year:ISOTHA 3.68e-03 NaN NaN NaN
## year:ISOTJK 1.83e-02 1.17e-03 15.60 <2e-16 ***
## year:ISOTKM 2.07e-02 5.32e-04 38.88 <2e-16 ***
## year:ISOTON -9.25e-03 NaN NaN NaN
## year:ISOTTO 5.18e-03 NaN NaN NaN
## year:ISOTUN 2.67e-02 1.19e-03 22.49 <2e-16 ***
## year:ISOTUR 1.31e-02 NaN NaN NaN
## year:ISOTZA -4.68e-03 NaN NaN NaN
## year:ISOUGA 9.07e-03 NaN NaN NaN
## year:ISOUKR 2.10e-02 4.29e-04 48.80 <2e-16 ***
## year:ISOURY 1.48e-02 4.18e-04 35.47 <2e-16 ***
## year:ISOUSA 4.05e-03 NaN NaN NaN
## year:ISOUZB 1.29e-02 4.58e-04 28.22 <2e-16 ***
## year:ISOVNM -9.25e-03 NaN NaN NaN
## year:ISOVUT 4.03e-03 NaN NaN NaN
## year:ISOYEM 1.49e-03 NaN NaN NaN
## year:ISOZAF 7.80e-03 7.68e-04 10.15 <2e-16 ***
## year:ISOZMB 2.83e-03 NaN NaN NaN
##
## Phi coefficients (precision model with identity link):
## Estimate Std. Error z value Pr(>|z|)
## (phi) 78.682 0.365 216 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 3.1e+03 on 285 Df
## Pseudo R-squared: 0.94
## Number of iterations in BFGS optimization: 572#add model predictions
model_visual$pred2 = predict(
beta_fit_2,
newdata = model_visual
)
model_visual %>%
ggplot(aes(year, color = ISO)) +
geom_point(aes(y = obs)) +
geom_line(aes(y = pred2)) +
ah_y## Warning: Removed 131 rows containing missing values (geom_point).
#nonlinear year effects but no interaction
betareg::betareg(
value ~ year + I(year^2) + ISO,
data = ah_final_long %>% filter(ISO %in% most_data_ISO),
# link = "cloglog",
control = betareg:::betareg.control(
maxit = 999999,
)
) -> beta_fit_3
#summary
beta_fit_3 %>% summary()##
## Call:
## betareg::betareg(formula = value ~ year + I(year^2) + ISO, data = ah_final_long %>%
## filter(ISO %in% most_data_ISO), control = betareg:::betareg.control(maxit = 999999,
## ))
##
## Standardized weighted residuals 2:
## Min 1Q Median 3Q Max
## -3.963 -0.600 -0.089 0.571 8.533
##
## Coefficients (mean model with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.30e+01 4.93e+00 12.78 < 2e-16 ***
## year -8.16e-02 5.52e-03 -14.80 < 2e-16 ***
## I(year^2) 2.67e-05 1.54e-06 17.30 < 2e-16 ***
## ISOAFG -5.04e+00 2.00e-01 -25.25 < 2e-16 ***
## ISOARG -1.89e+00 1.59e-01 -11.84 < 2e-16 ***
## ISOARM -2.95e+00 1.78e-01 -16.58 < 2e-16 ***
## ISOAUS -4.23e-01 2.93e-01 -1.45 0.14774
## ISOAUT 2.55e-01 1.95e-01 1.31 0.19074
## ISOAZE -3.05e+00 1.77e-01 -17.19 < 2e-16 ***
## ISOBDI -3.82e+00 2.14e-01 -17.86 < 2e-16 ***
## ISOBEL 5.23e-01 2.32e-01 2.26 0.02391 *
## ISOBEN -4.24e+00 2.02e-01 -21.03 < 2e-16 ***
## ISOBFA -3.91e+00 2.07e-01 -18.94 < 2e-16 ***
## ISOBGD -5.59e+00 2.17e-01 -25.78 < 2e-16 ***
## ISOBGR -9.23e-01 3.57e-01 -2.59 0.00973 **
## ISOBHR -4.32e+00 1.87e-01 -23.16 < 2e-16 ***
## ISOBHS -1.26e+00 3.63e-01 -3.46 0.00055 ***
## ISOBIH -1.09e+00 3.74e-01 -2.92 0.00356 **
## ISOBLR -2.09e+00 1.92e-01 -10.91 < 2e-16 ***
## ISOBLZ -1.29e+00 2.85e-01 -4.52 6.3e-06 ***
## ISOBOL -3.32e+00 1.90e-01 -17.48 < 2e-16 ***
## ISOBRA -1.95e+00 1.74e-01 -11.20 < 2e-16 ***
## ISOBRB -1.05e+00 3.03e-01 -3.45 0.00055 ***
## ISOBRN -7.47e-01 3.51e-01 -2.13 0.03310 *
## ISOBWA -2.19e+00 2.35e-01 -9.35 < 2e-16 ***
## ISOCAF -3.65e+00 2.20e-01 -16.60 < 2e-16 ***
## ISOCAN -3.10e-01 2.20e-01 -1.41 0.15880
## ISOCHE 3.40e-01 2.12e-01 1.60 0.10896
## ISOCHL -2.32e+00 2.58e-01 -8.99 < 2e-16 ***
## ISOCHN 3.59e-01 2.25e-01 1.60 0.10945
## ISOCIV -2.77e+00 2.57e-01 -10.77 < 2e-16 ***
## ISOCMR -3.96e+00 2.02e-01 -19.60 < 2e-16 ***
## ISOCOG -1.05e+00 3.63e-01 -2.89 0.00385 **
## ISOCOL -2.69e+00 1.83e-01 -14.72 < 2e-16 ***
## ISOCPV -1.83e+00 3.22e-01 -5.68 1.3e-08 ***
## ISOCRI -2.31e+00 1.95e-01 -11.85 < 2e-16 ***
## ISOCUB -9.56e-01 3.69e-01 -2.59 0.00959 **
## ISOCZE -2.38e-01 2.06e-01 -1.15 0.24839
## ISODEU 1.44e-02 1.80e-01 0.08 0.93615
## ISODNK 1.09e+00 2.29e-01 4.76 2.0e-06 ***
## ISODOM -3.27e+00 1.96e-01 -16.71 < 2e-16 ***
## ISODZA -2.73e+00 2.46e-01 -11.08 < 2e-16 ***
## ISOECU -2.94e+00 1.71e-01 -17.20 < 2e-16 ***
## ISOEGY -5.37e+00 1.89e-01 -28.36 < 2e-16 ***
## ISOESP -7.64e-01 1.82e-01 -4.20 2.7e-05 ***
## ISOEST -3.38e-01 2.59e-01 -1.31 0.19108
## ISOETH -3.18e+00 3.30e-01 -9.64 < 2e-16 ***
## ISOFIN 2.92e-01 2.48e-01 1.18 0.23931
## ISOFJI -2.52e+00 2.03e-01 -12.39 < 2e-16 ***
## ISOFRA -2.28e-01 1.81e-01 -1.26 0.20710
## ISOFSM -1.75e+00 5.06e-01 -3.46 0.00053 ***
## ISOGAB -1.82e+00 4.99e-01 -3.64 0.00027 ***
## ISOGEO -2.78e+00 1.80e-01 -15.43 < 2e-16 ***
## ISOGHA -3.51e+00 2.00e-01 -17.60 < 2e-16 ***
## ISOGMB -4.30e+00 1.99e-01 -21.64 < 2e-16 ***
## ISOGNB -3.88e+00 1.98e-01 -19.63 < 2e-16 ***
## ISOGRC -2.28e+00 2.38e-01 -9.61 < 2e-16 ***
## ISOGRD -1.84e+00 2.50e-01 -7.37 1.8e-13 ***
## ISOGTM -3.36e+00 2.04e-01 -16.53 < 2e-16 ***
## ISOGUY -1.47e+00 2.73e-01 -5.40 6.6e-08 ***
## ISOHND -2.10e+00 2.71e-01 -7.75 9.1e-15 ***
## ISOHTI -3.60e+00 1.91e-01 -18.88 < 2e-16 ***
## ISOHUN -6.81e-01 1.77e-01 -3.84 0.00012 ***
## ISOIDN -3.98e+00 1.80e-01 -22.12 < 2e-16 ***
## ISOIND -4.57e+00 1.75e-01 -26.08 < 2e-16 ***
## ISOIRL -1.18e+00 1.99e-01 -5.90 3.6e-09 ***
## ISOIRN -4.25e+00 1.97e-01 -21.52 < 2e-16 ***
## ISOIRQ -3.79e+00 2.01e-01 -18.86 < 2e-16 ***
## ISOISL 1.57e-01 3.01e-01 0.52 0.60318
## ISOISR -1.54e+00 2.82e-01 -5.45 5.0e-08 ***
## ISOITA -5.17e-01 1.79e-01 -2.89 0.00386 **
## ISOJAM -1.44e+00 3.51e-01 -4.10 4.1e-05 ***
## ISOJPN 2.50e-01 2.45e-01 1.02 0.30715
## ISOKAZ -1.72e+00 1.99e-01 -8.64 < 2e-16 ***
## ISOKEN -3.66e+00 1.91e-01 -19.20 < 2e-16 ***
## ISOKGZ -2.10e+00 1.89e-01 -11.07 < 2e-16 ***
## ISOKHM -1.63e+00 2.89e-01 -5.62 2.0e-08 ***
## ISOKOR -5.48e-01 3.04e-01 -1.80 0.07125 .
## ISOKWT -4.74e+00 1.99e-01 -23.88 < 2e-16 ***
## ISOLBR -3.47e+00 2.11e-01 -16.44 < 2e-16 ***
## ISOLBY -3.19e+00 2.19e-01 -14.59 < 2e-16 ***
## ISOLCA -2.07e+00 3.03e-01 -6.83 8.4e-12 ***
## ISOLKA -2.56e+00 2.02e-01 -12.68 < 2e-16 ***
## ISOLTU -1.90e+00 1.94e-01 -9.78 < 2e-16 ***
## ISOLUX 1.11e-01 3.10e-01 0.36 0.71898
## ISOLVA -1.08e+00 2.25e-01 -4.81 1.5e-06 ***
## ISOMAR -5.11e+00 2.03e-01 -25.20 < 2e-16 ***
## ISOMDA -1.88e+00 1.94e-01 -9.68 < 2e-16 ***
## ISOMDG -3.23e+00 2.18e-01 -14.84 < 2e-16 ***
## ISOMDV -2.33e+00 2.33e-01 -10.01 < 2e-16 ***
## ISOMEX -1.86e+00 1.63e-01 -11.43 < 2e-16 ***
## ISOMHL -1.35e+00 3.57e-01 -3.78 0.00016 ***
## ISOMKD -1.09e+00 3.75e-01 -2.90 0.00372 **
## ISOMLI -3.97e+00 2.02e-01 -19.65 < 2e-16 ***
## ISOMLT -7.89e-01 3.37e-01 -2.34 0.01931 *
## ISOMMR -2.70e+00 2.24e-01 -12.02 < 2e-16 ***
## ISOMNG -1.09e+00 3.75e-01 -2.90 0.00372 **
## ISOMUS -1.02e+00 3.50e-01 -2.91 0.00365 **
## ISOMYS -2.57e+00 2.29e-01 -11.24 < 2e-16 ***
## ISONER -4.05e+00 2.71e-01 -14.95 < 2e-16 ***
## ISONGA -4.29e+00 1.97e-01 -21.74 < 2e-16 ***
## ISONIC -3.41e+00 1.93e-01 -17.64 < 2e-16 ***
## ISONLD 8.50e-01 2.48e-01 3.42 0.00062 ***
## ISONOR 1.97e-02 2.35e-01 0.08 0.93312
## ISONPL -4.62e+00 1.99e-01 -23.21 < 2e-16 ***
## ISONZL -6.84e-01 2.81e-01 -2.44 0.01479 *
## ISOPAK -5.72e+00 2.20e-01 -26.01 < 2e-16 ***
## ISOPAN -2.45e+00 2.16e-01 -11.36 < 2e-16 ***
## ISOPER -3.08e+00 2.05e-01 -15.01 < 2e-16 ***
## ISOPHL -3.27e+00 2.00e-01 -16.33 < 2e-16 ***
## ISOPOL -1.05e+00 2.01e-01 -5.24 1.6e-07 ***
## ISOPRT -9.77e-01 1.71e-01 -5.73 1.0e-08 ***
## ISOPRY -1.95e+00 2.68e-01 -7.27 3.6e-13 ***
## ISOROU -9.34e-01 2.18e-01 -4.29 1.8e-05 ***
## ISORUS -1.55e+00 1.70e-01 -9.17 < 2e-16 ***
## ISOSGP -1.36e+00 3.71e-01 -3.68 0.00024 ***
## ISOSLV -3.37e+00 1.94e-01 -17.39 < 2e-16 ***
## ISOSRB -2.89e+00 2.04e-01 -14.18 < 2e-16 ***
## ISOSUR -6.96e-01 3.59e-01 -1.94 0.05218 .
## ISOSVK -1.33e+00 2.34e-01 -5.69 1.3e-08 ***
## ISOSVN 2.98e-01 2.92e-01 1.02 0.30680
## ISOSWE 3.20e-01 2.29e-01 1.40 0.16236
## ISOSWZ -2.38e+00 2.30e-01 -10.35 < 2e-16 ***
## ISOSYR -3.87e+00 2.03e-01 -19.07 < 2e-16 ***
## ISOTCD -4.16e+00 2.65e-01 -15.68 < 2e-16 ***
## ISOTGO -3.26e+00 2.05e-01 -15.90 < 2e-16 ***
## ISOTHA -1.09e+00 2.79e-01 -3.92 8.7e-05 ***
## ISOTJK -2.71e+00 1.81e-01 -14.99 < 2e-16 ***
## ISOTKM -2.44e+00 1.85e-01 -13.20 < 2e-16 ***
## ISOTON -8.77e-01 3.74e-01 -2.34 0.01906 *
## ISOTTO -2.20e+00 2.33e-01 -9.46 < 2e-16 ***
## ISOTUN -2.68e+00 2.29e-01 -11.71 < 2e-16 ***
## ISOTUR -3.83e+00 1.70e-01 -22.45 < 2e-16 ***
## ISOTZA -3.34e+00 2.09e-01 -15.99 < 2e-16 ***
## ISOUGA -3.76e+00 2.03e-01 -18.48 < 2e-16 ***
## ISOUKR -1.71e+00 1.89e-01 -9.07 < 2e-16 ***
## ISOURY -1.35e+00 1.91e-01 -7.07 1.5e-12 ***
## ISOUSA -7.63e-01 1.92e-01 -3.97 7.2e-05 ***
## ISOUZB -2.25e+00 1.88e-01 -11.92 < 2e-16 ***
## ISOVNM -1.09e+00 3.75e-01 -2.90 0.00372 **
## ISOVUT -1.64e+00 3.36e-01 -4.89 1.0e-06 ***
## ISOYEM -4.83e+00 2.45e-01 -19.66 < 2e-16 ***
## ISOZAF -2.63e+00 2.00e-01 -13.18 < 2e-16 ***
## ISOZMB -2.30e+00 2.59e-01 -8.88 < 2e-16 ***
##
## Phi coefficients (precision model with identity link):
## Estimate Std. Error z value Pr(>|z|)
## (phi) 36.58 1.45 25.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 2.54e+03 on 145 Df
## Pseudo R-squared: 0.866
## Number of iterations: 4971 (BFGS) + 16 (Fisher scoring)#add model predictions
model_visual$pred3 = predict(
beta_fit_3,
newdata = model_visual
)
model_visual %>%
ggplot(aes(year, color = ISO)) +
geom_point(aes(y = obs)) +
geom_line(aes(y = pred3)) +
ah_y## Warning: Removed 131 rows containing missing values (geom_point).
GG_save("figs/ah model 3.png")## Warning: Removed 131 rows containing missing values (geom_point).#save coefs
ISO_coefs = tibble(
ISO = ah_final$ISO,
n = ah_final[years_with_data] %>% miss_by_case(reverse = T),
mean_z_1800s = mean_zs_1800s,
) %>% left_join(
tibble(ISO = beta_fit %>% get_ISO_coefs() %>% names(),
model_1_coef = beta_fit %>% get_ISO_coefs())
) %>% left_join(
tibble(ISO = beta_fit_2 %>% get_ISO_coefs() %>% names(),
model_2_coef = beta_fit_2 %>% get_ISO_coefs())
) %>% left_join(
tibble(ISO = beta_fit_3 %>% get_ISO_coefs() %>% names(),
model_3_coef = beta_fit_3 %>% get_ISO_coefs())
)## Joining, by = "ISO"
## Joining, by = "ISO"
## Joining, by = "ISO"Bayesian beta regression
# install.packages("betaBayes")
# library(betaBayes)
#
# betaBayes::beta4reg(
# value ~ year + ISO,
# data = ah_final_long %>% filter(ISO %in% most_data_ISO),
# # model = "mode",
# # link = "probit",
# # prior = list(th1a0 = 0, th1b0 = 0, th2a0 = 1, th2b0 = 1),
# # na.action = na.exclude
# ) -> beta_bayes_fit_1
#
# beta_bayes_fit_1 %>% summary()
# devtools::install_github("pjoao266/bayesbr")
bayesbr::bayesbr(
value ~ year + ISO,
data = ah_final_long
) -> beta_bayes_fit_1##
## SAMPLING FOR MODEL 'bayesbr' NOW (CHAIN 1).
## Chain 1: Rejecting initial value:
## Chain 1: Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: First shape parameter is -nan, but must be > 0! (in 'model_bayesbr' at line 77)
##
## Chain 1: Rejecting initial value:
## Chain 1: Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: First shape parameter is -nan, but must be > 0! (in 'model_bayesbr' at line 77)
##
## Chain 1: Rejecting initial value:
## Chain 1: Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: First shape parameter is -nan, but must be > 0! (in 'model_bayesbr' at line 77)
##
## Chain 1: Rejecting initial value:
## Chain 1: Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: First shape parameter is -nan, but must be > 0! (in 'model_bayesbr' at line 77)
##
## Chain 1: Rejecting initial value:
## Chain 1: Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: First shape parameter is 0, but must be > 0! (in 'model_bayesbr' at line 77)
##
## Chain 1: Rejecting initial value:
## Chain 1: Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: First shape parameter is -nan, but must be > 0! (in 'model_bayesbr' at line 77)
##
## Chain 1: Rejecting initial value:
## Chain 1: Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: First shape parameter is 0, but must be > 0! (in 'model_bayesbr' at line 77)
##
## Chain 1: Rejecting initial value:
## Chain 1: Error evaluating the log probability at the initial value.
## Chain 1: Exception: beta_lpdf: First shape parameter is 0, but must be > 0! (in 'model_bayesbr' at line 77)
##
## Chain 1:
## Chain 1: Gradient evaluation took 0.001226 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 12.26 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 806.82 seconds (Warm-up)
## Chain 1: 326.957 seconds (Sampling)
## Chain 1: 1133.78 seconds (Total)
## Chain 1:## Warning: There were 5000 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is 2.12, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-essbeta_bayes_fit_1 %>% summary()##
## Call:
## bayesbr::bayesbr(formula = value ~ year + ISO, data = ah_final_long)
##
## Number of Chains: 1
## Iter: 10000
## Warmup: 5000L
##
## Quantile residuals:
## Min 1Q Median 3Q Max
## -4.186 -0.475 0.034 0.576 5.349
##
## Coefficients (mean model):
## Mean Median Std. Dev. HPD_inf HPD_sup
## (Intercept) -21.152 -21.169 0.12179 -21.327 -20.939
## year 0.013 0.013 0.00007 0.013 0.013
## ISOAFG -4.595 -4.592 0.07514 -4.724 -4.477
## ISOARG -1.749 -1.746 0.04906 -1.833 -1.673
## ISOARM -2.705 -2.745 0.11560 -2.860 -2.475
## ISOAUS -0.291 -0.354 0.28993 -0.759 0.212
## ISOAUT 0.347 0.339 0.07034 0.238 0.513
## ISOAZE -2.764 -2.777 0.10255 -2.911 -2.542
## ISOBDI -3.191 -3.192 0.16269 -3.488 -2.912
## ISOBEL 0.357 0.317 0.11681 0.204 0.580
## ISOBEN -3.733 -3.751 0.15986 -4.055 -3.429
## ISOBFA -3.494 -3.514 0.08647 -3.643 -3.340
## ISOBGD -5.150 -5.138 0.22451 -5.589 -4.763
## ISOBGR -0.552 -0.567 0.28409 -1.051 0.033
## ISOBHR -3.595 -3.584 0.10468 -3.762 -3.398
## ISOBHS -0.698 -0.662 0.15611 -1.009 -0.471
## ISOBIH -0.930 -0.943 0.10647 -1.095 -0.750
## ISOBLR -1.864 -1.869 0.13286 -2.094 -1.605
## ISOBLZ -1.125 -1.116 0.14964 -1.362 -0.902
## ISOBOL -2.833 -2.828 0.09063 -3.015 -2.682
## ISOBRA -1.825 -1.822 0.07087 -1.948 -1.709
## ISOBRB -1.066 -1.097 0.12006 -1.232 -0.840
## ISOBRN -0.695 -0.702 0.06918 -0.820 -0.543
## ISOBWA -1.736 -1.738 0.11610 -1.972 -1.541
## ISOCAF -3.228 -3.216 0.19396 -3.617 -2.828
## ISOCAN -0.179 -0.175 0.11563 -0.413 0.014
## ISOCHE 0.152 0.164 0.12413 -0.122 0.341
## ISOCHL -1.779 -1.833 0.24230 -2.103 -1.353
## ISOCHN 0.211 0.200 0.08330 0.066 0.392
## ISOCIV -1.954 -1.950 0.07054 -2.083 -1.835
## ISOCMR -3.579 -3.565 0.10203 -3.789 -3.397
## ISOCOG -1.052 -1.051 0.19400 -1.359 -0.692
## ISOCOL -2.512 -2.523 0.14994 -2.711 -2.171
## ISOCPV -1.701 -1.709 0.09767 -1.902 -1.511
## ISOCRI -1.917 -1.921 0.04942 -2.005 -1.836
## ISOCUB -0.212 -0.179 0.11461 -0.465 -0.042
## ISOCZE -0.122 -0.139 0.14085 -0.351 0.176
## ISODEU 0.161 0.157 0.07604 0.042 0.271
## ISODNK 0.974 0.977 0.05326 0.875 1.083
## ISODOM -2.970 -2.974 0.10798 -3.158 -2.754
## ISODZA -2.593 -2.630 0.15587 -2.805 -2.340
## ISOECU -2.711 -2.713 0.04041 -2.776 -2.622
## ISOEGY -5.023 -5.009 0.05401 -5.162 -4.949
## ISOESP -0.044 -0.028 0.14302 -0.297 0.174
## ISOEST -0.310 -0.307 0.06362 -0.434 -0.199
## ISOETH -2.143 -2.120 0.17216 -2.461 -1.829
## ISOFIN -0.030 -0.030 0.09596 -0.180 0.132
## ISOFJI -2.118 -2.115 0.10631 -2.320 -1.908
## ISOFRA -0.128 -0.137 0.09449 -0.294 0.075
## ISOFSM -1.594 -1.633 0.22649 -1.972 -1.069
## ISOGAB -0.949 -0.941 0.17757 -1.305 -0.637
## ISOGEO -2.570 -2.582 0.08343 -2.734 -2.416
## ISOGHA -3.198 -3.200 0.06187 -3.328 -3.085
## ISOGMB -3.877 -3.890 0.14254 -4.128 -3.618
## ISOGNB -3.700 -3.704 0.06714 -3.809 -3.597
## ISOGRC -1.839 -1.846 0.20883 -2.193 -1.503
## ISOGRD -1.485 -1.432 0.22242 -1.883 -1.161
## ISOGTM -2.987 -3.013 0.08605 -3.114 -2.814
## ISOGUY -1.667 -1.677 0.08121 -1.794 -1.488
## ISOHND -1.837 -1.829 0.17998 -2.150 -1.552
## ISOHRV -0.090 -0.147 0.15960 -0.315 0.167
## ISOHTI -3.325 -3.313 0.12433 -3.566 -3.120
## ISOHUN -0.670 -0.658 0.09855 -0.886 -0.512
## ISOIDN -3.713 -3.708 0.05956 -3.826 -3.606
## ISOIND -4.225 -4.197 0.11727 -4.441 -4.024
## ISOIRL -0.991 -0.986 0.10468 -1.218 -0.833
## ISOIRN -4.046 -4.066 0.07230 -4.150 -3.919
## ISOIRQ -3.418 -3.432 0.05081 -3.489 -3.300
## ISOISL -0.050 -0.024 0.13231 -0.374 0.163
## ISOISR -1.344 -1.367 0.08170 -1.478 -1.172
## ISOITA -0.373 -0.373 0.10114 -0.581 -0.196
## ISOJAM -0.981 -0.970 0.14701 -1.332 -0.701
## ISOJPN 0.250 0.263 0.11672 0.038 0.436
## ISOKAZ -1.369 -1.373 0.15904 -1.679 -1.082
## ISOKEN -3.160 -3.145 0.07921 -3.355 -3.053
## ISOKGZ -1.707 -1.718 0.07035 -1.812 -1.585
## ISOKHM -1.307 -1.299 0.06696 -1.463 -1.200
## ISOKOR -0.288 -0.264 0.13315 -0.536 -0.050
## ISOKWT -4.246 -4.222 0.08902 -4.472 -4.133
## ISOLBR -2.957 -2.957 0.07610 -3.131 -2.841
## ISOLBY -2.915 -2.930 0.13923 -3.123 -2.661
## ISOLCA -1.835 -1.823 0.15969 -2.068 -1.535
## ISOLKA -2.336 -2.328 0.09190 -2.513 -2.208
## ISOLTU -1.681 -1.691 0.15721 -1.908 -1.370
## ISOLUX 0.521 0.537 0.11228 0.257 0.710
## ISOLVA -1.177 -1.194 0.15037 -1.420 -0.933
## ISOMAR -4.650 -4.636 0.14991 -4.910 -4.326
## ISOMDA -1.656 -1.650 0.04281 -1.755 -1.583
## ISOMDG -2.778 -2.752 0.14481 -2.987 -2.581
## ISOMDV -1.938 -1.931 0.22584 -2.332 -1.533
## ISOMEX -1.836 -1.839 0.06483 -1.948 -1.717
## ISOMHL -1.230 -1.210 0.14496 -1.560 -0.955
## ISOMKD -0.598 -0.595 0.08401 -0.766 -0.425
## ISOMLI -3.445 -3.402 0.10550 -3.650 -3.273
## ISOMLT -0.604 -0.578 0.09549 -0.804 -0.478
## ISOMMR -2.719 -2.724 0.08491 -2.869 -2.497
## ISOMNG -0.657 -0.660 0.18007 -1.098 -0.291
## ISOMUS -0.559 -0.639 0.44250 -1.272 0.185
## ISOMYS -2.232 -2.218 0.06367 -2.340 -2.119
## ISONER -3.203 -3.202 0.24192 -3.705 -2.736
## ISONGA -3.794 -3.792 0.06832 -3.950 -3.682
## ISONIC -3.071 -3.080 0.08973 -3.219 -2.873
## ISONLD 1.076 1.082 0.04753 0.987 1.155
## ISONOR 0.074 0.081 0.06379 -0.046 0.198
## ISONPL -4.125 -4.130 0.08690 -4.284 -3.977
## ISONZL -0.248 -0.246 0.05444 -0.344 -0.121
## ISOPAK -5.078 -5.076 0.05660 -5.188 -4.966
## ISOPAN -2.122 -2.076 0.15417 -2.450 -1.888
## ISOPER -2.833 -2.830 0.10291 -2.981 -2.642
## ISOPHL -3.039 -3.033 0.05830 -3.165 -2.933
## ISOPOL -1.065 -1.058 0.11949 -1.279 -0.870
## ISOPRT -0.541 -0.530 0.06543 -0.681 -0.431
## ISOPRY -1.575 -1.564 0.06496 -1.719 -1.477
## ISOROU -0.798 -0.800 0.06563 -0.960 -0.677
## ISORUS -1.337 -1.338 0.09609 -1.483 -1.177
## ISOSGP -0.324 -0.298 0.12719 -0.556 -0.121
## ISOSLV -2.909 -2.914 0.05984 -2.995 -2.806
## ISOSRB -3.082 -3.079 0.08718 -3.255 -2.950
## ISOSUR 0.020 0.018 0.06151 -0.109 0.134
## ISOSVK -1.470 -1.476 0.04954 -1.536 -1.348
## ISOSVN 0.078 0.077 0.11750 -0.107 0.355
## ISOSWE 0.197 0.194 0.06945 0.091 0.354
## ISOSWZ -2.142 -2.160 0.06806 -2.259 -2.010
## ISOSYR -3.313 -3.328 0.15378 -3.533 -2.921
## ISOTCD -3.059 -3.063 0.06169 -3.186 -2.956
## ISOTGO -2.704 -2.695 0.11868 -2.971 -2.515
## ISOTHA -0.902 -0.888 0.16362 -1.214 -0.598
## ISOTJK -2.339 -2.329 0.05209 -2.452 -2.253
## ISOTKM -2.128 -2.137 0.08302 -2.270 -1.938
## ISOTON -0.465 -0.480 0.20251 -0.865 -0.127
## ISOTTO -1.976 -1.987 0.22162 -2.331 -1.515
## ISOTUN -1.998 -1.998 0.09116 -2.146 -1.824
## ISOTUR -3.535 -3.561 0.07769 -3.650 -3.396
## ISOTZA -2.834 -2.824 0.14277 -3.105 -2.575
## ISOUGA -3.180 -3.171 0.14678 -3.423 -2.917
## ISOUKR -1.399 -1.381 0.07538 -1.518 -1.274
## ISOURY -1.202 -1.230 0.15398 -1.448 -0.874
## ISOUSA -0.860 -0.873 0.10825 -1.015 -0.640
## ISOUZB -2.045 -2.041 0.04046 -2.124 -1.961
## ISOVNM -0.177 -0.085 0.41308 -0.960 0.450
## ISOVUT -1.219 -1.245 0.14309 -1.462 -0.910
## ISOYEM -4.088 -4.074 0.25774 -4.519 -3.645
## ISOZAF -2.414 -2.426 0.08266 -2.551 -2.259
## ISOZMB -1.977 -1.984 0.16848 -2.245 -1.606
## ISOZWE -1.174 -1.144 0.19737 -1.604 -0.852
##
## Coefficients (precision model):
## Mean Median Std. Dev. HPD_inf HPD_sup
## (phi) 26 26 0.66 25 27
##
##
## Pseudo R-squared: 0.75288
## AIC: -4702.8
## BIC: -3659
## DIC: -4647.8
## WAIC (SE): -4636.9 (101.35)
## LOOIC (SE): -4640.6 (100.81)#add interaction
# bayesbr::bayesbr(
# value ~ year * ISO,
# data = ah_final_long
# ) -> beta_bayes_fit_2
#
# beta_bayes_fit_2 %>% summary()
ISO_coefs = ISO_coefs %>%
left_join(
tibble(
model_bayes_1_coef = beta_bayes_fit_1 %>% coef %>% .$mean %>% .[-c(1:2)],
ISO = beta_bayes_fit_1 %>% coef %>% .$mean %>% names() %>% .[-c(1:2)] %>% str_replace("^ISO", "")
)
)## Joining, by = "ISO"Compare
#join the beta reg coefs
d = full_join(
d,
ISO_coefs,
by = "ISO"
)
#pairwise
ggpairs(d, columns = c("L_and_V12plusGEO", "R", "QNWplusSASplusGEO", "datapoints", "mean_z_1800s", "model_1_coef", "model_3_coef", "model_bayes_1_coef"), columnLabels = c("Lynn 2012", "Rindermann", "Becker", "AH datapoints", "z-score 1800s", "beta reg 1", "beta reg 3", "beta reg bayes"))## Warning: Removed 49 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 49 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 51 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 134 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning: Removed 49 rows containing missing values (geom_point).## Warning: Removed 49 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 51 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 134 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning: Removed 51 rows containing missing values (geom_point).
## Removed 51 rows containing missing values (geom_point).## Warning: Removed 47 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 107 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning: Removed 104 rows containing missing values (geom_point).
## Removed 104 rows containing missing values (geom_point).## Warning: Removed 105 rows containing missing values (geom_point).## Warning: Removed 104 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 134 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning: Removed 134 rows containing missing values (geom_point).
## Removed 134 rows containing missing values (geom_point).## Warning: Removed 135 rows containing missing values (geom_point).## Warning: Removed 134 rows containing missing values (geom_point).## Warning: Removed 134 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 134 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values## Warning: Removed 104 rows containing missing values (geom_point).
## Removed 104 rows containing missing values (geom_point).## Warning: Removed 105 rows containing missing values (geom_point).## Warning: Removed 104 rows containing missing values (geom_point).## Warning: Removed 134 rows containing missing values (geom_point).## Warning: Removed 104 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning: Removed 106 rows containing missing values (geom_point).
## Removed 106 rows containing missing values (geom_point).## Warning: Removed 107 rows containing missing values (geom_point).## Warning: Removed 106 rows containing missing values (geom_point).## Warning: Removed 135 rows containing missing values (geom_point).## Warning: Removed 106 rows containing missing values (geom_point).## Warning: Removed 106 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 107 rows containing missing values## Warning: Removed 105 rows containing missing values (geom_point).
## Removed 105 rows containing missing values (geom_point).## Warning: Removed 106 rows containing missing values (geom_point).## Warning: Removed 105 rows containing missing values (geom_point).## Warning: Removed 135 rows containing missing values (geom_point).## Warning: Removed 105 rows containing missing values (geom_point).## Warning: Removed 107 rows containing missing values (geom_point).## Warning: Removed 105 rows containing non-finite values (stat_density).
GG_save("figs/ggpairs.png")## Warning: Removed 49 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 49 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 51 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 134 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning: Removed 49 rows containing missing values (geom_point).## Warning: Removed 49 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 51 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 134 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning: Removed 51 rows containing missing values (geom_point).
## Removed 51 rows containing missing values (geom_point).## Warning: Removed 47 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 107 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning: Removed 104 rows containing missing values (geom_point).
## Removed 104 rows containing missing values (geom_point).## Warning: Removed 105 rows containing missing values (geom_point).## Warning: Removed 104 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 134 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 104 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning: Removed 134 rows containing missing values (geom_point).
## Removed 134 rows containing missing values (geom_point).## Warning: Removed 135 rows containing missing values (geom_point).## Warning: Removed 134 rows containing missing values (geom_point).## Warning: Removed 134 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 134 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values## Warning: Removed 104 rows containing missing values (geom_point).
## Removed 104 rows containing missing values (geom_point).## Warning: Removed 105 rows containing missing values (geom_point).## Warning: Removed 104 rows containing missing values (geom_point).## Warning: Removed 134 rows containing missing values (geom_point).## Warning: Removed 104 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 105 rows containing missing values## Warning: Removed 106 rows containing missing values (geom_point).
## Removed 106 rows containing missing values (geom_point).## Warning: Removed 107 rows containing missing values (geom_point).## Warning: Removed 106 rows containing missing values (geom_point).## Warning: Removed 135 rows containing missing values (geom_point).## Warning: Removed 106 rows containing missing values (geom_point).## Warning: Removed 106 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 107 rows containing missing values## Warning: Removed 105 rows containing missing values (geom_point).
## Removed 105 rows containing missing values (geom_point).## Warning: Removed 106 rows containing missing values (geom_point).## Warning: Removed 105 rows containing missing values (geom_point).## Warning: Removed 135 rows containing missing values (geom_point).## Warning: Removed 105 rows containing missing values (geom_point).## Warning: Removed 107 rows containing missing values (geom_point).## Warning: Removed 105 rows containing non-finite values (stat_density).#rescale age heaping to IQ scale
#use Rindermann as basis, find overlap set
R_ah_overlap = d %>% select(R, model_3_coef, ISO) %>% na.omit()
(R_ah_overlap_sumstats = R_ah_overlap %>% describe2())#rescale
d$ah_IQ = d$model_3_coef - R_ah_overlap_sumstats$mean[2]
d$ah_IQ = d$ah_IQ / R_ah_overlap_sumstats$sd[2]
d$ah_IQ = R_ah_overlap_sumstats$mean[1] + d$ah_IQ * R_ah_overlap_sumstats$sd[1]
#best estimates details
GG_scatter(d, "model_3_coef", "R", case_names = "ISO") +
scale_y_continuous("Rindermann's national IQs") +
scale_x_continuous("Age heaping IQs (beta reg)")## `geom_smooth()` using formula 'y ~ x'
GG_save("figs/ah R.png")## `geom_smooth()` using formula 'y ~ x'#world region summary stats
(ah_IQs = describe2(d %>% select(ah_IQ), group = d$UN_macroregion))## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning Inf## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning -
## Inf(R_IQs = describe2(d %>% select(R), group = d$UN_macroregion))#plot regional means
tibble(
region = ah_IQs$group,
ah_IQ = ah_IQs$mean,
R_IQ = R_IQs$mean
) %>%
GG_scatter("ah_IQ", "R_IQ", case_names = "region") +
scale_x_continuous("Age heaping derived IQ") +
scale_y_continuous("Rindermann IQ")## `geom_smooth()` using formula 'y ~ x'
GG_save("figs/regional_means.png")## `geom_smooth()` using formula 'y ~ x'#missing data and IQsMeta
#session info
write_sessioninfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Linux Mint 19.3
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_DK.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_DK.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.1 bayesbr_0.0.1.0 GGally_2.1.2
## [4] betareg_3.1-4 readxl_1.4.0 kirkegaard_2022-08-31
## [7] psych_2.2.3 assertthat_0.2.1 weights_1.0.4
## [10] Hmisc_4.7-0 Formula_1.2-4 survival_3.4-0
## [13] lattice_0.20-45 magrittr_2.0.3 forcats_0.5.1
## [16] stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4
## [19] readr_2.1.2 tidyr_1.2.0 tibble_3.1.7
## [22] ggplot2_3.3.6 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_2.0-3 ellipsis_0.3.2
## [4] modeltools_0.2-23 htmlTable_2.4.0 base64enc_0.1-3
## [7] fs_1.5.2 rstudioapi_0.13 mice_3.14.0
## [10] farver_2.1.0 rstan_2.21.5 flexmix_2.3-18
## [13] fansi_1.0.3 lubridate_1.8.0 xml2_1.3.3
## [16] codetools_0.2-18 splines_4.2.1 mnormt_2.0.2
## [19] knitr_1.39 jsonlite_1.8.0 nloptr_2.0.1
## [22] rematch_1.0.1 broom_0.8.0 cluster_2.1.4
## [25] dbplyr_2.1.1 png_0.1-7 compiler_4.2.1
## [28] httr_1.4.3 backports_1.4.1 Matrix_1.5-1
## [31] fastmap_1.1.0 cli_3.3.0 prettyunits_1.1.1
## [34] htmltools_0.5.2 tools_4.2.1 coda_0.19-4
## [37] gtable_0.3.0 glue_1.6.2 Rcpp_1.0.8.3
## [40] cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.4.1
## [43] gdata_2.18.0 nlme_3.1-159 lmtest_0.9-40
## [46] xfun_0.30 ps_1.7.0 lme4_1.1-29
## [49] rvest_1.0.2 lifecycle_1.0.1 gtools_3.9.2
## [52] stringdist_0.9.8 MASS_7.3-57 zoo_1.8-10
## [55] scales_1.2.0 hms_1.1.1 parallel_4.2.1
## [58] sandwich_3.0-1 inline_0.3.19 RColorBrewer_1.1-3
## [61] yaml_2.3.5 gridExtra_2.3 StanHeaders_2.21.0-7
## [64] loo_2.5.1 sass_0.4.1 rpart_4.1.16
## [67] reshape_0.8.9 latticeExtra_0.6-29 stringi_1.7.6
## [70] highr_0.9 checkmate_2.1.0 pkgbuild_1.3.1
## [73] boot_1.3-28 matrixStats_0.62.0 rlang_1.0.2
## [76] pkgconfig_2.0.3 evaluate_0.15 labeling_0.4.2
## [79] htmlwidgets_1.5.4 processx_3.5.3 tidyselect_1.1.2
## [82] plyr_1.8.7 R6_2.5.1 generics_0.1.2
## [85] DBI_1.1.2 mgcv_1.8-40 pillar_1.7.0
## [88] haven_2.5.0 foreign_0.8-82 withr_2.5.0
## [91] nnet_7.3-17 modelr_0.1.8 crayon_1.5.1
## [94] fdrtool_1.2.17 utf8_1.2.2 tmvnsim_1.0-2
## [97] tzdb_0.3.0 rmarkdown_2.14 jpeg_0.1-9
## [100] grid_4.2.1 data.table_1.14.2 callr_3.7.0
## [103] reprex_2.0.1 digest_0.6.29 RcppParallel_5.1.5
## [106] stats4_4.2.1 munsell_0.5.0 bslib_0.3.1#write data
d %>%
select(ISO, Country, UN_macroregion, ah_IQ, R, L_and_V12plusGEO, QNWplusSASplusGEO, mean_z_1800s, model_1_coef, model_3_coef, model_bayes_1_coef) %>%
write_csv("data/national_data.csv")