Guns
Ritchie & Hart
June 26, 2015
On Thu, Jun 25, 2015 at 11:51 PM, Tyler Ritchie wrote:
require(ggplot2)
di_rank<-c(1:25)
country<-c("Norway", "Australia", "Switzerland", "Netherlands", "United States",
"Germany", "New Zealand", "Canada", "Singapore", "Denmark",
"Ireland", "Sweden", "Iceland", "England", "South Korea",
"Japan", "Israel", "France", "Austria", "Belgium",
"Luxembourg", "Finland", "Slovenia", "Italy", "Spain")
guns_100<-c(31.3, 15.0, 45.7, 3.9, 88.8,
30.3, 22.6, 30.8, 0.5, 12.0,
8.6, 31.6, 30.3, 6.2, 1.1,
0.6, 7.3, 31.2, 30.4, 17.2,
15.3, 45.3, 13.5, 11.9, 10.4)
gun_homicides_100k<-c(0.05, 0.14, 0.77, 0.33, 3.2,
0.19, 0.16, 0.51, 0.02, 0.27,
0.48, 0.41, 0.0, 0.07, 0.03,
0.01, 0.09, 0.06, 0.22, 0.68,
0.62, 0.45, 0.10, 0.71, 0.20)
guns<-data.frame(rank=di_rank, country=country, guns_100=guns_100, gun_murder_100k=gun_homicides_100k)
p<-ggplot(guns, aes(guns_100, gun_murder_100k, label=country))
p + geom_point() + theme_bw() + labs(x= "Guns per 100 people", y="Homicides by gun per 100k people")
NEEDS MORE PATRIOTISM
Tweak the country list for one-word country names:
country[5]<-"USA"
country[7]<-"New-Zealand"
country[14]<-"United-Kingdom"
country[15]<-"South-Korea"
guns<-data.frame(rank=di_rank, country=country, guns_100=guns_100, gun_murder_100k=gun_homicides_100k,row.names=country)Get some flags:
require(RCurl)
require(png)
flags <- list()
for(i in country){
#I expext there's a better way to look up the flag unicodes than this....
html <- getURLContent(paste0("http://emojipedia.org/flag-for-",i))
pt1 <- substr(html,401,401)
pt2 <- substr(html,402,402)
extract_short_name<-function(multibyte_char1,multibyte_char2){
shorten<-function(mb){
out<-character_string<-eval(parse(text=gsub("\\", "", deparse(mb), fixed=TRUE)))
short<-substr(out,5,10)
short}
paste0(shorten(multibyte_char1),"-",shorten(multibyte_char2))}
name<-extract_short_name(pt1,pt2)
#but it works
#so use them to harvest some flags
URL<-paste0("https://raw.githubusercontent.com/Ranks/emojione/master/assets/png/",toupper(name),".png")
flags[[i]] <- readPNG(getURLContent(URL))
}Your plot:
p<-ggplot(guns, aes(guns_100, gun_murder_100k, label=country))
p<-p+geom_point() + theme_bw() + labs(x= "Guns per 100 people", y="Homicides by gun per 100k people")
p
but patrioticker:
for(i in country){
xy <- guns[i, 3:4]
p<-p+annotation_raster(flags[[i]],xmin=xy[[1]], xmax=xy[[1]]+4, ymin=xy[[2]], ymax=xy[[2]]+0.2,interpolate=T)
}
p
WORLDWIDE
Grab data from Guardian
require(dplyr)
require(countrycode)
worldwide<-read.table("~/Downloads/World firearms murders and ownership.csv",h=T,sep=",")
worldwide %>% mutate(country=Country.Territory,
gun_murder_100k=Homicide.by.firearm.rate.per.100.000.pop,
murder_100k=gun_murder_100k/(X..of.homicides.by.firearm/100),
guns_100=Average.firearms.per.100.people,
code=ISO.code,
region=countrycode(code,"iso2c","region")
) %>%
select(region,country,code,gun_murder_100k,murder_100k,guns_100) -> global
rownames(global)<-global$country
require(gapminder)
gapminder %>% filter(year==2007) %>% select(country,continent,gdpPercap,pop) -> gapdata
global <- merge(global,gapdata,by="country",all.x=T)Now we have all sorts of data!
Total Homicide rate
summary(lm(murder_100k~guns_100,data=global)) #ns##
## Call:
## lm(formula = murder_100k ~ guns_100, data = global)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.016 -7.392 -5.124 -1.822 72.551
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.4299 1.8955 5.502 2.89e-07 ***
## guns_100 -0.1500 0.1074 -1.397 0.166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.91 on 100 degrees of freedom
## (83 observations deleted due to missingness)
## Multiple R-squared: 0.01914, Adjusted R-squared: 0.009327
## F-statistic: 1.951 on 1 and 100 DF, p-value: 0.1656p<-ggplot(global, aes(guns_100, murder_100k))
p+geom_point() + theme_bw() + labs(x= "Guns per 100 people", y="Homicides total per 100k people")
Total weighted by population
summary(lm(murder_100k~guns_100,data=global,weights=pop)) #ns##
## Call:
## lm(formula = murder_100k ~ guns_100, data = global, weights = pop)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -129063 -20001 -12285 651 252968
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.38362 1.40041 5.272 1.35e-06 ***
## guns_100 -0.02129 0.04875 -0.437 0.664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 69000 on 72 degrees of freedom
## (111 observations deleted due to missingness)
## Multiple R-squared: 0.002642, Adjusted R-squared: -0.01121
## F-statistic: 0.1907 on 1 and 72 DF, p-value: 0.6636p<-ggplot(global, aes(guns_100, murder_100k))
p+geom_point(aes(size=pop)) + theme_bw() + labs(x= "Guns per 100 people", y="Homicides total per 100k people")
Gun Homicide rate
summary(lm(gun_murder_100k~guns_100,data=global)) #ns##
## Call:
## lm(formula = gun_murder_100k ~ guns_100, data = global)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.684 -4.840 -3.882 -1.932 63.240
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.71826 1.43472 3.986 0.000125 ***
## guns_100 -0.08527 0.08115 -1.051 0.295806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.74 on 105 degrees of freedom
## (78 observations deleted due to missingness)
## Multiple R-squared: 0.0104, Adjusted R-squared: 0.00098
## F-statistic: 1.104 on 1 and 105 DF, p-value: 0.2958p<-ggplot(global, aes(guns_100, gun_murder_100k))
p+geom_point() + theme_bw() + labs(x= "Guns per 100 people", y="Gun homicides total per 100k people")
Gun homicide weighted by population
summary(lm(gun_murder_100k~guns_100,data=global,weights=pop)) #ns##
## Call:
## lm(formula = gun_murder_100k ~ guns_100, data = global, weights = pop)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -114720 -14081 -7790 -1876 198612
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.71496 1.04705 3.548 0.000678 ***
## guns_100 -0.00292 0.03645 -0.080 0.936362
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 51590 on 74 degrees of freedom
## (109 observations deleted due to missingness)
## Multiple R-squared: 8.673e-05, Adjusted R-squared: -0.01343
## F-statistic: 0.006419 on 1 and 74 DF, p-value: 0.9364p<-ggplot(global, aes(guns_100, gun_murder_100k))
p+geom_point(aes(size=pop)) + theme_bw() + labs(x= "Guns per 100 people", y="Gun homicides total per 100k people")
Gun homicides as a % of total homicides
summary(lm((gun_murder_100k/murder_100k)~guns_100,data=global)) #still ns!!##
## Call:
## lm(formula = (gun_murder_100k/murder_100k) ~ guns_100, data = global)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31929 -0.22398 -0.04491 0.20523 0.56772
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.308008 0.034926 8.819 3.82e-14 ***
## guns_100 0.002128 0.001979 1.075 0.285
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2564 on 100 degrees of freedom
## (83 observations deleted due to missingness)
## Multiple R-squared: 0.01143, Adjusted R-squared: 0.001547
## F-statistic: 1.156 on 1 and 100 DF, p-value: 0.2848p<-ggplot(global, aes(guns_100, (gun_murder_100k/murder_100k)))
p+geom_point() + theme_bw() + labs(x= "Guns per 100 people", y="Gun homicides as a % of total homicides")
Weighted by population gun homicides as a % of total homicides
summary(lm((gun_murder_100k/murder_100k)~guns_100,data=global,weights=pop)) #slightly positive##
## Call:
## lm(formula = (gun_murder_100k/murder_100k) ~ guns_100, data = global,
## weights = pop)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -6369.3 -511.2 53.4 928.8 5866.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.250223 0.032783 7.633 7.42e-11 ***
## guns_100 0.004028 0.001141 3.530 0.000729 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1615 on 72 degrees of freedom
## (111 observations deleted due to missingness)
## Multiple R-squared: 0.1475, Adjusted R-squared: 0.1357
## F-statistic: 12.46 on 1 and 72 DF, p-value: 0.0007294summary(lm((gun_murder_100k/murder_100k)~guns_100,data=subset(global,guns_100<50),weights=pop)) #but entirely driven by the US as an outlier##
## Call:
## lm(formula = (gun_murder_100k/murder_100k) ~ guns_100, data = subset(global,
## guns_100 < 50), weights = pop)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -6282.9 -512.9 36.3 913.7 5851.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.243596 0.040618 5.997 7.61e-08 ***
## guns_100 0.004988 0.003619 1.378 0.172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1626 on 71 degrees of freedom
## (101 observations deleted due to missingness)
## Multiple R-squared: 0.02606, Adjusted R-squared: 0.01234
## F-statistic: 1.9 on 1 and 71 DF, p-value: 0.1724p<-ggplot(global, aes(guns_100, (gun_murder_100k/murder_100k)))
p+geom_point(aes(size=pop)) + theme_bw() + labs(x= "Guns per 100 people", y="Gun homicides as a % of total homicides")+geom_smooth(method="lm")
GDP?
p<-ggplot(global, aes(guns_100, murder_100k))
p+geom_point(aes(colour=gdpPercap)) + theme_bw() + labs(x= "Guns per 100 people", y="Homicides total per 100k people")
summary(lm(murder_100k~guns_100,data=subset(global,gdpPercap>20000))) # 0.03 murders/gun by country##
## Call:
## lm(formula = murder_100k ~ guns_100, data = subset(global, gdpPercap >
## 20000))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.97848 -0.40184 -0.09228 0.24864 2.74923
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.638168 0.256568 2.487 0.02056 *
## guns_100 0.030594 0.008985 3.405 0.00243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8378 on 23 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.3351, Adjusted R-squared: 0.3062
## F-statistic: 11.59 on 1 and 23 DF, p-value: 0.00243summary(lm(murder_100k~guns_100,data=subset(global,gdpPercap>20000),weights=pop)) #0.05 murders/gun by population##
## Call:
## lm(formula = murder_100k ~ guns_100, data = subset(global, gdpPercap >
## 20000), weights = pop)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -9826.0 -1984.9 53.5 1627.4 14504.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.288906 0.237937 1.214 0.237
## guns_100 0.050032 0.004327 11.562 4.63e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4591 on 23 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.8532, Adjusted R-squared: 0.8468
## F-statistic: 133.7 on 1 and 23 DF, p-value: 4.626e-11p<-ggplot(subset(global,gdpPercap>20000), aes(guns_100, murder_100k))
p+geom_point(aes(colour=gdpPercap)) + theme_bw() + labs(x= "Guns per 100 people", y="Homicides total per 100k people")+geom_smooth(method="lm")
summary(lm(murder_100k~guns_100,data=subset(global,gdpPercap>20000&guns_100<50))) #but still drivein entirely by outliers status of USA##
## Call:
## lm(formula = murder_100k ~ guns_100, data = subset(global, gdpPercap >
## 20000 & guns_100 < 50))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76035 -0.43817 -0.10847 0.08302 2.41240
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.098011 0.247976 4.428 0.000212 ***
## guns_100 0.002637 0.010867 0.243 0.810547
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6866 on 22 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.002669, Adjusted R-squared: -0.04266
## F-statistic: 0.05886 on 1 and 22 DF, p-value: 0.8105p<-ggplot(subset(global,gdpPercap>20000&guns_100<50), aes(guns_100, murder_100k))
p+geom_point(aes(colour=gdpPercap)) + theme_bw() + labs(x= "Guns per 100 people", y="Homicides total per 100k people")+geom_smooth(method="lm")
Continet?
subset(global,gdpPercap>20000)->global
p<-ggplot(global, aes(guns_100, murder_100k,colour=continent))
p<-p+geom_text(aes(label=code)) + theme_bw() + labs(x= "Guns per 100 people", y="Homicides total per 100k people")
p + geom_smooth(method="lm",se=F) ## Warning: Removed 4 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 5 rows containing missing values (geom_text).
Get flags
names<-as.character(global$country)
names<-gsub(" ", "-",names)
names[30]<-"USA"
row.names(global)<-names
allflags <- list()
for(i in names){
#I expext there's a better way to look up the flag unicodes than this....
html <- getURLContent(paste0("http://emojipedia.org/flag-for-",i))
pt1 <- substr(html,401,401)
pt2 <- substr(html,402,402)
extract_short_name<-function(multibyte_char1,multibyte_char2){
shorten<-function(mb){
out<-character_string<-eval(parse(text=gsub("\\", "", deparse(mb), fixed=TRUE)))
short<-substr(out,5,10)
short}
paste0(shorten(multibyte_char1),"-",shorten(multibyte_char2))}
name<-extract_short_name(pt1,pt2)
#but it works
#so use them to harvest some flags
URL<-paste0("https://raw.githubusercontent.com/Ranks/emojione/master/assets/png/",toupper(name),".png")
allflags[[i]] <- readPNG(getURLContent(URL))
}Total homicides
p<-ggplot(global, aes(guns_100, murder_100k),ylim=c(0,5.5),xlim=c(0,90))
p<-p+geom_point() + theme_bw() + labs(x= "Guns per 100 people", y="Homicides total per 100k people")
for(i in names){
xy <- global[i, 6:5]
p<-p+annotation_raster(allflags[[i]],xmin=xy[[1]], xmax=xy[[1]]+4, ymin=xy[[2]], ymax=xy[[2]]+0.2,interpolate=T)
}
p
Gun homicides
p<-ggplot(global, aes(guns_100, gun_murder_100k),ylim=c(0,5.5),xlim=c(0,90))
p<-p+geom_point() + theme_bw() + labs(x= "Guns per 100 people", y="Homicides by gun per 100k people")
for(i in names){
xy <- global[i, c(6,4)]
p<-p+annotation_raster(allflags[[i]],xmin=xy[[1]], xmax=xy[[1]]+4, ymin=xy[[2]], ymax=xy[[2]]+0.2,interpolate=T)
}
p
(Emoji set designed and offered free by Emoji One)