Skip to content

Commit

Permalink
fix bug bootstrap scheme for tvvar
Browse files Browse the repository at this point in the history
  • Loading branch information
jmbh committed Mar 26, 2019
1 parent 4183dfd commit f93137a
Show file tree
Hide file tree
Showing 4 changed files with 22 additions and 15 deletions.
7 changes: 5 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: mgm
Type: Package
Title: Estimating Time-Varying k-Order Mixed Graphical Models
Version: 1.2-6
Date: 2019-07-01
Version: 1.2-7
Date: 2019-26-03
Author: Jonas Haslbeck
Maintainer: Jonas Haslbeck <jonashaslbeck@gmail.com>
Description: Estimation of k-Order time-varying Mixed Graphical Models and mixed VAR(p) models via elastic-net regularized neighborhood regression. For details see linked paper.
Expand All @@ -12,3 +12,6 @@ License: GPL (>= 2)
Imports: matrixcalc, glmnet, stringr, Hmisc, qgraph, gtools
LazyData: true
Packaged: 2019-07-01 00:00:01 CET; jo



3 changes: 3 additions & 0 deletions NEWS
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
Changes in Version 1.2-7
o Fixed bug in bootstrap scheme for time-varying VAR models

Changes in Version 1.2-6
o Fix aliasing of internal functions

Expand Down
9 changes: 5 additions & 4 deletions R/mvar.R
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,8 @@ mvar <- function(data, # n x p data matrix
ind_Gauss <- which(type == 'g')
if(scale) for(i in ind_Gauss) data[, i] <- scale(data[, i])



# ----- Basic Input Checks -----

# Checks on Data
Expand Down Expand Up @@ -224,6 +226,8 @@ mvar <- function(data, # n x p data matrix

# ----- Use bootstrap instead of original data (called from resample()) -----

n_lags <- length(lags)

if(!is.null(args$bootstrap)) {
if(args$bootstrap) {

Expand All @@ -238,7 +242,7 @@ mvar <- function(data, # n x p data matrix
# overwrite data with bootstrap sample, passed on from resample()
data_response <- data_response[args$boot_ind, ] # args$boot_ind is defined such that it only selects valid rows
l_data_lags <- lapply(l_data_lags, function(x) x[args$boot_ind, ])

# overwrite weights
weights_design <- weights[args$boot_ind]

Expand All @@ -249,9 +253,6 @@ mvar <- function(data, # n x p data matrix

# This checks glmnet requirements wrt categorical predictors, for each lag set


n_lags <- length(lags)

for(lag in 1:n_lags) {

glmnetRequirements(data = l_data_lags[[lag]],
Expand Down
18 changes: 9 additions & 9 deletions R/resample.R
Original file line number Diff line number Diff line change
Expand Up @@ -320,21 +320,19 @@ resample <- function(object, # one of the four mgm model objects (mgm, mvar, tvm
lags = o_call$lags,
consec = o_call$consec)


# add false for excluded 1:max_lags time points at beginning of time series
# max_lags <- max(o_call$lags)
# full_included_vec <- c(rep(FALSE, max_lags), data_lagged$included)

# Make use of time-vector to be able to have block-bootstrap on actual time scale
timepoints_design <- o_call$timepoints[data_lagged$included]

# Break data into blocks of equal time-duration (unsing timepoints vector)
Qt <- quantile(timepoints_design, probs = seq(0, 1, length = blocks + 1))
ind_blocks <- cut(x = timepoints_design, # important: indices in the design matrix
breaks = Qt,
labels = FALSE)

ind_blocks[1] <- 1 # for some reason the first entry is NA

# Rows included in the design matrix
ind_valid_rows <- (1:nrow(data))[data_lagged$included]

# Storage
l_ind <- list()

Expand All @@ -346,8 +344,8 @@ resample <- function(object, # one of the four mgm model objects (mgm, mvar, tvm

for(b2 in 1:blocks) {

# Take bootstrap samples within blocks
block_b2 <- which(ind_blocks == b2)
# Take bootstrap samples within blocks defined above from rows that are included (ind_valid_rows)
block_b2 <- ind_valid_rows[which(ind_blocks == b2)]
l_ind_blocks[[b2]] <- sample(x = block_b2,
size = length(block_b2),
replace = TRUE)
Expand All @@ -367,6 +365,8 @@ resample <- function(object, # one of the four mgm model objects (mgm, mvar, tvm

for(b in 1:nB) {

# browser()

set.seed(seeds[b])
if(verbatim) print(paste0("Seed = ", seeds[b]))
tt <- proc.time()[3]
Expand Down

0 comments on commit f93137a

Please sign in to comment.