Russian history: frequentist versus Bayesian case
Alexander Levakov, Senior Research Fellow, Ph.D
December, 2016
library(rjags)## Loading required package: coda## Linked to JAGS 4.1.0## Loaded modules: basemod,bugslibrary(mcmc)
library(stringr)
library(BEST)Research goal
The Tsardom of Russia, also known as the Tsardom of Muscovy, was the name of the centralized Russian state from Ivan IV’s assumption of the title of Tsar in 1547 until Peter the Great’s foundation of the Russian Empire in 1721. The House of Romanov was the second dynasty (race), after the Rurik dynasty, to rule over Russia, and reigned from 1613 until the abdication of Czar Nicholas II on March 15, 1917, as a result of the February Revolution See: https://en.wikipedia.org/wiki/Tsardom_of_Russia, https://en.wikipedia.org/wiki/House_of_Romanov and http://www.spsl.nsc.ru/history/descr/main_p.htm.
We want to estimate the parameters of ruling duration for both races, using both classical (frequentist’s) and Bayesian approach (http://www.indiana.edu/~kruschke/BEST/). Let’s do it!
Data
Tsardom<-c(1462-1505, 1505-1533, 1533-1584, 1584-1598, 1598-1605, 1605-1605, 1605-1606, 1606-1610, 1610-1613, 1613-1645, 1645-1676, 1676-1682, 1682-1696)
Tsardom<-abs(Tsardom)
Tsardom## [1] 43 28 51 14 7 0 1 4 3 32 31 6 14hist(Tsardom)
Romanov<-c(1696-1725,1725-1727,1727-1730,1730-1740,1740-1741,1741-1761,1761-1762,1762-1796,1796-1801,1801-1825,1825-1855,1855-1881,1881-1894,1894-1917)
Romanov<-abs(Romanov)
Romanov## [1] 29 2 3 10 1 20 1 34 5 24 30 26 13 23hist(Romanov)
Tsars<-c(Tsardom,Romanov)
Tsars## [1] 43 28 51 14 7 0 1 4 3 32 31 6 14 29 2 3 10 1 20 1 34 5 24
## [24] 30 26 13 23hist(Tsars)
Summary of the data
summary(Tsardom)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 4 14 18 31 51summary(Romanov)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 3.50 16.50 15.79 25.50 34.00summary(Tsars)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 3.50 14.00 16.85 28.50 51.00qc<-c(0.1,0.25,0.5,0.6, 0.75,0.9,0.95)
quantile(Tsardom,qc)## 10% 25% 50% 60% 75% 90% 95%
## 1.4 4.0 14.0 16.8 31.0 40.8 46.2quantile(Romanov,qc)## 10% 25% 50% 60% 75% 90% 95%
## 1.3 3.5 16.5 22.4 25.5 29.7 31.4quantile(Tsars,qc)## 10% 25% 50% 60% 75% 90% 95%
## 1.0 3.5 14.0 21.8 28.5 32.8 40.3For the last Romanov we see \(T=23\) that is Nicholas II ruled over Russia until the abdication being in the 60% percentile of the sample while for the Tsardom this percentile is \(16.8\).
Null hypothesis
We postulate our null hypothesis of equality of two samples (Romanov, Tsardom) means as \(H_0: T_{Romanov} = T_{Tsardom}\). Consequently our alternative hypothesis will be \(H_1: {T_{Romanov}} \neq {T_{Tsardom}}\).
Frequentist estimation
#Student t-test for two means
t.test(x = Tsardom,y = Romanov)##
## Welch Two Sample t-test
##
## data: Tsardom and Romanov
## t = 0.38586, df = 21.45, p-value = 0.7034
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -9.70439 14.13296
## sample estimates:
## mean of x mean of y
## 18.00000 15.78571#Test for variances equality
var.test(x = Tsardom, y = Romanov)##
## F test to compare two variances
##
## data: Tsardom and Romanov
## F = 2.003, num df = 12, denom df = 13, p-value = 0.2285
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.6352478 6.4883991
## sample estimates:
## ratio of variances
## 2.003048#Test for normality
shapiro.test(Tsardom)##
## Shapiro-Wilk normality test
##
## data: Tsardom
## W = 0.8838, p-value = 0.0803shapiro.test(Romanov)##
## Shapiro-Wilk normality test
##
## data: Romanov
## W = 0.89445, p-value = 0.09371#wilcoxon test for medians
wilcox.test(x = Tsardom,y = Romanov,exact = F)##
## Wilcoxon rank sum test with continuity correction
##
## data: Tsardom and Romanov
## W = 97.5, p-value = 0.7707
## alternative hypothesis: true location shift is not equal to 0All of the tests prove null hypothesis \(H_0: T_{Romanov} = T_{Tsardom}\).
Bayesian estimation
\[ P(\Theta|D) = P(D|\Theta)P(\Theta)/P(D)\]
Here we produce by MCMC (Markov chain Monte Carlo) procedure the posterior distribution of the duration values parameters \(\Theta\) (means, difference of means and the effect size) given the data \(D\). Then we verify our results using diagnostics.
#Random number generator seed
set.seed(12345)
#Bayesian alternative to the t-test
BESTout <- BESTmcmc(Tsardom, Romanov)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 27
## Unobserved stochastic nodes: 5
## Total graph size: 77
##
## Initializing model
##
## Adaptive phase, 500 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 33334 iterations x 3 chains
##
##
## MCMC took 1.107 minutes.summary(BESTout)## mean median mode HDI% HDIlo HDIup compVal %>compVal
## mu1 17.556 17.53 17.241 95 6.568 28.517
## mu2 15.800 15.80 15.800 95 8.350 23.145
## muDiff 1.755 1.75 1.763 95 -11.663 14.704 0 60.8
## sigma1 18.748 17.98 16.403 95 11.081 28.042
## sigma2 13.266 12.77 12.166 95 8.190 19.388
## sigmaDiff 5.482 5.11 4.500 95 -5.106 16.925 0 86.3
## nu 38.799 30.10 13.950 95 2.094 99.679
## log10nu 1.460 1.48 1.530 95 0.776 2.114
## effSz 0.108 0.11 0.132 95 -0.680 0.902 0 60.8
For the last Romanov Nicholas II we see \(T_{Romanov} = 23.2\) as HIDup for the Bayesian estimate of the mean i.e. 95% probability to fit expected value.
Conclusions
- Both races of Russian rulers (Tsardom, Romanov) have very close means of ruling duration - 17.6 and 15.8 years in the 95% HDI (highest density interval) region \((6.59 < T_{Tsardom} < 28.6, 8.36 < T_{Romanov} < 23.2)\).
- The median of the difference of means \(1.79\) is in the 95% HDI \((-11.7 < T_{Tsardom} - T_{Romanov} < 14.8)\).
- The effect size with the median \(0.12\) is in the 95% HDI region of \((-.697 < Effect < 0.893)\) wthat proves equality for ruling duration means.
- Bayesian approach produces more robust estimation than classical frequentest one.