####################
## Intro to R     ##
## Lecture 5      ##
## Dave Armstrong ##
## ICPSR Sum Prog ##
####################


library(lattice)
library(car)
data(Duncan)

xyplot(prestige ~ income, data=Duncan)

xyplot(prestige ~ income, data=Duncan, pch=16, col="black")


xyplot(prestige ~ income | type, data=Duncan, pch=1, col="black", 
    as.table=T, layout=c(3,1))

xyplot(interlocks ~ log(assets) | nation, groups=sector, data=Ornstein)


xyplot(prestige ~ income, groups=type, data=Duncan)


xyplot(prestige ~ income | type, data=Duncan, pch=1, col="black", 
    as.table=T, panel=function(x,y){
        panel.points(x,y,col="black")
        panel.abline(lm(y ~ x))})


##Page 7
pch.vec <- c(1,2,4)
col.vec <- c("black", "red", "blue")

xyplot(prestige ~ income, data=Duncan, 
    pch=pch.vec[as.numeric(Duncan$type)], 
    col=col.vec[as.numeric(Duncan$type)])

xyplot(prestige ~ income, data=Duncan, groups=type, 
	pch=pch.vec, col=col.vec)
	
	
    pch=pch.vec[as.numeric(Duncan$type)], 
    col=col.vec[as.numeric(Duncan$type)])



xyplot(prestige ~ income, data=Duncan, panel=function(x,y){
    panel.points(x,y, pch=pch.vec[as.numeric(Duncan$type)], 
        col=col.vec[as.numeric(Duncan$type)])
        }
)


mod2 <- lm(prestige~type + income, data=Duncan)
b <- mod2$coef
a_bc <- b[1]
a_prof <- b[1] + b[2]
a_wc <- b[1] + b[3]


## Page 8
xyplot(prestige ~ income, data=Duncan, panel=function(x,y){
    panel.points(x,y, pch=pch.vec[as.numeric(Duncan$type)], 
        col=col.vec[as.numeric(Duncan$type)])
    panel.abline(a=a_bc, b=b[4], lty=1, col="black")
    panel.abline(a=a_prof, b=b[4], lty=2, col="red")
    panel.abline(a=a_wc, b=b[4], lty=3, col="blue")}
)




xyplot(prestige ~ income, data=Duncan, panel=function(x,y){
    panel.points(x,y, pch=pch.vec[as.numeric(Duncan$type)], 
        col=col.vec[as.numeric(Duncan$type)])
    panel.abline(a=a_bc, b=b[4], lty=1, col="black")
    panel.abline(a=a_prof, b=b[4], lty=2, col="red")
    panel.abline(a=a_wc, b=b[4], lty=3, col="blue")
        },
    key=list(space="top", 
        text=list(c("Blue Collar", "Professional", "White Collar")), 
        points=list(pch=pch.vec, col=col.vec), 
        lines=list(lty=c(1,2,3), col=col.vec))
)

xyplot(prestige ~ income, data=Duncan, panel=function(x,y){
    panel.points(x,y, pch=pch.vec[as.numeric(Duncan$type)], 
        col=col.vec[as.numeric(Duncan$type)])
    panel.abline(a=a_bc, b=b[4], lty=1, col="black")
    panel.abline(a=a_prof, b=b[4], lty=2, col="red")
    panel.abline(a=a_wc, b=b[4], lty=3, col="blue")
        },
    key=list(space="top", 
        points=list(pch=pch.vec, col=col.vec), 
        text=list(c("Blue Collar", "Professional", "White Collar")), 
        lines=list(lty=c(1,2,3), col=col.vec))
)


xyplot(prestige ~ income | type, data=Duncan, 
	panel=function(x,y,subscripts){
		panel.points(x, y, pch=1, col="black")
		panel.points(Duncan$education[subscripts], y, pch=2, col="red")})



## Page 9
levels(Duncan$type) <- c("Blue Collar", "Professional", "White Collar")

xyplot(prestige ~ income | type, data=Duncan, as.table=T, 
    panel=function(x,y,subscripts){    
        panel.text(50,50, paste("packet # = ", packet.number(), sep=""), cex=2)
}
)



mod <- lm(prestige ~ income*type, data=Duncan)
s <- seq(min(Duncan$income), max(Duncan$income), length=25)
newdat<- data.frame(
	income = s, 
	type = rep(c("bc", "prof", "wc"), each=25),
	prestige = 50)
	
