R/loglikLOOCVVAR1fused.r

In ragt2ridges: Ridge Estimation of Vector Auto-Regressive (VAR) Processes

loglikLOOCVVAR1fused <- function(lambdas, 
                                 Y, 
                                 id, 
                                 unbalanced=matrix(nrow=0, ncol=2), 
                                 ...){

	########################################################################
	# 
	# DESCRIPTION:
	# Evaluation of the (minus) leave-one-out cross-validated log-likelihood 
	# of the VAR(1) model for given choices of the ridge penalty parameters 
	# (lambdaA and lambdaO). The functions also works with a (possibly) 
	# unbalanced experimental set-up. The VAR(1)-process is assumed to have 
	# mean zero.
	#
	# ARGUMENTS:
	# -> lambdas      : Numeric of length two. It contains the ridge penalty 
	#                   parameters to be used in the estimation of A and the 
	#                   Omega, the precision matrix of the errors (also 
	#                   called innovations). 
	# -> Y            : Three-dimensional array containing the data. The 
	#                   first, second and third dimensions correspond to 
	#                   covariates, time and samples, respectively. The data 
	#                   are assumed to centered covariate-wise. 
	# -> unbalanced   : A matrix with two columns, indicating the unbalances 
	#                   in the design. Each row represents a missing design 
	#                   point in the (time x individual)-layout. The first 
	#                   and second column indicate the time and individual 
	#                   (respectively) specifics of the missing design point.
	# -> ...          : Other arguments to be passed to ridgeVAR1.
	#
	# DEPENDENCIES:
	# require("rags2ridges")          # functions from package : ridgeVAR1
	#
	# NOTES:
	# ...
	#
	########################################################################

	# input checks
	if (!is(Y, "array")){ 
		stop("Input (Y) is of wrong class.") 
	}
	if (length(dim(Y)) != 3){ 
		stop("Input (Y) is of wrong dimensions: either covariate, time or sample dimension is missing.") 
	}
	if (!is(id, "numeric") & !is(id, "integer")){ 
		stop("Input (id) is of wrong class.") 
	}
	if (length(id) != dim(Y)[3]){ 
		stop("Input (id) is of wrong length: should equal sample dimension of Y.") 
	}		
	if (!is(lambdas, "numeric")){ 
		stop("Input (lambdas) is of wrong class.") 
	}
	if (length(lambdas) != 3){ 
		stop("Input (lambdas) is of wrong length.") 
	}
	if (any(is.na(lambdas))){ 
		stop("Input (lambdas) is not a vector of non-negative numbers.") 
	}
	if (any(lambdas < 0)){ 
		stop("Input (lambdas) is not a vector of non-negative numbers.") 
	}
	if (!is(unbalanced, "matrix")){ 
		stop("Input (unbalanced) is of wrong class.") 
	}    
	if (ncol(unbalanced) != 2){ 
		stop("Wrong dimensions of the matrix unbalanced.") 
	} 

	# determine leave-one-out scheme
	LOOscheme <- cbind(rep(2:dim(Y)[2], dim(Y)[3]), 
	                   sort(rep(1:dim(Y)[3], dim(Y)[2]-1)))	
	if (nrow(unbalanced) > 0){
		LOO2unbalanced <- numeric()
		for (k in 1:nrow(unbalanced)){
			LOO2unbalanced <- c(LOO2unbalanced, 
			                    which(apply(LOOscheme, 1, function(Y, Z){ all(Y == Z) }, Z=unbalanced[k,])))
		}
		LOOscheme <- LOOscheme[-LOO2unbalanced,]
	}

	loglik <- 0
	for (k in 1:nrow(LOOscheme)){
		# evaluate LOOCV estimates of VAR(1) parameters A and Se
		VAR1hats <- ridgeVAR1fused(Y, 
		                           id, 
		                           lambdas[1], 
		                           lambdas[2], 
		                           lambdas[3], 
		                           unbalanced=rbind(unbalanced, LOOscheme[k,,drop=FALSE]), 
		                           ...)
  
		# evaluate LOOCV loglikelihood
		idGroup <- id[LOOscheme[k,2]]
		loglik  <- loglik + 
			   .armaVAR1_loglik_LOOCVinternal(Y[,LOOscheme[k,1]  , LOOscheme[k,2]], 
 				                          Y[,LOOscheme[k,1]-1, LOOscheme[k,2]], 
	                                                  VAR1hats$As[c((nrow(Y)*idGroup+1):(nrow(Y)*(idGroup+1))),], 
						          VAR1hats$P)
	}
	# return minus LOOCV loglikelihood
	return(-loglik)
}