library(bootstrap)

# estimating prediction error through bootstrap and cross validation for hormone data

# data

z <- c(99, 152, 293, 155, 196, 53, 184, 171, 52, 376, 385, 402, 29, 76, 296, 151, 177, 209, 119, 188, 115, 88, 58, 49, 150, 107, 125)
l <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0)

length(z)
length(l)

p <- t(rbind(z,l))		# rbind combinds two vectors into rows
						# t give the transpose.
dim(p)
                                      
y <- c(25.8, 20.5, 14.3, 23.2, 20.6, 31.1, 20.9, 20.9, 30.4, 16.3, 11.6, 11.8, 32.5, 32.0, 18.0, 24.1, 26.5, 25.8, 28.8, 22.0, 29.7, 28.9, 32.8, 32.5, 25.4, 31.7, 28.5)
length(y)

# functions

theta.fit <- function(x, y){
	lsfit(x,y)
}

theta.predict <- function (fit, x){
	cbind(1, x)%*%fit$coef					# There was an error here x was z
}

sq.err <- function(y, yhat){
	(y-yhat)^2
}

# Regressions with one x variable, z

fit.out <- theta.fit(z,y)
fit.out

#bootstrap

result1<- bootpred(z,y,20,theta.fit,theta.predict,err.meas=sq.err)

miss.clas <- function(y,yhat){
	1*(yhat!=y)
}

#cross validation

result2<- crossval(z,y,theta.fit,theta.predict, ngroup = 27)

# Regressions with two x variables, p

fit.out <- theta.fit(p,y)
fit.out

#bootstrap

result1<- bootpred(p,y,20,theta.fit,theta.predict,err.meas=sq.err)

miss.clas <- function(y,yhat){
	1*(yhat!=y)
}

#cross validation

result2<- crossval(p,y,theta.fit,theta.predict,ngroup=10)		# ngroup can't be very large for this dataset and multiple regression.