Skip to content

Commit 5478b3c

Browse files
committed
fRegress.double()
1 parent fcffd4a commit 5478b3c

File tree

1 file changed

+269
-0
lines changed

1 file changed

+269
-0
lines changed

R/fRegress.double.R

Lines changed: 269 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,269 @@
1+
fRegress.double <- function(y, xfdlist, betalist, wt=NULL,
2+
y2cMap=NULL, SigmaE=NULL, returnMatrix=FALSE, ...)
3+
{
4+
5+
# FREGRESS.DOUBLE Fits a scalar dependent variable using the concurrent
6+
# functional regression model using inner products
7+
# of functional covariates and functional regression
8+
# functions.
9+
#
10+
# Arguments:
11+
# Y ... an object for the dependent variable,
12+
# which be a numeric vector
13+
# XFDLIST ... a list object of length p with each list
14+
# containing an object for an independent variable.
15+
# the object may be:
16+
# a functional data object or
17+
# a vector
18+
# if XFDLIST is a functional data object or a vector,
19+
# it is converted to a list of length 1.
20+
# BETALIST ... a list object of length p with each list
21+
# containing a functional parameter object for
22+
# the corresponding regression function. If any of
23+
# these objects is a functional data object, it is
24+
# converted to the default functional parameter object.
25+
# if BETALIST is a functional parameter object
26+
# it is converted to a list of length 1.
27+
# WT ... a vector of nonnegative weights for observations
28+
# Y2CMAP ... the matrix mapping from the vector of observed values
29+
# to the coefficients for the dependent variable.
30+
# This is output by function SMOOTH_BASIS. If this is
31+
# supplied, confidence limits are computed, otherwise not.
32+
# SIGMAE ... Estimate of the covariances among the residuals. This
33+
# can only be estimated after a preliminary analysis
34+
# with .
35+
# RETURNMATRIX ... If False, a matrix in sparse storage model can be returned
36+
# from a call to function BsplineS. See this function for
37+
# enabling this option.
38+
#
39+
# Returns LIST ... A list containing seven members with names:
40+
# yfdPar ... first argument of
41+
# xfdlist ... second argument of
42+
# betalist ... third argument of
43+
# betaestlist ... estimated regression functions
44+
# yhatfdobj ... functional data object containing fitted functions
45+
# Cmatinv ... inverse of the coefficient matrix, needed for
46+
# function .STDERR that computes standard errors
47+
# wt ... weights for observations
48+
# df ... degrees of freedom for fit
49+
# This list object is converted to a class with the name ""
50+
# function predict. is an example of a method that can be called simply
51+
# as predict(List). In this call List can be any object of the
52+
# "".
53+
54+
# Last modified 25 August 2020 by Jim Ramsay
55+
56+
# check Y and compute sample size N
57+
58+
print("inside fRegress.numeric")
59+
60+
if (!inherits(y, "numeric")) stop("Y is not a numeric vector.")
61+
62+
# ----------------------------------------------------------------
63+
# YFDPAR is scalar or multivariate
64+
# ----------------------------------------------------------------
65+
66+
arglist <- fRegressArgCheck(y, xfdlist, betalist, wt)
67+
68+
yfdPar <- arglist$yfdPar
69+
xfdlist <- arglist$xfdlist
70+
betalist <- arglist$betalist
71+
wt <- arglist$wt
72+
73+
ymat <- as.matrix(y)
74+
N <- dim(ymat)[1]
75+
p <- length(xfdlist)
76+
77+
print("preliminaries")
78+
79+
Zmat <- NULL
80+
Rmat <- NULL
81+
pjvec <- rep(0,p)
82+
ncoef <- 0
83+
for (j in 1:p) {
84+
print(j)
85+
xfdj <- xfdlist[[j]]
86+
print(class(xfdj))
87+
xcoef <- xfdj$coefs
88+
xbasis <- xfdj$basis
89+
betafdParj <- betalist[[j]]
90+
bbasis <- betafdParj$fd$basis
91+
bnbasis <- bbasis$nbasis
92+
pjvec[j] <- bnbasis
93+
Jpsithetaj <- inprod(xbasis,bbasis)
94+
Zmat <- cbind(Zmat,crossprod(xcoef,Jpsithetaj))
95+
print("betafdParj$estimate")
96+
if (betafdParj$estimate) {
97+
lambdaj <- betafdParj$lambda
98+
if (lambdaj > 0) {
99+
Lfdj <- betafdParj$Lfd
100+
Rmatj <- lambdaj*eval.penalty(bbasis,Lfdj)
101+
} else {
102+
Rmatj <- matrix(0,bnbasis,bnbasis)
103+
}
104+
if (ncoef > 0) {
105+
zeromat <- matrix(0,ncoef,bnbasis)
106+
Rmat <- rbind(cbind(Rmat, zeromat),
107+
cbind(t(zeromat), Rmatj))
108+
} else {
109+
Rmat <- Rmatj
110+
ncoef <- ncoef + bnbasis
111+
}
112+
}
113+
}
114+
115+
# -----------------------------------------------------------
116+
# set up the linear equations for the solution
117+
# -----------------------------------------------------------
118+
119+
print("Cmat and Dmat")
120+
121+
# solve for coefficients defining BETA
122+
123+
if (any(wt != 1)) {
124+
rtwt <- sqrt(wt)
125+
Zmatwt <- Zmat*rtwt
126+
ymatwt <- ymat*rtwt
127+
Cmat <- t(Zmatwt) %*% Zmatwt + Rmat
128+
Dmat <- t(Zmatwt) %*% ymatwt
129+
} else {
130+
Cmat <- t(Zmat) %*% Zmat + Rmat
131+
Dmat <- t(Zmat) %*% ymat
132+
}
133+
134+
# print("solving")
135+
136+
eigchk(Cmat)
137+
138+
Cmatinv <- solve(Cmat)
139+
140+
betacoef <- Cmatinv %*% Dmat
141+
142+
143+
# compute and print degrees of freedom measure
144+
145+
df <- sum(diag(Zmat %*% Cmatinv %*% t(Zmat)))
146+
147+
print("betaestlist")
148+
149+
# set up fdPar object for BETAESTFDPAR
150+
151+
betaestlist <- betalist
152+
mj2 <- 0
153+
for (j in 1:p) {
154+
betafdParj <- betalist[[j]]
155+
betafdj <- betafdParj$fd
156+
ncoefj <- betafdj$basis$nbasis
157+
mj1 <- mj2 + 1
158+
mj2 <- mj2 + ncoefj
159+
indexj <- mj1:mj2
160+
betacoefj <- betacoef[indexj]
161+
betaestfdj <- betafdj
162+
betaestfdj$coefs <- as.matrix(betacoefj)
163+
betaestfdParj <- betafdParj
164+
betaestfdParj$fd <- betaestfdj
165+
betaestlist[[j]] <- betaestfdParj
166+
}
167+
168+
# set up fd object for predicted values
169+
170+
print("yhatmat")
171+
172+
yhatmat <- matrix(0,N,1)
173+
for (j in 1:p) {
174+
xfdj <- xfdlist[[j]]
175+
if (inherits(xfdj, "fd")) {
176+
xbasis <- xfdj$basis
177+
xnbasis <- xbasis$nbasis
178+
xrng <- xbasis$rangeval
179+
nfine <- max(501,10*xnbasis+1)
180+
tfine <- seq(xrng[1], xrng[2], len=nfine)
181+
deltat <- tfine[2]-tfine[1]
182+
xmat <- eval.fd(tfine, xfdj, 0, returnMatrix)
183+
betafdParj <- betaestlist[[j]]
184+
betafdj <- betafdParj$fd
185+
betamat <- eval.fd(tfine, betafdj, 0, returnMatrix)
186+
fitj <- deltat*(crossprod(xmat,betamat) -
187+
0.5*(outer(xmat[1, ],betamat[1, ]) +
188+
outer(xmat[nfine,],betamat[nfine,])))
189+
yhatmat <- yhatmat + fitj
190+
} else{
191+
betaestfdParj <- betaestlist[[j]]
192+
betavecj <- betaestfdParj$fd$coefs
193+
yhatmat <- yhatmat + xfdj %*% t(betavecj)
194+
}
195+
}
196+
yhatfdobj <- yhatmat
197+
198+
# -----------------------------------------------------------------------
199+
# Compute pointwise standard errors of regression coefficients
200+
# if both y2cMap and SigmaE are supplied.
201+
# -----------------------------------------------------------------------
202+
203+
print("standard errors")
204+
205+
if (!(is.null(y2cMap) || is.null(SigmaE))) {
206+
207+
# check dimensions of y2cMap and SigmaE
208+
209+
y2cdim <- dim(y2cMap)
210+
if (y2cdim[2] != dim(SigmaE)[1]) stop(
211+
"Dimensions of Y2CMAP not correct.")
212+
213+
214+
# compute linear mapping c2bMap takinging coefficients for
215+
# response into coefficients for regression functions
216+
217+
c2bMap <- Cmatinv %*% t(Zmat)
218+
y2bmap <- c2bMap
219+
bvar <- y2bmap %*% as.matrix(SigmaE) %*% t(y2bmap)
220+
betastderrlist <- vector("list",p)
221+
mj2 <- 0
222+
for (j in 1:p) {
223+
betafdParj <- betalist[[j]]
224+
betabasisj <- betafdParj$fd$basis
225+
ncoefj <- betabasisj$nbasis
226+
mj1 <- mj2 + 1
227+
mj2 <- mj2 + ncoefj
228+
indexj <- mj1:mj2
229+
bvarj <- bvar[indexj,indexj]
230+
betarng <- betabasisj$rangeval
231+
nfine <- max(c(501,10*ncoefj+1))
232+
tfine <- seq(betarng[1], betarng[2], len=nfine)
233+
bbasismat <- eval.basis(tfine, betabasisj, 0, returnMatrix)
234+
bstderrj <- sqrt(diag(bbasismat %*% bvarj %*% t(bbasismat)))
235+
bstderrfdj <- smooth.basis(tfine, bstderrj, betabasisj)$fd
236+
betastderrlist[[j]] <- bstderrfdj
237+
}
238+
} else {
239+
betastderrlist = NULL
240+
bvar = NULL
241+
c2bMap = NULL
242+
}
243+
244+
# -----------------------------------------------------------------------
245+
# Set up output list object
246+
# -----------------------------------------------------------------------
247+
248+
fRegressList <-
249+
list(yfdPar = y,
250+
xfdlist = xfdlist,
251+
betalist = betalist,
252+
betaestlist = betaestlist,
253+
yhatfdobj = yhatfdobj,
254+
Cmat = Cmat,
255+
Dmat = Dmat,
256+
Cmatinv = Cmatinv,
257+
wt = wt,
258+
df = df,
259+
y2cMap = y2cMap,
260+
SigmaE = SigmaE,
261+
betastderrlist = betastderrlist,
262+
bvar = bvar,
263+
c2bMap = c2bMap)
264+
265+
class(fRegressList) <- "fRegress"
266+
267+
return(fRegressList)
268+
269+
}

0 commit comments

Comments
 (0)