bwSelPredict <- function(data,
type,
level,
obj,
test,
modeltype,
# k = NULL,
# consec = NULL,
... ) # is lags passed on?
{
# -------------------- Input checks & calc Aux Vars -------------------
p <- ncol(data)
# Fetch arguments
args <- list(...)
# if(is.null(args$k)) args$k <- 2
args$d <- args$k - 1
n_test <- length(test)
consec <- args$consec
# Scale Gaussians
ind_Gauss <- which(type == 'g')
for(i in ind_Gauss) data[, i] <- scale(data[, i])
data_df <- as.data.frame(data)
# Categoricals into factors (Needed to use formula to construct design matrix)
for(i in which(type=='c')) data_df[, i] <- as.factor(data_df[, i])
# -------------------- Get VAR data structure -------------------
if(modeltype == 'mvar') {
# Glue all lagged variables together
n_lags <- length(args$lags)
# Prepare Data (already cuts max(lags) first observations to compute design matrix)
data_lagged <- lagData(data = data,
lags = args$lags,
consec = consec)
# out of list
data_response <- data_lagged$data_response
l_data_lags <- data_lagged$l_data_lags
data_response <- apply(data_response, 2, as.numeric) # to avoid confusion with labels of categories if there are factors
# Detele cases
data_response <- data_response[data_lagged$included, ]
l_data_lags <- lapply(l_data_lags, function(z) z <- z[data_lagged$included, ])
data_response <- apply(data_response, 2, as.numeric) # to avoid confusion with labels of categories if there are factors
}
# -------------------- Compute Predictions -------------------
m_pred <- matrix(NA, nrow = n_test, ncol = p)
for(v in 1:p) {
# ----- Create Design Matrices -----
if(modeltype == 'mvar') {
# append response with predictors
y <- data_response[,v] # response variable v
data_v <- cbind(do.call(cbind, l_data_lags)) # combine
data_v <- as.data.frame(data_v) # because model.matrix() below requries a data.frame
# Dummy coding
form <- as.formula('y ~ (.)')
# Construct standard design matrix (to get number of parameters for tau threshold)
X_standard <- model.matrix(form, data = data_v)[, -1] # delete intercept (added by glmnet later)
if(args$overparameterize) {
# # Compute augmented type and level vectors
type_aug <- rep(type, n_lags)
level_aug <- rep(level, n_lags)
# Construct over-parameterized design matrix
X_over <- ModelMatrix(data = data_v,
type = type_aug,
level = level_aug,
labels = colnames(data_v),
d = 1,
v = NULL)
X <- X_over
} else {
X <- X_standard
}
}
if(modeltype == 'mgm') {
# Create HOI design matrix (copied from mgm)
if(args$d > (p - 1)) {
stop("Order of interactions cannot be larger than the number of predictors.")
} else if (args$d == 1){ form <- as.formula(paste(colnames(data_df)[v],"~ (.)"))
} else { form <- as.formula(paste(colnames(data_df)[v],"~ (.)^", args$d)) }
# Construct standard design matrix (to get number of parameters for tau threshold)
X_standard <- model.matrix(form, data = data_df)[, -1] # delete intercept (added by glmnet later)
if(args$overparameterize) {
# Construct over-parameterized design matrix
X_over <- ModelMatrix(data = data,
type = type,
level = level,
labels = colnames(data),
d = args$d,
v = v) # defined above from k
X <- X_over
} else {
X <- X_standard
}
y <- data[, v]
} # end if: modeltype mgm
# ----- Make Predictions for test points -----
for(i in 1:n_test) {
coefs <- obj$tvmodels[[i]]$nodemodels[[v]]$model
row <- test[i]
if(type[v] == 'c') {
n_cat <- level[v]
# Compute exp(potentials)
Potentials <- rep(NA, n_cat)
for(k in 1:n_cat) Potentials[k] <- exp(coefs[[k]][1] + X[row,] %*% coefs[[k]][-1])
# compute category-probabilities
Probabilities <- Potentials / sum(Potentials)
# Classify
m_pred[i, v] <- sort(unique(y))[which.max(Probabilities)]
} else {
## Prediction Continuous
coefs <- as.numeric(coefs) # get coefficients
m_pred[i, v] <- coefs[1] + X[row,] %*% coefs[-1] # predict
}
} # end for: i (over test locations)
} # end for: v
# -------------------- Compute Errors -------------------
if(modeltype == 'mvar') {
if(n_test==1) {
m_true <- matrix(data_response[test, ], nrow = n_test) # for the case n_test=1
} else {
m_true <- data_response[test, ]
}
}
if(modeltype == 'mgm') {
if(n_test==1) {
m_true <- matrix(data[test, ], nrow = n_test) # for the case n_test=1
} else {
m_true <- data[test, ]
}
}
m_error <- matrix(NA, nrow = n_test, ncol = p)
# 0/1 loss for categorical, RMSE for continuous (which here is just absolute error, because 1 element)
for(i in 1:p) {
if(type[v] == 'c') {
m_error[, i] <- abs(m_true[, i] - m_pred[, i])
} else {
m_error[, i] <- abs(m_true[, i] - m_pred[, i]) # RMSE for 1 case = abs error
}
}
# -------------------- Output -------------------
predOutlist <- list()
predOutlist$errors <- m_error
predOutlist$error_time <- rowMeans(m_error)
predOutlist$error_mean <- mean(m_error)
return(predOutlist)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.