fRegress <- function(y, xfdlist, betalist, wt=NULL, y2cMap=NULL, SigmaE=NULL) { # FREGRESS Fits a functional linear model using multiple # functional independent variables with the dependency being # pointwise or concurrent. # The case of a scalar independent variable is included by treating # it as a functional independent variable with a constant basis # and a unit coefficient. # # Arguments: # Y ... an object for the dependent variable, # which may be: # a functional data object, # a functional parameter (fdPar) object, or # a vector # XFDLIST ... a list object of length p with each list # containing an object for an independent variable. # the object may be: # a functional data object or # a vector # if XFDLIST is a functional data object or a vector, # it is converted to a list of length 1. # BETALIST ... a list object of length p with each list # containing a functional parameter object for # the corresponding regression function. If any of # these objects is a functional data object, it is # converted to the default functional parameter object. # if BETALIST is a functional parameter object # it is converted to a list of length 1. # WT ... a vector of nonnegative weights for observations # Y2CMAP ... the matrix mapping from the vector of observed values # to the coefficients for the dependent variable. # This is output by function SMOOTH_BASIS. If this is # supplied, confidence limits are computed, otherwise not. # SIGMAE ... Estimate of the covariances among the residuals. This # can only be estimated after a preliminary analysis # with FREGRESS. # # Returns FREGRESSLIST ... A list containing seven members with names: # yfdPar ... first argument of FREGRESS # xfdlist ... second argument of FREGRESS # betalist ... third argument of FREGRESS # betaestlist ... estimated regression functions # yhatfdobj ... functional data object containing fitted functions # Cmat ... coefficient matrix for the linear system defining # the regression coefficient basis vector # Dmat ... right side vector for the linear system defining # the regression coefficient basis vector # Cmatinv ... inverse of the coefficient matrix, needed for # function FREGRESS.STDERR that computes standard errors # wt ... weights for observations # df ... degrees of freedom for fit # Last modified 29 October 2009 by Jim Ramsay arglist <- fRegressArgCheck(y, xfdlist, betalist, wt) yfdPar <- arglist$yfdPar xfdlist <- arglist$xfdlist betalist <- arglist$betalist wt <- arglist$wt p <- length(xfdlist) N <- dim(xfdlist[[1]]$coef)[1] wtconst <- var(wt) == 0 # -------------------------------------------------------------------------- # branch depending on whether the dependent variable is functional or scalar # -------------------------------------------------------------------------- if (inherits(yfdPar, "fdPar")) { # ---------------------------------------------------------------- # YFDPAR is a functional parameter object # ---------------------------------------------------------------- # extract dependent variable information yfdobj <- yfdPar$fd ylambda <- yfdPar$lambda yLfdobj <- yfdPar$Lfd ycoef <- yfdobj$coefs ycoefdim <- dim(ycoef) N <- ycoefdim[2] ybasisobj <- yfdobj$basis rangeval <- ybasisobj$rangeval ynbasis <- ybasisobj$nbasis if (length(ycoefdim) > 2) stop("YFDOBJ from YFDPAR is not univariate.") # -------- set up the linear equations for the solution ----------- # compute the total number of coefficients to be estimated ncoef <- 0 for (j in 1:p) { betafdParj <- betalist[[j]] if (betafdParj$estimate) { ncoefj <- betafdParj$fd$basis$nbasis ncoef <- ncoef + ncoefj } } Cmat <- matrix(0,ncoef,ncoef) Dmat <- rep(0,ncoef) # loop through rows of CMAT mj2 <- 0 for (j in 1:p) { betafdParj <- betalist[[j]] if (betafdParj$estimate) { betafdj <- betafdParj$fd betabasisj <- betafdj$basis ncoefj <- betabasisj$nbasis # row indices of CMAT and DMAT to fill mj1 <- mj2 + 1 mj2 <- mj2 + ncoefj indexj <- mj1:mj2 # compute right side of equation DMAT xfdj <- xfdlist[[j]] if (wtconstant) { xyfdj <- xfdj*yfdobj } else { xyfdj <- (xfdj*wt)*yfdobj } wtfdj <- sum(xyfdj) Dmatj <- inprod(betabasisj,onesfd,0,0,rangeval,wtfdj) Dmat[indexj] <- Dmatj # loop through columns of CMAT mk2 <- 0 for (k in 1:j) { betafdPark <- betalist[[k]] if (betafdPark$estimate) { betafdk <- betafdPark$fd betabasisk <- betafdk$basis ncoefk <- betabasisk$nbasis # column indices of CMAT to fill mk1 <- mk2 + 1 mk2 <- mk2 + ncoefk indexk <- mk1:mk2 # set up two weight functions xfdk <- xfdlist[[k]] if (wtconstant) { xxfdjk <- xfdj*xfdk } else { xxfdjk <- (xfdj*wt)*xfdk } wtfdjk <- sum(xxfdjk) Cmatjk <- inprod(betabasisj, betabasisk, 0, 0, rangeval, wtfdjk) Cmat[indexj,indexk] <- Cmatjk Cmat[indexk,indexj] <- t(Cmatjk) } } # attach penalty term to diagonal block lambda <- betafdParj$lambda if (lambda > 0) { Lfdj <- betafdParj$Lfd Rmatj <- eval.penalty(betafdParj$fd$basis, Lfdj) Cmat[indexj,indexj] <- Cmat[indexj,indexj] + lambda*Rmatj } } } Cmat <- (Cmat+t(Cmat))/2 # check Cmat for singularity eigchk(Cmat) # solve for coefficients defining BETA Lmat <- chol(Cmat) Lmatinv <- solve(Lmat) Cmatinv <- Lmatinv %*% t(Lmatinv) betacoef <- Cmatinv %*% Dmat # set up fdPar objects for reg. fns. in BETAESTLIST betaestlist <- betalist mj2 <- 0 for (j in 1:p) { betafdParj <- betalist[[j]] if (betafdParj$estimate) { betafdj <- betafdParj$fd ncoefj <- betafdj$basis$nbasis mj1 <- mj2 + 1 mj2 <- mj2 + ncoefj indexj <- mj1:mj2 coefj <- betacoef[indexj] betafdj$coefs <- as.matrix(coefj) betafdParj$fd <- betafdj } betaestlist[[j]] <- betafdParj } # set up fd objects for predicted values in YHATFDOBJ nfine <- max(501,10*ynbasis+1) tfine <- seq(rangeval[1], rangeval[2], len=nfine) yhatmat <- matrix(0,nfine,N) for (j in 1:p) { xfdj <- xfdlist[[j]] xmatj <- eval.fd(tfine, xfdj) betafdParj <- betaestlist[[j]] betafdj <- betafdParj$fd betavecj <- eval.fd(tfine, betafdj) yhatmat <- yhatmat + xmatj*as.vector(betavecj) } yhatfdobj <- smooth.basis(tfine, yhatmat, ybasisobj)$fd df <- NA # ----------------------------------------------------------------------- # Compute pointwise standard errors of regression coefficients # if both y2cMap and SigmaE are supplied. # ----------------------------------------------------------------------- # ------------------------------------------------------------------- # Set up output list object # ------------------------------------------------------------------- fRegressList <- list(yfdPar = yfdPar, xfdlist = xfdlist, betalist = betalist, betaestlist = betaestlist, yhatfdobj = yhatfdobj, Cmat = Cmat, Dmat = Dmat, Cmatinv = Cmatinv, wt = wt, df = df, y2cMap = y2cMap, SigmaE = SigmaE, betastderrlist = betastderrlist, bvar = bvar, c2bMap = c2bMap) } # ------------------------------------------------------------------- if (inherits(yfdPar,"numeric")) { # ---------------------------------------------------------------- # YFDPAR is scalar or multivariate # ---------------------------------------------------------------- ymat <- as.matrix(yfdPar) N <- dim(ymat)[1] Zmat <- NULL Rmat <- NULL pjvec <- rep(0,p) ncoef <- 0 for (j in 1:p) { xfdj <- xfdlist[[j]] xcoef <- xfdj$coefs xbasis <- xfdj$basis betafdParj <- betalist[[j]] bbasis <- betafdParj$fd$basis bnbasis <- bbasis$nbasis pjvec[j] <- bnbasis Jpsithetaj <- inprod(xbasis,bbasis) Zmat <- cbind(Zmat,crossprod(xcoef,Jpsithetaj)) if (betafdParj$estimate) { lambdaj <- betafdParj$lambda if (lambdaj > 0) { Lfdj <- betafdParj$Lfd Rmatj <- lambdaj*eval.penalty(bbasis,Lfdj) } else { Rmatj <- matrix(0,bnbasis,bnbasis) } if (ncoef > 0) { zeromat <- matrix(0,ncoef,bnbasis) Rmat <- rbind(cbind(Rmat, zeromat), cbind(t(zeromat), Rmatj)) } else { Rmat <- Rmatj ncoef <- ncoef + bnbasis } } } # ----------------------------------------------------------- # set up the linear equations for the solution # ----------------------------------------------------------- # solve for coefficients defining BETA if (any(wt != 1)) { rtwt <- sqrt(wt) Zmatwt <- Zmat*rtwt ymatwt <- ymat*rtwt Cmat <- t(Zmatwt) %*% Zmatwt + Rmat Dmat <- t(Zmatwt) %*% ymatwt } else { Cmat <- t(Zmat) %*% Zmat + Rmat Dmat <- t(Zmat) %*% ymat } eigchk(Cmat) Cmatinv <- solve(Cmat) betacoef <- Cmatinv %*% Dmat # compute and print degrees of freedom measure df <- sum(diag(Zmat %*% Cmatinv %*% t(Zmat))) # set up fdPar object for BETAESTFDPAR betaestlist <- betalist mj2 <- 0 for (j in 1:p) { betafdParj <- betalist[[j]] betafdj <- betafdParj$fd ncoefj <- betafdj$basis$nbasis mj1 <- mj2 + 1 mj2 <- mj2 + ncoefj indexj <- mj1:mj2 betacoefj <- betacoef[indexj] betaestfdj <- betafdj betaestfdj$coefs <- as.matrix(betacoefj) betaestfdParj <- betafdParj betaestfdParj$fd <- betaestfdj betaestlist[[j]] <- betaestfdParj } # set up fd object for predicted values yhatmat <- matrix(0,N,1) for (j in 1:p) { xfdj <- xfdlist[[j]] if (inherits(xfdj, "fd")) { xbasis <- xfdj$basis xnbasis <- xbasis$nbasis xrng <- xbasis$rangeval nfine <- max(501,10*xnbasis+1) tfine <- seq(xrng[1], xrng[2], len=nfine) deltat <- tfine[2]-tfine[1] xmat <- eval.fd(tfine, xfdj) betafdParj <- betaestlist[[j]] betafdj <- betafdParj$fd betamat <- eval.fd(tfine, betafdj) fitj <- deltat*(crossprod(xmat,betamat) - 0.5*(outer(xmat[1, ],betamat[1, ]) + outer(xmat[nfine,],betamat[nfine,]))) yhatmat <- yhatmat + fitj } else{ betaestfdParj <- betaestlist[[j]] betavecj <- betaestfdParj$fd$coefs yhatmat <- yhatmat + xfdj %*% t(betavecj) } } yhatfdobj <- yhatmat # ----------------------------------------------------------------------- # Compute pointwise standard errors of regression coefficients # if both y2cMap and SigmaE are supplied. # ----------------------------------------------------------------------- if (!(is.null(y2cMap) || is.null(SigmaE))) { # check dimensions of y2cMap and SigmaE y2cdim <- dim(y2cMap) if (y2cdim[2] != dim(SigmaE)[1]) stop( "Dimensions of Y2CMAP not correct.") # compute linear mapping c2bMap takinging coefficients for # response into coefficients for regression functions c2bMap <- Cmatinv %*% t(Zmat) y2bmap <- c2bMap bvar <- y2bmap %*% as.matrix(SigmaE) %*% t(y2bmap) betastderrlist <- vector("list",p) mj2 <- 0 for (j in 1:p) { betafdParj <- betalist[[j]] betabasisj <- betafdParj$fd$basis ncoefj <- betabasisj$nbasis mj1 <- mj2 + 1 mj2 <- mj2 + ncoefj indexj <- mj1:mj2 bvarj <- bvar[indexj,indexj] betarng <- betabasisj$rangeval nfine <- max(c(501,10*ncoefj+1)) tfine <- seq(betarng[1], betarng[2], len=nfine) bbasismat <- eval.basis(tfine, betabasisj) bstderrj <- sqrt(diag(bbasismat %*% bvarj %*% t(bbasismat))) bstderrfdj <- smooth.basis(tfine, bstderrj, betabasisj)$fd betastderrlist[[j]] <- bstderrfdj } } else { betastderrlist = NULL bvar = NULL c2bMap = NULL } # ----------------------------------------------------------------------- # Set up output list object # ----------------------------------------------------------------------- fRegressList <- list(yfdPar = yfdPar, xfdlist = xfdlist, betalist = betalist, betaestlist = betaestlist, yhatfdobj = yhatfdobj, Cmat = Cmat, Dmat = Dmat, Cmatinv = Cmatinv, wt = wt, df = df, y2cMap = y2cMap, SigmaE = SigmaE, betastderrlist = betastderrlist, bvar = bvar, c2bMap = c2bMap) } return(fRegressList) } # ------------------------------------------------------------------------ eigchk <- function(Cmat) { # check Cmat for singularity eigval <- eigen(Cmat)$values ncoef <- length(eigval) if (eigval[ncoef] < 0) { neig <- min(length(eigval),10) cat("\nSmallest eigenvalues:\n") print(eigval[(ncoef-neig+1):ncoef]) cat("\nLargest eigenvalues:\n") print(eigval[1:neig]) stop("Negative eigenvalue of coefficient matrix.") } if (eigval[ncoef] == 0) stop("Zero eigenvalue of coefficient matrix.") logcondition <- log10(eigval[1]) - log10(eigval[ncoef]) if (logcondition > 12) { warning("Near singularity in coefficient matrix.") cat(paste("\nLog10 Eigenvalues range from\n", log10(eigval[ncoef])," to ",log10(eigval[1]),"\n")) } }