* R code  for Robust Statistics homework problems;

Copy and paste code one problem part at a time.

The easiest way to get rpack.txt and robdata.txt is to copy and paste the following two commands into R.

source("http://parker.ad.siu.edu/Olive/rpack.txt")
source("http://parker.ad.siu.edu/Olive/robdata.txt")

#The following commands are sometimes used.

library(glmnet)
library(leaps)

You will sometimes need packages and may need to install packages
the first time you use a given computer.

install.packages("glmnet")
install.packages("leaps")


If rpack.txt and robdata.txt are on G drive, get rpack and robdata with
the following commands:

source("G:/rpack.txt") 
source("G:/robdata.txt")


#Problem 1.2 R commands 

#a)

X <- matrix(rnorm(300),nrow=100,ncol=3)
Y <- (X %*% 1:3)^3 + rnorm(100)
bols <- lsfit(X,Y)$coef[-1]
plot(X %*% bols, Y)

#c)

plot(10*X %*% bols, Y)
title("Scaled OLS View")

#Problem 1.3 R commands 

#a)

cltsim <- function(n=100, nruns=100){
ybar <- 1:nruns
for(i in 1:nruns){
  ybar[i] <- mean(rexp(n))}
list(ybar=ybar)}

#b)

z1 <- cltsim(n=1)
z5 <- cltsim(n=5)
z25 <- cltsim(n=25)
z200 <- cltsim(n=200)
par(mfrow=c(2,2))
hist(z1$ybar)
hist(z5$ybar)
hist(z25$ybar)
hist(z200$ybar)

#d) 

cltsim <- function(n=100, nruns=100){
ybar <- 1:nruns
for(i in 1:nruns){
  ybar[i] <- mean(rnorm(n))}
list(ybar=ybar)}

z1 <- cltsim(n=1)
z5 <- cltsim(n=5)
z25 <- cltsim(n=25)
z200 <- cltsim(n=200)
par(mfrow=c(2,2))
hist(z1$ybar)
hist(z5$ybar)
hist(z25$ybar)
hist(z200$ybar)

par(mfrow=c(1,1))


#Problem 2.22 R commands

height <- rnorm(87, mean=1692, sd = 65)
height[61:65] <- 19.0 
 
par(mfrow=c(2,2))
plot(height)
title("a) Dot plot of heights")
hist(height)
length(height)
val <- quantile(height)[4] - quantile(height)[2]
val
wid <- 4*1.06*min(sqrt(var(height)),val/1.34)*(87^(-1/5))
wid
dens<- density(height,width=wid)
plot(dens$x,dens$y)
lines(dens$x,dens$y)
title("c) Density of heights")
boxplot(height)
title("d) Boxplot of heights")

par(mfrow=c(1,1))

#Problem 2.23 R commands

#a)

y <- rnorm(10000) 
mad(y, constant=1)

#b)

y <- rexp(10000) 
mad(y, constant=1)


#Problem 2.24 R commands 

tmn <-function(x, tp = 0.25){
mean(x, trim = tp)}

#a)

y <- rnorm(10000) 
tmn(y)

#b)  

y <- rexp(10000) 
mean(y)
tmn(y)

#Problem 2.25 R commands 

metmn <-function(x, k = 6){
madd <- mad(x, constant = 1)
med <- median(x)
mean(x[(x >= med - k * madd) & (x <= med + k * madd)])}

y <- rnorm(10000) 
metmn(y)

#Problem 2.26 R commands 

y <- rnorm(10000) 
ratmn(y)

#Problem 2.27 R commands 

y <- rnorm(10000) 
rstmn(y)

#Problem 2.28 R commands 

cci(height)

#Problem 2.29 R commands 

medci(height)


#Problem 2.30 R commands 

tmci(height)

#Problem 2.36 R commands 

qplot<-function(y){	
plot(sort(y), ppoints(y))
title("QPLOT")}

#a)

