Model Selection

Model selection using AIC and WAIC.

library(rethinking)
data("Howell1")
myData<-Howell1
head(myData)

plot(myData$height,myData$weight)

Using R’s built in lm() to fit a linear and quadratic model.

lmFit<-lm(weight~height,data=myData)
lmFitQuad<-lm(weight~height + I(height^2),data=myData)

lmFit

Call:
lm(formula = weight ~ height, data = myData)

Coefficients:
(Intercept)       height  
   -33.7561       0.5017  

lmFitQuad

Call:
lm(formula = weight ~ height + I(height^2), data = myData)

Coefficients:
(Intercept)       height  I(height^2)  
  20.844760    -0.472953     0.004033  

broom::glance(lmFitQuad)$AIC

[1] 3041.213

Plotting the fit

myData2<-cbind(myData,predict(lmFit),predict(lmFitQuad))
head(myData2)

myData2ordered<-myData2[order(myData2$height),]
plot(myData2$height,myData2$weight)
lines(myData2ordered$height,myData2ordered[['predict(lmFit)']])

lines(myData2ordered$height,myData2ordered[['predict(lmFitQuad)']])

Extracting AIC from the lm() fits. Lower values of AIC are favored.

broom::glance(lmFit)$AIC - broom::glance(lmFitQuad)$AIC

[1] 256.1191

library(stats)
AIC(lmFit)

[1] 3297.332

AIC(lmFitQuad)

[1] 3041.213

AIC(lmFit)-AIC(lmFitQuad)

[1] 256.1191

extractAIC(lmFit)

[1]    2.000 1751.527

extractAIC(lmFitQuad)

[1]    3.000 1495.408

extractAIC(lmFit)-extractAIC(lmFitQuad)

[1]  -1.0000 256.1191

Thus, we would select the curved line over the straight line.

Next, writing our own likelihood functions to fit a linear and quadratic models.

nllRegressionLinear<-function(parmVector,weightDat,heightDat,print=F){

  intercept<-parmVector[1]
  B<-parmVector[2]
  mySD<-parmVector[3]
  
  w<-weightDat
  h<-heightDat
  mu<-intercept + B*h
  
  nllik<- -sum(dnorm(x=w,mean=mu,sd=mySD,log=TRUE))

  if (print) cat("nllik= ",nllik,sep=" ",fill=T)
  return(nllik)}

parVecLinear<-c(0.1,0.1,10) # Initial parameter values 
outRegLinear<-optim(par=parVecLinear,fn=nllRegressionLinear,method="L-BFGS-B",lower=c(-Inf,-Inf, 0.01),upper=c(Inf,Inf, Inf), weightDat=myData$weight, heightDat=myData$height,print=F)
outRegLinear$par 

[1] -33.7561543   0.5016994   4.9837499

outRegLinear$value

[1] 1645.666

Calculating AIC where AIC = 2k - 2ln(L) where L is the maximized likelihood. L is given by outRegLinear\(value from above. k is the number of parameters fit, which is equal to the length of parVecLinear. For our linear model, its an intercept, slope, and sd. Note since we calculated the negative log likelihood, we calculate AIC using 2*k + 2*outRegLinear\)value

lmAIC = 2*length(parVecLinear) + 2*outRegLinear$value
lmAIC

[1] 3297.332

broom::glance(lmFit)

nllRegressionQuad<-function(parmVector,weightDat,heightDat,print=F){

  intercept<-parmVector[1]
  B1<-parmVector[2]
  B2<-parmVector[3]
  mySD<-parmVector[4]
  
  w<-weightDat
  h<-heightDat
  mu<-intercept + B1*h + B2*h^2
  
  nllik<- -sum(dnorm(x=w,mean=mu,sd=mySD,log=TRUE))

  if (print) cat("nllik= ",nllik,sep=" ",fill=T)
  return(nllik)}

parVecQuad<-c(20.1,0.1,0.1,10) 
outRegQuad<-optim(par=parVecQuad,fn=nllRegressionQuad,method="L-BFGS-B",lower=c(-Inf,-Inf,-Inf, 0.01),upper=c(Inf,Inf,Inf,Inf), weightDat=myData$weight, heightDat=myData$height,print=F)
outRegQuad$par 

[1] 20.092051849 -0.460139727  0.003981269  3.931391664

outRegQuad$value

[1] 1516.635

broom::glance(lmFitQuad)

Calculating AIC for the quadratic model. Note that we have one additional parameter compared to the linear model.

quadAIC = 2*length(parVecQuad) + 2*outRegQuad$value
quadAIC

[1] 3041.271

DeltaAIC<-lmAIC-quadAIC
DeltaAIC

[1] 256.0611

Plotting the fit

plot(myData$height,myData$weight)
abline(a=outRegLinear$par[1],b=outRegLinear$par[2])

curve(outRegQuad$par[1] + outRegQuad$par[2]*x + outRegQuad$par[3]*x^2 , 55, 180,add=TRUE)

Finally, using the QUAP function in the Rethinking package.

library(rethinking)
lm.qa <- quap(
    alist(
        weight ~ dnorm(mu,sd) ,  # normal likelihood
        sd ~ dexp(1),     # exponential prior
        mu<-intercept + beta*height,
        intercept~dnorm(0,10),
        beta~dnorm(0,10)
    ) ,
    data=list(weight=myData$weight,height=myData$height) )

precis(lm.qa)

WAIC(lm.qa)

library(rethinking)
quad.qa <- quap(
    alist(
        weight ~ dnorm(mu,sd) ,  # normal likelihood
        sd ~ dexp(1),     # exponential prior
        mu<-intercept + beta1*height + beta2*height2,
        intercept~dnorm(0,10),
        beta1~dnorm(0,10),
        beta2~dnorm(0,10)
    ) ,
    data=list(weight=myData$weight,height=myData$height, height2=myData$height^2),
    start=list(sd=5))

precis(quad.qa,digits=4)

WAIC(quad.qa)

WAIC(lm.qa)$WAIC - WAIC(quad.qa)$WAIC

[1] 256.1963

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

plot(cars)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.