# R source file for the solutions to problems given in Chapter 3
# Time-stamp: <2008-02-18 08:24:28 chl>

##############################################################################
n <- 100
x <- rnorm(n)
y <- 2*x+rexp(n,1)
plot(x,y)
points(x[x<1 & y>5],y[x<1 & y>5],pch=19)
ols <- lm(y~x)
abline(ols)
summary(ols)
library(quantreg)
lad <- rq(y~x)
abline(lad$coefficients[1],lad$coefficients[2],col=2)
summary(lad)

##############################################################################
library(MASS)
hm <- rlm(y~x)
abline(hm$coefficients[1],hm$coefficients[2],col=3)
summary(hm)
#library(lqs)  # library lqs is now included in the MASS package
lts <- ltsreg(y~x)
abline(lts$coefficients[1],lts$coefficients[2],col=4)
legend.text <-
c(paste("OLS",round(ols$coefficients[2],3)),paste("LAD",round(lad$coefficients[2],3)),
paste("Huber",round(hm$coefficients[2],3)),paste("LTS",round(lts$coefficients[2],3)))
legend("topleft",legend.text,lty=1,col=1:4)
title(main="Regression lines and Slopes estimates")

##############################################################################
n <- 100
true.sdx <- 2.5
true.slope <- 3
x <- rnorm(n,23,true.sdx)
y <- true.slope*x+rnorm(n)
sdx <- sd(x)
sdy <- sd(y)
cor.xy <- cor(x,y)
cor.xy*((n-1)/n*sdy)/((n-1)/n*sdx)
cor.xy*sdy/sdx
summary(lm(y~x))$coefficients[[2]]
plot(x,y)

##############################################################################
n <- 100
x <- runif(n,-10,10)
expected.y <- x^2
y <- x^2+rnorm(n)
cor(x,y)
plot(x,y)

##############################################################################
sim <- function(n=100) {
    x <- runif(n,-10,10)
    y <- x^2+rnorm(n)
    return(cor(x,y))
}

##############################################################################
set.seed(12345)
res <- replicate(100,sim())
summary(res)
hist(res)

##############################################################################
data(anscombe)
var.names <- colnames(anscombe)
par(mfrow=c(2,2))
for (i in 1:4)
	plot(anscombe[,i],anscombe[,i+4],xlab=var.names[i],ylab=var.names[i+4])

##############################################################################
HIV <- matrix(c(3,2,4,22),nrow=2)
dimnames(HIV) <- list(c("Cases","Noncases"),c("Exposed","Unexposed"))
HIV
apply(HIV,1,sum)
apply(HIV,2,sum)

##############################################################################
er <- (HIV[1,1]/sum(HIV[,1]))-(HIV[1,2]/sum(HIV[,2]))
rr <- (HIV[1,1]/sum(HIV[,1]))/(HIV[1,2]/sum(HIV[,2]))
or <- (HIV[1,1]*HIV[2,2])/(HIV[2,1]*HIV[1,2])

##############################################################################
ESOP <- matrix(c(255,9,520,191),nrow=2)
dimnames(ESOP) <- list(c("Control","Case"),c("0-9 g/day","10+ g/day"))
ESOP

##############################################################################
(ESOP[1,1]*ESOP[2,2])/(ESOP[1,2]*ESOP[2,1])

##############################################################################
library(vcd)
or <- oddsratio(ESOP,log=FALSE)
or
confint(or)
plot(or)
fourfold(or)

##############################################################################
library(epicalc)
cci(191,520,9,255,design="case-control")

##############################################################################
LEUK <- read.table("leuk.txt",header=TRUE)
library(survival)
LEUK.surv <- Surv(LEUK$time,LEUK$cens)
LEUK.km <- survfit(LEUK.surv~group,data=LEUK)
plot(LEUK.km[1])
lines(LEUK.km[2],col="gray60")
legend("topright",c("6-MP","Placebo"),lty=1,col=c("black","gray60"),bty="n")
title(xlab="Weeks since randomization",ylab="Proportion in remission")