#height from Problem 2.22 commands:
height <- rnorm(87, mean=1692, sd = 65)
height[61:65] <- 19.0 
 
qplot(height)

#b)

Y <- rnorm(1000)
qplot(y)


#Problem 2.37 R commands

#a)  Copy and paste the rpack source command from the top of the file.


#bi)

rcisim(n=10,type=1)    #Don't forget the source command!

rcisim(n=50,type=1) 

rcisim(n=100,type=1) 

#b ii)

rcisim(n=10,type=2) 

rcisim(n=50,type=2) 

rcisim(n=100,type=2) 


#######


##Problem 3.26, needs rpack

#b)

simx2 <- matrix(rnorm(200),nrow=100,ncol=2)
outx2 <- matrix(10 + rnorm(80),nrow=40,ncol=2)
outx2 <- rbind(outx2,simx2)
maha(outx2)


##Problem 3.27, needs rpack and get outx2 from Problem 3.26

library(MASS)
rmaha(outx2)

#Problem 3.28, needs rpack

rcovsim(100)

rcovsim(100)

rcovsim(100)


###R commands for 3.29, 3.30, and 3.31

##copy and paste the two source commands from near the top of this file.

#source("J:/rpack.txt")  # assumes flash drive is on drive J
%source("J:/robdata.txt") 

library(MASS)


##Problem 3.29, needs rpack and robdata

#a)

library(MASS)
ddcomp(buxx)  #left click on the outliers and right click Stop 4 times 

#b)

ddcomp(cbrainx)  #left click on the outliers and right click Stop 4 times 

#c)

ddcomp(museum[,-1])  #left click on the outliers and right click Stop 4 times 


##Problem 3.30 R commands, needs rpack

#repeat the following command several times

concmv() #right click Stop 5 times


##Problem 3.31, needs rpack

ddmv(p=2,gam=.4)     #right click Stop 5 times
ddmv(p=4,gam=.35)
ddmv(p=4,gam=.30)


##Problem 3.32, needs rpack, robdata

cmba2(buxx) #right click Stop on the plot 7 times

cmba2(cbrainx) #right click Stop 7 times

cmba2(museum[,-1]) # right click Stop 7 times


##Problem 3.33, needs rpack

#a)

covesim(n=100,p=4)

#b) 

covesim(n=100,p=4,outliers=1,pm=15)


##Problem 3.34, needs robdata

library(MASS)

z <- cbind(cbrainx,cbrainy)
z<-z[,-c(2,4,10)]
dim(z)
#Get the log determinant of the FMCD estimator
#applied to a permutation of the data. If the MCD
#estimator was computed, the log determinant would
#be a minimum and would not change under a permutation
#of the rows of the data matrix.
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit


###HW commands for 3.35-3.42

#source("G:/rpack.txt")

library(MASS)


##Problem 3.35, needs rpack

#b)

library(MASS)
ddsim(n=20,p=2)
ddsim(n=30,p=2)

#if the plots are too fast use
ddsim(n=20,p=2,ident=T) #rightclick stop 10 times
ddsim(n=30,p=2,ident=T) #rightclick stop 10 times



##Problem 3.36, needs rpack

library(MASS)
corrsim(n=20,p=2,nruns=10)
corrsim(n=30,p=2,nruns=10)


##Problem 3.37, needs rpack

#b)

library(MASS)
n <- 400
p <- 3
eps <- 0.4
x <- matrix(rnorm(n * p), ncol = p, nrow = n)
zu <- runif(n)
x[zu < eps,] <- x[zu < eps,]*5

#c)

ddplot(x) #right click Stop once on the plot


##Problem 3.38, needs rpack

#b)

#plot for classical covering ellipsoid
simx2 <- matrix(rnorm(200),nrow=100,ncol=2)
outx2 <- matrix(10 + rnorm(80),nrow=40,ncol=2)
outx2 <- rbind(outx2,simx2)
ellipse(outx2)

