##### Script for Slidedeck e8
###   Correlation Test


### Preamble

source("http://rfs.kvasaheim.com/stat200.R")

dt = read.csv("http://rfs.kvasaheim.com/data/crime.csv")
attach(dt)


### First model

par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,0,0)+1)
par(font.lab=2,cex.lab=1.1)

plot.new()
plot.window( xlim=c(0,2500), ylim=c(0,2500))

points(vcrime90,vcrime00, pch=21, bg="salmon")

axis(1); axis(2)
title(xlab="Violent Crime Rate (1990)", line=2.25)
title(ylab="Violent Crime Rate (2000)", line=3.00)




xbar=mean(vcrime90)
ybar=mean(vcrime00)
ito = c("green","red")


par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,0,0)+1)
par(font.lab=2,cex.lab=1.1)

plot.new()
plot.window( xlim=c(0,2500), ylim=c(0,2500))

abline(v=xbar, lty=3)
abline(h=ybar, lty=3)

for(i in 1:length(vcrime90)) {
	x = vcrime90[i]
	y = vcrime00[i]
 points(x,y, pch=21, bg=ito[((x>xbar)*2-1)*((y>ybar)*2-1)])
}

axis(1); axis(2)
title(xlab="Violent Crime Rate (1990)", line=2.25)
title(ylab="Violent Crime Rate (2000)", line=3.00)





cov(vcrime90,vcrime00)
cor(vcrime90,vcrime00)






# Three Correlations

par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(0,0,0,0)+1)
par(font.lab=2,cex.lab=1.1)
par(xaxt="n", yaxt="n")
par(mfrow=c(1,3))

nn=2e3
xx = runif(nn, min=-3,max=3)
y1 = +0.8*xx+rnorm(nn,m=0,s=0.5)
y2 = +0.6*xx^2+rnorm(nn,m=0,s=0.5)-2
y3 = -0.8*xx+rnorm(nn,m=0,s=0.5)


plot.new()
plot.window( xlim=c(-5,5), ylim=c(-5,5))
points(xx,y1, pch=20, col=rgb(0,0,0,0.5))
mtext(side=1, at=0, text="Positive")



plot.new()
plot.window( xlim=c(-5,5), ylim=c(-5,5))
points(xx,y2, pch=20, col=rgb(0,0,0,0.5))
mtext(side=1, at=0, text="None")



plot.new()
plot.window( xlim=c(-5,5), ylim=c(-5,5))
points(xx,y3, pch=20, col=rgb(0,0,0,0.5))
mtext(side=1, at=0, text="Negative")










### Example 1: 
#   Property Crime Rates


cor.test(pcrime90,pcrime00)


# graphic
par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,0,0)+1)
par(font.lab=2,cex.lab=1.1)

plot.new()
plot.window( xlim=c(0,8000), ylim=c(0,6000))

for(i in 1:length(pcrime00)) {
  x = pcrime90[i]
  y = pcrime00[i]
  ito = rgb(x/15000+y/15000,0,0)
  points(x,y, pch=20, col=ito)
}

axis(1, at=0:8*1000); axis(2, at=0:3*2000)
title(xlab="Property Crime Rate (1990)", line=2.25)
title(ylab="Property Crime Rate (2000)", line=3.00)






### Example 2:
#   Wealth and Professionalism

cor.test(gspcap00,profleg)


par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,0,0)+1)
par(font.lab=2,cex.lab=1.1)

plot.new()
plot.window( xlim=c(0,80), ylim=c(0,1))

text(0,0.9,label="Size by population\nColor by legislative professionalism", pos=4, cex=0.8)

for(i in 1:length(gspcap00)) {
  x = gspcap00[i]/1000
  y = profleg[i]
  if(is.na(y) || is.na(x)) next
  p = (pop00[i]-pop90[i])/pop90[i]*10
  ito = rgb(1-y,1,1-y)
  points(x,y, pch=21, bg=ito, cex=p)
}

axis(1, at=0:8*10); axis(2, at=seq(0,1,0.2))
title(xlab="GSP per capita (2000) [$1000]", line=2.25)
title(ylab="Legislature Professionalism", line=3.00)







### Example 3:
#   Wealth and Education

cor.test(gspcap00,enroll00)


par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,0,0)+1)
par(font.lab=2,cex.lab=1.1)

plot.new()
plot.window( xlim=c(0,80), ylim=c(0.8,1.2))

text(0,1.15,label="Size by population\nColor by violent crime rate (2000)", pos=4, cex=0.8)

for(i in 1:length(gspcap00)) {
  x = gspcap00[i]/1000
  y = enroll00[i]/100
  if(is.na(y) || is.na(x)) next
  p = pop00[i]/max(pop00)*10
  c = vcrime00[i]/max(vcrime00)
  ito = rgb(1,c,c)
  points(x,y, pch=21, bg=ito, cex=p)
}

axis(1, at=0:14*10); axis(2, at=seq(0,1.6,0.1))
title(xlab="GSP per capita (2000) [$1000]", line=2.25)
title(ylab="Enrollment Percentage", line=3.00)






### Example 4:
#   Wealth and Crime

cor.test(gspcap00,vcrime00)


par(xaxs="i",yaxs="i")
par(family="serif",las=1)
par(mar=c(3,3,0,0)+1)
par(font.lab=2,cex.lab=1.1)

plot.new()
plot.window( xlim=c(0,80), ylim=c(0,1550))

text(0,1500,label="Size by population", pos=4, cex=0.8)

for(i in 1:length(gspcap00)) {
  x = gspcap00[i]/1000
  y = vcrime00[i]
  if(is.na(y) || is.na(x)) next
  p = pop00[i]/max(pop00)*20
  c = c(6,3,"red",4)[census4[i]]
  points(x,y, pch=21, bg="dodgerblue", cex=sqrt(p))
}

axis(1, at=0:14*10); axis(2, at=seq(0,1700,250))
title(xlab="GSP per capita (2000) [$1000]", line=2.25)
title(ylab="Violent Crime Rate (2000)", line=3.00)



### End of File