fit_pspt3d_sar <- function(y,vary_init,sp1,sp2,time,Xfull,Zfull,Wsp = NULL,
nvarpar,nvarnopar,nvarspt,
#weights = NULL, GLAM = FALSE,
cspt,dsptlist,bdegspt,pordspt,nknotsspt,
cnopar,dnoparlist,bdegnopar,pordnopar,nknotsnopar,
names_varnopar,names_varpar,
psanova = TRUE,
f1_main = TRUE, f2_main = TRUE, ft_main = TRUE,
f12_int = FALSE, f1t_int = FALSE,
f2t_int = FALSE, f12t_int = FALSE,
sar = FALSE, rho_init = 0, rho_fixed = FALSE,
ar1 = FALSE, phi_init = 0, phi_fixed = FALSE,
bold = NULL, maxit = 20, thr = 1e-2,
trace = FALSE, var_num = FALSE)
{
nsp <- length(sp1); ntime <- length(time)
if (!is.null(Wsp)) Wsp <- Matrix::Matrix(Wsp)
In <- Matrix::Diagonal(nsp); It <- Matrix::Diagonal(ntime)
X <- Matrix::Matrix(Xfull)
Z <- Matrix::Matrix(Zfull)
if (!is.null(nknotsnopar)) nvarnopar <- length(nknotsnopar)
if (!sar){
rho_fixed <- TRUE
rho_init <- 0
} else if (is.null(rho_fixed)) rho_fixed <- FALSE
if (!ar1){
phi_fixed <- TRUE
phi_init <- 0
} else if (is.null(phi_fixed)) phi_fixed <- FALSE
# Build vector and matrices for variance components in mixed model
var_comp <- par_var_comp3d(la = var(as.numeric(y)), np_fixed = ncol(X),
pordspt = pordspt, cspt = cspt, dsptlist = dsptlist,
nvarnopar = nvarnopar,cnopar = cnopar,
pordnopar = pordnopar, dnoparlist = dnoparlist,
psanova = psanova, f1_main = f1_main,
f2_main = f2_main, ft_main = ft_main,
f12_int = f12_int, f1t_int = f1t_int,
f2t_int = f2t_int, f12t_int = f12t_int)
# Number of parameters
np <- var_comp$np; np_eff <- var_comp$np_eff
# Vector of parameters
la <- var_comp$la
# Other vectors for SAP
g1u <- var_comp$g1u; g2u <- var_comp$g2u; g3u <- var_comp$g3u
g11b <- var_comp$g11b; g21b <- var_comp$g21b
g12b <- var_comp$g12b; g31b <- var_comp$g31b
g22b <- var_comp$g22b; g32b <- var_comp$g32b
g1t <- var_comp$g1t; g2t <- var_comp$g2t; g3t <- var_comp$g3t
G1inv.n <- var_comp$G1inv.n; G2inv.n <- var_comp$G2inv.n
G3inv.n <- var_comp$G3inv.n
if (psanova) {
g12u <- var_comp$g12u; g21u <- var_comp$g21u
g12b <- var_comp$g12b; g21b <- var_comp$g21b
g13u <- var_comp$g13u; g31u <- var_comp$g31u
g13b <- var_comp$g13b; g31b <- var_comp$g31b
g23u <- var_comp$g23u; g32u <- var_comp$g32u
g23b <- var_comp$g23b; g32b <- var_comp$g32b
g123u <- var_comp$g123u; g213u <- var_comp$g213u
g321u <- var_comp$g321u
g123b <- var_comp$g123b; g213b <- var_comp$g213b
g132b <- var_comp$g132b; g312b <- var_comp$g312b
g231b <- var_comp$g231b; g321b <- var_comp$g321b
G4inv.n <- var_comp$G4inv.n; G5inv.n <- var_comp$G5inv.n
G6inv.n <- var_comp$G6inv.n; G7inv.n <- var_comp$G7inv.n
G8inv.n <- var_comp$G8inv.n; G9inv.n <- var_comp$G9inv.n
G10inv.n <- var_comp$G10inv.n; G11inv.n <- var_comp$G11inv.n
G12inv.n <- var_comp$G12inv.n
}
# Do not remove var_comp, it is used in the next loop...
# 0 0
# 0 I
D <- Matrix::Matrix(diag(c(rep(0,np_eff[1]), rep(1,sum(np_eff[-1])))))
ed <- rep(0,length(la)-1)
# add rho and phi to la vector of parameters
la <- c(la,rho_init,phi_init)
# Initialise the parameters
if (is.null(bold)) bold = rep(0,sum(np_eff))
eta <- X %*% bold[1:np_eff[1]] + Z %*% bold[-(1:np_eff[1])] #+ offset
start <- proc.time()[3]
for (iq in 1:maxit) { # Nested loops for SAP and phi (AR1) and rho (SAR)
for (it in 1:maxit) {
rho <- la[length(la)-1]
phi <- la[length(la)]
if (is.null(Wsp)) {
A <- In
} else {
A <- In - rho*Wsp
}
# fast index ntime, slow index nsp
A_It_y <- matrix(RH(A,RH(It, array(y,dim=c(ncol(It),ncol(A))))),ncol=1)
# Es lo mismo que esta operación (lo he comprobado):
# A_It_y <- matrix(kronecker(A,I_t)%*%y,ncol=1)
# Build covariance G matrix for random effects: block diagonal matrix
lG <- build_G3d(la = la, lg = var_comp,
nvarnopar = nvarnopar, dnoparlist = dnoparlist,
psanova = psanova, f1_main = f1_main, f2_main = f2_main,
ft_main = ft_main, f12_int = f12_int, f1t_int = f1t_int,
f2t_int = f2t_int, f12t_int = f12t_int)
G <- lG$G; Ginv <- lG$Ginv;
G_eff <- lG$G_eff; Ginv_eff <- lG$Ginv_eff
rm(lG)
sig2u <- la[1]
if (phi_fixed & phi_init==0) { # (Omega=I)
Xstar <- X
Zstar <- Z
A_It_ystar <- A_It_y } else {
# CONSTRUIMOS MATRIZ OMEGA (SI phi=0 ENTONCES OMEGA=I)
call.Omega <- build_Omega_ar1(phi,ntime)
Omega <- call.Omega$Omega
Omegainv <- call.Omega$Omegainv
chol.Omegainv <- lltA(Omegainv)
Xstar <- apply(aperm(RH(In,RH(chol.Omegainv,
array(X,dim=c(ncol(chol.Omegainv),ncol(In),ncol(X))))),
perm=c(2,1,3)),2,rbind)
Zstar <- apply(aperm(RH(In,RH(chol.Omegainv,
array(Z,dim=c(ncol(chol.Omegainv),ncol(In),ncol(Z))))),
perm=c(2,1,3)),2,rbind)
A_It_ystar <- matrix(RH(In,RH(chol.Omegainv,
array(A_It_y,dim=c(ncol(chol.Omegainv),ncol(In))))),ncol=1)
rm(call.Omega,Omegainv,chol.Omegainv)
}
if (is.null(weights)) {
w <- as.vector(matrix(1,nrow=length(y))) } else { w <- weights }
mat <- construct_matrices(Xstar,Zstar,A_It_ystar)
C <- Matrix::Matrix( construct_block(mat$XtX, Matrix::t(mat$ZtX*G_eff),
mat$ZtX, Matrix::t(mat$ZtZ*G_eff)) )
# ES LA MATRIZ C DE COEFICIENTES DEL SISTEMA (12) EN PAPER SAP
# NO SE PUEDE RESOLVER SI SE INCLUYEN LOS ZEROS EN MATRICES X,Z,G
Hinv <- try(Matrix::solve((1/la[1])*C + D,tol = 1e-40))
if (class(Hinv) == "try-error")
Hinv <- MASS::ginv((1/la[1])*as.matrix(C) + as.matrix(D))
b <- as.vector( (1/la[1])*Hinv %*% mat$u )
bfixed <- b[1:np_eff[1]]
names(bfixed) <- gsub("X_","",colnames(X))
names(bfixed) <- paste("fixed_",names(bfixed),sep="")
brandom <- G_eff*b[-(1:np_eff[1])]
names(brandom) <- gsub("Z_","",colnames(Z))
names(brandom) <- paste("random_",names(brandom),sep="")
# Compute effective dimensions and variances
# Only the diagonal of ZtPZ
dZtPZ <- 1/la[1]*apply((Matrix::t(Hinv[-(1:np_eff[1]),])*mat$ZtXtZ),2,sum)
dZtPZ_wide <- rep(0,length(G))
brandom_wide <- rep(0,length(G))
index.zeros.G <- G==0
index <- 1
for (ide in 1:length(G)) {
if (!index.zeros.G[ide]) {
brandom_wide[ide] <- brandom[index]
dZtPZ_wide[ide] <- dZtPZ[index]
index <- index + 1
}
}
# MIRAR AÑADIR tau_init, tau_fixed, taunopar_fixed y taunoparinit
ltau_edf <- update_tau3d(la = la, lg = var_comp, G = G,
dZtPZ_wide = dZtPZ_wide, brandom_wide = brandom_wide,
psanova = psanova, f1_main = f1_main, f2_main = f2_main,
ft_main = ft_main, f12_int = f12_int, f1t_int = f1t_int,
f2t_int = f2t_int, f12t_int = f12t_int,
nvarnopar = nvarnopar, dnoparlist = dnoparlist)
tau1 <- ltau_edf$tau1; tau2 <- ltau_edf$tau2; tau3 <- ltau_edf$tau3
ed1 <- ltau_edf$ed1; ed2 <- ltau_edf$ed2; ed3 <- ltau_edf$ed3
tauspt <- c(tau1,tau2,tau3); edfspt <- c(ed1,ed2,ed3)
names(tauspt) <- names(edfspt) <- c("sp1","sp2","time")
taunopar <- ltau_edf$taunopar; edfnopar <- ltau_edf$edfnopar
if (psanova){
tau4 <- ltau_edf$tau4; tau5 <- ltau_edf$tau5; tau6 <- ltau_edf$tau6;
tau7 <- ltau_edf$tau7; tau8 <- ltau_edf$tau8; tau9 <- ltau_edf$tau9;
tau10 <- ltau_edf$tau10; tau11 <- ltau_edf$tau11; tau12 <- ltau_edf$tau12
ed4 <- ltau_edf$ed4; ed5 <- ltau_edf$ed5; ed6 <- ltau_edf$ed6;
ed7 <- ltau_edf$ed7; ed8 <- ltau_edf$ed8; ed9 <- ltau_edf$ed9;
ed10 <- ltau_edf$ed10; ed11 <- ltau_edf$ed11; ed12 <- ltau_edf$ed12
tauspt <- c(tauspt,tau4,tau5,tau6,tau7,tau8,tau9,tau10,tau11,tau12)
edfspt <- c(edfspt,ed4,ed5,ed6,ed7,ed8,ed9,ed10,ed11,ed12)
names(tauspt) <- names(edfspt) <- c("f1_main","f2_main","ft_main",
"f12.1","f12.2","f1t.1","f1t.2","f2t.1","f2t.2",
"f12t.1","f12t.2","f12t.3")
}
# OJO: USA LAS VARIABLES STAR PARA CALCULAR SSR Y EDF'S
# Regression (Fahrmeir et al.) pp. 180
ssr <- as.numeric(mat$yty - t(c(bfixed, brandom)) %*% (2*mat$u - C %*% b))
# New variance components
lanew <- c(sig2u, tauspt)
edftot <- sum(edfspt) + ncol(Xstar)
if (nvarnopar>0){
lanew <- c(lanew,taunopar)
edftot <- edftot + sum(edfnopar)
}
if ((!rho_fixed) || !(rho==0)) edftot <- edftot + 1
if ((!phi_fixed) || !(phi==0)) edftot <- edftot + 1
sig2u <- as.numeric((ssr/(length(y) - edftot)))
# Update first component of la
lanew[1] <- sig2u
# Update last two components of la
lanew <- c(lanew,rho,phi)
dla <- mean(abs(la - lanew))
la <- lanew
if (trace) {
cat('\n Iteration SAP: ',it)
cat('\n sig2u ',la[1])
cat('\n edfspt:',round(edfspt,2))
if (!is.null(edfnopar)) cat('\n edfnopar: ',round(edfnopar,2),'\n')
}
# convergence check
if (dla < thr) break
} # end for (it in 1:maxit)
if(!rho_fixed & !phi_fixed) {
par_init <- c(rho,phi)
lb <- c(-Inf,-1)
ub <- c(1,1)
} else if (!rho_fixed & phi_fixed) {
par_init <- c(rho,phi_init)
lb <- c(-Inf,phi_init)
ub <- c(1,phi_init)
} else if (rho_fixed & !phi_fixed) {
par_init <- c(rho_init,phi)
lb <- c(rho_init,-1)
ub <- c(rho_init,1)
} else {
rho <- rho_init
phi <- phi_init
}
if (!rho_fixed | !phi_fixed){
if (trace) print_level = 1 else print_level = 0
rho_phi_optim <- nloptr::nloptr(x0=par_init,
eval_f=llik_reml_fn3d,
#eval_grad=score_phi,
lb=lb, ub=ub,
opts=list("algorithm"="NLOPT_LN_BOBYQA",
"xtol_rel"=thr,
"ftol_abs"=thr,
"maxtime"=100,
"print_level"=print_level),
sig2u=sig2u,nsp=nsp,ntime=ntime,
Wsp=Wsp,y=y,X=X,Z=Z,G_eff=G_eff,
np_eff=np_eff,
bfixed=bfixed,
rho_fixed=rho_fixed,
phi_fixed=phi_fixed)
rho <- rho_phi_optim$solution[1]
phi <- rho_phi_optim$solution[2]
}
# Check Convergence in rho and phi parameters
rho_old <- la[length(la)-1]
phi_old <- la[length(la)]
drhophi <- mean(abs(c(rho_old,phi_old)-c(rho,phi)))
la[length(la)-1] <- rho
la[length(la)] <- phi
if (trace & any(c(sar,ar1))) {
cat("\n Iteration SAR-AR1: ",iq)
if (sar) cat("\n rho ",rho)
if (ar1) cat("\n phi",phi)
}
if (drhophi < thr) break
} # End loop for SAR-AR1
end <- proc.time()[3]
cat("\n Time to fit the model: ", (end-start), "seconds")
eta <- Xstar%*%bfixed + Zstar%*%brandom #+ offset
if (is.null(names_varnopar)) { edfnopar <- NULL; taunopar<-NULL }
# FINAL ESTIMATES OF PARAMETERS
sig2u <- la[1]
# CREO QUE HABRÍA QUE ACTUALIZAR TAMBIÉN LOS EDF... REPASAR...
if(!psanova) tauspt <- la[2:4] else tauspt <- la[2:13]
rho <- la[length(la)-1]
phi <- la[length(la)]
# Valor de log.lik y log.lik.reml en el óptimo
foptim <- llik_reml_var3d(c(rho,phi),
sig2u=sig2u,nsp=nsp,ntime=ntime,
Wsp=Wsp,y=y,X=X,Z=Z,G_eff=G_eff,
np_eff=np_eff,bfixed=bfixed,
rho_fixed=rho_fixed,phi_fixed=phi_fixed)
llik_reml <- foptim$llik_reml
llik <- foptim$llik
se_rho_an <- foptim$se_rho
se_phi_an <- foptim$se_phi
cov_rho_phi_an <- foptim$cov_rho_phi
rm(foptim)
if (var_num){
hessian_optim <- numDeriv::hessian(llik_reml_fn3d,c(rho,phi),
sig2u=sig2u,nsp=nsp,ntime=ntime,
Wsp=Wsp,y=y,X=X,Z=Z,G_eff=G_eff,
np_eff=np_eff,
bfixed=bfixed,
rho_fixed=rho_fixed,
phi_fixed=phi_fixed)
var_rho_phi_num <- solve(hessian_optim)
rownames(var_rho_phi_num) <- colnames(var_rho_phi_num) <- c("rho","phi")
} else var_rho_phi_num <- NULL
########################## 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. OJO: Tiene en cuenta autocorrelación...
if (phi_fixed & phi_init == 0) { # (Omega=I)
Omega <- Omegainv <- It
} else {
# CONSTRUIMOS MATRIZ OMEGA (SI phi=0 ENTONCES OMEGA=I)
call.Omega <- build_Omega_ar1(phi,ntime)
Omega <- Matrix::Matrix(call.Omega$Omega)
Omegainv <- Matrix::Matrix(call.Omega$Omegainv)
rm(call.Omega)
}
Rinv <- kronecker(In,Omegainv)
# #chol.Rinv <- chol(Rinv)
Uchol_Rinv <- Matrix::Matrix(lltA(as(Rinv,"matrix")))
#chol_Rinv <- Matrix::Cholesky(Rinv)
#factors_chol_Rinv <- Matrix::expand(chol_Rinv)
#Lchol_Rinv <- Matrix::t(factors_chol_Rinv$P) %*% factors_chol_Rinv$L
#Uchol_Rinv <- Matrix::t(factors_chol_Rinv$L) %*% factors_chol_Rinv$P
#rm(chol_Rinv,factors_chol_Rinv)
# Comprobación factor Choleski
#Rinv2 <- Matrix::tcrossprod(Lchol_Rinv)
#Rinv3 <- Matrix::crossprod(Uchol_Rinv)
#range(Rinv - Rinv2); range(Rinv - Rinv3)
#Xstar <- chol.Rinv %*% X
#Zstar <- chol.Rinv %*% Z
Xstar <- Uchol_Rinv %*% X
Zstar <- Uchol_Rinv %*% Z
# #Astar <- as.matrix(chol.Rinv)%*%kronecker(A,It)
start <- proc.time()[3]
A_cov1 <- rbind(cbind(Matrix::crossprod(Xstar),
Matrix::t(Xstar) %*% Zstar),
cbind(Matrix::t(Zstar) %*% Xstar,
Matrix::crossprod(Zstar) +
sig2u*Matrix::Matrix(diag(Ginv_eff)) ))
cov1_eff <- try( sig2u*Matrix::solve(A_cov1) )
if (class(cov1_eff) == "try-error") {
cov1_eff <- try( sig2u*solve(as.matrix(A_cov1),tol=1e-40) )
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 <- Matrix::Matrix(cov1_eff)
rownames(cov1_eff) <- colnames(cov1_eff) <- c(names(bfixed),names(brandom))
se_bfixed <- sqrt(diag(as.matrix(cov1_eff[names(bfixed),names(bfixed)])))
names(se_bfixed) <- names(bfixed)
se_brandom <- sqrt(diag(as.matrix(cov1_eff[names(brandom),names(brandom)])))
names(se_brandom) <- names(brandom)
end <- proc.time()[3]
cat("\n Time to compute covariances: ", (end-start), "seconds")
# # Matriz Covarianzas Frequentist tipo Sandwich (Algo menor las varianzas)
# # C.Rinv.C <- (1/sig2u)*rbind(cbind(t(X)%*%kronecker(Diagonal(nsp),Omegainv)%*%X,
# # t(X)%*%kronecker(Diagonal(nsp),Omegainv)%*%Z),
# # cbind(t(Z)%*%kronecker(Diagonal(nsp),Omegainv)%*%X,
# # t(Z)%*%kronecker(Diagonal(nsp),Omegainv)%*%Z))
# C.Rinv.C <- (1/sig2u)*rbind(cbind(Matrix::crossprod(Xstar),
# Matrix::t(Xstar) %*% Zstar),
# cbind(Matrix::t(Zstar) %*% Xstar,
# Matrix::crossprod(Zstar)))
# cov2.eff <- cov1_eff %*% C.Rinv.C %*% cov1_eff
# CALCULAR VALORES ESTIMADOS Y RESIDUOS
fit_Ay <- X %*% bfixed + Z %*% brandom # + offset
fit <- try((Matrix::solve(A) %x% It) %*% fit_Ay)
if (class(fit) == "try-error") {
fit <- (MASS::ginv(as(A,"matrix")) %x% It) %*% eta }
###### Commented Code: same result, much slower
# var_fit_Ay <- cbind(X,Z) %*% cov1_eff %*% Matrix::t(cbind(X,Z))
# var_fit <- try((Matrix::solve(A) %x% It) %*% var_fit_Ay %*%
# (Matrix::t(Matrix::solve(A)) %x% It))
# if (class(var_fit) == "try-error") { (MASS::ginv(as(A,"matrix")) %x% It) %*%
# var_fit_Ay %*% (t(MASS::ginv(as(A,"matrix"))) %x% It) }
# se_fit_Ay <- sqrt(diag(as.matrix(var_fit_Ay)))
# se_fit <- sqrt(diag(as.matrix(var_fit)))
se_fit_Ay <- Matrix::rowSums((cbind(X,Z) %*% cov1_eff) * cbind(X,Z))^0.5
se_fit <- Matrix::rowSums(
(((Matrix::solve(A) %x% It) %*% cbind(X,Z))
%*% cov1_eff) *
((Matrix::solve(A) %x% It) %*% cbind(X,Z)) )^0.5
# resids N(0,sigma_u^2*Omega) and resids.norm N(0,sigma_eps*I)
residuals <- as.vector((A %x% It) %*% y - fit_Ay)
# REPASAR LOS RESIDS_NORM COMPARADO CON RESIDS.NORM
#residuals_norm <- sqrt(1-phi^2)*(chol(as.matrix(Rinv)) %*% residuals)
residuals_norm <- sqrt(1-phi^2)*(Uchol_Rinv %*% residuals)
# Compute AIC y BIC based on loglik functions (Fahrmeir, pp. 664 and 677)
aic <- -2*llik + 2*edftot
bic <- -2*llik + log(length(y))*edftot
res <- list(edfspt = edfspt, edfnopar = edfnopar, edftot = edftot,
tauspt = tauspt, taunopar = taunopar,
psanova = psanova, sar = sar, ar1 = ar1,
fitted.values = as.vector(fit),
fit_Ay = as.vector(fit_Ay),
se_fitted.values = as.vector(se_fit),
se_fit_Ay = as.vector(se_fit_Ay),
residuals = as.vector(residuals),
residuals_norm = as.vector(residuals_norm),
sig2 = sig2u,
rho = rho, phi = phi,
se_rho = se_rho_an, se_phi = se_phi_an,
cov_rho_phi_an = cov_rho_phi_an,
var_rho_phi_num = var_rho_phi_num,
bfixed = bfixed, brandom = brandom,
se_bfixed = se_bfixed, se_brandom = se_brandom,
llik = llik, llik_reml = llik_reml,
aic = aic, bic = bic,
vcov_b = cov1_eff,
sp1 = sp1,sp2 = sp2, time = time)
} # end of function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.