preds <- predict(mod, newdata=newdat, interval="confidence")
preds <- as.data.frame(preds) 
preds$income <- newdat$income
preds$type <- newdat$type
xyplot(fit ~ income | type, data=preds, type="l", 
	panel = function(x, y, subscripts){
		panel.lines(x,y, col="black")
		panel.lines(x, preds$lwr[subscripts], lty=2, col="black")
		panel.lines(x, preds$upr[subscripts], lty=2, col="black")})


vals <- by(Duncan$income, list(Duncan$type), function(x)x)

xyplot(fit ~ income | type, data=preds, type="l", ylim=c(0,100),
	panel = function(x, y, subscripts){
		panel.lines(x,y, col="black")
		panel.lines(x, preds$lwr[subscripts], lty=2, col="black")
		panel.lines(x, preds$upr[subscripts], lty=2, col="black")
		panel.rug(vals[[packet.number()]], col="black")})



ints <- c(a_bc, a_prof, a_wc)
slope <- b[4]

## Page 10
xyplot(prestige ~ income | type, data=Duncan, as.table=T, 
    panel=function(x,y){    
        panel.points(x,y, pch=1, col="black")
        panel.abline(a=ints[packet.number()], b=slope)
}
)


xyplot(prestige ~ income | type, data=Duncan, as.table=T, 
    panel=function(x,y,subscripts){    
        panel.points(x,y, pch=1, col="black")
        panel.points(Duncan$education[subscripts], y, col="red", pch=2)
} 
)

xyplot(prestige ~ income | type, data=Duncan, as.table=T, 
    xlim=c(0,100), 
    panel=function(x,y,subscripts){    
        panel.points(x,y, pch=1, col="black")
        panel.points(Duncan$education[subscripts], y, col="red", pch=2)
}
)

xyplot(prestige ~ income | type, data=Duncan, as.table=T, 
    xlim=c(0,100), 
    panel=function(x,y,subscripts){    
        panel.points(x,y, pch=1, col="black")
        panel.points(Duncan$education[subscripts], y, col="red", pch=2)
}, 
    key=list(space="top", points=list(pch=c(1,2), col=c("black", "red")), 
        text=list(c("Income", "Education")))
)

hist(Duncan$income)
histogram(~income, 	
	data=Duncan, 
	breaks=seq(0,90,10), 
	col="white")
histogram(~income, 	
	data=Duncan, 
	breaks=h$breaks, 
	col="white")




dens <- density(Duncan$income, from=0, to=100)
plot(dens, type="l")

densityplot(~income, 
    data=Duncan, 
    col = "black")

hist(Duncan$income, 
	freq=F)
lines(dens, lwd=1.5)

xyplot(~income, data=Duncan,
    ylim=c(0,.02), xlim=c(-5,95),
    panel = function(x,y){
        panel.histogram(x, 
           breaks=seq(0,90,10), 
           col="white")
        panel.densityplot(x, 
           col="black", 
           lwd=1.5, 
           plot.points=F)
    })

histogram(~weight | Diet, data=ChickWeight, type="density", 
	ylim=c(0,.01),
	panel = function(x){
		panel.histogram(x, breaks=seq(0,400, 50), type="density")
		panel.densityplot(x, plot.points=F, col="black", lwd=1.5)
})

histogram(~log(assets) | nation, data=Ornstein, type="density", 
	panel = function(x){
	panel.histogram(x, breaks=seq(0,max(log(Ornstein$assets)), length=10), type="density")
	panel.densityplot(x, plot.points=F, col="black", lwd=1.5)
})

## Page 11
library(Rcmdr)

scatter3d(Duncan$income, Duncan$prestige, Duncan$education, surface=F)
scatter3d(Duncan$income, Duncan$prestige, Duncan$education, surface=T, 
    fit="linear")


scatter3d(Duncan$income, Duncan$prestige, Duncan$education, surface=F,
    grou=Duncan$type)

## Page 12
mod3d <- lm(prestige ~ income + education + type, data=Duncan)
summary(mod3d)

scatter3d(Duncan$income, Duncan$prestige, Duncan$education, surface=T, 
    groups=Duncan$type, fit="linear", parallel=T)

scatter3d(Duncan$income, Duncan$prestige, Duncan$education, surface=T, 
    grou=Duncan$type, fit="linear", parallel=F)