#plot for RMVN covering ellipsoid
zout <- covrmvn(outx2)
ellipse(outx2,center=zout$center,cov=zout$cov)


##Problem 3.39, needs rpack 

#c)
# get x using Problem 3.37


mplot(x)  #right click Stop  once


##Problem 3.40, needs rpack

#c)

# get x using Problem 3.37

wddplot(x)  


##Problem 3.41, needs rpack

#source("G:/robdata.txt")
ddplot4(buxx,alpha=0.2) #right click Stop once

## Problem 3.42 a)

#Copy and paste two source commands and the three library commands
#from near the top of this file into R.
#Need glmnet library

out <- cv.glmnet(buxx,buxy,alpha=1) #lasso
                   lam <- out$lambda.min
                   fit <- predict(out,s=lam,newx=buxx)
res <- buxy - fit
vars <- 1:4
lcoef <- predict(out,type="coefficients",s=lam)
                lcoef<-as.vector(lcoef)[-1]
                vin <- vars[lcoef!=0]
                pp <- length(vin)+1
AERplot2(yhat=fit,y=buxy,res=res,d=pp)
title("lasso")

#Problem 3.42 b) 

ind <-getB(buxx)$indx
yb <- buxy[ind]
xb <- buxx[ind,]
xnotinb <- buxx[-ind,]
ynotinb <- buxy[-ind]
out <- cv.glmnet(xb,yb,alpha=1) #lasso
                   lam <- out$lambda.min
                   fit <- predict(out,s=lam,newx=xb)
       fit2 <- predict(out,s=lam,newx=xnotinb)
yy <- c(yb,ynotinb)
yhat <- c(fit,fit2)
res <- yy - yhat
vars <- 1:4
lcoef <- predict(out,type="coefficients",s=lam)
                lcoef<-as.vector(lcoef)[-1]
                vin <- vars[lcoef!=0]
                pp <- length(vin)+1
AERplot2(yhat=yhat,y=yy,res=res,d=pp)
title("lasso using covmb set B")

#Problem 3.42 c)

ddplot5(buxx)


##Problem 3.43, needs rpack

## Problem 3.43 a) 

mldsim6(n=100,p=10,outliers=1,gam=0.25,pm=20, osteps=9)

mldsim6(n=100,p=45,outliers=1,gam=0.25,pm=60, osteps=9)

#Problem 1.15 b) #takes a few minutes since p = 1000

mldsim7(n=100,p=1000,outliers=2,gam=0.4,pm=900,osteps=0)

mldsim7(n=100,p=1000,outliers=2,gam=0.4,pm=900,osteps=9)

#Problem 3.43 c)

mldsim7(n=100,p=100,outliers=3,gam=0.25,pm=8,osteps=0)

mldsim7(n=100,p=100,outliers=3,gam=0.25,pm=8,osteps=9)


##Problem 3.44, needs rpack


#a) diagnostic for sphericity

x <- matrix(rnorm(1000), nrow=100, ncol=10)
x <- 10*x
MDSQ <- mahalanobis(x,mean(x),var(x))
DSQ <- mahalanobis(x,mean(x),diag(10))
plot(MDSQ,DSQ)
#abline(0,100)


##Problem 4.1, needs rpack

##do not forget the two source commands near the top of this file


ddplot4(buxx,alpha=0.2) #right click Stop once


##Problem 4.2, needs rpack

predsim()


##do not forget the two source commands near the top of this file


##Problem 4.3, needs rpack

#a) 

predsim2(n=200,nv=100,p=10,nruns=5000,xtype=3,dtype=1)
qchisq(0.95,df=10)

#b)

predsim2(n=200,nv=100,p=10,nruns=5000,xtype=1,dtype=1)


#c) Takes a little while.

predsim2(n=50,nv=20,p=100,nruns=5000,xtype=1,dtype=1)


#Problem 5.8 R commands, needs rpack 

#b) 

nx <- matrix(rnorm(300),nrow=100,ncol=3)
y <- exp( 4 + nx%*%c(1,1,1) + 0.5*rnorm(100) )

