########## R script: glmmPQL-tst ##########

# For illustrating how glmmPQL() can be used
# to do semiparametric regresssion (at least
# for the simplest case for now).


# Set the type of model to be fit.
  
do.scatter.smooth <- T

# Specify high level settings

family.val <- "poisson"
#family.val <- "binomial"

#basis <- "truncated lines"
basis <- "thin plate splines"

# Load the MASS library to access glmmPQL()

library(MASS)

if (do.scatter.smooth)
{

   # Specify the settings

   seed <- 49032
   set.seed(seed)
   n <- 1000

   etafun <- function(x)
      return(cos(4*pi*x))

   num.knots <- 15

   # Generate data

   x <- runif(n)

   if (family.val=="binomial")
   {
      mu <- 1/(1+exp(-(etafun(x))))
      y <- rbinom(n,1,mu)
   }

   if (family.val=="poisson")
   {
      mu <- exp(etafun(x))
      y <- rpois(n,mu)
   }

   # Set knots

   knots <- seq(0,1,length=(num.knots+2))[-c(1,(num.knots+2))]
   knots <- quantile(unique(x),knots)

   # Construct X and Z matrices
     
   X <- cbind(rep(1,n),x)

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

   if (basis=="thin plate splines")
   {
      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))))
   }

   # Obtain the fit using glmmPQL()

   g <-as.factor(rep(1,length(y)))
   dat <-data.frame(y=y,x=x)

   if (family.val=="binomial")
      fit <- glmmPQL(y~-1+X,dat,random=list(g=pdIdent(~Z-1)),family=binomial)

   if (family.val=="poisson")
      fit <- glmmPQL(y~-1+X,dat,random=list(g=pdIdent(~Z-1)),family=poisson)

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

   eta.hat <- as.vector(X%*%beta.hat + Z%*%u.hat)

   if (family.val=="binomial")
      mu.hat <-  1/(1+exp(-(eta.hat)))

   if (family.val=="poisson")
      mu.hat <-  exp(eta.hat)

   # Plot the results

   plot(x,y) 
   lines(x[order(x)],mu.hat[order(x)],lwd=2,col="red")
   lines(x[order(x)],mu[order(x)],lwd=2,col="blue")
}