## ----echo=FALSE,eval=FALSE----------------------------------------------------
#  options(width=80)

## ----eval=FALSE---------------------------------------------------------------
#  library(catdata)
#  data(unemployment, package="catdata")
#  attach(unemployment)

## ----eval=FALSE---------------------------------------------------------------
#  library(mgcv)

## ----eval=FALSE---------------------------------------------------------------
#  durbin[durbin==2] <- 0
#  
#  gamage <- gam(durbin ~ s(age, bs="ps", m=c(2,1), k=15), family=binomial())

## ----eval=FALSE---------------------------------------------------------------
#  minage <- min(age)
#  maxage <- max(age)
#  ageindex <- seq(from=minage, to=maxage, by=0.01)
#  
#  pred <- predict(gamage, newdata=data.frame("age"=ageindex), type="response")

## ----eval=FALSE---------------------------------------------------------------
#  bspline<-function(x,k,i,m=2)
#  {if (m==-1)
#  {res<-as.numeric(x<k[i+1]&x>=k[i])}
#  else{
#  z0<-(x-k[i])/(k[i+m+1]-k[i])
#  z1<-(k[i+m+2]-x)/(k[i+m+2]-k[i+1])
#  res<- z0*bspline(x,k,i,m-1)+z1*bspline(x,k,i+1,m-1)}
#  res}

## ----eval=FALSE---------------------------------------------------------------
#  k <- gamage$smooth[[1]]$knots
#  
#  meanage <- c()
#  
#  for (i in minage:maxage){
#  meanage[i] <- sum(durbin[age==i])
#  if(sum(durbin[age==i])!=0){
#  meanage[i] <- mean(durbin[age==i])
#  }
#  }

## ----fig.width=8,eval=FALSE---------------------------------------------------
#  par(cex=1.3, lwd=1.5)
#  plot(ageindex, pred, type="l", ylim=c(0,0.8), col="gray", xlab="age/years",
#       ylab="short-term unemployment")
#  for(i in 1:15)
#  {lines(ageindex,(0.5*gamage$coefficients[i+1]+0.7)*(bspline(x=ageindex,k=k,
#  i=i+1,m=2)),type="l", col="gray")}
#  points(minage:maxage, meanage[minage:maxage], pch=20)

## ----eval=FALSE---------------------------------------------------------------
#  gamage2 <- gam(durbin ~ s(age, bs="ps", fx=TRUE, m=c(2,1),k=15), family=binomial())
#  
#  pred2 <- predict(gamage2, newdata=data.frame("age"=ageindex), type="response")

## ----fig.width=8,eval=FALSE---------------------------------------------------
#  par(cex=1.3, lwd=1.5)
#  plot(ageindex, pred2, type="l", ylim=c(0,0.8), col="gray", xlab="age/years",
#       ylab="short-term unemployment")
#  for(i in 1:15)
#  {lines(ageindex, (0.1*gamage2$coefficients[i+1]+0.3)*
#    bspline(x=ageindex,k=k,i=i+1,m=2),type="l", col="gray")}
#  points(minage:maxage, meanage[minage:maxage],pch=20)

## ----echo=FALSE,results='hide',eval=FALSE-------------------------------------
#  detach(unemployment)
#  rm(unemployment)