# NOT RUN {
library(gamair); library(mgcv)
## 4.2.1
data(engine); attach(engine)
plot(size,wear,xlab="Engine capacity",ylab="Wear index")
tf <- function(x,xj,j) {
## generate jth tent function from set defined by knots xj
dj <- xj*0;dj[j] <- 1
approx(xj,dj,x)$y
}
tf.X <- function(x,xj) {
## tent function basis matrix given data x
## and knot sequence xj
nk <- length(xj); n <- length(x)
X <- matrix(NA,n,nk)
for (j in 1:nk) X[,j] <- tf(x,xj,j)
X
}
sj <- seq(min(size),max(size),length=6) ## generate knots
X <- tf.X(size,sj) ## get model matrix
b <- lm(wear~X-1) ## fit model
s <- seq(min(size),max(size),length=200)## prediction data
Xp <- tf.X(s,sj) ## prediction matrix
plot(size,wear) ## plot data
lines(s,Xp%*%coef(b)) ## overlay estimated f
## 4.2.2
prs.fit <- function(y,x,xj,sp) {
X <- tf.X(x,xj) ## model matrix
D <- diff(diag(length(xj)),differences=2) ## sqrt penalty
X <- rbind(X,sqrt(sp)*D) ## augmented model matrix
y <- c(y,rep(0,nrow(D))) ## augmented data
lm(y~X-1) ## penalized least squares fit
}
sj <- seq(min(size),max(size),length=20) ## knots
b <- prs.fit(wear,size,sj,2) ## penalized fit
plot(size,wear) ## plot data
Xp <- tf.X(s,sj) ## prediction matrix
lines(s,Xp%*%coef(b)) ## plot the smooth
## 4.2.3
rho = seq(-9,11,length=90)
n <- length(wear)
V <- rep(NA,90)
for (i in 1:90) { ## loop through smoothing params
b <- prs.fit(wear,size,sj,exp(rho[i])) ## fit model
trF <- sum(influence(b)$hat[1:n]) ## extract EDF
rss <- sum((wear-fitted(b)[1:n])^2) ## residual SS
V[i] <- n*rss/(n-trF)^2 ## GCV score
}
plot(rho,V,type="l",xlab=expression(log(lambda)),
main="GCV score")
sp <- exp(rho[V==min(V)]) ## extract optimal sp
b <- prs.fit(wear,size,sj,sp) ## re-fit
plot(size,wear,main="GCV optimal fit")
lines(s,Xp%*%coef(b))
## 4.2.3 mixed model connection
## copy of llm from 2.2.4...
llm <- function(theta,X,Z,y) {
## untransform parameters...
sigma.b <- exp(theta[1])
sigma <- exp(theta[2])
## extract dimensions...
n <- length(y); pr <- ncol(Z); pf <- ncol(X)
## obtain \hat \beta, \hat b...
X1 <- cbind(X,Z)
ipsi <- c(rep(0,pf),rep(1/sigma.b^2,pr))
b1 <- solve(crossprod(X1)/sigma^2+diag(ipsi),
t(X1)%*%y/sigma^2)
## compute log|Z'Z/sigma^2 + I/sigma.b^2|...
ldet <- sum(log(diag(chol(crossprod(Z)/sigma^2 +
diag(ipsi[-(1:pf)])))))
## compute log profile likelihood...
l <- (-sum((y-X1%*%b1)^2)/sigma^2 - sum(b1^2*ipsi) -
n*log(sigma^2) - pr*log(sigma.b^2) - 2*ldet - n*log(2*pi))/2
attr(l,"b") <- as.numeric(b1) ## return \hat beta and \hat b
-l
}
X0 <- tf.X(size,sj) ## X in original parameterization
D <- rbind(0,0,diff(diag(20),difference=2))
diag(D) <- 1 ## augmented D
X <- t(backsolve(t(D),t(X0))) ## re-parameterized X
Z <- X[,-c(1,2)]; X <- X[,1:2] ## mixed model matrices
## estimate smoothing and variance parameters...
m <- optim(c(0,0),llm,method="BFGS",X=X,Z=Z,y=wear)
b <- attr(llm(m$par,X,Z,wear),"b") ## extract coefficients
## plot results...
plot(size,wear)
Xp1 <- t(backsolve(t(D),t(Xp))) ## re-parameterized pred. mat.
lines(s,Xp1%*%as.numeric(b),col="grey",lwd=2)
library(nlme)
g <- factor(rep(1,nrow(X))) ## dummy factor
m <- lme(wear~X-1,random=list(g=pdIdent(~Z-1)))
lines(s,Xp1%*%as.numeric(coef(m))) ## and to plot
## 4.3.1 Additive
tf.XD <- function(x,xk,cmx=NULL,m=2) {
## get X and D subject to constraint
nk <- length(xk)
X <- tf.X(x,xk)[,-nk] ## basis matrix
D <- diff(diag(nk),differences=m)[,-nk] ## root penalty
if (is.null(cmx)) cmx <- colMeans(X)
X <- sweep(X,2,cmx) ## subtract cmx from columns
list(X=X,D=D,cmx=cmx)
} ## tf.XD
am.fit <- function(y,x,v,sp,k=10) {
## setup bases and penalties...
xk <- seq(min(x),max(x),length=k)
xdx <- tf.XD(x,xk)
vk <- seq(min(v),max(v),length=k)
xdv <- tf.XD(v,vk)
## create augmented model matrix and response...
nD <- nrow(xdx$D)*2
sp <- sqrt(sp)
X <- cbind(c(rep(1,nrow(xdx$X)),rep(0,nD)),
rbind(xdx$X,sp[1]*xdx$D,xdv$D*0),
rbind(xdv$X,xdx$D*0,sp[2]*xdv$D))
y1 <- c(y,rep(0,nD))
## fit model..
b <- lm(y1~X-1)
## compute some useful quantities...
n <- length(y)
trA <- sum(influence(b)$hat[1:n]) ## EDF
rsd <- y-fitted(b)[1:n] ## residuals
rss <- sum(rsd^2) ## residual SS
sig.hat <- rss/(n-trA) ## residual variance
gcv <- sig.hat*n/(n-trA) ## GCV score
Vb <- vcov(b)*sig.hat/summary(b)$sigma^2 ## coeff cov matrix
## return fitted model...
list(b=coef(b),Vb=Vb,edf=trA,gcv=gcv,fitted=fitted(b)[1:n],
rsd=rsd,xk=list(xk,vk),cmx=list(xdx$cmx,xdv$cmx))
} ## am.fit
am.gcv <- function(lsp,y,x,v,k) {
## function suitable for GCV optimization by optim
am.fit(y,x,v,exp(lsp),k)$gcv
}
## find GCV optimal smoothing parameters...
fit <- optim(c(0,0), am.gcv, y=trees$Volume, x=trees$Girth,
v=trees$Height,k=10)
sp <- exp(fit$par) ## best fit smoothing parameters
## Get fit at GCV optimal smoothing parameters...
fit <- am.fit(trees$Volume,trees$Girth,trees$Height,sp,k=10)
am.plot <- function(fit,xlab,ylab) {
## produces effect plots for simple 2 term
## additive model
start <- 2 ## where smooth coeffs start in beta
for (i in 1:2) {
## sequence of values at which to predict...
x <- seq(min(fit$xk[[i]]),max(fit$xk[[i]]),length=200)
## get prediction matrix for this smooth...
Xp <- tf.XD(x,fit$xk[[i]],fit$cmx[[i]])$X
## extract coefficients and cov matrix for this smooth
stop <- start + ncol(Xp)-1; ind <- start:stop
b <- fit$b[ind];Vb <- fit$Vb[ind,ind]
## values for smooth at x...
fv <- Xp%*%b
## standard errors of smooth at x....
se <- rowSums((Xp%*%Vb)*Xp)^.5
## 2 s.e. limits for smooth...
ul <- fv + 2*se;ll <- fv - 2 * se
## plot smooth and limits...
plot(x,fv,type="l",ylim=range(c(ul,ll)),xlab=xlab[i],
ylab=ylab[i])
lines(x,ul,lty=2);lines(x,ll,lty=2)
start <- stop + 1
}
} ## am.plot
par(mfrow=c(1,3))
plot(fit$fitted,trees$Vol,xlab="fitted volume ",
ylab="observed volume")
am.plot(fit,xlab=c("Girth","Height"),
ylab=c("s(Girth)","s(Height)"))
## 4.4 Generalized additive
gam.fit <- function(y,x,v,sp,k=10) {
## gamma error log link 2 term gam fit...
eta <- log(y) ## get initial eta
not.converged <- TRUE
old.gcv <- -100 ## don't converge immediately
while (not.converged) {
mu <- exp(eta) ## current mu estimate
z <- (y - mu)/mu + eta ## pseudodata
fit <- am.fit(z,x,v,sp,k) ## penalized least squares
if (abs(fit$gcv-old.gcv)<1e-5*fit$gcv) {
not.converged <- FALSE
}
old.gcv <- fit$gcv
eta <- fit$fitted ## updated eta
}
fit$fitted <- exp(fit$fitted) ## mu
fit
} ## gam.fit
gam.gcv <- function(lsp,y,x,v,k=10) {
gam.fit(y,x,v,exp(lsp),k=k)$gcv
}
fit <- optim(c(0,0),gam.gcv,y=trees$Volume,x=trees$Girth,
v=trees$Height,k=10)
sp <- exp(fit$par)
fit <- gam.fit(trees$Volume,trees$Girth,trees$Height,sp)
par(mfrow=c(1,3))
plot(fit$fitted,trees$Vol,xlab="fitted volume ",
ylab="observed volume")
am.plot(fit,xlab=c("Girth","Height"),
ylab=c("s(Girth)","s(Height)"))
## 4.6 mgcv
library(mgcv) ## load the package
library(gamair) ## load the data package
data(trees)
ct1 <- gam(Volume~s(Height)+s(Girth),
family=Gamma(link=log),data=trees)
ct1
plot(ct1,residuals=TRUE)
## 4.6.1
ct2 <- gam(Volume~s(Height,bs="cr")+s(Girth,bs="cr"),
family=Gamma(link=log),data=trees)
ct2
ct3 <- gam(Volume ~ s(Height) + s(Girth,bs="cr",k=20),
family=Gamma(link=log),data=trees)
ct3
ct4 <- gam(Volume ~ s(Height) + s(Girth),
family=Gamma(link=log),data=trees,gamma=1.4)
ct4
plot(ct4,residuals=TRUE)
## 4.6.2
ct5 <- gam(Volume ~ s(Height,Girth,k=25),
family=Gamma(link=log),data=trees)
ct5
plot(ct5,too.far=0.15)
ct6 <- gam(Volume ~ te(Height,Girth,k=5),
family=Gamma(link=log),data=trees)
ct6
plot(ct6,too.far=0.15)
## 4.6.3
gam(Volume~Height+s(Girth),family=Gamma(link=log),data=trees)
trees$Hclass <- factor(floor(trees$Height/10)-5,
labels=c("small","medium","large"))
ct7 <- gam(Volume ~ Hclass+s(Girth),
family=Gamma(link=log),data=trees)
par(mfrow=c(1,2))
plot(ct7,all.terms=TRUE)
anova(ct7)
AIC(ct7)
summary(ct7)
# }
Run the code above in your browser using DataCamp Workspace