########## 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") }