#################################################################################
#################################################################################
# R code associated with 
# "Functional Generalized Linear Models with Images as Predictors"
# by Philip T. Reiss and R. Todd Ogden, submitted to Biometrics
# November 20, 2008
#
# This R function can be used to fit the functional principal component
# regression model described in the paper.
#
# This code is distributed freely with no warranties of any kind.  It has not
# been checked for bugs as thoroughly as we would have liked.  If you use the
# code, we'd very much appreciate hearing from you about it.  Please contact the
# first author at phil.reiss@nyumc.org for further information or code, to 
# provide feedback, or to report problems.

#################################################################################
# fpcr.mgcv
#################################################################################
# Fit functional principal component regression using the gam function from the 
# mgcv package
# 
# Arguments:
# y			scalar outcome vector
# sigmat	matrix of signals: n x N, where n is the length of y and N is the
#				number of sites at which each signal is defined
# basis		N x K matrix of values of K basis functions at the N sites
# halfpens  a list, each component of which is a matrix M with K columns such 
#				that M'M defines a penalty on the basis coefficients
# nc		number of principal components
# covt		covariates: an n-row matrix, or a vector of length n 
# mean.signal.term	logical: should the mean of each subject's signal be 
#					included as a covariate?
# The remaining arguments are passed to the mgcv function gam; see the gam 
# documentation for details.  Please contact the authors if further information
# is required regarding the "basis" and "halfpen" arguments.
# 
# Value:
# An object of class "gam" (see gamObject in the mgcv documentation), with two 
# additional components: $nc (same as argument nc) and $fhat (coefficient 
# function estimate).  

fpcr.mgcv <- function(y, sigmat, basis, halfpens, nc, covt = NULL, mean.signal.term = FALSE, 
                 family = gaussian(), gamma = 1, sp=NULL, min.sp=NULL) {
    require(mgcv)

    # Design matrix
    X0 <- if (mean.signal.term) cbind(1, rowMeans(sigmat)) 
          else matrix(1, length(y), 1)    
    if (!is.null(covt)) X0 <- cbind(X0, as.matrix(covt))
    n.unpen.cols <- 1 + mean.signal.term + if (!is.null(covt)) NCOL(covt) else 0
    # Decorrelate the signals from the other columns of the design matrix
    sigs.decor <- if (ncol(X0) == 1) scale(sigmat, center = T, scale = F) 
                  else lm(sigmat ~ X0 - 1)$resid
    SB <- sigs.decor %*% basis
    V.nc <- svd(SB)$v[ , 1 : nc]
    X <- cbind(X0, SB %*% V.nc)

    # Define list of penalties for paraPen
    S <- list()
    for (i in 1:length(halfpens)) {
        S[[i]] <- matrix(0, ncol(X), ncol(X))
        S[[i]][-(1:n.unpen.cols), -(1:n.unpen.cols)] <- 
             crossprod(halfpens[[i]] %*% V.nc) 
    }

    obje = gam(y~X-1, paraPen=list(X=S), family=family, gamma=gamma, sp=sp,    
               min.sp=min.sp)
    obje$fhat = basis %*% V.nc %*% obje$coef[-(1:n.unpen.cols)]
    obje$nc = nc
    obje
}