#c) 

tplot(nx,y)

#d)

out <- lsfit(nx,log(y))
ls.print(out)

#Problem 5.9 R commands, needs rpack and robdata

#b)

pisim(n=100, type = 1)

#c)

pisim(n=100, type = 3)

#d)

piplot(cbrainx,cbrainy)


#Problem 5.10 R commands, needs rpack and robdata

#a)

cmar <- par("mar")
par(mfrow = c(2, 1))
par(mar=c(4.0,4.0,2.0,0.5))
#main command is below

MLRplot(buxx,buxy)

#restore default

par(mfrow=c(1,1))
par(mar=cmar)

#c)

library(MASS)
rrplot2(buxx,buxy)

#d)

ffplot2(buxx,buxy)


#Problem 5.11 R commands, needs rpack and robdata

#b)

diagplot(buxx,buxy)


#Problem 5.12 needs rpack

##do not forget the two source commands near the top of this file

par(mfrow=c(2,2))
olskout(n=10)
title("a) n=10, k=1")
olskout(n=100)
title("b) n=100, k=1")
olskout(n=1000)
title("c) n=1000, k=1")
olskout(n=10000,k=40)
title("d) n=10000, k=40")
par(mfrow=c(1,1))

#olskout(n=10,k=4)
#olskout(n=100,k=4)
#olskout(n=1000,k=4)
#olskout(n=10000,k=4)
#olskout(n=10000,k=400)
#olskout(n=100000,k=400)
#olskout(n=10,pm=40,k=6)
#olskout(n=10,pm=400,k=1)
#olskout(n=10,pm=100,k=1)
#olskout(n=100,pm=400,k=1)

######


#Problem 6.1 R commands, needs rpack

b)

nltv(0.5)
nltv(0.75) 
nltv(0.9)
nltv(0.9999)

#Problem 6.2 R commands, needs rpack

b)

deltv(0.5)
deltv(0.75) 
deltv(0.9)
deltv(0.9999)

#Problem 6.3 R commands, needs rpack

b)

cltv(0.5)
cltv(0.75) 
cltv(0.9)
cltv(0.9999)

#Problem 6.4 R commands, needs rpack and robdata

b)

mbamv(belx,bely)

c)

#for better looking plots than the default

cmar <- par("mar")
par(mfrow = c(2, 1))
par(mar=c(4.0,4.0,2.0,0.5))
#main command is below

mbamv2(buxx,buxy)

#restore default

par(mfrow=c(1,1))
par(mar=cmar)

#Problem 6.5 R commands, needs rpack and robdata

b)

mlrplot2(belx,bely)

c)

mlrplot2(cbrainx,cbrainy)

d)

mlrplot2(museum[,3:11],museum[,2])

e)

mlrplot2(buxx,buxy)

f)

#use the following command several times
mlrplot2(hx,hy)

library(MASS) 
ffplot2(hx,hy)



##Problem 6.6 R commands, needs rpack and robdata

library(MASS)

tvreg(buxx,buxy,ii=1)

tvreg2(buxx,buxy, M = 0) #use the M that you find, M= 0 is bad


###########


##Problem 7.15 R command, needs rpack 

regbootsim3(nruns=500)


##Problem 7.16 R commands, needs rpack

#make sure you installed package leaps with command
#install.packages("leaps")

#Copy and paste two source commands from near the top of this file into R.

#a)

library(leaps)
n<- 1000
p<-4
q<-p-1
b <- 0 * 1:q
b[1] <- 1 #beta = (1,1,0,0)^T
x <- matrix(rnorm(n * q), nrow = n, ncol = q)
y <- 1 + x %*% b + rnorm(n)
outols <- lsfit(x,y)
outols$coef  #Tn = betahat

#b) 

tem<-regboot(x,y)
betas <- tem$betas
betas[1:5,]

#c) 

apply(betas,2,mean)

