🐷🦷 cursed data from Hodges chapter 4

ohno.R

# Hodges pig jawbone
# http://www.biostat.umn.edu/~hodges/RPLMBook/Datasets/06_Pig_jawbone/Ex6.html


library(mgcv)
library(ggplot2)

source("smooth.construct.tr.smooth.spec.R")

whitey <- read.csv("Pig_jawbone_Whitey_only.csv")
whitey$fTransect <- factor(whitey$Transect, ordered=TRUE)

# base plot
p <- ggplot(whitey, aes(x=distmm, y=ElMod)) +
  geom_point(colour="grey60") +
  facet_wrap(~fTransect) +
  theme_minimal()
print(p)


# prediction grid
predg <- expand.grid(distmm = seq(0, 1.5, length.out=150),
                     fTransect = factor(1:9, levels=levels(whitey$fTransect)))


## fit some models


# hodges (4.11)
b_411 <- gam(ElMod ~ fTransect + fTransect:distmm + fTransect:I(distmm^2) +
                     s(distmm, by=fTransect, bs="tr", k=25),
             method="REML", data=whitey)

summary(b_411)
# what does that look like?
p_411 <- cbind(predg, ElMod=predict(b_411, predg))
p + geom_line(data=p_411)

# hodges (4.12) - (4.14)
# this adds a "global" smooth of distmm
b_412 <- gam(ElMod ~ fTransect + fTransect:distmm + fTransect:I(distmm^2) +
                     s(distmm, bs="tr", k=25) +
                     s(distmm, by=fTransect, bs="tr", k=25),
             method="REML", data=whitey)
summary(b_412)

# compre these models
p_412 <- cbind(predg, ElMod=predict(b_412, predg))
p + geom_line(data=p_411, colour="red") + geom_line(data=p_412, colour="blue")



# a better version of 4.11, using the fs basis (fewer DF, smoopars)
b_fs11 <- gam(ElMod ~ fTransect + fTransect:distmm + fTransect:I(distmm^2) +
                     s(distmm, fTransect, bs="fs",
                       xt=list(bs="tr", k=25)),
             method="REML", data=whitey)
summary(b_fs11)
# get rid of those fixed effects?
b_fs11_2 <- gam(ElMod ~ fTransect +
                     s(distmm, fTransect, bs="fs",
                       xt=list(bs="tr", k=25)),
             method="REML", data=whitey)
summary(b_fs11_2)

# compare? v. similar but diff. extrapolative behaviour?
p_fs11 <- cbind(predg, ElMod=predict(b_fs11, predg))
p_fs11_2 <- cbind(predg, ElMod=predict(b_fs11_2, predg))
p + geom_line(data=p_fs11, colour="red") + geom_line(data=p_fs11_2, colour="blue")

# not much point in doing fs versions of (4.12)-(4.14) as the above fs models
# already do the global smooth + diffs thing?

# mainly to compare degrees of freedom...
AIC(b_411, b_412, b_fs11, b_fs11_2)


## what about setting-up the RW model?

# setup design matrix, this is awful code
pattern <- matrix(0, 9, 8)
pattern[lower.tri(pattern)] <- 1

Xrw <- apply(pattern, 2, rep, times=as.numeric(table(whitey$Transect)))

# construct penalty
Srw <- diag(8)

whitey$Xrw <- Xrw

b_rw <- gam(ElMod ~ distmm + I(distmm^2) + Xrw +
                    s(distmm, bs="tr", k=25),
            paraPen=list(Xrw=list(Srw)),
            method="REML", data=whitey)
summary(b_rw)


# model with just a random effect and overall smooth
b_re <- gam(ElMod ~ distmm + I(distmm^2) +
                    s(fTransect, bs="re") +
                    s(distmm, bs="tr", k=25),
            method="REML", data=whitey)
summary(b_re)

# compare fixed effects sizes
coef(b_re)[1:3]
coef(b_rw)[1:3]

# compare fits
p_re <- cbind(predg, ElMod=predict(b_re, predg))


predg$Xrw <- apply(pattern, 2, rep, times=as.numeric(table(predg$fTransect)))


p_rw <- cbind(predg, ElMod=predict(b_rw, predg))
p + geom_line(data=p_re, colour="red") + geom_line(data=p_rw, colour="blue")

cor(predict(b_re), predict(b_rw))

Pig_jawbone_Whitey_only.csv

smooth.construct.tr.smooth.spec.R

     ## Adding a penalized truncated power basis class and methods
     ## as favoured by Ruppert, Wand and Carroll (2003) 
     ## Semiparametric regression CUP. (No advantage to actually
     ## using this, since mgcv can happily handle non-identity 
     ## penalties.)
     
     smooth.construct.tr.smooth.spec<-function(object,data,knots) {
     ## a truncated power spline constructor method function
     ## object$p.order = null space dimension
       m <- object$p.order[1]
       if (is.na(m)) m <- 2 ## default 
       if (m<1) stop("silly m supplied")
       if (object$bs.dim<0) object$bs.dim <- 10 ## default
       nk<-object$bs.dim-m-1 ## number of knots
       if (nk<=0) stop("k too small for m")
       x <- data[[object$term]]  ## the data
       x.shift <- mean(x) # shift used to enhance stability
       k <- knots[[object$term]] ## will be NULL if none supplied
       if (is.null(k)) # space knots through data
       { n<-length(x)
         k<-quantile(x[2:(n-1)],seq(0,1,length=nk+2))[2:(nk+1)]
       }
       if (length(k)!=nk) # right number of knots?
       stop(paste("there should be ",nk," supplied knots"))
       x <- x - x.shift # basis stabilizing shift
       k <- k - x.shift # knots treated the same!
       X<-matrix(0,length(x),object$bs.dim)
       for (i in 1:(m+1)) X[,i] <- x^(i-1)
       for (i in 1:nk) X[,i+m+1]<-(x-k[i])^m*as.numeric(x>k[i])
       object$X<-X # the finished model matrix
       if (!object$fixed) # create the penalty matrix
       { object$S[[1]]<-diag(c(rep(0,m+1),rep(1,nk)))
       }
       object$rank<-nk  # penalty rank
       object$null.space.dim <- m+1  # dim. of unpenalized space
       ## store "tr" specific stuff ...
       object$knots<-k;object$m<-m;object$x.shift <- x.shift
      
       object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
      
       class(object)<-"tr.smooth"  # Give object a class
       object
     }
     
     Predict.matrix.tr.smooth<-function(object,data) {
     ## prediction method function for the `tr' smooth class
       x <- data[[object$term]]
       x <- x - object$x.shift # stabilizing shift
       m <- object$m;     # spline order (3=cubic)
       k<-object$knots    # knot locations
       nk<-length(k)      # number of knots
       X<-matrix(0,length(x),object$bs.dim)
       for (i in 1:(m+1)) X[,i] <- x^(i-1)
       for (i in 1:nk) X[,i+m+1] <- (x-k[i])^m*as.numeric(x>k[i])
       X # return the prediction matrix
     }