ME <- readxl::read_xls("~/Yoga/ViolentCrimeStatistics.xls")
names(ME) <- ME[3, ]
ME <- ME[c(4:23), ]
ME$Year %<>% substr(., 1, 4)
ME <- ME[, -c(21, 22)]
names(ME) %<>% gsub("\\\n", " ", .)
ME %<>% mutate_all(funs(as.numeric))
ME %<>% mutate(VCRate = `Violent crime rate`/100000) %>% select(Year, VCRate, everything())
ME$VCRate %>% diff %>% c(0, .) %>% cbind(ME$Year) %>% as.data.frame() %>% setNames(object = ., 
    c("VCd", "Year")) %>% ggplot(data = ., mapping = aes(x = Year, y = VCd)) + geom_line() + 
    labs(title = "Derivative of violent crime rate by Year", subtitle = "", caption = "", 
        x = "Year", y = "dViolent Crime / dYear") + theme(plot.title = element_text(hjust = 0.5), 
    plot.subtitle = element_text(hjust = 0.5))

The derivative of the rate of change shows that decreases in crime were actually more dramatic between 1991 and 1995 than between 1996 & 2010 when the Maharishi Study specifies.

ME$Fac <- factor(rep(1:4, each = 5))


ME_lm <- lm(VCRate ~ gam::s(Year):Fac, data = ME)
ME %>% mutate(PredVC = predict(ME_lm)) %>% ggplot(data = ., mapping = aes(x = Year, 
    y = VCRate)) + geom_line() + geom_line(aes(y = PredVC), color = "red") + geom_smooth(method = "lm", 
    aes(fill = Fac)) + labs(title = "Splined Regression of Violent Crimes by 5 Year Groupings", 
    subtitle = "Groupings According to the Study", caption = "Black is actual rates \n Red line is prediction from Generalized Additive Splines Model", 
    x = "Years", y = "Violent Crime Rate") + theme(plot.title = element_text(hjust = 0.5), 
    plot.subtitle = element_text(hjust = 0.5))

1996, when the TM group reached it’s critical point - as claimed by the study, did indeed show an inflection from a previous short lived rise in violent crime back towards a downward trend.

Two Way Anova to determine if there are significant differences between rates across year groupings.

car::Anova(lm(VCRate ~ Fac * Year, data = ME, contrasts = "contr.poly"), type = 2, 
    test.statistic = "F", digits = 3, icontrasts = "contr.poly") %T>% assign("anova_maha", 
    ., envir = .GlobalEnv) %>% stargazer::stargazer(summary = F, type = "html")
Sum Sq Df F value Pr(> F)
Fac 0.00000 3 76.513 0.00000
Year 0.00000 1 183.579 0
Fac:Year 0.00000 3 11.071 0.001
Residuals 0.00000 12

We can conclude that there are interactions between groups of years

Pairwise comparisons of Violent crime rates across year groupings

ME_lm %>% lsmeans::lsmeans(., pairwise ~ Fac, adjust = "tukey")
## $lsmeans
##  Fac  lsmean       SE df lower.CL upper.CL
##  1   0.00576 0.000215 15  0.00530  0.00622
##  2   0.00517 0.000101 15  0.00495  0.00538
##  3   0.00534 0.000101 15  0.00512  0.00555
##  4   0.00605 0.000213 15  0.00560  0.00651
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast  estimate       SE df t.ratio p.value
##  1 - 2     0.000592 0.000172 15  3.453  0.0167 
##  1 - 3     0.000425 0.000288 15  1.475  0.4754 
##  1 - 4    -0.000291 0.000414 15 -0.704  0.8941 
##  2 - 3    -0.000167 0.000171 15 -0.980  0.7629 
##  2 - 4    -0.000884 0.000287 15 -3.083  0.0342 
##  3 - 4    -0.000716 0.000170 15 -4.210  0.0038 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates

As indicated in the study, group 3 years prior to the critical mass of the TM group compared to group 4 years 2006-2010 during the critical mass of the TM group, there is a significant change in slope (p < .01). However, when comparing group 4 years to group 1 years there is no significant difference (p = .8941) showing that the downward trend was already precedent prior to the critical mass of the TM Group.