|
1 |
| -fRegress.CV <- function(y, xfdlist, betalist, wt=NULL, CVobs=1:N, |
2 |
| - returnMatrix=FALSE, ...) |
3 |
| -{ |
4 |
| - |
5 |
| -# FREGRESS.CV computes cross-validated error sum of squares |
6 |
| -# for scalar or functional responses. NOTE: ordinary and |
7 |
| -# generalized cross validation scores are now returned by fRegress |
8 |
| -# when scalar responses are used. |
9 |
| - |
10 |
| -# last modified 28 July 2012 by Jim Ramsay |
11 |
| - |
12 |
| -# check the arguments |
13 |
| - |
14 |
| -argList <- fRegressArgCheck(y, xfdlist, betalist, wt) |
15 |
| -yfdPar <- argList$yfdPar |
16 |
| -xfdlist <- argList$xfdlist |
17 |
| -betalist <- argList$betalist |
18 |
| -wt <- argList$wt |
19 |
| - |
20 |
| -# extract dimensions of the data and analysis |
21 |
| - |
22 |
| -p <- length(xfdlist) |
23 |
| -N <- dim(xfdlist[[1]]$coef)[2] |
24 |
| -M <- length(CVobs) |
25 |
| - |
26 |
| -# branch to either scalar or functional dependent variable |
27 |
| - |
28 |
| -if (inherits(yfdPar, "numeric")) { |
29 |
| - |
30 |
| - # scalar dependent variable case |
31 |
| - |
32 |
| - yvec <- yfdPar |
33 |
| - SSE.CV <- 0 |
34 |
| - errfd <- c() |
35 |
| - for (m in 1:M) { |
36 |
| - i <- CVobs[m] |
37 |
| - # eliminate case i from the weights |
38 |
| - wti <- wt[-i] |
39 |
| - xfdlisti <- vector("list",p) |
40 |
| - for (j in 1:p) { |
41 |
| - xfdj <- xfdlist[[j]] |
42 |
| - if (inherits(xfdj, "numeric")) { |
43 |
| - betafdParj <- betalist[[j]] |
44 |
| - betafdj <- betafdParj$fd |
45 |
| - basisj <- betafdj$basis |
46 |
| - betarangej <- basisj$rangeval |
47 |
| - conbasisj <- create.constant.basis(betarangej) |
48 |
| - xfdj <- fd(matrix(xfdj,1,N), conbasisj) |
49 |
| - } |
50 |
| - basisj <- xfdj$basis |
51 |
| - coefj <- xfdj$coefs |
52 |
| - if (dim(coefj)[1] == 1) coefj <- matrix(coefj[-i],1,N-1) |
53 |
| - else coefj <- as.matrix(coefj[,-i]) |
54 |
| - xfdlisti[[j]] <- fd(coefj,basisj) |
55 |
| - } |
56 |
| - yveci <- yvec[-i] |
57 |
| - fRegressListi <- fRegress(yveci, xfdlisti, betalist, wti) |
58 |
| - betaestlisti <- fRegressListi$betaestlist |
59 |
| - yhati <- 0 |
60 |
| - for (j in 1:p) { |
61 |
| - betafdParj <- betaestlisti[[j]] |
62 |
| - betafdj <- betafdParj$fd |
63 |
| - xfdj <- xfdlist[[j]] |
64 |
| - bbasisj <- betafdj$basis |
65 |
| - rangej <- bbasisj$rangeval |
66 |
| - nfine <- max(501, bbasisj$nbasis*10+1) |
67 |
| - tfine <- seq(rangej[1], rangej[2], len=nfine) |
68 |
| - delta <- tfine[2]-tfine[1] |
69 |
| - betavec <- eval.fd(tfine, betafdj, 0, returnMatrix) |
70 |
| - xveci <- eval.fd(tfine, xfdj[i], 0, returnMatrix) |
71 |
| - yhati <- yhati + delta*(sum(xveci*betavec) - |
72 |
| - 0.5*( xveci[1] *betavec[1] + |
73 |
| - xveci[nfine]*betavec[nfine] )) |
74 |
| - } |
75 |
| - errfd[i] = yvec[i] - yhati; |
76 |
| - SSE.CV <- SSE.CV + errfd[i]^2 |
77 |
| - } |
78 |
| - } else { |
79 |
| - |
80 |
| - # functional dependent variable case |
81 |
| - |
82 |
| - yfd <- yfdPar$fd |
83 |
| - SSE.CV <- 0 |
84 |
| - errcoefs <- c() |
85 |
| - for(m in 1:length(CVobs)){ |
86 |
| - # index of case to eliminate |
87 |
| - i <- CVobs[m] |
88 |
| - # eliminate case i from the weights |
89 |
| - wti <- wt[-i] |
90 |
| - # eliminate case i from covariates |
91 |
| - txfdlist <- xfdlist |
92 |
| - for(k in 1:p){ |
93 |
| - txfdlist[[k]] <- xfdlist[[k]][-i] |
94 |
| - } |
95 |
| - # eliminate case i from dependent variable |
96 |
| - yfdi <- yfd[-i] |
97 |
| - # carry out the functional regression analysis |
98 |
| - tres <- fRegress(yfdi,txfdlist,betalist,wti) |
99 |
| - # extract the regression coefficient functions |
100 |
| - betaestlisti <- tres$betaestlist |
101 |
| - # compute the fit to the data for case i |
102 |
| - yhatfdi <- 0 |
103 |
| - for(k in 1:p){ |
104 |
| - betafdPark = betaestlisti[[k]] |
105 |
| - betafdk = betafdPark$fd |
106 |
| - xfdk = xfdlist[[k]] |
107 |
| - xfdik = xfdk[i] |
108 |
| - tempfd = xfdik*betafdk |
109 |
| - yhatfdi <- yhatfdi + tempfd |
110 |
| - } |
111 |
| - # compute the residual function |
112 |
| - errfdi <- yfd[i] - yhatfdi |
113 |
| - # increment the error sum of squares by the integral of the |
114 |
| - # square of the residual function |
115 |
| - SSE.CV <- SSE.CV + inprod(errfdi,errfdi) |
116 |
| - # add the coefficients for the residual function |
117 |
| - errcoefs <- cbind(errcoefs,errfdi$coefs) |
118 |
| - } |
119 |
| - # set up the functional data object for the residual fns |
120 |
| - errfd <- fd(errcoefs,errfdi$basis) |
121 |
| - names(errfd$fdnames)[[3]] <- "Xval Errors" |
122 |
| -} |
123 |
| -return(list(SSE.CV=SSE.CV,errfd.cv=errfd)) |
124 |
| -} |
125 |
| - |
126 |
| - |
| 1 | +fRegress.CV <- function(y, xfdlist, betalist, wt=NULL, CVobs=1:N, |
| 2 | + returnMatrix=FALSE, ...) |
| 3 | +{ |
| 4 | + |
| 5 | +# FREGRESS.CV computes cross-validated error sum of squares |
| 6 | +# for scalar or functional responses. NOTE: ordinary and |
| 7 | +# generalized cross validation scores are now returned by fRegress |
| 8 | +# when scalar responses are used. |
| 9 | + |
| 10 | +# last modified 16 December 2020 by Jim Ramsay |
| 11 | + |
| 12 | +# check the arguments |
| 13 | + |
| 14 | +argList <- fRegressArgCheck(y, xfdlist, betalist, wt) |
| 15 | + |
| 16 | +yfdobj <- argList$yfd |
| 17 | +xfdlist <- argList$xfdlist |
| 18 | +betalist <- argList$betalist |
| 19 | +wt <- argList$wt |
| 20 | + |
| 21 | +# extract dimensions of the data and analysis |
| 22 | + |
| 23 | +p <- length(xfdlist) |
| 24 | +N <- dim(xfdlist[[1]]$coef)[2] |
| 25 | +M <- length(CVobs) |
| 26 | + |
| 27 | +# branch to either scalar or functional dependent variable |
| 28 | + |
| 29 | +if (inherits(yfdobj, "numeric")) { |
| 30 | + |
| 31 | + # scalar dependent variable case |
| 32 | + |
| 33 | + yvec <- yfdobj |
| 34 | + SSE.CV <- 0 |
| 35 | + errfd <- c() |
| 36 | + for (m in 1:M) { |
| 37 | + i <- CVobs[m] |
| 38 | + # eliminate case i from the weights |
| 39 | + wti <- wt[-i] |
| 40 | + xfdlisti <- vector("list",p) |
| 41 | + for (j in 1:p) { |
| 42 | + xfdj <- xfdlist[[j]] |
| 43 | + if (inherits(xfdj, "numeric")) { |
| 44 | + betafdParj <- betalist[[j]] |
| 45 | + betafdj <- betafdParj$fd |
| 46 | + basisj <- betafdj$basis |
| 47 | + betarangej <- basisj$rangeval |
| 48 | + conbasisj <- create.constant.basis(betarangej) |
| 49 | + xfdj <- fd(matrix(xfdj,1,N), conbasisj) |
| 50 | + } |
| 51 | + basisj <- xfdj$basis |
| 52 | + coefj <- xfdj$coefs |
| 53 | + if (dim(coefj)[1] == 1) coefj <- matrix(coefj[-i],1,N-1) |
| 54 | + else coefj <- as.matrix(coefj[,-i]) |
| 55 | + xfdlisti[[j]] <- fd(coefj,basisj) |
| 56 | + } |
| 57 | + yveci <- yvec[-i] |
| 58 | + fRegressListi <- fRegress(yveci, xfdlisti, betalist, wti) |
| 59 | + betaestlisti <- fRegressListi$betaestlist |
| 60 | + yhati <- 0 |
| 61 | + for (j in 1:p) { |
| 62 | + betafdParj <- betaestlisti[[j]] |
| 63 | + betafdj <- betafdParj$fd |
| 64 | + xfdj <- xfdlist[[j]] |
| 65 | + bbasisj <- betafdj$basis |
| 66 | + rangej <- bbasisj$rangeval |
| 67 | + nfine <- max(501, bbasisj$nbasis*10+1) |
| 68 | + tfine <- seq(rangej[1], rangej[2], len=nfine) |
| 69 | + delta <- tfine[2]-tfine[1] |
| 70 | + betavec <- eval.fd(tfine, betafdj, 0, returnMatrix) |
| 71 | + xveci <- eval.fd(tfine, xfdj[i], 0, returnMatrix) |
| 72 | + yhati <- yhati + delta*(sum(xveci*betavec) - |
| 73 | + 0.5*( xveci[1] *betavec[1] + |
| 74 | + xveci[nfine]*betavec[nfine] )) |
| 75 | + } |
| 76 | + errfd[i] = yvec[i] - yhati; |
| 77 | + SSE.CV <- SSE.CV + errfd[i]^2 |
| 78 | + } |
| 79 | + } else { |
| 80 | + |
| 81 | + # functional dependent variable case |
| 82 | + |
| 83 | + yfd <- yfdobj |
| 84 | + SSE.CV <- 0 |
| 85 | + errcoefs <- c() |
| 86 | + for(m in 1:length(CVobs)){ |
| 87 | + # index of case to eliminate |
| 88 | + i <- CVobs[m] |
| 89 | + # eliminate case i from the weights |
| 90 | + wti <- wt[-i] |
| 91 | + # eliminate case i from covariates |
| 92 | + txfdlist <- xfdlist |
| 93 | + for(k in 1:p){ |
| 94 | + txfdlist[[k]] <- xfdlist[[k]][-i] |
| 95 | + } |
| 96 | + # eliminate case i from dependent variable |
| 97 | + yfdi <- yfd[-i] |
| 98 | + # carry out the functional regression analysis |
| 99 | + tres <- fRegress(yfdi,txfdlist,betalist,wti) |
| 100 | + # extract the regression coefficient functions |
| 101 | + betaestlisti <- tres$betaestlist |
| 102 | + # compute the fit to the data for case i |
| 103 | + yhatfdi <- 0 |
| 104 | + for(k in 1:p){ |
| 105 | + betafdPark = betaestlisti[[k]] |
| 106 | + betafdk = betafdPark$fd |
| 107 | + xfdk = xfdlist[[k]] |
| 108 | + xfdik = xfdk[i] |
| 109 | + tempfd = xfdik*betafdk |
| 110 | + yhatfdi <- yhatfdi + tempfd |
| 111 | + } |
| 112 | + # compute the residual function |
| 113 | + errfdi <- yfd[i] - yhatfdi |
| 114 | + # increment the error sum of squares by the integral of the |
| 115 | + # square of the residual function |
| 116 | + SSE.CV <- SSE.CV + inprod(errfdi,errfdi) |
| 117 | + # add the coefficients for the residual function |
| 118 | + errcoefs <- cbind(errcoefs,errfdi$coefs) |
| 119 | + } |
| 120 | + # set up the functional data object for the residual fns |
| 121 | + errfd <- fd(errcoefs,errfdi$basis) |
| 122 | + names(errfd$fdnames)[[3]] <- "Xval Errors" |
| 123 | +} |
| 124 | +return(list(SSE.CV=SSE.CV, errfd.cv=errfd)) |
| 125 | +} |
| 126 | + |
| 127 | + |
0 commit comments