R/fregre.igls.r
In moviedo5/fda.usc: Functional Data Analysis and Utilities for Statistical Computing
#' @title Fit of Functional Generalized Least Squares Model Iteratively
#' @aliases fregre.igls
#'
#' @param formula A two-sided linear formula object describing the
#' model, with the response on the left of a \code{~} operator and the
#' terms, separated by \code{+} operators, on the right.
#' @param data An optional data frame containing the variables named in
#' \code{model}, \code{correlation}, \code{weights}, and
#' \code{subset}. By default the variables are taken from the environment
#' from which \code{gls} is called.
#' @param basis.x List of basis for functional explanatory data estimation.
#' @param basis.b List of basis for \eqn{\beta(t)} parameter estimation.
#' @param rn List of Ridge parameter.
#' @param lambda List of Roughness penalty parameter.
#' @param weights weights
#' @param correlation List describing the correlation structure. Defaults to
#' \code{NULL}, corresponding to uncorrelated errors.
#' See the following internal functions for a description and a code example in script file.
#' \itemize{
#' \item \code{corUnstruc(x)}, fit an unstrutured correlation.
#' \item \code{cor.AR(x, order.max = 8, p=1, method = "lm")} fit an Autoregressive Models to Time Series using \code{\link{ar}} function.
#' \item \code{cor.ARMA(x, p, d = 0, q = 0, method = "lm", order.max = 1)} Fit an ARIMA model to a univariate time series using \code{\link{arima}} function.
#' \item \code{corExpo(xy,range, method = "euclidean",p=2)} Fit an exponential correlation structure.
#' }
#' @param maxit Number of maximum of interactions.
# @param weights An optional \code{\link{varFunc}} object or one-sided formula
# describing the within-group heteroscedasticity structure. If given as
# a formula, it is used as the argument to \code{\link{varFixed}},
# corresponding to fixed variance weights. See the documentation on
# \code{\link{varClasses}} for a description of the available \code{\link{varFunc}}
# classes. Defaults to \code{NULL}, corresponding to homoscedastic errors.
#' @param control Control parameters.
#' @param \dots Further arguments passed to or from other methods.
#' @description This function fits iteratively a functional linear model using generalized
#' least squares. The errors are allowed to be correlated and/or have unequal variances.
#' \enumerate{
#' \item Begin with a preliminary estimation of \eqn{\hat{\theta}=\theta_0} (for instance,
#' \eqn{\theta_0=0}). Compute \eqn{\hat{W}}.
#' \item Estimate \eqn{b_\Sigma =(Z'\hat{W}Z)^{-1}Z'\hat{W}y}
#' \item Based on the residuals, \eqn{\hat{e}=\left(y-Zb_\Sigma \right)}, update
#' \eqn{\hat{\theta}=\rho\left({\hat{e}}\right)} where \eqn{\rho} depends on the
#' dependence structure chosen.
#' \item Repeats steps 2 and 3 until convergence (small changes in \eqn{b_\Sigma} and/or \eqn{\hat{\theta}}).
#' }
#' @return An object of class \code{fregre.igls} representing the functional linear model
#' fit with temporal dependence errors.
#' Beside, the class(z) is similar to "fregre.lm" plus the following objects:
#' \itemize{
#' \item \code{corStruct}: Fitted AR or ARIMA model.
# \item{formula.ini}{ formula in call.}
# \item{fdataob}{ }
# \item{rn}{ rn}
# \item{vs.list}{ }
# \item{correlation}{ }
#' }
#'
#' @references Oviedo de la Fuente, M., Febrero-Bande, M., Pilar Munoz, and
#' Dominguez, A. (2018). Predicting seasonal influenza transmission using
#' functional regression models with temporal dependence. PloS one, 13(4), e0194250.
#' \doi{10.1371/journal.pone.0194250}
#' @examples
#' \dontrun{
#' data(tecator)
#' x=tecator$absorp.fdata
#' x.d2<-fdata.deriv(x,nderiv=)
#' tt<-x[["argvals"]]
#' dataf=as.data.frame(tecator$y)
#' # plot the response
#' plot(ts(tecator$y$Fat))
#' ldata=list("df"=dataf,"x.d2"=x.d2)
#' res.gls=fregre.igls(Fat~x.d2,data=ldata,
#' correlation=list("cor.ARMA"=list()),
#' control=list("p"=1))
#' res.gls
#' res.gls$corStruct
#' }
#' @keywords models regression
#' @export
fregre.igls<-function (formula, data, basis.x = NULL, basis.b = NULL, correlation,
maxit = 100, rn, lambda, weights = rep(1, n), control, ...)
{
#print("fregre.igls")
tf <- terms.formula(formula)
terms <- attr(tf, "term.labels")
nt <- length(terms)
if (attr(tf, "response") > 0) {
response <- as.character(attr(tf, "variables")[2])
pf <- rf <- paste(response, "~", sep = "")
}
else pf <- rf <- "~"
vtab <- rownames(attr(tf, "factors"))
vnf = intersect(terms, names(data$df))
vnf2 = intersect(vtab[-1], names(data$df)[-1])
vfunc2 = setdiff(terms, vnf)
vint = setdiff(terms, vtab)
vfunc = setdiff(vfunc2, vint)
off <- attr(tf, "offset")
name.coef = nam = par.fregre = beta.l = list()
kterms = 1
n <- length(data[["df"]][, response])
XX = data.frame(data[["df"]][, c(response)], weights)
namxx = names(XX) = c(response, "weights")
aa <- NULL
if (length(vnf) > 0) {
for (i in 1:length(vnf)) {
if (kterms > 1)
pf <- paste(pf, "+", vnf[i], sep = "")
else pf <- paste(pf, vnf[i], sep = "")
kterms <- kterms + 1
}
if (attr(tf, "intercept") == 0) {
pf <- paste(pf, -1, sep = "")
}
mf <- as.data.frame(model.matrix(formula(pf), data$df))
vnf2 <- names(mf)[-1]
pf <- rf <- paste(response, "~", sep = "")
for (i in 1:length(vnf2)) pf <- paste(pf, "+", vnf2[i],
sep = "")
cname <- names(XX)
XX <- data.frame(XX, mf)
names(XX) <- c(cname, names(mf))
}
else {
XX <- data.frame(XX, model.matrix(formula(paste(pf, "1")),
data$df))
names(XX)[3] <- "(Intercept)"
}
if (missing(rn)) {
rn0 = FALSE
rn = list()
}
else rn0 <- TRUE
if (missing(lambda)) {
lambda0 = FALSE
lambda = list()
}
else lambda0 <- TRUE
mat <- rep(0, len = ncol(XX) - 2)
imat2 <- ncol(XX) - 2
mat2 <- diag(0, nrow = imat2)
if (length(vfunc) > 0) {
mean.list = vs.list = JJ = list()
bsp1 <- bsp2 <- TRUE
for (i in 1:length(vfunc)) {
if (is(data[[vfunc[i]]],"fdata")) {
tt <- data[[vfunc[i]]][["argvals"]]
rtt <- data[[vfunc[i]]][["rangeval"]]
np <- length(tt)
fdat <- data[[vfunc[i]]]
dat <- data[[vfunc[i]]]$data
if (is.null(basis.x[[vfunc[i]]]))
basis.x[[vfunc[i]]] <- create.fdata.basis(fdat,
l = 1:7)
else if (basis.x[[vfunc[i]]]$type == "pc" | basis.x[[vfunc[i]]]$type ==
"pls")
bsp1 = FALSE
if (is.null(basis.b[[vfunc[i]]]) & bsp1)
basis.b[[vfunc[i]]] <- create.fdata.basis(fdat)
else if (basis.x[[vfunc[i]]]$type == "pc" | basis.x[[vfunc[i]]]$type ==
"pls")
bsp2 = FALSE
if (bsp1 & bsp2) {
if (is.null(rownames(dat)))
rownames(fdat$data) <- 1:nrow(dat)
fdnames = list(time = tt, reps = rownames(fdat[["data"]]),
values = "values")
xcc <- fdata.cen(data[[vfunc[i]]])
mean.list[[vfunc[i]]] = xcc[[2]]
if (!is.null(basis.x[[vfunc[i]]]$dropind)) {
int <- setdiff(1:basis.x[[vfunc[i]]]$nbasis,
basis.x[[vfunc[i]]]$dropind)
basis.x[[vfunc[i]]]$nbasis <- length(int)
basis.x[[vfunc[i]]]$dropind <- NULL
basis.x[[vfunc[i]]]$names <- basis.x[[vfunc[i]]]$names[int]
}
if (!is.null(basis.b[[vfunc[i]]]$dropind)) {
int <- setdiff(1:basis.b[[vfunc[i]]]$nbasis,
basis.b[[vfunc[i]]]$dropind)
basis.b[[vfunc[i]]]$nbasis <- length(int)
basis.b[[vfunc[i]]]$dropind <- NULL
basis.b[[vfunc[i]]]$names <- basis.b[[vfunc[i]]]$names[int]
}
x.fd = Data2fd(argvals = tt, y = t(xcc[[1]]$data),
basisobj = basis.x[[vfunc[i]]], fdnames = fdnames)
r = x.fd[[2]][[3]]
J = inprod(basis.x[[vfunc[i]]], basis.b[[vfunc[i]]])
Z = t(x.fd$coefs) %*% J
colnames(J) = colnames(Z) = name.coef[[vfunc[i]]] = paste(vfunc[i],
".", basis.b[[vfunc[i]]]$names, sep = "")
XX = cbind(XX, Z)
for (j in 1:length(colnames(Z))) {
if (kterms >= 1)
pf <- paste(pf, "+", colnames(Z)[j], sep = "")
else pf <- paste(pf, colnames(Z)[j], sep = "")
kterms <- kterms + 1
}
JJ[[vfunc[i]]] <- J
nbasisb <- basis.b[[vfunc[i]]]$nbasis
R = diag(0, ncol = nbasisb, nrow = nbasisb)
R <- eval.penalty(basis.b[[vfunc[i]]], Lfdobj = int2Lfd(2))
lenm <- ncol(mat2)
if (rn0)
stop("Ridge regressions is only implemented for functional principal component basis")
MM <- matrix(0, nrow = lenm, ncol = nbasisb)
MM2 <- matrix(0, nrow = nbasisb, ncol = imat2 +
nbasisb)
mat2 <- cbind(mat2, MM)
mat2 <- rbind(mat2, MM2)
if (!is.null(lambda[[vfunc[i]]])) {
mat2[(imat2 + 1):(imat2 + nbasisb), (imat2 +
1):(imat2 + nbasisb)] <- lambda[[vfunc[i]]] *
R
}
imat2 <- imat2 + nbasisb
}
else {
basis <- basis.x[[vfunc[i]]]
l <- basis$l
vs <- t(basis$basis$data)
basis$coefs <- basis$coefs[, l, drop = FALSE]
Z <- basis$coefs
response = "y"
colnames(Z) = name.coef[[vfunc[i]]] = paste(vfunc[i],
".", rownames(basis$basis$data), sep = "")
XX = cbind(XX, Z)
vs.list[[vfunc[i]]] = basis$basis
mean.list[[vfunc[i]]] = basis$mean
for (j in 1:length(colnames(Z))) {
if (kterms >= 1)
pf <- paste(pf, "+", name.coef[[vfunc[i]]][j],
sep = "")
else pf <- paste(pf, name.coef[[vfunc[i]]][j],
sep = "")
kterms <- kterms + 1
}
if (!is.null(rn[[vfunc[i]]])) {
mat <- c(mat, rep(rn[[vfunc[i]]], len = length(l)))
}
else mat <- c(mat, rep(0, len = length(l)))
lenl <- length(l)
lenm <- ncol(mat2)
MM <- matrix(0, nrow = lenm, ncol = lenl)
MM2 <- matrix(0, nrow = lenl, ncol = imat2 + lenl)
mat2 <- cbind(mat2, MM)
mat2 <- rbind(mat2, MM2)
if (!is.null(lambda[[vfunc[i]]])) {
R <- P.penalty(1:length(l))
mat2[(imat2 + 1):(imat2 + lenl), (imat2 +
1):(imat2 + lenl)] <- lambda[[vfunc[i]]] *
R
}
imat2 <- imat2 + lenl
}
}
else {
if (is(data[[vfunc[i]]],"fd")) {
fdat <- data[[vfunc[i]]]
if (is.null(basis.x[[vfunc[i]]]))
basis.x[[vfunc[i]]] <- fdat$basis
else if (inherits(basis.x[[vfunc[i]]], "pca.fd"))
bsp1 = FALSE
if (is.null(basis.b[[vfunc[i]]]) & bsp1)
basis.b[[vfunc[i]]] <- create.fdata.basis(fdat,
l = 1:max(5, floor(basis.x[[vfunc[i]]]$nbasis/5)),
type.basis = basis.x[[vfunc[i]]]$type,
rangeval = fdat$basis$rangeval)
else if (inherits(basis.x[[vfunc[i]]], "pca.fd"))
bsp2 = FALSE
if (bsp1 & bsp2) {
r = fdat[[2]][[3]]
if (!is.null(basis.x[[vfunc[i]]]$dropind)) {
int <- setdiff(1:basis.x[[vfunc[i]]]$nbasis,
basis.x[[vfunc[i]]]$dropind)
basis.x[[vfunc[i]]]$nbasis <- length(int)
basis.x[[vfunc[i]]]$dropind <- NULL
basis.x[[vfunc[i]]]$names <- basis.x[[vfunc[i]]]$names[int]
}
if (!is.null(basis.b[[vfunc[i]]]$dropind)) {
int <- setdiff(1:basis.b[[vfunc[i]]]$nbasis,
basis.b[[vfunc[i]]]$dropind)
basis.b[[vfunc[i]]]$nbasis <- length(int)
basis.b[[vfunc[i]]]$dropind <- NULL
basis.b[[vfunc[i]]]$names <- basis.b[[vfunc[i]]]$names[int]
}
J = inprod(basis.x[[vfunc[i]]], basis.b[[vfunc[i]]])
mean.list[[vfunc[i]]] <- mean.fd(x.fd)
x.fd <- center.fd(x.fd)
Z = t(x.fd$coefs) %*% J
colnames(J) = colnames(Z) = name.coef[[vfunc[i]]] = paste(vfunc[i],
".", basis.b[[vfunc[i]]]$names, sep = "")
XX = cbind(XX, Z)
for (j in 1:length(colnames(Z))) {
if (kterms >= 1)
pf <- paste(pf, "+", colnames(Z)[j],
sep = "")
else pf <- paste(pf, colnames(Z)[j], sep = "")
kterms <- kterms + 1
}
JJ[[vfunc[i]]] <- J
}
else {
l <- ncol(basis.x[[vfunc[i]]]$scores)
vs <- basis.x[[vfunc[i]]]$harmonics$coefs
Z <- basis.x[[vfunc[i]]]$scores
response = "y"
colnames(Z) = name.coef[[vfunc[i]]] = paste(vfunc[i],
".", colnames(basis.x[[vfunc[i]]]$harmonics$coefs),
sep = "")
XX = cbind(XX, Z)
vs.list[[vfunc[i]]] = vs
mean.list[[vfunc[i]]] = basis.x[[vfunc[i]]]$meanfd
for (j in 1:length(colnames(Z))) {
if (kterms >= 1)
pf <- paste(pf, "+", name.coef[[vfunc[i]]][j],
sep = "")
else pf <- paste(pf, name.coef[[vfunc[i]]][j],
sep = "")
kterms <- kterms + 1
}
}
}
else stop("Please, enter functional covariate")
}
}
}
if (!is.data.frame(XX))
XX = data.frame(XX)
par.fregre$formula = pf
par.fregre$data = XX
y <- XX[, 1]
scores <- as.matrix(XX[, -(1:2)])
W <- diag(weights)
error <- FALSE
if (missing(correlation))
correlation0 <- FALSE
else {
if (is.matrix(correlation))
error = FALSE
else error = TRUE
correlation0 <- TRUE
}
if (!rn0 & !lambda0 & !correlation0) {
# print(1)
W0 <- FALSE
z = lm(formula = pf, data = XX, ...)
e <- z$residuals
S <- solve(t(scores) %*% W %*% scores)
class(z) <- c(class(z), "fregre.lm")
mat0 <- diag(0, ncol(scores))
}
else {
# print(11)
mat0 <- diag(0, ncol(scores))
if (lambda0) {
mat0 <- mat2
}
if (rn0) {
mat0 <- diag(mat)
}
if (rn0 & lambda0)
warning("Only ridge penalization is done by rn argument (lambda argument is ignored)")
W <- diag(weights) %*% W
S <- solve(t(scores) %*% W %*% scores + mat0)
Cinv <- S %*% t(scores) %*% W
ycen = y - mean(y)
coefs <- Cinv %*% XX[, 1]
z <- list()
z$fitted.values <- yp <- drop(scores %*% coefs)
e <- z$residuals <- XX[, 1] - z$fitted.values
it <- 1
eps = 0.001
eps2 = 1e-10
err2 = sqrt(sum(coefs^2))
MM <- matrix(1:n, ncol = 1)
corStruct <- list()
# print(error)
while (error) {
W <- W0 <- diag(n)
if (is.list(correlation)) {
ncor <- length(correlation)
name.cor <- names(correlation)
}
else {
ncor <- 1
name.cor <- correlation[[1]]
correlation = list(correlation = list())
names(correlation) <- name.cor
}
wei0 <- wei <- diag(weights)
ee <- e
for (icor in 1:ncor) {
# print(icor); print(correlation)
if (!is.list(correlation[[icor]])) {
name.cor <- correlation[[icor]]
correlation[[icor]] = list()
names(correlation)[icor] <- name.cor
}
if (is.null(correlation[[icor]][["group"]])) {
n <- length(e)
# print(name.cor)
# print(name.cor[icor])
par.cor <- switch(name.cor[icor], cor.ARMA = {
cor.name <- "cor.ARMA"
#print("cor.ARMA")
if (is.null(correlation[[icor]][["index"]])) index <- 1:n else {
index.df <- correlation[[icor]][["index"]]
index <- data[["df"]][, index.df]
if (is.vector(index)) e <- e[order(index)] else {
par.ord <- list()
for (i in 1:ncol(index)) par.ord[[i]] <- index[,
i]
index <- as.numeric(do.call("order",
par.ord))
e <- e[index]
}
}
list(x = e[index])
}, cor.AR = {
cor.name <- "cor.AR"
# print("cor.AR")
if (is.null(correlation[[icor]][["index"]])) index <- 1:n else {
index.df <- correlation[[icor]][["index"]]
index <- data[["df"]][, index.df]
if (is.vector(index)) e <- e[order(index)]
else {
par.ord <- list()
for (i in 1:ncol(index)) par.ord[[i]] <- index[,i]
index <- as.numeric(do.call("order", par.ord))
e <- e[index]
}
}
list(x = e[index])
}, corExpo = {
index.df <- correlation[[icor]][["index"]]
dxy <- data[["df"]][, index.df]
index <- 1:length(e)
cor.name <- "corExpo"
list(xy = dxy)
}, corVgm = {
index.df <- correlation[[icor]][["index"]]
xy <- data[["df"]][, index.df]
dxy0 <- as.matrix(dist(xy, diag = TRUE, upper = TRUE))
index <- 1:n
cor.name <- "corVgm"
residual <- e[index]
list(dxy = dxy0, xy = xy, df = data.frame(residual))
}, cor.Exp = {
index.df <- correlation[[icor]][["index"]]
xy <- data[["df"]][, index.df]
index <- 1:n
cor.name <- "cor.Exp"
residual <- e[index]
dff <- data.frame(data[["df"]][, index.df],
residual)
names(dff) <- c(index.df, "residual")
list(df = dff)
}, corUnstruc = {
cor.name <- "corUnstruc"
list(correlation[[icor]][["index"]])
})
e <- e[index]
# print(names(par.cor))
# print(names(correlation))
# print("1")
if (!missing(control)) {
npar <- length(par.cor)
#print(names(par.cor))
#print(names(par.ord))
# print("2")
for (i in 1:length(control)) {
par.cor[npar + i] <- control[i]
names(par.cor)[npar + i] <- names(control)[i]
}
}
# print("antes do.call")
# print(cor.name)
# print(names(par.cor))
aa <- do.call(cor.name, par.cor)
W0 <- aa[[1]]
corStruct <- aa[-1]
}
else {
#print("hay grupo")
n <- nrow(data$df)
W00 <- matrix(0, ncol = n, nrow = n)
name.group <- correlation[[icor]]$group
ncomb <- length(name.group)
# print(name.group)
if (ncomb > 1) {
gr <- (data[["df"]][, name.group[1]])
for (ngr in 2:ncomb) gr <- paste(gr, (data[["df"]][,
name.group[ngr]]), sep = "")
gr <- factor(gr, levels = unique(gr))
}
#else gr <- data[["df"]][[name.group]]
else gr <- data[["df"]][,name.group]
if (!is.factor(gr))
gr <- factor(gr)
lev <- levels(gr)
nlev <- length(lev)
#print(name.cor[icor])
par.cor <- switch(name.cor[icor],
cor.ARMA = {
#print("entra cor.ARMAAAAAA")
index.df <- correlation[[icor]][["index"]]
#print(correlation)
#print(index.df)
index <- data[["df"]][, index.df]
e2 <- NULL
for (i in 1:nlev) {
jj <- gr == lev[i]
e0 <- e[jj]
e2 <- cbind(e2, e0[order(index[jj])])
}
cor.name <- "cor.ARMA"
#print(correlation)
#print("okk")
if (is.null(correlation[[icor]][["p"]])) p0<-1
else p0<-correlation[[icor]][["p"]]
#print(p0)
list(x = e2,p=p0,order.max=p0)
},
corCloud = {
index.df <- correlation[[icor]][["index"]]
index <- data[["df"]][, index.df]
gr <- data[["df"]][, correlation[[icor]][["group"]]]
xy <- data[["df"]][, index.df]
cor.name <- "corCloud"
residual <- e
dff <- data.frame(data[["df"]][, index.df],
residual)
names(dff) <- c(index.df, "residual")
list(df = dff, gr = gr)
},
cor.Exp = {
index.df <- correlation[[icor]][["index"]]
index <- data[["df"]][, index.df]
gr <- data[["df"]][, correlation[[icor]][["group"]]]
xy <- data[["df"]][, index.df]
cor.name <- "cor.Exp"
residual <- e
dff <- data.frame(data[["df"]][, index.df],
residual)
names(dff) <- c(index.df, "residual")
list(df = dff, gr = gr)
},
corUnstruc = {
cor.name <- "corUnstruc"
index <- 1:nrow(dff)
par.cor <- list(correlation[[icor]][["index"]])
})
if (!missing(control)) {
npar <- length(par.cor)
for (k in 1:length(control)) {
par.cor[npar + k] <- control[k]
names(par.cor)[npar + k] <- names(control)[k]
}
}
#print("calll")
#print(cor.name)
#print(par.cor)
#print("antes call")
aa <- do.call(cor.name, par.cor)
W0 <- aa[[1]]
corStruct <- aa[-1]
for (j in 1:nlev) {
n <- table(gr)[j]
ilev <- gr == lev[j]
e <- ee[ilev]
if (name.cor[icor] == "cor.Exp" | name.cor[icor] ==
"corCloud")
W00[ilev, ilev] <- W0
else W00[ilev, ilev] <- W0[order(index[ilev]),
order(index[ilev])]
}
W0 <- W00
}
wei0 <- wei0 %*% W0
}
W0 <- wei0
W <- try(solve(W0), silent = TRUE)
# print(class(W) )
if (is(W,"try-error")) {
sv <- svd(W0)
W <- drop((sv$v %*% diag(1/sv$d) %*% t(sv$u)))
}
mat0[1, 1] <- 0
W2 <- t(scores) %*% W %*% scores + mat0
S <- try(solve(W2), silent = TRUE)
if (is(S,"try-error")) {
sv <- svd(W2)
S <- drop((sv$v %*% diag(1/sv$d) %*% t(sv$u)))
}
Cinv <- S %*% t(scores) %*% W
coefs <- Cinv %*% matrix(y, ncol = 1)
yp <- drop(scores %*% coefs)
coefs <- drop(coefs)
e <- y - yp
err3 <- coefs
# err4 <- max(abs((err3 - err2)/err2))
# riq <- diff(summary(e)[c(2,5)])
riq <- diff(quantile(e,probs=c(.25,.75)))
err4 <- max(abs(err3 - err2)*riq)/max(abs(err2)+eps2)
if (err4 < eps)
error <- FALSE
else {
if (it == maxit)
error <- FALSE
else it <- it + 1
err2 <- err3
}
}
H <- scores %*% Cinv
df <- traza(H)
coefs <- drop(coefs)
z$coefficients <- coefs
z$mean.list <- mean.list
z$df.residual <- n - df
z$H <- H
z$r2 <- 1 - sum(z$residuals^2)/sum(ycen^2)
if (is(basis.x[[vfunc[1]]],"basisfd")) {
z$call[[1]] = "fregre.basis"
}
else {
if (basis.x[[vfunc[1]]]$type == "pc")
z$call[[1]] = "fregre.pc"
if (basis.x[[vfunc[1]]]$type == "pls")
z$call[[1]] = "fregre.pls"
}
class(z) <- c("fregre.fd", "fregre.lm")
rdf <- n - df
sr2 <- sum(e^2)/rdf
r2 <- 1 - sum(e^2)/sum(ycen^2)
r2.adj <- 1 - (1 - r2) * ((n - 1)/rdf)
GCV <- sum(e^2)/(n - df)^2
z$terms <- terms
z$residuals <- drop(e)
z$fitted.values <- yp
z$y <- y
z$rank <- df
z$df.residual <- rdf
Z = cbind(rep(1, len = n), Z)
colnames(Z)[1] = "(Intercept)"
std.error = sqrt(diag(S) * sr2)
t.value = coefs/std.error
p.value = 2 * pt(abs(t.value), n - df, lower.tail = FALSE)
coefficients <- cbind(coefs, std.error, t.value, p.value)
colnames(coefficients) <- c("Estimate", "Std. Error",
"t value", "Pr(>|t|)")
z$coefs <- coefficients
class(z) <- "lm"
z$it <- it
}
for (i in 1:length(vfunc)) {
if (bsp1)
beta.l[[vfunc[i]]] = fd(z$coefficients[name.coef[[vfunc[i]]]],
basis.b[[vfunc[i]]])
else {
if (inherits(data[[vfunc[i]]], "fdata")) {
beta.est <- z$coefficients[name.coef[[vfunc[i]]]] *
vs.list[[vfunc[i]]]
beta.est$data <- colSums(beta.est$data)
beta.est$names$main <- "beta.est"
beta.est$data <- matrix(as.numeric(beta.est$data),
nrow = 1)
beta.est$names$main <- "beta.est"
beta.est$data <- matrix(as.numeric(beta.est$data),
nrow = 1)
if (basis.x[[vfunc[i]]]$type == "pls") {
if (basis.x[[vfunc[i]]]$norm) {
sd.X <- sqrt(apply(data[[vfunc[i]]]$data,
2, var))
beta.est$data <- beta.est$data/sd.X
}
}
beta.l[[vfunc[i]]] <- beta.est
}
else {
beta.est <- z$coefficients[name.coef[[vfunc[i]]]] *
t(vs.list[[vfunc[i]]])
beta.est <- colSums(beta.est)
beta.l[[vfunc[i]]] <- fd(beta.est, basis.x[[vfunc[i]]]$harmonics$basis)
}
}
}
z$sr2 <- sum(e^2)/z$df.residual
z$Vp = z$sr2 * S
z$beta.l = beta.l
z$formula = pf
z$mean = mean.list
z$formula.ini = formula
z$basis.x = basis.x
z$basis.b = basis.b
z$JJ <- JJ
z$data = z$data
z$XX = XX
z$data <- data
z$fdataobj <- data[[vfunc[1]]]
if (correlation0) {
rn0 <- TRUE
z$corStruct <- corStruct
z$it <- it
}
z$rn <- rn0
z$lambda <- lambda0
z$W <- W
z$W0 <- W0
z$correlation <- correlation
z$correl <- correlation0
z$vs.list = vs.list
class(z) <- c( "fregre.igls","lm")
z
}
moviedo5/fda.usc documentation built on Nov. 6, 2024, 1:21 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' @title Fit of Functional Generalized Least Squares Model Iteratively
#' @aliases fregre.igls
#'
#' @param formula A two-sided linear formula object describing the
#' model, with the response on the left of a \code{~} operator and the
#' terms, separated by \code{+} operators, on the right.
#' @param data An optional data frame containing the variables named in
#' \code{model}, \code{correlation}, \code{weights}, and
#' \code{subset}. By default the variables are taken from the environment
#' from which \code{gls} is called.
#' @param basis.x List of basis for functional explanatory data estimation.
#' @param basis.b List of basis for \eqn{\beta(t)} parameter estimation.
#' @param rn List of Ridge parameter.
#' @param lambda List of Roughness penalty parameter.
#' @param weights weights
#' @param correlation List describing the correlation structure. Defaults to
#' \code{NULL}, corresponding to uncorrelated errors.
#' See the following internal functions for a description and a code example in script file.
#' \itemize{
#' \item \code{corUnstruc(x)}, fit an unstrutured correlation.
#' \item \code{cor.AR(x, order.max = 8, p=1, method = "lm")} fit an Autoregressive Models to Time Series using \code{\link{ar}} function.
#' \item \code{cor.ARMA(x, p, d = 0, q = 0, method = "lm", order.max = 1)} Fit an ARIMA model to a univariate time series using \code{\link{arima}} function.
#' \item \code{corExpo(xy,range, method = "euclidean",p=2)} Fit an exponential correlation structure.
#' }
#' @param maxit Number of maximum of interactions.
# @param weights An optional \code{\link{varFunc}} object or one-sided formula
# describing the within-group heteroscedasticity structure. If given as
# a formula, it is used as the argument to \code{\link{varFixed}},
# corresponding to fixed variance weights. See the documentation on
# \code{\link{varClasses}} for a description of the available \code{\link{varFunc}}
# classes. Defaults to \code{NULL}, corresponding to homoscedastic errors.
#' @param control Control parameters.
#' @param \dots Further arguments passed to or from other methods.
#' @description This function fits iteratively a functional linear model using generalized
#' least squares. The errors are allowed to be correlated and/or have unequal variances.
#' \enumerate{
#' \item Begin with a preliminary estimation of \eqn{\hat{\theta}=\theta_0} (for instance,
#' \eqn{\theta_0=0}). Compute \eqn{\hat{W}}.
#' \item Estimate \eqn{b_\Sigma =(Z'\hat{W}Z)^{-1}Z'\hat{W}y}
#' \item Based on the residuals, \eqn{\hat{e}=\left(y-Zb_\Sigma \right)}, update
#' \eqn{\hat{\theta}=\rho\left({\hat{e}}\right)} where \eqn{\rho} depends on the
#' dependence structure chosen.
#' \item Repeats steps 2 and 3 until convergence (small changes in \eqn{b_\Sigma} and/or \eqn{\hat{\theta}}).
#' }
#' @return An object of class \code{fregre.igls} representing the functional linear model
#' fit with temporal dependence errors.
#' Beside, the class(z) is similar to "fregre.lm" plus the following objects:
#' \itemize{
#' \item \code{corStruct}: Fitted AR or ARIMA model.
# \item{formula.ini}{ formula in call.}
# \item{fdataob}{ }
# \item{rn}{ rn}
# \item{vs.list}{ }
# \item{correlation}{ }
#' }
#'
#' @references Oviedo de la Fuente, M., Febrero-Bande, M., Pilar Munoz, and
#' Dominguez, A. (2018). Predicting seasonal influenza transmission using
#' functional regression models with temporal dependence. PloS one, 13(4), e0194250.
#' \doi{10.1371/journal.pone.0194250}
#' @examples
#' \dontrun{
#' data(tecator)
#' x=tecator$absorp.fdata
#' x.d2<-fdata.deriv(x,nderiv=)
#' tt<-x[["argvals"]]
#' dataf=as.data.frame(tecator$y)
#' # plot the response
#' plot(ts(tecator$y$Fat))
#' ldata=list("df"=dataf,"x.d2"=x.d2)
#' res.gls=fregre.igls(Fat~x.d2,data=ldata,
#' correlation=list("cor.ARMA"=list()),
#' control=list("p"=1))
#' res.gls
#' res.gls$corStruct
#' }
#' @keywords models regression
#' @export
fregre.igls<-function (formula, data, basis.x = NULL, basis.b = NULL, correlation,
maxit = 100, rn, lambda, weights = rep(1, n), control, ...)
{
#print("fregre.igls")
tf <- terms.formula(formula)
terms <- attr(tf, "term.labels")
nt <- length(terms)
if (attr(tf, "response") > 0) {
response <- as.character(attr(tf, "variables")[2])
pf <- rf <- paste(response, "~", sep = "")
}
else pf <- rf <- "~"
vtab <- rownames(attr(tf, "factors"))
vnf = intersect(terms, names(data$df))
vnf2 = intersect(vtab[-1], names(data$df)[-1])
vfunc2 = setdiff(terms, vnf)
vint = setdiff(terms, vtab)
vfunc = setdiff(vfunc2, vint)
off <- attr(tf, "offset")
name.coef = nam = par.fregre = beta.l = list()
kterms = 1
n <- length(data[["df"]][, response])
XX = data.frame(data[["df"]][, c(response)], weights)
namxx = names(XX) = c(response, "weights")
aa <- NULL
if (length(vnf) > 0) {
for (i in 1:length(vnf)) {
if (kterms > 1)
pf <- paste(pf, "+", vnf[i], sep = "")
else pf <- paste(pf, vnf[i], sep = "")
kterms <- kterms + 1
}
if (attr(tf, "intercept") == 0) {
pf <- paste(pf, -1, sep = "")
}
mf <- as.data.frame(model.matrix(formula(pf), data$df))
vnf2 <- names(mf)[-1]
pf <- rf <- paste(response, "~", sep = "")
for (i in 1:length(vnf2)) pf <- paste(pf, "+", vnf2[i],
sep = "")
cname <- names(XX)
XX <- data.frame(XX, mf)
names(XX) <- c(cname, names(mf))
}
else {
XX <- data.frame(XX, model.matrix(formula(paste(pf, "1")),
data$df))
names(XX)[3] <- "(Intercept)"
}
if (missing(rn)) {
rn0 = FALSE
rn = list()
}
else rn0 <- TRUE
if (missing(lambda)) {
lambda0 = FALSE
lambda = list()
}
else lambda0 <- TRUE
mat <- rep(0, len = ncol(XX) - 2)
imat2 <- ncol(XX) - 2
mat2 <- diag(0, nrow = imat2)
if (length(vfunc) > 0) {
mean.list = vs.list = JJ = list()
bsp1 <- bsp2 <- TRUE
for (i in 1:length(vfunc)) {
if (is(data[[vfunc[i]]],"fdata")) {
tt <- data[[vfunc[i]]][["argvals"]]
rtt <- data[[vfunc[i]]][["rangeval"]]
np <- length(tt)
fdat <- data[[vfunc[i]]]
dat <- data[[vfunc[i]]]$data
if (is.null(basis.x[[vfunc[i]]]))
basis.x[[vfunc[i]]] <- create.fdata.basis(fdat,
l = 1:7)
else if (basis.x[[vfunc[i]]]$type == "pc" | basis.x[[vfunc[i]]]$type ==
"pls")
bsp1 = FALSE
if (is.null(basis.b[[vfunc[i]]]) & bsp1)
basis.b[[vfunc[i]]] <- create.fdata.basis(fdat)
else if (basis.x[[vfunc[i]]]$type == "pc" | basis.x[[vfunc[i]]]$type ==
"pls")
bsp2 = FALSE
if (bsp1 & bsp2) {
if (is.null(rownames(dat)))
rownames(fdat$data) <- 1:nrow(dat)
fdnames = list(time = tt, reps = rownames(fdat[["data"]]),
values = "values")
xcc <- fdata.cen(data[[vfunc[i]]])
mean.list[[vfunc[i]]] = xcc[[2]]
if (!is.null(basis.x[[vfunc[i]]]$dropind)) {
int <- setdiff(1:basis.x[[vfunc[i]]]$nbasis,
basis.x[[vfunc[i]]]$dropind)
basis.x[[vfunc[i]]]$nbasis <- length(int)
basis.x[[vfunc[i]]]$dropind <- NULL
basis.x[[vfunc[i]]]$names <- basis.x[[vfunc[i]]]$names[int]
}
if (!is.null(basis.b[[vfunc[i]]]$dropind)) {
int <- setdiff(1:basis.b[[vfunc[i]]]$nbasis,
basis.b[[vfunc[i]]]$dropind)
basis.b[[vfunc[i]]]$nbasis <- length(int)
basis.b[[vfunc[i]]]$dropind <- NULL
basis.b[[vfunc[i]]]$names <- basis.b[[vfunc[i]]]$names[int]
}
x.fd = Data2fd(argvals = tt, y = t(xcc[[1]]$data),
basisobj = basis.x[[vfunc[i]]], fdnames = fdnames)
r = x.fd[[2]][[3]]
J = inprod(basis.x[[vfunc[i]]], basis.b[[vfunc[i]]])
Z = t(x.fd$coefs) %*% J
colnames(J) = colnames(Z) = name.coef[[vfunc[i]]] = paste(vfunc[i],
".", basis.b[[vfunc[i]]]$names, sep = "")
XX = cbind(XX, Z)
for (j in 1:length(colnames(Z))) {
if (kterms >= 1)
pf <- paste(pf, "+", colnames(Z)[j], sep = "")
else pf <- paste(pf, colnames(Z)[j], sep = "")
kterms <- kterms + 1
}
JJ[[vfunc[i]]] <- J
nbasisb <- basis.b[[vfunc[i]]]$nbasis
R = diag(0, ncol = nbasisb, nrow = nbasisb)
R <- eval.penalty(basis.b[[vfunc[i]]], Lfdobj = int2Lfd(2))
lenm <- ncol(mat2)
if (rn0)
stop("Ridge regressions is only implemented for functional principal component basis")
MM <- matrix(0, nrow = lenm, ncol = nbasisb)
MM2 <- matrix(0, nrow = nbasisb, ncol = imat2 +
nbasisb)
mat2 <- cbind(mat2, MM)
mat2 <- rbind(mat2, MM2)
if (!is.null(lambda[[vfunc[i]]])) {
mat2[(imat2 + 1):(imat2 + nbasisb), (imat2 +
1):(imat2 + nbasisb)] <- lambda[[vfunc[i]]] *
R
}
imat2 <- imat2 + nbasisb
}
else {
basis <- basis.x[[vfunc[i]]]
l <- basis$l
vs <- t(basis$basis$data)
basis$coefs <- basis$coefs[, l, drop = FALSE]
Z <- basis$coefs
response = "y"
colnames(Z) = name.coef[[vfunc[i]]] = paste(vfunc[i],
".", rownames(basis$basis$data), sep = "")
XX = cbind(XX, Z)
vs.list[[vfunc[i]]] = basis$basis
mean.list[[vfunc[i]]] = basis$mean
for (j in 1:length(colnames(Z))) {
if (kterms >= 1)
pf <- paste(pf, "+", name.coef[[vfunc[i]]][j],
sep = "")
else pf <- paste(pf, name.coef[[vfunc[i]]][j],
sep = "")
kterms <- kterms + 1
}
if (!is.null(rn[[vfunc[i]]])) {
mat <- c(mat, rep(rn[[vfunc[i]]], len = length(l)))
}
else mat <- c(mat, rep(0, len = length(l)))
lenl <- length(l)
lenm <- ncol(mat2)
MM <- matrix(0, nrow = lenm, ncol = lenl)
MM2 <- matrix(0, nrow = lenl, ncol = imat2 + lenl)
mat2 <- cbind(mat2, MM)
mat2 <- rbind(mat2, MM2)
if (!is.null(lambda[[vfunc[i]]])) {
R <- P.penalty(1:length(l))
mat2[(imat2 + 1):(imat2 + lenl), (imat2 +
1):(imat2 + lenl)] <- lambda[[vfunc[i]]] *
R
}
imat2 <- imat2 + lenl
}
}
else {
if (is(data[[vfunc[i]]],"fd")) {
fdat <- data[[vfunc[i]]]
if (is.null(basis.x[[vfunc[i]]]))
basis.x[[vfunc[i]]] <- fdat$basis
else if (inherits(basis.x[[vfunc[i]]], "pca.fd"))
bsp1 = FALSE
if (is.null(basis.b[[vfunc[i]]]) & bsp1)
basis.b[[vfunc[i]]] <- create.fdata.basis(fdat,
l = 1:max(5, floor(basis.x[[vfunc[i]]]$nbasis/5)),
type.basis = basis.x[[vfunc[i]]]$type,
rangeval = fdat$basis$rangeval)
else if (inherits(basis.x[[vfunc[i]]], "pca.fd"))
bsp2 = FALSE
if (bsp1 & bsp2) {
r = fdat[[2]][[3]]
if (!is.null(basis.x[[vfunc[i]]]$dropind)) {
int <- setdiff(1:basis.x[[vfunc[i]]]$nbasis,
basis.x[[vfunc[i]]]$dropind)
basis.x[[vfunc[i]]]$nbasis <- length(int)
basis.x[[vfunc[i]]]$dropind <- NULL
basis.x[[vfunc[i]]]$names <- basis.x[[vfunc[i]]]$names[int]
}
if (!is.null(basis.b[[vfunc[i]]]$dropind)) {
int <- setdiff(1:basis.b[[vfunc[i]]]$nbasis,
basis.b[[vfunc[i]]]$dropind)
basis.b[[vfunc[i]]]$nbasis <- length(int)
basis.b[[vfunc[i]]]$dropind <- NULL
basis.b[[vfunc[i]]]$names <- basis.b[[vfunc[i]]]$names[int]
}
J = inprod(basis.x[[vfunc[i]]], basis.b[[vfunc[i]]])
mean.list[[vfunc[i]]] <- mean.fd(x.fd)
x.fd <- center.fd(x.fd)
Z = t(x.fd$coefs) %*% J
colnames(J) = colnames(Z) = name.coef[[vfunc[i]]] = paste(vfunc[i],
".", basis.b[[vfunc[i]]]$names, sep = "")
XX = cbind(XX, Z)
for (j in 1:length(colnames(Z))) {
if (kterms >= 1)
pf <- paste(pf, "+", colnames(Z)[j],
sep = "")
else pf <- paste(pf, colnames(Z)[j], sep = "")
kterms <- kterms + 1
}
JJ[[vfunc[i]]] <- J
}
else {
l <- ncol(basis.x[[vfunc[i]]]$scores)
vs <- basis.x[[vfunc[i]]]$harmonics$coefs
Z <- basis.x[[vfunc[i]]]$scores
response = "y"
colnames(Z) = name.coef[[vfunc[i]]] = paste(vfunc[i],
".", colnames(basis.x[[vfunc[i]]]$harmonics$coefs),
sep = "")
XX = cbind(XX, Z)
vs.list[[vfunc[i]]] = vs
mean.list[[vfunc[i]]] = basis.x[[vfunc[i]]]$meanfd
for (j in 1:length(colnames(Z))) {
if (kterms >= 1)
pf <- paste(pf, "+", name.coef[[vfunc[i]]][j],
sep = "")
else pf <- paste(pf, name.coef[[vfunc[i]]][j],
sep = "")
kterms <- kterms + 1
}
}
}
else stop("Please, enter functional covariate")
}
}
}
if (!is.data.frame(XX))
XX = data.frame(XX)
par.fregre$formula = pf
par.fregre$data = XX
y <- XX[, 1]
scores <- as.matrix(XX[, -(1:2)])
W <- diag(weights)
error <- FALSE
if (missing(correlation))
correlation0 <- FALSE
else {
if (is.matrix(correlation))
error = FALSE
else error = TRUE
correlation0 <- TRUE
}
if (!rn0 & !lambda0 & !correlation0) {
# print(1)
W0 <- FALSE
z = lm(formula = pf, data = XX, ...)
e <- z$residuals
S <- solve(t(scores) %*% W %*% scores)
class(z) <- c(class(z), "fregre.lm")
mat0 <- diag(0, ncol(scores))
}
else {
# print(11)
mat0 <- diag(0, ncol(scores))
if (lambda0) {
mat0 <- mat2
}
if (rn0) {
mat0 <- diag(mat)
}
if (rn0 & lambda0)
warning("Only ridge penalization is done by rn argument (lambda argument is ignored)")
W <- diag(weights) %*% W
S <- solve(t(scores) %*% W %*% scores + mat0)
Cinv <- S %*% t(scores) %*% W
ycen = y - mean(y)
coefs <- Cinv %*% XX[, 1]
z <- list()
z$fitted.values <- yp <- drop(scores %*% coefs)
e <- z$residuals <- XX[, 1] - z$fitted.values
it <- 1
eps = 0.001
eps2 = 1e-10
err2 = sqrt(sum(coefs^2))
MM <- matrix(1:n, ncol = 1)
corStruct <- list()
# print(error)
while (error) {
W <- W0 <- diag(n)
if (is.list(correlation)) {
ncor <- length(correlation)
name.cor <- names(correlation)
}
else {
ncor <- 1
name.cor <- correlation[[1]]
correlation = list(correlation = list())
names(correlation) <- name.cor
}
wei0 <- wei <- diag(weights)
ee <- e
for (icor in 1:ncor) {
# print(icor); print(correlation)
if (!is.list(correlation[[icor]])) {
name.cor <- correlation[[icor]]
correlation[[icor]] = list()
names(correlation)[icor] <- name.cor
}
if (is.null(correlation[[icor]][["group"]])) {
n <- length(e)
# print(name.cor)
# print(name.cor[icor])
par.cor <- switch(name.cor[icor], cor.ARMA = {
cor.name <- "cor.ARMA"
#print("cor.ARMA")
if (is.null(correlation[[icor]][["index"]])) index <- 1:n else {
index.df <- correlation[[icor]][["index"]]
index <- data[["df"]][, index.df]
if (is.vector(index)) e <- e[order(index)] else {
par.ord <- list()
for (i in 1:ncol(index)) par.ord[[i]] <- index[,
i]
index <- as.numeric(do.call("order",
par.ord))
e <- e[index]
}
}
list(x = e[index])
}, cor.AR = {
cor.name <- "cor.AR"
# print("cor.AR")
if (is.null(correlation[[icor]][["index"]])) index <- 1:n else {
index.df <- correlation[[icor]][["index"]]
index <- data[["df"]][, index.df]
if (is.vector(index)) e <- e[order(index)]
else {
par.ord <- list()
for (i in 1:ncol(index)) par.ord[[i]] <- index[,i]
index <- as.numeric(do.call("order", par.ord))
e <- e[index]
}
}
list(x = e[index])
}, corExpo = {
index.df <- correlation[[icor]][["index"]]
dxy <- data[["df"]][, index.df]
index <- 1:length(e)
cor.name <- "corExpo"
list(xy = dxy)
}, corVgm = {
index.df <- correlation[[icor]][["index"]]
xy <- data[["df"]][, index.df]
dxy0 <- as.matrix(dist(xy, diag = TRUE, upper = TRUE))
index <- 1:n
cor.name <- "corVgm"
residual <- e[index]
list(dxy = dxy0, xy = xy, df = data.frame(residual))
}, cor.Exp = {
index.df <- correlation[[icor]][["index"]]
xy <- data[["df"]][, index.df]
index <- 1:n
cor.name <- "cor.Exp"
residual <- e[index]
dff <- data.frame(data[["df"]][, index.df],
residual)
names(dff) <- c(index.df, "residual")
list(df = dff)
}, corUnstruc = {
cor.name <- "corUnstruc"
list(correlation[[icor]][["index"]])
})
e <- e[index]
# print(names(par.cor))
# print(names(correlation))
# print("1")
if (!missing(control)) {
npar <- length(par.cor)
#print(names(par.cor))
#print(names(par.ord))
# print("2")
for (i in 1:length(control)) {
par.cor[npar + i] <- control[i]
names(par.cor)[npar + i] <- names(control)[i]
}
}
# print("antes do.call")
# print(cor.name)
# print(names(par.cor))
aa <- do.call(cor.name, par.cor)
W0 <- aa[[1]]
corStruct <- aa[-1]
}
else {
#print("hay grupo")
n <- nrow(data$df)
W00 <- matrix(0, ncol = n, nrow = n)
name.group <- correlation[[icor]]$group
ncomb <- length(name.group)
# print(name.group)
if (ncomb > 1) {
gr <- (data[["df"]][, name.group[1]])
for (ngr in 2:ncomb) gr <- paste(gr, (data[["df"]][,
name.group[ngr]]), sep = "")
gr <- factor(gr, levels = unique(gr))
}
#else gr <- data[["df"]][[name.group]]
else gr <- data[["df"]][,name.group]
if (!is.factor(gr))
gr <- factor(gr)
lev <- levels(gr)
nlev <- length(lev)
#print(name.cor[icor])
par.cor <- switch(name.cor[icor],
cor.ARMA = {
#print("entra cor.ARMAAAAAA")
index.df <- correlation[[icor]][["index"]]
#print(correlation)
#print(index.df)
index <- data[["df"]][, index.df]
e2 <- NULL
for (i in 1:nlev) {
jj <- gr == lev[i]
e0 <- e[jj]
e2 <- cbind(e2, e0[order(index[jj])])
}
cor.name <- "cor.ARMA"
#print(correlation)
#print("okk")
if (is.null(correlation[[icor]][["p"]])) p0<-1
else p0<-correlation[[icor]][["p"]]
#print(p0)
list(x = e2,p=p0,order.max=p0)
},
corCloud = {
index.df <- correlation[[icor]][["index"]]
index <- data[["df"]][, index.df]
gr <- data[["df"]][, correlation[[icor]][["group"]]]
xy <- data[["df"]][, index.df]
cor.name <- "corCloud"
residual <- e
dff <- data.frame(data[["df"]][, index.df],
residual)
names(dff) <- c(index.df, "residual")
list(df = dff, gr = gr)
},
cor.Exp = {
index.df <- correlation[[icor]][["index"]]
index <- data[["df"]][, index.df]
gr <- data[["df"]][, correlation[[icor]][["group"]]]
xy <- data[["df"]][, index.df]
cor.name <- "cor.Exp"
residual <- e
dff <- data.frame(data[["df"]][, index.df],
residual)
names(dff) <- c(index.df, "residual")
list(df = dff, gr = gr)
},
corUnstruc = {
cor.name <- "corUnstruc"
index <- 1:nrow(dff)
par.cor <- list(correlation[[icor]][["index"]])
})
if (!missing(control)) {
npar <- length(par.cor)
for (k in 1:length(control)) {
par.cor[npar + k] <- control[k]
names(par.cor)[npar + k] <- names(control)[k]
}
}
#print("calll")
#print(cor.name)
#print(par.cor)
#print("antes call")
aa <- do.call(cor.name, par.cor)
W0 <- aa[[1]]
corStruct <- aa[-1]
for (j in 1:nlev) {
n <- table(gr)[j]
ilev <- gr == lev[j]
e <- ee[ilev]
if (name.cor[icor] == "cor.Exp" | name.cor[icor] ==
"corCloud")
W00[ilev, ilev] <- W0
else W00[ilev, ilev] <- W0[order(index[ilev]),
order(index[ilev])]
}
W0 <- W00
}
wei0 <- wei0 %*% W0
}
W0 <- wei0
W <- try(solve(W0), silent = TRUE)
# print(class(W) )
if (is(W,"try-error")) {
sv <- svd(W0)
W <- drop((sv$v %*% diag(1/sv$d) %*% t(sv$u)))
}
mat0[1, 1] <- 0
W2 <- t(scores) %*% W %*% scores + mat0
S <- try(solve(W2), silent = TRUE)
if (is(S,"try-error")) {
sv <- svd(W2)
S <- drop((sv$v %*% diag(1/sv$d) %*% t(sv$u)))
}
Cinv <- S %*% t(scores) %*% W
coefs <- Cinv %*% matrix(y, ncol = 1)
yp <- drop(scores %*% coefs)
coefs <- drop(coefs)
e <- y - yp
err3 <- coefs
# err4 <- max(abs((err3 - err2)/err2))
# riq <- diff(summary(e)[c(2,5)])
riq <- diff(quantile(e,probs=c(.25,.75)))
err4 <- max(abs(err3 - err2)*riq)/max(abs(err2)+eps2)
if (err4 < eps)
error <- FALSE
else {
if (it == maxit)
error <- FALSE
else it <- it + 1
err2 <- err3
}
}
H <- scores %*% Cinv
df <- traza(H)
coefs <- drop(coefs)
z$coefficients <- coefs
z$mean.list <- mean.list
z$df.residual <- n - df
z$H <- H
z$r2 <- 1 - sum(z$residuals^2)/sum(ycen^2)
if (is(basis.x[[vfunc[1]]],"basisfd")) {
z$call[[1]] = "fregre.basis"
}
else {
if (basis.x[[vfunc[1]]]$type == "pc")
z$call[[1]] = "fregre.pc"
if (basis.x[[vfunc[1]]]$type == "pls")
z$call[[1]] = "fregre.pls"
}
class(z) <- c("fregre.fd", "fregre.lm")
rdf <- n - df
sr2 <- sum(e^2)/rdf
r2 <- 1 - sum(e^2)/sum(ycen^2)
r2.adj <- 1 - (1 - r2) * ((n - 1)/rdf)
GCV <- sum(e^2)/(n - df)^2
z$terms <- terms
z$residuals <- drop(e)
z$fitted.values <- yp
z$y <- y
z$rank <- df
z$df.residual <- rdf
Z = cbind(rep(1, len = n), Z)
colnames(Z)[1] = "(Intercept)"
std.error = sqrt(diag(S) * sr2)
t.value = coefs/std.error
p.value = 2 * pt(abs(t.value), n - df, lower.tail = FALSE)
coefficients <- cbind(coefs, std.error, t.value, p.value)
colnames(coefficients) <- c("Estimate", "Std. Error",
"t value", "Pr(>|t|)")
z$coefs <- coefficients
class(z) <- "lm"
z$it <- it
}
for (i in 1:length(vfunc)) {
if (bsp1)
beta.l[[vfunc[i]]] = fd(z$coefficients[name.coef[[vfunc[i]]]],
basis.b[[vfunc[i]]])
else {
if (inherits(data[[vfunc[i]]], "fdata")) {
beta.est <- z$coefficients[name.coef[[vfunc[i]]]] *
vs.list[[vfunc[i]]]
beta.est$data <- colSums(beta.est$data)
beta.est$names$main <- "beta.est"
beta.est$data <- matrix(as.numeric(beta.est$data),
nrow = 1)
beta.est$names$main <- "beta.est"
beta.est$data <- matrix(as.numeric(beta.est$data),
nrow = 1)
if (basis.x[[vfunc[i]]]$type == "pls") {
if (basis.x[[vfunc[i]]]$norm) {
sd.X <- sqrt(apply(data[[vfunc[i]]]$data,
2, var))
beta.est$data <- beta.est$data/sd.X
}
}
beta.l[[vfunc[i]]] <- beta.est
}
else {
beta.est <- z$coefficients[name.coef[[vfunc[i]]]] *
t(vs.list[[vfunc[i]]])
beta.est <- colSums(beta.est)
beta.l[[vfunc[i]]] <- fd(beta.est, basis.x[[vfunc[i]]]$harmonics$basis)
}
}
}
z$sr2 <- sum(e^2)/z$df.residual
z$Vp = z$sr2 * S
z$beta.l = beta.l
z$formula = pf
z$mean = mean.list
z$formula.ini = formula
z$basis.x = basis.x
z$basis.b = basis.b
z$JJ <- JJ
z$data = z$data
z$XX = XX
z$data <- data
z$fdataobj <- data[[vfunc[1]]]
if (correlation0) {
rn0 <- TRUE
z$corStruct <- corStruct
z$it <- it
}
z$rn <- rn0
z$lambda <- lambda0
z$W <- W
z$W0 <- W0
z$correlation <- correlation
z$correl <- correlation0
z$vs.list = vs.list
class(z) <- c( "fregre.igls","lm")
z
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.