## Plotting parameter values
cexval <- 1.9
cex.axisval <- cexval
cex.labval <- cexval
pchval <- 19
youngCol <- "blue"
oldCol <- "red"
plotWidth <- 720
plotHeight <- 720
plotBG <- "transparent"
marvals <- c(5,5,2,2)
lwdval <- 3

### set up a cubic function for f(X)
beta3 <- -0.0004  ## cubic coefficient
beta2 <- 0.025 ## quadratic coefficient
beta1 <- 2.0 ## linear coefficient
beta0 <- 50.0 ## constant intercept term
fx <- function(x) beta3*(x^3) + beta2*(x^2) + beta1*x + beta0  ## True function f(X)

Nsim <- 1000 ## number of data points to simulate
smallN <- 50 ## size for smaller dataset 
noiseSD <- 10.0 ## SD of additive noise
minX <- 40  ## minimum value for predictor
maxX <- 80  ## maximum value for predictor
Xvec <- minX:maxX

simfn1 <- function(nsim,f){ ## function to generate nsim points
  ## from another function f
  x <- sample(x=Xvec, size=nsim, replace=TRUE) ## Simulate N points uniformly along minX to maxX at integers
  truefx <- f(x)  ## The true values of f(X) for the simulated x
  epsilon <- rnorm(n=nsim, mean=0, sd=noiseSD) ## simulate zero mean noise for each data point of sd = noise SD
  y <- truefx + epsilon ## Generate observed data y as f(X) + epsilon
  datset <- data.frame(x,y,truefx,epsilon)
  names(datset) <- c("age","weight","trueFx","noise")
  datset
}

set.seed(1)
dat <- simfn1(nsim=Nsim, f=fx)