#d) 

temp<-vselboot(x,y)
vsbetas <- temp$betas
vsbetas[1:5,]

#e)

apply(vsbetas,2,mean)




##Problem 7.19 R commands 


## Problem 7.19 a)

#Copy and paste two source commands from near the top of this file into R.

library(glmnet)

out <- cv.glmnet(buxx,buxy,alpha=1) #lasso
                   lam <- out$lambda.min
                   fit <- predict(out,s=lam,newx=buxx)
res <- buxy - fit
vars <- 1:4
lcoef <- predict(out,type="coefficients",s=lam)
                lcoef<-as.vector(lcoef)[-1]
                vin <- vars[lcoef!=0]
                pp <- length(vin)+1
AERplot2(yhat=fit,y=buxy,res=res,d=pp)
title("lasso")

#Problem 7.19 b) 

ind <-getB(buxx)$indx
yb <- buxy[ind]
xb <- buxx[ind,]
xnotinb <- buxx[-ind,]
ynotinb <- buxy[-ind]
out <- cv.glmnet(xb,yb,alpha=1) #lasso
                   lam <- out$lambda.min
                   fit <- predict(out,s=lam,newx=xb)
       fit2 <- predict(out,s=lam,newx=xnotinb)
yy <- c(yb,ynotinb)
yhat <- c(fit,fit2)
res <- yy - yhat
vars <- 1:4
lcoef <- predict(out,type="coefficients",s=lam)
                lcoef<-as.vector(lcoef)[-1]
                vin <- vars[lcoef!=0]
                pp <- length(vin)+1
AERplot2(yhat=yhat,y=yy,res=res,d=pp)
title("lasso using covmb set B")

#Problem 7.19 c)

ddplot5(buxx)




##Problem 7.20  elastic net with outliers

#a)


library(glmnet)
out <- enet2(buxx,buxy) #elastic net
                   lam <- out$lam
                   tem <- out$out
                   fit <- predict(tem,s=lam,newx=buxx)
res <- buxy - fit
vars <- 1:4
lcoef <- predict(tem,type="coefficients",s=lam)
                lcoef<-as.vector(lcoef)[-1]
                vin <- vars[lcoef!=0]
                pp <- length(vin)+1
AERplot2(yhat=fit,y=buxy,res=res,d=pp)
title("a) elastic net")

#7.20 b)

#resistant elastic net

ind <-getB(buxx)$indx
yb <- buxy[ind]
xb <- buxx[ind,]
xnotinb <- buxx[-ind,]
ynotinb <- buxy[-ind]
out <- enet2(xb,yb) #elastic net
                   lam <- out$lam
                   tem <- out$out
                   fit <- predict(tem,s=lam,newx=xb)
       fit2 <- predict(tem,s=lam,newx=xnotinb)
yy <- c(yb,ynotinb)
yhat <- c(fit,fit2)
res <- yy - yhat
vars <- 1:4
lcoef <- predict(tem,type="coefficients",s=lam)
                lcoef<-as.vector(lcoef)[-1]
                vin <- vars[lcoef!=0]
                pp <- length(vin)+1
AERplot2(yhat=yhat,y=yy,res=res,d=pp)
title("b) elastic net using covmb set B")


##Problem 7.21 a)

out <- cv.glmnet(cbrainx,cbrainy,alpha=1) #lasso
                   lam <- out$lambda.min
                   fit <- predict(out,s=lam,newx=cbrainx)
res <- cbrainy - fit
vars <- 1:11
lcoef <- predict(out,type="coefficients",s=lam)
                lcoef<-as.vector(lcoef)[-1]
                vin <- vars[lcoef!=0]
                pp <- length(vin)+1
AERplot2(yhat=fit,y=cbrainy,res=res,d=pp)
title("lasso")

#Problem 7.21 b)

