Regression and Other Stories: Human Development Index

Human Development Index - Looking at data in different ways. See Chapter 2 in Regression and Other Stories.

---

Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("foreign")
library("maps")

Load data

hdi <- read.table(root("HDI/data","hdi.dat"), header=TRUE)
head(hdi)

  rank            state   hdi canada.dist
1    1      Connecticut 0.962           2
2    2    Massachusetts 0.961           2
3    3       New Jersey 0.961           2
4    4 Washington, D.C. 0.960           4
5    5         Maryland 0.960           3
6    6           Hawaii 0.959           2

votes <- read.dta(root("HDI/data","state vote and income, 68-00.dta"))
head(votes)

  st_fips st_year   st_state st_stateabb st_total st_dem st_rep st_repshare
1       1    1968    Alabama          AL     1050    197    147   0.4273256
2       2    1968     Alaska          AK       83     35     38   0.5205479
3       4    1968    Arizona          AZ      487    171    267   0.6095890
4       5    1968   Arkansas          AR      610    185    189   0.5053476
5       6    1968 California          CA     7252   3244   3468   0.5166866
6       8    1968   Colorado          CO      811    335    409   0.5497312
  st_demshare st_income st_inc10k
1   0.5726744  11662.76  1.166276
2   0.4794520  20388.45  2.038845
3   0.3904110  14864.88  1.486488
4   0.4946524  11124.08  1.112408
5   0.4833135  19362.40  1.936240
6   0.4502688  15724.20  1.572420

Pre-process

income2000 <- votes[votes[,"st_year"]==2000, "st_income"]
state.income <- c(income2000[1:8],NA,income2000[9:50])
state.abb.long <- c(state.abb[1:8],"DC",state.abb[9:50])
state.name.long <- c(state.name[1:8],"Washington, D.C.",state.name[9:50])
hdi.ordered <- rep(NA, 51)
can <- rep(NA, 51)
for (i in 1:51){
  ok <- hdi[,"state"]==state.name.long[i]
  hdi.ordered[i] <- hdi[ok,"hdi"]
  can[i] <- hdi[ok,"canada.dist"]
}
no.dc <- state.abb.long != "DC"

Plot average state income and Human Development Index

par(mar=c(3,3,2.5,1), mgp=c(1.5,.2,0), tck=-.01, pty="s")
plot(state.income, hdi.ordered,
     xlab="Average state income in 2000", ylab="Human Development Index", type="n")
text(state.income, hdi.ordered, state.abb.long)

Plot rank of average state income and Human Development Index

par(mar=c(3,3,2.5,1), mgp=c(1.5,.2,0), tck=-.01, pty="s")
plot(rank(state.income[no.dc]), rank(hdi.ordered[no.dc]),
     xlab="Rank of average state income in 2000", ylab="Rank of Human Development Index", type="n")
text(rank(state.income[no.dc]), rank(hdi.ordered[no.dc]), state.abb)

print(cor(rank(hdi.ordered[no.dc]),rank(state.income[no.dc])), digits=2)

[1] 0.86

Plot a map of Human Devlopment index

statemaps <- function(a, grayscale=FALSE, ...){
  if (length(a)==51){
    no.dc <- c(1:8,10:51)
    a <- a[no.dc]
  }
  if (length(a)==50){
    lower48 <- state.abb!="AK" & state.abb!="HI"
    a <- a[lower48]
  }
  else if (length(a)!=48) stop("wrong number of states")
  mapping <- list(1,2,3,4,5,6,7,9,10,11,12,13,14,15,16,17,18,19,20:22,23:24,25,26,27,28,29,30,31,32,33,34:37,38:40,41,42,43,44,45,46,47,48,49,50,51,52,53:55,56:60,61,62,63)
  # for (i in 1:length(mapping)) print(regions[mapping[[i]]])
  a.long <- rep(NA, 63)
  projection <- "bonne"
  for (i in 1:48){
    a.long[mapping[[i]]] <- a[i]
  }
 if (grayscale){
   a.long.scaled <- .95*(a.long-min(a,na.rm=TRUE))/(max(a,na.rm=TRUE)-min(a,na.rm=TRUE))
   shades <- a.long.scaled
   not.dc <- !is.na(a.long.scaled)
   shades[not.dc] <- gray(shades[not.dc])
   map('state', proj=projection, param=25, lty=0, ...)
   m <- map('state', proj=projection, param=25, fill=TRUE, plot=FALSE)
   polygon(m$x,m$y, col=shades, lwd=0.5, border="gray30")
 }
 else {
   map('state', proj=projection, param=25, lty=0, ...)
   m <- map('state', proj=projection, param=25, fill=TRUE, plot=FALSE)
   polygon(m$x,m$y, col=a.long, lwd=0.5, border="gray30")
 }
}

par(mar=c(0,0,0,0))
statemaps(ifelse (can==1, "darkgreen", ifelse (can==2, "green4",
                     ifelse (can==3, "green3", ifelse (can==4, "green2",
                     ifelse (can==5, "palegreen2", "yellowgreen"))))))
mtext("Human Development Index by State", line=-2)
mtext("(Colors indicate number of states you need to drive through to reach the Canadian border.)",
      side=1, cex=.75, line=-2)