## Plotting first figure
png(filename="Lecture1plots/plot1.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$age,dat$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval/5, cex.axis=cex.axisval, cex.lab=cex.labval)
dev.off()

## Fitting the model with the true form
lm1 <- lm(weight ~ age + I(age^2) + I(age^3),data=dat)

coeflm <- as.numeric(coef(lm1))
fxx <- function(x) coeflm[4]*(x^3) + coeflm[3]*(x^2) + coeflm[2]*x + coeflm[1]

## Plotting second figure
png(filename="Lecture1plots/plot2.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$age,dat$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval/5, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(fxx,add=TRUE,lwd=lwdval,col="red") ## the polynomial 
                       ## regression estimated cubic curve
## curve(fx,add=TRUE,lwd=lwdval,col="blue") ## the true 
                                          ## cubic curve
dev.off()

pred60 <- fxx(60) ## fitted polynomial prediction at age 60

## Plotting second figure for a prediction
png(filename="Lecture1plots/plot2_5.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$age,dat$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval/5, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(fxx,add=TRUE,lwd=lwdval,col="red") ## the polynomial 
## regression estimated cubic curve
segments(x0=60,y0=115,x1=60,y1=pred60,col="red",lty=2,lwd=1.5)
segments(x0=35,y0=pred60,x1=60,y1=pred60,col="red",lty=2,lwd=1.5)
## curve(fx,add=TRUE,lwd=lwdval,col="blue") ## the true 
## cubic curve
dev.off()


#### calculating mean weight at age 60
dat60 <- subset(dat,subset=age==60)
mean60 <- mean(dat60$weight)

## Plotting third figure
png(filename="Lecture1plots/plot3.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$age,dat$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval/5, cex.axis=cex.axisval, cex.lab=cex.labval)
points(60,mean60,pch=pchval,cex=cexval,col="red")
points(60,mean60,pch=4,cex=5,col="red")
segments(x0=60,y0=115,x1=60,y1=mean60,col="red",lty=2,lwd=1.5)
segments(x0=35,y0=mean60,x1=60,y1=mean60,col="red",lty=2,lwd=1.5)
dev.off()


### Calculating mean weight at each age
meanWeight <- vector(length = length(Xvec))
for(i in Xvec){
  datsub <- dat[dat$age==i, ]
  meanWeight[i-minX+1] <- mean(datsub$weight)
}

## Plotting fourth figure
png(filename="Lecture1plots/plot4.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(dat$age,dat$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval/5, cex.axis=cex.axisval, cex.lab=cex.labval)
points(Xvec, meanWeight, cex=cexval/3, col="red")
lines(Xvec, meanWeight, lwd=lwdval, col="red")
dev.off()
  
### Use a smaller dataset
datSmall <- dat[1:smallN, ]

## Plotting fifth figure
png(filename="Lecture1plots/plot5.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
dev.off()

### Using kernel smoothing
ksmooth1 <- ksmooth(datSmall$age, datSmall$weight, kernel="box", bandwidth = 10)

## Plotting sixth figure
png(filename="Lecture1plots/plot6.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
lines(ksmooth1$x,ksmooth1$y,lwd=lwdval,col="red")
abline(v=55,lty=2,lwd=lwdval, col="red")
abline(v=60,lwd=lwdval, col="red")
abline(v=65,lty=2,lwd=lwdval, col="red")
dev.off()

### Using kernel smoothing
ksmooth2 <- ksmooth(datSmall$age, datSmall$weight, kernel="box", bandwidth = 20)

## Plotting seventh figure
png(filename="Lecture1plots/plot7.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
lines(ksmooth2$x,ksmooth2$y,lwd=lwdval,col="red")
abline(v=50,lty=2,lwd=lwdval, col="red")
abline(v=60,lwd=lwdval, col="red")
abline(v=70,lty=2,lwd=lwdval, col="red")
dev.off()


## Sort the data on age because ksmooth sorts the data on age
## when it is applied
datSmallSort <- datSmall[order(datSmall$age), ]

kernels <- c(1:50)
mse <- vector(length=length(kernels))
for(i in 1:length(kernels)){
  ksmoothTemp <- ksmooth(datSmallSort$age, datSmallSort$weight, x.points = datSmallSort$age, kernel="box", bandwidth = kernels[i])
  mse[i] <- mean((datSmallSort$weight - ksmoothTemp$y)^2)
}

## Plotting eighth figure
png(filename="Lecture1plots/plot8.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(kernels,mse, pch=".", xlab="Bandwidth", ylab="MSE", cex=0.5, cex.axis=cex.axisval, cex.lab=cex.labval)
lines(kernels,mse,lwd=lwdval,col="red")
dev.off()


set.seed(2)
datSmall2 <- simfn1(nsim=smallN, f=fx)
datSmallFULL <- rbind(datSmall,datSmall2)

## Plotting ninth figure
png(filename="Lecture1plots/plot9.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmallFULL$age,datSmallFULL$weight, pch=".", xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
points(datSmall2$age,datSmall2$weight, pch=19, cex=cexval, col="blue")
dev.off()

## Plotting tenth figure
png(filename="Lecture1plots/plot10.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmallFULL$age,datSmallFULL$weight, pch=".", xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
points(datSmall2$age,datSmall2$weight, pch=19, cex=cexval, col="blue")
points(datSmall$age,datSmall$weight, pch=19, cex=cexval, col="red")
lines(ksmooth1$x,ksmooth1$y,lwd=lwdval,col="red")
legend(x="bottomright",legend=c("Training","Test"),pch=rep(pchval,2),col=c("red","blue"),cex=cexval)
dev.off()


### Generating mean square error vs. test data fit

datSmallSort2 <- datSmall2[order(datSmall2$age), ]

kernels <- c(1:50)
mse2 <- vector(length=length(kernels))
for(i in 1:length(kernels)){
  ksmoothTemp <- ksmooth(datSmallSort$age, datSmallSort$weight, x.points = datSmallSort2$age, kernel="box", bandwidth = kernels[i])
  mse2[i] <- mean((datSmallSort2$weight - ksmoothTemp$y)^2)
}

minIndex <- which.min(mse2)
kernels[minIndex]
mse2[minIndex]

## Plotting eleventh figure
png(filename="Lecture1plots/plot11.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(rep(kernels,2),c(mse,mse2), pch=".", xlab="Bandwidth", ylab="MSE", cex=0.5, cex.axis=cex.axisval, cex.lab=cex.labval)
lines(kernels,mse,lwd=lwdval,col="red")
lines(kernels,mse2,lwd=lwdval,col="blue")
legend(x="bottomright",legend=c("Training","Test"),lwd=rep(lwdval,2),col=c("red","blue"),cex=cexval)
dev.off()


## Plotting twelfth figure
png(filename="Lecture1plots/plot12.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(rep(kernels,2),c(mse,mse2), pch=".", xlab="Bandwidth", ylab="MSE", cex=0.5, cex.axis=cex.axisval, cex.lab=cex.labval)
lines(kernels,mse,lwd=lwdval,col="red")
lines(kernels,mse2,lwd=lwdval,col="blue")
segments(x0=kernels[minIndex],y0=35,x1=kernels[minIndex],y1=mse2[minIndex],col="blue",lty=2,lwd=lwdval)
legend(x="bottomright",legend=c("Training","Test"),lwd=rep(lwdval,2),col=c("red","blue"),cex=cexval)
dev.off()

#### Simulating some normalized 2 predictor data

set.seed(3)
bp <- rnorm(100)
bil <- rnorm(100)

sum(abs(bp)<0.18) # margin to have 10 center points 1D
sum(sqrt(bp^2+bil^2)<0.553) # margin to have 10 center points 2D

## 13th figure
png(filename="Lecture1plots/plot13.png", width=150, height=720, bg="transparent")
stripchart(bp,vertical=TRUE,pch=19,cex=0.3,col="blue",cex.axis=cex.axisval, cex.lab=cex.labval)
points(-0.18,pch=21,col="black",bg="red",cex=1.3)
points(0.18,pch=21,col="black",bg="red",cex=1.3)
dev.off()

library(plotrix)  ## to get draw.circle function
## 14th figure
png(filename="Lecture1plots/plot14.png", width=720, height=720, bg="transparent")
par(mar=marvals,pty="s") ## setting margin width
plot(bil,bp, pch=19, xlab="billirubin", ylab="blood pressure", cex=0.5, col="blue",cex.axis=cex.axisval, cex.lab=cex.labval)
## draw.circle(0,0,radius=0.553,border="red",col="transparent",lwd=lwdval)
draw.circle(0,0,radius=0.6,border="red",col="transparent",lwd=lwdval)
dev.off()


################# Linear modeling

lmfit <- list()
coeflmfit <- list()
fxxfit <- list()

lmfit[[1]] <- lm(weight ~ age,data=datSmall)
lmfit[[2]] <- lm(weight ~ age + I(age^2),data=datSmall)
lmfit[[3]] <- lm(weight ~ age + I(age^2) + I(age^3),data=datSmall)
lmfit[[4]] <- lm(weight ~ age + I(age^2) + I(age^3) + I(age^4),data=datSmall)
lmfit[[5]] <- lm(weight ~ age + I(age^2) + I(age^3) + I(age^4) + I(age^5),data=datSmall)

for(i in 1:5) coeflmfit[[i]] <- as.numeric(coef(lmfit[[i]]))

fxxfit1 <- function(x) coeflmfit[[1]][2]*x + coeflmfit[[1]][1]
fxxfit2 <- function(x) coeflmfit[[2]][3]*x^2 + coeflmfit[[2]][2]*x + coeflmfit[[2]][1]
fxxfit3 <- function(x) coeflmfit[[3]][4]*x^3 + coeflmfit[[3]][3]*x^2 + coeflmfit[[3]][2]*x + coeflmfit[[3]][1]
fxxfit4 <- function(x) coeflmfit[[4]][5]*x^4 + coeflmfit[[4]][4]*x^3 + coeflmfit[[4]][3]*x^2 + coeflmfit[[4]][2]*x + coeflmfit[[4]][1]
fxxfit5 <- function(x) coeflmfit[[5]][6]*x^5 + coeflmfit[[5]][5]*x^4 + coeflmfit[[5]][4]*x^3 + coeflmfit[[5]][3]*x^2 + coeflmfit[[5]][2]*x + coeflmfit[[5]][1]


## 15th figure
png(filename="Lecture1plots/plot15.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(fxxfit1,add=TRUE,lwd=lwdval,col="red") ## linear fit
legend("bottomright",legend="linear fit",col="red",lwd=lwdval)
dev.off()


## 16th figure
png(filename="Lecture1plots/plot16.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(fxxfit1,add=TRUE,lwd=lwdval,col="red") ## linear fit
curve(fxxfit2,add=TRUE,lwd=lwdval,col="dark green") ## quadratic fit 
legend("bottomright",legend=c("linear fit","quadratic fit"),col=c("red","dark green"),lwd=lwdval)
dev.off()

## 17th figure
png(filename="Lecture1plots/plot17.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(fxxfit1,add=TRUE,lwd=lwdval,col="red") ## linear fit
curve(fxxfit2,add=TRUE,lwd=lwdval,col="dark green") ## quadratic fit
curve(fxxfit3,add=TRUE,lwd=lwdval,col="blue") ## cubic fit 
legend("bottomright",legend=c("linear fit","quadratic fit","cubic fit"),col=c("red","dark green","blue"),lwd=lwdval)
dev.off()


## 18th figure
png(filename="Lecture1plots/plot18.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(fxxfit1,add=TRUE,lwd=lwdval,col="red") ## linear fit
curve(fxxfit2,add=TRUE,lwd=lwdval,col="dark green") ## quadratic fit
curve(fxxfit3,add=TRUE,lwd=lwdval,col="blue") ## cubic fit 
curve(fxxfit4,add=TRUE,lwd=lwdval,col="dark orange") ## quartic fit 
legend("bottomright",legend=c("linear fit","quadratic fit","cubic fit","quartic fit"),col=c("red","dark green","blue","dark orange"),lwd=lwdval)
## curve(fx,add=TRUE,lwd=lwdval,col="blue") ## the true curve
dev.off()


## 19th figure
png(filename="Lecture1plots/plot19.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(fxxfit1,add=TRUE,lwd=lwdval,col="red") ## linear fit
curve(fxxfit2,add=TRUE,lwd=lwdval,col="dark green") ## quadratic fit
curve(fxxfit3,add=TRUE,lwd=lwdval,col="blue") ## cubic fit 
curve(fxxfit4,add=TRUE,lwd=lwdval,col="dark orange") ## quartic fit 
curve(fxxfit5,add=TRUE,lwd=lwdval,col="purple") ## quintic fit 
legend("bottomright",legend=c("linear fit","quadratic fit","cubic fit","quartic fit","quintic fit"),col=c("red","dark green","blue","dark orange","purple"),lwd=lwdval)
## curve(fx,add=TRUE,lwd=lwdval,col="blue") ## the true curve
dev.off()

### Generating mean square error vs. training andtest data fit

orders <- c(1:5)
mseLin <- mseLin2 <- vector(length=length(orders))

for(i in 1:5) mseLin[i] <- mean((datSmall$weight - predict(lmfit[[i]],newdata=datSmall))^2)
for(i in 1:5) mseLin2[i] <- mean((datSmall2$weight - predict(lmfit[[i]],newdata=datSmall2))^2)


minLinIndex <- which.min(mseLin2)
orders[minLinIndex]
mseLin2[minIndex]

## 20th figure
png(filename="Lecture1plots/plot20.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(rep(orders,2),c(mseLin,mseLin2), pch=".", xlab="Polynomial order", ylab="MSE", cex=0.5, cex.axis=cex.axisval, cex.lab=cex.labval)
lines(orders,mseLin,lwd=lwdval,col="red")
lines(orders,mseLin2,lwd=lwdval,col="blue")
legend(x="topright",legend=c("Training","Test"),lwd=rep(lwdval,2),col=c("red","blue"),cex=cexval)
dev.off()


## 21st figure
png(filename="Lecture1plots/plot21.png", width=870, height=650, bg="transparent")
par(mar=marvals) ## setting margin width
plot(datSmall$age,datSmall$weight, pch=pchval, xlab="X, age (years)", ylab="f(X), weight (lb)", cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
curve(fxxfit1,add=TRUE,lwd=lwdval,col="red") ## linear fit
curve(fxxfit2,add=TRUE,lwd=lwdval,col="dark green") ## quadratic fit
curve(fxxfit3,add=TRUE,lwd=lwdval,col="blue") ## cubic fit 
curve(fxxfit4,add=TRUE,lwd=lwdval,col="dark orange") ## quartic fit 
curve(fxxfit5,add=TRUE,lwd=lwdval,col="purple") ## quintic fit 
curve(fx,add=TRUE,lwd=5,lty=2,col="dark blue") ## the true curve
legend("bottomright",legend=c("linear fit","quadratic fit","cubic fit","quartic fit","quintic fit","true model"),col=c("red","dark green","blue","dark orange","purple","dark blue"),lwd=c(rep(lwdval,5),5),lty=c(rep(1,4),2))
dev.off()


#### finding turning points of quartic
root1 <- (-2*beta2 + sqrt(4*beta2^2 - 12*beta3*beta1))/(6*beta3)
root2 <- (-2*beta2 - sqrt(4*beta2^2 - 12*beta3*beta1))/(6*beta3)
root1
root2
fx(root1)
fx(root2)

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

#### 2D polynomial function + noise
b0 <- 40.0
b1 <- 1.5
b2 <- 0.7
b3 <- -0.4
b4 <- 0.2
b5 <- 0.08

fxx2D <- function(x1,x2) b0 + b1*x1 + b2*x2 + b3*x1^2 + b4*x2^2 + b5*x1*x2

fxx2Dmatrix <- function(x) b0 + b1*x[1] + b2*x[2] + b3*x[1]^2 + b4*x[2]^2 + b5*x[1]*x[2]

## Generate the true surface plot
rangevals <- 2
x1lim <- x2lim <- seq(-rangevals,rangevals,length.out=100)
surfvals <- outer(x1lim,x2lim,fxx2D)
nrz <- nrow(surfvals)
ncz <- ncol(surfvals)

## Generate 50 random locations in the surface.
set.seed(6)
npoints <- 50
x1vals <- runif(n=npoints, min=(-rangevals), max=rangevals)
x2vals <- runif(n=npoints, min=(-rangevals), max=rangevals)
noisevals <- rnorm(n=npoints,mean=0,sd=0.3)
respvals <- apply(cbind(x1vals,x2vals),MARGIN=1,FUN=fxx2Dmatrix)
datvals <- respvals + noisevals
polydat <- data.frame(x1vals,x2vals,respvals,noisevals,datvals)

# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("blue", "green") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the surfvals at the facet centres
zfacet <- surfvals[-1, -1] + surfvals[-1, -ncz] + surfvals[-nrz, -1] + surfvals[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

png(filename="Lecture1plots/plot22.png", width=640, height=640, bg="transparent")
p  <- persp(x1lim, x2lim, surfvals, theta = -50, phi = 20, col=color[facetcol], shade=0.1, border=NA, ticktype = "detailed", xlab="executive function score", ylab="memory function score", zlab="frontal lobe volume",cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
obs <- trans3d(polydat$x1vals,polydat$x2vals,polydat$datvals, p)
pred <- trans3d(polydat$x1vals,polydat$x2vals,polydat$respvals, p)
points(obs, col="red",pch=19)
segments(obs$x, obs$y, pred$x, pred$y, col="black")
dev.off()

### linear fit

lm2Da <- lm(datvals ~ x1vals + x2vals, data = polydat)
coef2Da <- as.numeric(coef(lm2Da))
fxx2Dafit <- function(x1,x2) coef2Da[1] + coef2Da[2]*x1 + coef2Da[3]*x2
surfvals2Da <- outer(x1lim,x2lim,fxx2Dafit)
polydat$linvals <- fxx2Dafit(polydat$x1vals,polydat$x2vals)

png(filename="Lecture1plots/plot23.png", width=640, height=640, bg="transparent")
p  <- persp(x1lim, x2lim, surfvals2Da, theta = -50, phi = 20, col=color[facetcol], shade=0.1, border=NA, ticktype = "detailed", xlab="executive function score", ylab="memory function score", zlab="frontal lobe volume",cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
obs <- trans3d(polydat$x1vals,polydat$x2vals,polydat$datvals, p)
pred <- trans3d(polydat$x1vals,polydat$x2vals,polydat$linvals, p)
points(obs, col="red",pch=19)
segments(obs$x, obs$y, pred$x, pred$y, col="black")
dev.off()

### linear fit with interaction

lm2Db <- lm(datvals ~ x1vals * x2vals, data = polydat)
coef2Db <- as.numeric(coef(lm2Db))
fxx2Dbfit <- function(x1,x2) coef2Db[1] + coef2Db[2]*x1 + coef2Db[3]*x2 + coef2Db[4]*x1*x2
surfvals2Db <- outer(x1lim,x2lim,fxx2Dbfit)
polydat$intvals <- fxx2Dbfit(polydat$x1vals,polydat$x2vals)

png(filename="Lecture1plots/plot24.png", width=640, height=640, bg="transparent")
p  <- persp(x1lim, x2lim, surfvals2Db, theta = -50, phi = 20, col=color[facetcol], shade=0.1, border=NA, ticktype = "detailed", xlab="executive function score", ylab="memory function score", zlab="frontal lobe volume",cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
obs <- trans3d(polydat$x1vals,polydat$x2vals,polydat$datvals, p)
pred <- trans3d(polydat$x1vals,polydat$x2vals,polydat$intvals, p)
points(obs, col="red",pch=19)
segments(obs$x, obs$y, pred$x, pred$y, col="black")
dev.off()

### quadratic fit 

lm2Dc <- lm(datvals ~ x1vals * x2vals + I(x1vals^2) + I(x2vals^2), data = polydat)
coef2Dc <- as.numeric(coef(lm2Dc))
fxx2Dcfit <- function(x1,x2) coef2Dc[1] + coef2Dc[2]*x1 + coef2Dc[3]*x2 + coef2Dc[4]*x1^2 + coef2Dc[5]*x2^2 + coef2Dc[6]*x1*x2
surfvals2Dc <- outer(x1lim,x2lim,fxx2Dcfit)
polydat$quadvals <- fxx2Dcfit(polydat$x1vals,polydat$x2vals)

png(filename="Lecture1plots/plot25.png", width=640, height=640, bg="transparent")
p  <- persp(x1lim, x2lim, surfvals2Dc, theta = -50, phi = 20, col=color[facetcol], shade=0.1, border=NA, ticktype = "detailed", xlab="executive function score", ylab="memory function score", zlab="frontal lobe volume",cex=cexval, cex.axis=cex.axisval, cex.lab=cex.labval)
obs <- trans3d(polydat$x1vals,polydat$x2vals,polydat$datvals, p)
pred <- trans3d(polydat$x1vals,polydat$x2vals,polydat$quadvals, p)
points(obs, col="red",pch=19)
segments(obs$x, obs$y, pred$x, pred$y, col="black")
dev.off()