xcontin <- cbrainx[,-c(2,4,10)] #continuous predictors
ind <-getB(xcontin)$indx
yb <- cbrainy[ind]
xb <- cbrainx[ind,]
xnotinb <- cbrainx[-ind,]
ynotinb <- cbrainy[-ind]
out <- cv.glmnet(xb,yb,alpha=1) #lasso
                   lam <- out$lambda.min
                   fit <- predict(out,s=lam,newx=xb)
       fit2 <- predict(out,s=lam,newx=xnotinb)
yy <- c(yb,ynotinb)
yhat <- c(fit,fit2)
res <- yy - yhat
vars <- 1:4
lcoef <- predict(out,type="coefficients",s=lam)
                lcoef<-as.vector(lcoef)[-1]
                vin <- vars[lcoef!=0]
                pp <- length(vin)+1
AERplot2(yhat=yhat,y=yy,res=res,d=pp)
title("lasso using covmb set B")

#length(cbrainy) - length(ind) 
#8 cases with the largest distances are omitted from covmb2 set B

#Problem 7.21 c)

ddplot5(xcontin)


##Problem 7.22   needs rpack and library(leaps)

library(leaps)  #may take 20 minutes

regbootsim2(n=100,p=5,k=1,nruns=1000,type=1,psi=0)
vsbootsim3(n=100,p=5,k=1,nruns=1000,type=1,psi=0)


#########


##Problem 8.6

#8.6 a)

n<-100000; q<-7
b <- 0 * 1:q + 1
x <- matrix(rnorm(n * q), nrow = n, ncol = q)
y <- 1 + x %*% b + rnorm(n)
out<-lsfit(x,y)
res<-out$res
yhat<-y-res
dd<-length(out$coef)
AERplot(yhat,y,res=res,d=dd,alph=1) #usual response plot

# 8.6 b) 

AERplot(yhat,y,res=res,d=dd,alph=0.01)
#plots data outside the 99% pointwise PIs


##Problem 8.7  

#8.7 a) 

belx <- 50:73
 
bely <-c(0.44,0.47,0.47,0.59,0.66,0.73,0.81,0.88,1.06,1.2,1.35,1.49,
1.61,2.12,11.9,12.4,14.2,15.9,18.2,21.2,4.3,2.4,2.7,2.9)

out <- lm(bely~belx)
ESP <- predict(out)
Y <- bely
plot(ESP,Y)
abline(0,1)

#8.7 b) 

library(mgcv)
outgam <- gam(bely~s(belx))
EAP <- predict(outgam)
plot(EAP,Y)
abline(0,1)

##############



#Problem 9.6 R commands, needs rpack

n3x <- matrix(rnorm(300),nrow=100,ncol=3)
ln3x <- exp(n3x)
library(MASS)

a)

pairs(n3x)
pairs(ln3x)

b)

ncy <- (n3x%*%1:3)^3 + 0.1*rnorm(100)
plot(n3x%*%(1:3),ncy)

c)

trviews(n3x, ncy)

d)

#change betahat coefficients below

plot(n3x%*%c(22.469,61.242,75.285),n3x%*%1:3)

e)

lncy <- (ln3x%*%1:3)^3 + 0.1*rnorm(100)
trviews(ln3x,lncy)

f)

#change betahat coefficients below
plot(ln3x%*%c(94.848,216.719,328.444),ln3x%*%1:3)

#Problem 9.7 R commands, needs rpack

a)

library(MASS)
nx <- matrix(rnorm(300),nrow=100,ncol=3)
lnx <- exp(nx)
SP <- lnx%*%1:3
lnsincy <- sin(SP)/SP + 0.01*rnorm(100)

b)

trviews(lnx,lnsincy)
essp(lnx,lnsincy,M=40)

c)

ctrviews(lnx,lnsincy)

d)

lmsviews(lnx,lnsincy)

#Problem 9.8 R commands, needs rpack

a)

