########## R script: union-wage.Rs ##########

# For illustrating penalized spline fitting
# of ordinary scatterplot data. Note that the
# R packages "nlme" and "MASS" are required as 
# well as the file default.knots.r containing 
# the function default.knots().

# Last changed: 02 APR 2004 by M.P. Wand

create.pdf <- F

basis <- "truncated lines"
#basis <- "radial cubics"

library(MASS)
library(nlme)

trade.union <- read.table("union-wage.dat",header=T)

x <- trade.union$wage
y <- trade.union$union.member
source("default.knots.r")
knots <- default.knots(x)
n <- length(y)

plot(x,y,type="n",bty="l",xlab="wage",ylab="P(union membership)")
rug(x[y==0])
rug(x[y==1],side=3)

X <- cbind(rep(1,n),x)

if (basis=="truncated lines")
{
   Z <- outer(x,knots,"-")    
   Z <- Z*(Z>0) 
}

if (basis=="radial cubics")
{
   svd.Omega <- svd(abs(outer(knots,knots,"-"))^3)
   matrix.sqrt.Omega <- t(svd.Omega$v %*% (t(svd.Omega$u)
                                     *sqrt(svd.Omega$d)))
   Z <- t(solve(matrix.sqrt.Omega,t(abs(outer(x,knots,"-")^3))))
}

group <- as.factor(rep(1,length(y)))
data.fr <- groupedData(y~x|group,data=data.frame(x,y))

fit <- glmmPQL(y~-1+X,random=list(group=pdIdent(~-1+Z)),
               data=data.fr,family=binomial)

beta.hat <- fit$coef$fixed
u.hat <- unlist(fit$coef$random)

n.g <- 101
x.g <- seq(min(x),max(x),length=n.g)
X.g <- cbind(rep(1,n.g),x.g)

if (basis=="truncated lines")
{
   Z.g <- outer(x.g,knots,"-")
   Z.g <- Z.g*(Z.g>0)
}

if (basis=="radial cubics")
   Z.g <- t(solve(matrix.sqrt.Omega,t(abs(outer(x.g,knots,"-")^3))))

eta.g <- X.g%*%beta.hat + Z.g%*%u.hat
mu.g <-  1/(1+exp(-(eta.g)))

if (create.pdf)
   pdf("union-wage.pdf",width=8)

plot(x,y,type="n",bty="l",xlab="wage",ylab="P(union membership)")
rug(x[y==0])
rug(x[y==1],side=3)

lines(x.g,mu.g,lwd=4,col="magenta")

########## End of union-wage ##########