fitpsar2d <-
function(y,Xsp,Xpar=NULL,Xnopar=NULL,
knots=c(10,10),pord=c(2,2),bdeg=c(3,3),
x1lim=NULL,x2lim=NULL,
vary_init=var(y), tau_init = rep(1,2), tau_fixed = rep(FALSE,2),
rho_init = 0, rho_fixed = FALSE,
onlygam=FALSE,trace=FALSE,thr=1e-4,maxit=200,
offset=0,bold=NULL,weights=NULL,Wsp=NULL,post_estim=TRUE,
knots_nopar=NULL,bdeg_nopar=NULL,pord_nopar=NULL) {
# if(!exists("spline.des")) require(splines)
# if(!exists("ginv")) require(MASS)
# if(!exists("multiroot")) require(rootSolve)
nsp <- length(y)
x1 <- Xsp[,1]
x2 <- Xsp[,2]
if(trace) start <- proc.time()[3]
if(!onlygam)
{
if(is.null(x1lim)) x1lim <- c(min(x1)-0.01, max(x1)+0.01)
if(is.null(x2lim)) x2lim <- c(min(x2)-0.01, max(x2)+0.01)
MM1 <- MM_basis(x1, x1lim[1], x1lim[2], knots[1], bdeg[1], pord[1])
MM2 <- MM_basis(x2, x2lim[1], x2lim[2], knots[2], bdeg[2], pord[2])
X1 <- MM1$X; Z1<-MM1$Z; d1<-MM1$d; B1<-MM1$B
X2 <- MM2$X; Z2<-MM2$Z; d2<-MM2$d; B2<-MM2$B
c1 <- ncol(B1); c2 <- ncol(B2); rm(MM1,MM2,B1,B2)
g1u <- rep(1,pord[2])%x%d1
g2u <- d2%x%rep(1,pord[1])
g1b <- rep(1,c2-pord[2])%x%d1
g2b <- d2%x%rep(1,c1-pord[1])
X_spt <- Rten2(X2, X1)
Z_spt <- cbind(Rten2(X2, Z1), Rten2(Z2, X1),Rten2(Z2, Z1))
rm(X1,X2,Z1,Z2)
ncolX_spt<-ncol(X_spt)
ncolZ_spt<-ncol(Z_spt)
} else {
X_spt <- 1; Z_spt <- NULL
ncolX_spt <- 1; ncolZ_spt <- 0
}
X <- X_spt
Z <- Z_spt
if (is.null(Xpar))
{
X_par <- NULL
nvar_par <- 0
} else {
X_par <- Xpar
nvar_par <- ncol(X_par)
X <- cbind(X_par,X)
}
# Number of parameters in each part
if(!onlygam)
{
np <- c(prod(pord)+nvar_par, (c1-pord[1])*pord[2],
(c2-pord[2])*pord[1],(c1-pord[1])*(c2-pord[2]))
} else { np <- ncol(X) }
if(!is.null(Xnopar))
{
nvar_nopar <- ncol(Xnopar)
X_nopar <- NULL; Z_nopar <- NULL; d_nopar <- NULL; c_nopar <- NULL
if(is.null(knots_nopar)) knots_nopar <- rep(10,nvar_nopar)
if(is.null(bdeg_nopar)) bdeg_nopar <- rep(3,nvar_nopar)
if(is.null(pord_nopar)) pord_nopar <- rep(2,nvar_nopar)
for(i in 1:nvar_nopar)
{
var_np<-Xnopar[,i]
MM.var <- MM_basis(var_np,min(var_np)-0.01,max(var_np)+0.01,
knots_nopar[i],bdeg_nopar[i],pord_nopar[i])
X.var <- MM.var$X; Z.var <- MM.var$Z; d.var <- MM.var$d; B.var <- MM.var$B
X_nopar <- cbind(X_nopar,X.var[,c(-1)]) # Quita columna intercepto
Z_nopar <- cbind(Z_nopar,Z.var)
d_nopar <- cbind(d_nopar,d.var)
c_nopar <- c(c_nopar,ncol(B.var))
rm(MM.var,B.var)
}
X <- cbind(X,X_nopar)
Z <- cbind(Z,Z_nopar)
np[1] <- np[1]+ncol(X_nopar)
np <- c(np,c_nopar-pord_nopar)
} else {
nvar_nopar <- 0
X_nopar <- NULL; Z_nopar <- NULL; d_nopar <- NULL; c_nopar <- NULL
}
# 0 0
# 0 I
D <- diag(c(rep(0,np[1]), rep(1,sum(np[-1]))))
if(!onlygam)
{
G1inv.n <- c(g1u, rep(0, np[3]), g1b)
G2inv.n <- c(rep(0,np[2]), g2u, g2b)
# Initialize variance components
la <- c(vary_init,tau_init[1:2],rho_init)
ed <- c(0,0)
} else {
la <- c(vary_init,rho_init)
ed <- NULL
}
# AÑADE PARÁMETROS SMOOTHING COVARIABLES NO PARAMÉTRICAS
if(!is.null(Xnopar))
{
la <- c(la,rep(1,nvar_nopar))
ed <- c(ed,rep(0,nvar_nopar))
}
# Initialise the parameters
if(is.null(bold)) bold <- rep(0,sum(np))
if(!is.null(Z)){ eta <- X%*%bold[1:np[1]] + Z%*%bold[-(1:np[1])] + offset } else {
eta <- X%*%bold[1:np[1]] + offset }
for (iq in 1:maxit) { # Nested loops for SAP and Scores of rho
w <- rep(1,length(y))
if(trace) start1 <- proc.time()[3]
for (it in 1:maxit) { # Loop for SAP
if(trace) start2 <- proc.time()[3]
# CÁLCULO MATRIZ A PARA MODELO PS-SAR
if(!onlygam) rho <- la[4] else rho <- la[2]
if(!is.null(Wsp))
{
A <- Matrix::Diagonal(nsp)-rho*Wsp
} else {
A <- Matrix::Diagonal(nsp)
}
Ay <- A%*%y
# CASO GAUSSIANO
mat <- construct_matrices(as.matrix(X),as.matrix(Z),
as.matrix(Ay),w=1,GLAM=FALSE)
# Build penalty matrix: block diagonal matrix
if(!onlygam)
{
Ginv <- c((1/la[2])*g1u, (1/la[3])*g2u, (1/la[2])*g1b + (1/la[3])*g2b)
} else { Ginv <- NULL }
if(!is.null(Xnopar))
{
for(k in 1:nvar_nopar)
{
if(!onlygam){ Ginv <- c(Ginv,(1/la[4+k])*d_nopar[,k]) } else {
Ginv <- c(Ginv,(1/la[2+k])*d_nopar[,k]) }
}
}
G <- 1/Ginv
# X'WX X'WZG
# Z'WX Z'WZG
C <- construct_block(mat$XtX, t(mat$ZtX*G), mat$ZtX, t(mat$ZtZ*G))
# ES LA MATRIZ C DE COEFICIENTES DEL SISTEMA (12) EN PAPER SAP
Hinv <- try(solve((1/la[1])*C + D))
if(class(Hinv) == "try-error") {
Hinv <- MASS::ginv((1/la[1])*C + D)
}
b <- (1/la[1])*Hinv%*%mat$u
b_fixed <- b[1:np[1]]
b_random <- G*b[-(1:np[1])]
if(nvar_par==0) {b_fixed_par <- NULL} else {b_fixed_par <- b_fixed[1:nvar_par]}
sig2u<-la[1]
# Compute effective dimensions and variances
# Only the diagonal of ZtPZ
dZtPZ <- 1/la[1]*apply((t(Hinv[-(1:np[1]),])*mat$ZtXtZ),2,sum)
if(!onlygam)
{
# Tau 1
if(is.null(Xnopar)) {G1inv.d <- (1/la[2])*G1inv.n} else{
G1inv.d <- c((1/la[2])*G1inv.n,rep(0,sum(np[-c(1:4)])))}
ed1 <- sum(dZtPZ*(G1inv.d*G^2))
ed1 <- ifelse(ed1 == 0, 1e-50,ed1)
if(!tau_fixed[1])
{
tau1 <- ifelse(is.null(Xnopar),sum(b_random^2*G1inv.n)/ed1,
sum(b_random^2*c(G1inv.n,rep(0,sum(np[-c(1:4)]))))/ed1)
tau1 <- ifelse(tau1 == 0, 1e-50, tau1)
} else {tau1 <- tau_init[1]}
# Tau 2
if(is.null(Xnopar)) {G2inv.d <- (1/la[3])*G2inv.n} else{
G2inv.d <- c((1/la[3])*G2inv.n,rep(0,sum(np[-c(1:4)])))}
ed2 <- sum(dZtPZ*(G2inv.d*G^2))
ed2 <- ifelse(ed2 == 0, 1e-50, ed2)
if(!tau_fixed[2])
{
tau2 <- ifelse(is.null(Xnopar),sum(b_random^2*G2inv.n)/ed2,
sum(b_random^2*c(G2inv.n,rep(0,sum(np[-c(1:4)]))))/ed2)
tau2 <- ifelse(tau2 == 0, 1e-50, tau2)
} else {tau2 <- tau_init[2]}
}
if(!is.null(Xnopar))
{
Ginv.d_nopar <- matrix(0,nrow=sum(np[2:length(np)]),ncol=ncol(Xnopar))
tau_nopar <- rep(0,nvar_nopar); edf_nopar <- rep(0,nvar_nopar)
for (k in 1:nvar_nopar)
{
if(!onlygam)
{
Ginv.d_nopar[(sum(np[2:(4+k-1)])+1):sum(np[2:(4+k)]),k] <- (1/la[4+k])*d_nopar[,k]
} else {
cumnp_nopar <- cumsum(np[-1])
if(k==1)
{Ginv.d_nopar[(1:cumnp_nopar[1]),k] <- (1/la[2+k])*d_nopar[,k]} else {
Ginv.d_nopar[((cumnp_nopar[k-1]+1):cumnp_nopar[k]),k] <- (1/la[2+k])*d_nopar[,k]}
}
edf_nopar[k] <- sum(dZtPZ*(Ginv.d_nopar[,k]*G^2))
edf_nopar[k] <- ifelse(edf_nopar[k] == 0, 1e-50, edf_nopar[k])
tau_nopar[k] <- sum(b_random^2*Ginv.d_nopar[,k])/edf_nopar[k]
tau_nopar[k] <- ifelse(tau_nopar[k] == 0, 1e-50, tau_nopar[k])
}
}
P.Hv <- (1/sig2u)*(Matrix::Diagonal(nsp)-
(1/sig2u)*( cbind(X, Z %*% Matrix::Matrix(diag(G))) ) %*%
( Hinv %*% Matrix::t(cbind(X,Z))) )
# Sigma^2
ssr <- mat$yty - t(c(b_fixed, b_random)) %*% (2*mat$u - C %*% b)
# OJO: SE DISMINUYE EN UN GRADO DE LIBERTAD POR LA ESTIMACIÓN DE RHO
if(!onlygam) edf_tot <- ed1+ed2+length(b_fixed) else edf_tot <- length(b_fixed)
if(!rho_fixed) edf_tot <- edf_tot+1
if(!is.null(Xnopar)){edf_tot <- edf_tot+sum(edf_nopar)}
sig2u <- as.numeric((ssr / (length(y) - edf_tot)))
# New variance components and convergence check
if(!onlygam){lanew = c(sig2u, tau1, tau2,rho)} else {lanew = c(sig2u,rho)}
if(!is.null(Xnopar)) {lanew = c(lanew,tau_nopar)}
dla = mean(abs(la - lanew))
la = lanew
if(trace) {
cat(sprintf("%1$3d %2$10.6f", it, dla))
if(!onlygam) cat(sprintf("%8.3f", c(ed1, ed2)), '\n')
cat('la ',la,"\n")
}
if (dla < thr) break
} # end for (it in 1:maxit) SAP loop
if (trace) {
end2 <- proc.time()[3]
cat("Timings:\nSAP", (end2-start2), "seconds\n")
}
# Se calcula el score analítico de L_{R} respecto a rho.
if(!rho_fixed)
{
score.REML_rho <- function(rho)
{
A <- as(Matrix::Diagonal(nsp)-rho*Wsp,"CsparseMatrix")
score.REML_rho <- t(as.matrix(P.Hv)%*%as.matrix(A%*%y))%*%as.matrix(Wsp%*%y) -
sum(diag(as.matrix(solve(as.matrix(A),as.matrix(Wsp)))))
as.numeric(score.REML_rho)
}
# Se iguala el score a 0 y se resuelve numéricamente la ecuación unidimensional
rho_roots <- rootSolve::multiroot(score.REML_rho,start=rho)
if(abs(rho_roots$f.root[1])<1.0E-6 &
rho_roots$root[1]<1 & rho_roots$root[1]> (-1))
{rho <- rho_roots$root[1]}
} else {rho <- rho_init}
if(!onlygam) rho_old <- la[4] else rho_old <- la[2]
drho <- abs(rho_old-rho)
if(!onlygam) la[4] <- rho else la[2] <- rho
if(trace)
{
cat(sprintf("%1$3d %2$10.6f", iq, drho), '\n')
cat("la ",la,'\n')
}
if (drho < thr) break
} # End loop for SAP and Scores of phi and rho
if (trace)
{
end <- proc.time()[3]
cat("All process", end - start, "seconds\n")
}
eta <- X%*%b_fixed + Z%*%b_random + offset
if(is.null(Xnopar)){ edf_nopar <- NULL; tau_nopar<-NULL }
# ESTIMATES OF PARAMETERS
sig2u <- la[1]
if(!onlygam)
{ tau_sp <- la[2:3]; edf_sp <- c(ed1,ed2); rho <- la[4]
} else {
tau_sp<-NULL; edf_sp <- NULL; rho <- la[2]
}
G.eff <- G # NO ESTAMOS TRABAJANDO CON PS-ANOVA
Ginv.eff <- Ginv
if(!is.null(Wsp)) {
A <- as(Matrix::Diagonal(nsp)-rho*Wsp,"CsparseMatrix") } else {
A <- Matrix::Diagonal(nsp) }
# Valor del log.lik.reml en el óptimo
V <- Z %*% (Matrix::Matrix(diag(G.eff)) %*% Matrix::t(Z)) +
sig2u*Matrix::Diagonal(nsp)
Vinv <- try(Matrix::solve(V))
if(class(V) == "try-error")
{
Vinv <- Matrix::Matrix(MASS::ginv(as.matrix(V)))
}
Vinv.X <- Vinv %*% X
Xt.Vinv.X.inv <- try(Matrix::solve(Matrix::t(X) %*% (Vinv.X)))
if(class(Xt.Vinv.X.inv)=="try-error")
{
Xt.Vinv.X.inv <- Matrix::Matrix(
MASS::ginv(as.matrix(t(as.matrix(X)) %*% (Vinv.X))) )
}
P <- try(Vinv-Vinv.X %*% Xt.Vinv.X.inv %*% Matrix::t(Vinv.X))
if(class(P) == "try-error")
{
P <- Matrix::Matrix( Vinv - (Vinv.X %*%
MASS::ginv( as.matrix((Matrix::t(X) %*% Vinv.X) %*%
Matrix::t(Vinv.X)) ) ) )
}
log.lik.reml <- as.numeric( -0.5*(Matrix::determinant(V)$modulus+
Matrix::determinant(Matrix::t(X) %*% Vinv.X)$modulus +
Matrix::t(A %*% y) %*% (P %*% (A %*%y))) +
Matrix::determinant(A)$modulus )
log.lik <- as.numeric( -0.5*(Matrix::determinant(V)$modulus +
Matrix::t( A%*% y) %*% (P %*% (A %*% y))) +
Matrix::determinant(A)$modulus )
# Asignación Efectos Fijos y Aleatorios según tipo variables
b_fixed_par <- NULL; b_fixed_nopar <- NULL; b_random_nopar <- NULL
if(nvar_par>0) b_fixed_par <- b_fixed[1:ncol(X_par)]
if(nvar_nopar>0)
{
b_fixed_nopar <- b_fixed[(ncol(X)-ncol(X_nopar)+1):ncol(X)]
if(!onlygam) b_random_spt <- b_random[1:(ncol(Z)-ncol(Z_nopar))] else b_random_spt <- NULL
b_random_nopar <- b_random[(ncol(Z)-ncol(Z_nopar)+1):ncol(Z)]
if(nvar_par>0)
{
b_fixed_spt <- b_fixed[(ncol(X_par)+1):(ncol(X)-ncol(X_nopar))]
} else { b_fixed_spt <- b_fixed[1:(ncol(X)-ncol(X_nopar))] }
} else {
if(!onlygam) b_random_spt <- b_random else b_random_spt <- NULL
if(nvar_par>0) { b_fixed_spt <- b_fixed[(ncol(X_par)+1):ncol(X)]
} else { b_fixed_spt <- b_fixed[1:ncol(X)]}
}
if(post_estim==TRUE)
{
# # CÁLCULO MATRIZ COVARIANZAS EFECTOS FIJOS Y ALEATORIOS
# # pp.375 Fahrmeir et al.
# # Se reescala A multiplicándola por sig2u toda la matriz
# # Matriz Covarianzas Bayesiana.
A_cov1 <- rbind( cbind(Matrix::crossprod(X),Matrix::t(X) %*% Z),
cbind(Matrix::t(Z) %*% X,
Matrix::crossprod(Z) +
sig2u*Matrix::Matrix(diag(Ginv.eff))) )
cov1_eff <- try(sig2u*Matrix::solve(A_cov1))
if(class(cov1_eff) == "try-error") {
cov1_eff <- sig2u*MASS::ginv(as.matrix(A_cov1))
cov1_eff.ginv <- TRUE
} else { cov1_eff.ginv <- FALSE }
cov1_eff <- as.matrix(cov1_eff)
var1_b_fixed <- diag(cov1_eff)[1:(ncol(X))]
var1_b_random <- diag(cov1_eff)[(ncol(X)+1):(ncol(X)+ncol(Z))]
var1_b_fixed_par <- NULL; cov1_b_fixed_par <- NULL
var1_b_fixed_nopar <- NULL; var1_b_random_nopar <- NULL
if(nvar_par>0)
{
var1_b_fixed_par <- diag(cov1_eff)[1:ncol(X_par)]
cov1_b_fixed_par <- cov1_eff[1:ncol(X_par),1:ncol(X_par)]
}
if(nvar_nopar>0)
{
if(!onlygam)
{
var1_b_random_spt <- diag(cov1_eff)[(ncol(X)+1):
(ncol(X)+ncol(Z)-ncol(Z_nopar))]
var1_b_random_nopar <- diag(cov1_eff)[(ncol(X) +
length(var1_b_random_spt)+1):
(ncol(X)+ncol(Z))]
} else {
var1_b_random_spt <- NULL
var1_b_random_nopar <- diag(cov1_eff)[(ncol(X)+1):(ncol(X)+ncol(Z))]
}
if(nvar_par>0){
var1_b_fixed_spt <- diag(cov1_eff)[(ncol(X_par)+1):
(ncol(X)-ncol(X_nopar))]
var1_b_fixed_nopar <-
diag(cov1_eff)[(ncol(X_par)+length(b_fixed_spt)+1):
ncol(X)]
} else {
var1_b_fixed_spt <- diag(cov1_eff)[1:ncol(X_spt)]
var1_b_fixed_nopar <- diag(cov1_eff)[(ncol(X_spt)+1):ncol(X)]}
} else {
if(!onlygam)
{
var1_b_random_spt <- diag(cov1_eff)[(ncol(X)+1):
(ncol(X)+ncol(Z))]
} else var1_b_random_spt <- NULL
if(nvar_par>0){
var1_b_fixed_spt <- diag(cov1_eff)[(ncol(X_par)+1):ncol(X)]
} else {
var1_b_fixed_spt <- diag(cov1_eff)[1:ncol(X)]}
}
# # CALCULAR VALORES ESTIMADOS Y RESIDUOS
fit_Ay <- as.vector(X %*% b_fixed + Z %*% b_random + offset)
fit <- try(Matrix::solve(A) %*% fit_Ay)
if(class(fit) == "try-error")
{
fit <- MASS::ginv(as.matrix(A)) %*% eta
}
fit <- as.vector(fit)
# mH = kron(Ainv,It)*(X Z)*VAR(beta,alpha)*t(X Z)*Rinv*kron(A,It)
mH <- (1/sig2u)*(Matrix::solve(A) %*% cbind(X,Z)) %*%
(cov1_eff %*% (Matrix::t(cbind(X,Z)) %*% A))
var_fit_Ay <- as.matrix( cbind(X,Z) %*% cov1_eff %*%
Matrix::t(cbind(X,Z)) )
var_fit <- try(Matrix::solve(A) %*% (var_fit_Ay %*%
Matrix::t(Matrix::solve(A))) )
if(class(var_fit) == "try-error")
{
MASS::ginv(as.matrix(A)) %*% (var_fit_Ay %*%
t(MASS::ginv(as.matrix(A))))
}
var_fit <- as.vector(var_fit)
resids <- as.vector( (A %*% y) - fit_Ay )
# Compute AIC y BIC based on loglik functions (Fahrmeir, pp. 664 and 677)
aic <- -2*log.lik + 2*edf_tot
bic <- -2*log.lik + log(length(y))*edf_tot
# CÁLCULO AIC and BIC based on ssr: LIBRO FAHRMEIR pp. 564 (STAR MODELS WITH GAUSSIAN ERRORS)
# VIP: SSR CON RESIDUOS NORMALIZADOS
ssr <- sum(resids^2)
aic.sig2 <- length(y)*log(ssr/length(y))+2*(edf_tot+1)
bic.sig2 <- length(y)*log(ssr/length(y))+log(length(y))*(edf_tot+1)
# # COMPUTE FITS EFFECTS BY VARIABLE
fit_cov_par <- NULL; fit_cov_nopar <- NULL
fit_cov_nopar_fixed <- NULL; fit_cov_nopar_random <- NULL
var_fit_cov_nopar <- NULL; k.Z_nopar <- NULL; edf_nopar_tot <- NULL
if(nvar_par>0)
{
for(i in 1:nvar_par)
{
fit_cov_par <- cbind(fit_cov_par,X_par[,i]*b_fixed_par[i])
}
}
if(nvar_nopar>0)
{
edf_nopar_tot <- rep(0,nvar_nopar)
k.Z_nopar <- rep(0,nvar_nopar)
for(i in 1:nvar_nopar)
{
fit_cov_nopar_fixed <- cbind(fit_cov_nopar_fixed,
X_nopar[,i]*b_fixed_nopar[i])
var_np<-Xnopar[,i]
# Matriz selección para varianza f(x_j) Libro Ruppert et al, pp. 175
E_var_np <- matrix(0,ncol=ncol(cbind(X,Z)),nrow=ncol(cbind(X,Z)))
vindex.var_np_fixed <- (ncol(X)-ncol(X_nopar)+i)
E_var_np[vindex.var_np_fixed,vindex.var_np_fixed] <- 1
MM.var_np <- MM_basis(var_np,min(var_np)-0.01,max(var_np)+0.01,
knots_nopar[i],bdeg_nopar[i],pord_nopar[i])
X.var_np <- MM.var_np$X; Z.var_np <- MM.var_np$Z
k.Z_nopar[i] <- ncol(Z.var_np)
rm(MM.var_np)
if(i==1){
fit_cov_nopar_random <- cbind(fit_cov_nopar_random,
Z.var_np%*%b_random_nopar[1:k.Z_nopar[1]])
vindex.var_np_random <- c((ncol(X)+ncolZ_spt+1):
(ncol(X)+ncolZ_spt+k.Z_nopar[1]))
} else {
fit_cov_nopar_random <- cbind(fit_cov_nopar_random,
Z.var_np%*%b_random_nopar[(sum(k.Z_nopar[1:(i-1)])+1):
sum(k.Z_nopar[1:i])])
vindex.var_np_random<- c((ncol(X)+ncolZ_spt+
sum(k.Z_nopar[1:(i-1)])+1):
(ncol(X)+ncolZ_spt+sum(k.Z_nopar[1:i])))
}
# Se cambia la diagonal de la matriz de Selección
for(j in 1:length(vindex.var_np_random))
{
E_var_np[vindex.var_np_random[j],vindex.var_np_random[j]] <- 1
}
mH_var_np <- as.matrix( (1/sig2u)*(Matrix::solve(A) %*% cbind(X,Z) %*%
E_var_np) %*% (cov1_eff %*%
(Matrix::t(cbind(X,Z))%*%A)) )
edf_nopar_tot[i] <- sum(diag(mH_var_np))
var_fit_np <- rep(0,length(fit))
for(j in 1:length(fit))
{
var_fit_np[j] <- sig2u*sum(mH_var_np[j,]^2)
}
var_fit_cov_nopar <- cbind(var_fit_cov_nopar,var_fit_np)
}
fit_cov_nopar <- fit_cov_nopar_fixed+fit_cov_nopar_random
}
# # COMPUTE SPATIO-TEMPORAL TREND.
if(!onlygam)
{
spt_trend_fixed <- as.vector(X_spt%*%b_fixed_spt)
spt_trend_random <- as.vector(Z_spt%*%b_random_spt)
spt_trend <- spt_trend_fixed+spt_trend_random
E_spt_trend <- matrix(0,ncol=ncol(cbind(X,Z)),nrow=ncol(cbind(X,Z)))
vindex_spt_trend_fixed <- c((nvar_par+1):(ncol(X)-nvar_nopar))
for(i in 1:length(vindex_spt_trend_fixed))
{
E_spt_trend[vindex_spt_trend_fixed[i],
vindex_spt_trend_fixed[i]] <- 1
}
vindex_spt_trend_random <- c((ncol(X)+1):(ncol(X)+ncolZ_spt))
for(i in 1:length(vindex_spt_trend_random))
{
E_spt_trend[vindex_spt_trend_random[i],
vindex_spt_trend_random[i]] <- 1
}
mH_spt_trend <- as.matrix( (1/sig2u)*(Matrix::solve(A) %*% cbind(X,Z) %*%
E_spt_trend) %*%
(cov1_eff %*% (Matrix::t(cbind(X,Z)) %*% A)) )
edf_spt_tot <- sum(diag(mH_spt_trend))
var_spt_trend <- rep(0,length(spt_trend))
for(j in 1:length(spt_trend))
{
var_spt_trend[j] <- sig2u*sum(mH_spt_trend[j,]^2)
}
spt_trend_fixed <- matrix(spt_trend_fixed,ncol=nsp,byrow=FALSE)
spt_trend_random <- matrix(spt_trend_random,ncol=nsp,byrow=FALSE)
spt_trend <- matrix(spt_trend,ncol=nsp,byrow=FALSE)
var_spt_trend <- matrix(var_spt_trend,ncol=nsp,byrow=FALSE)
} else {
spt_trend_fixed <- as.vector(X_spt%*%b_fixed_spt)
spt_trend <- spt_trend_fixed
spt_trend_random <- NULL
var_spt_trend <- NULL
edf_spt_tot <- 1
}
} # end if(post_estim==TRUE)
res <- list (call = match.call(),
# Parámetros estimados
knots = knots, bdeg = bdeg, pord = pord,
knots_nopar=knots_nopar,bdeg_nopar=bdeg_nopar,pord_nopar=pord_nopar,
tau_sp = tau_sp,tau_nopar=tau_nopar,
edf_sp = edf_sp,edf_nopar=edf_nopar,edf_tot=edf_tot,
edf_spt_tot=edf_spt_tot,
sig2u=sig2u,rho=rho,
ssr=ssr,aic=aic,bic=bic,
aic.sig2=aic.sig2,bic.sig2=bic.sig2,
log.lik.reml=log.lik.reml,log.lik=log.lik,
b_fixed=b_fixed,b_random=b_random,
b_fixed_spt=b_fixed_spt,b_random_spt=b_random_spt,
b_fixed_par=b_fixed_par,b_fixed_nopar=b_fixed_nopar,
b_random_nopar=b_random_nopar,
var1_b_fixed=var1_b_fixed,
var1_b_random=var1_b_random,
var1_b_fixed_spt=var1_b_fixed_spt,
var1_b_random_spt=var1_b_random_spt,
var1_b_fixed_par=var1_b_fixed_par,
cov1_b_fixed_par=cov1_b_fixed_par,
var1_b_fixed_nopar=var1_b_fixed_nopar,
var1_b_random_nopar=var1_b_random_nopar,
# fit, resids y var(fit)
fit=fit,fit_Ay=fit_Ay,resids=resids,
var_fit=var_fit,var_fit_Ay=var_fit_Ay,fit_cov_par=fit_cov_par,
fit_cov_nopar=fit_cov_nopar,
fit_cov_nopar_fixed=fit_cov_nopar_fixed,
fit_cov_nopar_random=fit_cov_nopar_random,
var_fit_cov_nopar=var_fit_cov_nopar,
# tendencias espacio-temporales
spt_trend=spt_trend,
spt_trend_fixed=spt_trend_fixed,spt_trend_random=spt_trend_random,
var_spt_trend=var_spt_trend,eta=eta,maxit=it)
res
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.