nx <- matrix(rnorm(300),nrow=100,ncol=3)
SP <- nx%*%1:3
ncuby <- SP^3 + rnorm(100)
nexpy <- exp(SP) + rnorm(100)
nlinsy <- SP + 4*sin(SP) + 0.1*rnorm(100)
nsincy <- sin(SP)/SP + 0.01*rnorm(100)
nsiny <- sin(SP) + 0.1*rnorm(100)
nsqrty <- sqrt(abs(SP)) + 0.1*rnorm(100)
nsqy <- SP^2 + rnorm(100)

b)

plot(SP,ncuby)
plot(-SP,ncuby)

c)

library(MASS)
trviews(nx,ncuby)
essp(nx,ncuby, M=40)

d)

tem <- c(12.60514, 25.06613, 37.25504)
ESP <- nx%*%tem
plot(ESP,SP)

e)

lnx <- exp(nx)
SP <- lnx%*%1:3
lncuby <- (SP/3)^3 + rnorm(100)
lnlinsy <- SP + 10*sin(SP) + 0.1*rnorm(100)
lnsincy <- sin(SP)/SP + 0.01*rnorm(100)
lnsiny <- sin(SP/3) + 0.1*rnorm(100)
ESP <- lnx%*%tem



################

#Problem 11.19 R commands, needs rpack

out <- lrdata()
x <- out$x
y <- out$y
lressp(x,y)

#Problem 11.20 R commands, needs rpack

a)

out <- prdata()
x <- out$x
y <- out$y
pressp(x,y)


b)

prplot(x,y)


#Problem 11.21

a) 

belx <- 50:73
 
bely <-c(0.44,0.47,0.47,0.59,0.66,0.73,0.81,0.88,1.06,1.2,1.35,1.49,
1.61,2.12,11.9,12.4,14.2,15.9,18.2,21.2,4.3,2.4,2.7,2.9)

out <- lm(bely~belx)
ESP <- predict(out)
Y <- bely
plot(ESP,Y)
abline(0,1)

b) 

library(mgcv)
outgam <- gam(bely~s(belx))
EAP <- predict(outgam)
plot(EAP,Y)
abline(0,1)

#Problem 11.22

#ha = number patients with heart attack for each ck group
#ok = number of patients without a heart attack

a)

ck <- c(20,60,100,140,180,220,260,300,340,380,420,460)
ha <- c(2,13,30,30,21,19,18,13,19,15,7,8)
ok <- c(88,26,8,5,0,1,1,1,1,0,0,0)
heart <- as.data.frame(cbind(ck,ha,ok))
outw <- glm(cbind(ha,ok)~ck+I(ck^2)+I(ck^3),family=binomial, data=heart)
out <- glm(cbind(ha,ok)~log(ck),family=binomial, data=heart)
par(mfrow=c(2,2))

library(mgcv)
outgam <- gam(cbind(ha,ok)~s(ck),family=binomial, data=heart)
plot(outgam)
title("a) Shat")

outgam2 <- gam(cbind(ha,ok)~s(log(ck)),family=binomial, data=heart)
plot(outgam2)
title("b) Shat for log(ck) GAM")

EAP <- predict(outgam)
ESP <- predict(out)
plot(EAP,ESP)
abline(0,1)
title("c) EE Plot: log(ck) GLM")

Z <- ha/(ha+ok)
fit<-Z
plot(ESP,Z)
fit <- exp(ESP)/(1 + exp(ESP))
indx <- sort.list(ESP)
lines(ESP[indx], fit[indx])
title("d) Response Plot ")
par(mfrow=c(1,1))

b)

AIC(outw)

AIC(out) 



#HW 11.23;
#make sure cbrain.dat.txt is on j drive or the infile command won't work;
options ls = 70;
data cbrain;
infile "j:cbrain.dat.txt";
input obs x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 sex x12;

proc rsquare;
  model sex = x2 x3 x4 x5 x6 x7 x9 x10 / cp;
  
proc reg;
  model sex = x2 x3 x4 x5 x6 x7 x9 x10;
  output out = a p = pred r = resid;
  run;