https://github.com/cran/fda
Tip revision: ca5f7615f07fc7130085233d37e769f62fdb9f73 authored by J. O. Ramsay on 03 October 2008, 00:00:00 UTC
version 2.0.3
version 2.0.3
Tip revision: ca5f761
linmod.R
linmod <- function(xfdobj, yfdobj, wtvec=rep(1,nrep),
xLfdobj=int2Lfd(2), yLfdobj=int2Lfd(2),
xlambda=0, ylambda=0)
{
# This function fits an unrestricted functional linear model for a
# a functional dependent variable and a single function covariate.
# Smoothing is controlled by two parameters XLAMBDA and YLAMBDA,
# corresponding to the independent and dependent functional
# variables, respectively.
# Argument:
# XFDOBJ ... If the independent variable is multivariate, a design matrix.
# If the independent variable is functional, a "fd" object.
# YFDOBJ ... If the dependent variable is multivariate, a design matrix.
# If the dependent variable is functional, a "fd" object.
# WTVEC ... a vector of weights
# XLFDOBJ ... A linear differential operator object for the independent
# variable.
# YLFDOBJ ... A linear differential operator object for the dependent
# variable.
# XLAMBDA ... a smoothing parameter for the independent variable
# YLAMBDA ... a smoothing parameter for the dependent variable
# ZMATRNK ... actual rank of independent variable matrix for the
# functional DV/multivariate IV case
# Returns: a list containing
# ALPHA ... a vector of intercept values
# REGFD ... a functional data object for the regression function
# Last modified 2007.11.28 by Spencer Graves
# Previously modified: 4 December 2005
# check XLfdobj and YLfdobj
xLfdobj <- int2Lfd(xLfdobj)
yLfdobj <- int2Lfd(yLfdobj)
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# The multivariate IV and functional DV case
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (inherits(yfdobj, "fd") && !inherits(xfdobj, "fd")) {
stop ("Use function fRegress for this model.")
}
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# The functional IV and multivariate DV case
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (inherits(xfdobj, "fd") && !(inherits(yfdobj, "fd"))) {
stop ("Use function fRegress for this model.")
}
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# The functional IV and functional DV case
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (inherits(xfdobj, "fd") && inherits(yfdobj, "fd")) {
coefx <- xfdobj$coefs
coefy <- yfdobj$coefs
coefdx <- dim(coefx)
coefdy <- dim(coefy)
ndimx <- length(coefdx)
ndimy <- length(coefdy)
if (ndimx < 2) stop(
"Linear modeling impossible with 1 replication")
if (ndimx == 3) stop(
"This version cannot accommodate multiple functional IVs")
nrep <- coefdx[2]
if (coefdy[2] != nrep) stop (
"Numbers of observations in the first two arguments do not match.")
rangewt <- range(wtvec)
if (rangewt[1] < 0) stop("WTVEC must not contain negative values.")
if (rangewt[1] == rangewt[2]) wtvar <- FALSE else wtvar <- TRUE
xbasisobj <- xfdobj$basis
nbasisx <- xbasisobj$nbasis
typex <- xbasisobj$type
rangevx <- xbasisobj$rangeval
ybasisobj <- yfdobj$basis
nbasisy <- ybasisobj$nbasis
typey <- ybasisobj$type
rangevy <- ybasisobj$rangeval
if (length(wtvec) != nrep) stop("WTVEC of wrong length")
if (min(wtvec) <= 0) stop("All values of WTVEC must be positive.")
if (xlambda < 0) warning (
"Value of LAMBDA was negative, and 0 used instead.")
jmatx <- inprod(xbasisobj, xbasisobj)
penmatx <- inprod(xbasisobj, xbasisobj, xLfdobj, xLfdobj)
if (ndimx == 2) {
zmatx <- t(rbind(matrix(1,1,nrep),jmatx %*% coefx))
} else {
zmatx <- t(rbind(matrix(1,1,nrep),jmatx %*% coefx[,,1]))
}
jmaty <- inprod(ybasisobj, ybasisobj)
penmaty <- inprod(ybasisobj, ybasisobj, yLfdobj, yLfdobj)
alpha <- rep(0,nbasisy)
index <- 2:(nbasisx+1)
kmatx <- matrix(0,nbasisx+1,nbasisx+1)
kmatx[index,index] <- penmatx
if (wtvar) {
# zmatw <- sweep(zmatx,2,wtmat,"*")
zmatw <- sweep(zmatx,2,wtvec,"*")
tempx <- solve(crossprod(zmatw,zmatx) + xlambda*kmatx)
tempy <- solve(jmaty + ylambda*penmaty)
gmat <- tempx %*% crossprod(zmatw,t(coefy[,,1])) %*% jmaty %*% tempy
} else {
tempx <- solve(crossprod(zmatx) + xlambda*kmatx)
tempy <- solve(jmaty + ylambda*penmaty)
if (ndimy == 2) {
gmat <- tempx %*% crossprod(zmatx,t(coefy)) %*% jmaty %*% tempy
} else {
gmat <- tempx %*% crossprod(zmatx,t(coefy[,,1])) %*% jmaty %*% tempy
}
}
yhatcoef <- t(zmatx %*% gmat)
# bcoef <- matrix(0,nbasisx,nbasisy)
bcoef <- gmat[index,]
alpha <- as.matrix(gmat[1,])
fdnames <- yfdobj$fdnames
alphafdnames <- fdnames
alphafdnames[[2]] <- "Intercept"
alphafd <- fd(alpha, ybasisobj, alphafdnames)
regfd <- bifd(bcoef, xbasisobj, ybasisobj)
yhatfd <- fd(yhatcoef,ybasisobj, fdnames)
linmodlist <- list(alphafd=alphafd, regfd=regfd, yhatfd=yhatfd)
return( linmodlist )
}
}