tests/testthat/helper-zoo.R

dimnames_lm_x <- function(dimnames, n_cols_x, intercept) {
  
  if (intercept && (length(dimnames) > 1)) {
    
    result <- list(dimnames[[1]], c("(Intercept)", dimnames[[2]]))
    
  } else if (!intercept && (length(dimnames) > 1)) {
    
    result <- list(dimnames[[1]], dimnames[[2]])
    
  } else if (intercept) {
    
    result <- list(NULL, c("(Intercept)", paste0("x", rep(1:(n_cols_x - 1)))))
    
  } else {
    
    result <- list(NULL, paste0("x", rep(1:n_cols_x)))
    
  }
  
  return(result)
  
}

scale_z <- function(x, center = TRUE, scale = TRUE) {
  
  n_rows_x <- length(x)
  result <- as.numeric(NA)
  
  temp <- scale(x, center = center, scale = scale)
  
  if ((scale && (n_rows_x > 1)) || !scale) {
    result <- temp[length(temp)]
  }
  
  return(result)
  
}

rollapplyr_cube <- function(f, x, y, width) {
  
  if (is.matrix(x) || is.matrix(y)) {
    
    if (!is.matrix(x)) {
      
      temp_attr <- attributes(x)
      x <- as.matrix(zoo::coredata(x))
      attr(x, "dimnames") <- NULL
      attr(x, "index") <- temp_attr[["index"]]
      attr(x, "class") <- temp_attr[["class"]]
      
    }
    
    if (!is.matrix(y)) {
      
      temp_attr <- attributes(y)
      y <- as.matrix(zoo::coredata(y))
      attr(y, "dimnames") <- NULL
      attr(y, "index") <- temp_attr[["index"]]
      attr(y, "class") <- temp_attr[["class"]]
      
    }
    
    n_rows_xy <- nrow(x)
    n_cols_x <- ncol(x)
    n_cols_y <- ncol(y)
    result <- array(as.numeric(NA), c(n_cols_x, n_cols_y, n_rows_xy))
    
    for (i in 1:n_rows_xy) {
      
      result[ , , i] <- f(x[max(1, i - width + 1):i, , drop = FALSE],
                          y[max(1, i - width + 1):i, , drop = FALSE])
      
    }
    
    attr(result, "dim") <- c(n_cols_x, n_cols_y, n_rows_xy)
    
    x_dimnames <- dimnames(x)
    y_dimnames <- dimnames(y)
    attr(result, "dimnames") <- list(x_dimnames[[2]], y_dimnames[[2]], NULL)
    
  } else {
    
    n_rows_xy <- length(x)
    result <- array(as.numeric(NA), n_rows_xy)
    
    for (i in 1:n_rows_xy) {
      
      result[i] <- f(x[max(1, i - width + 1):i],
                     y[max(1, i - width + 1):i])
      
    }
    
    result <- as.numeric(result)
    
  }
  
  return(result)
  
}

crossprod_scale <- function(x, y, center = FALSE, scale = FALSE) {
  
  result <- crossprod(scale(x, center = center, scale = scale),
                      scale(y, center = center, scale = scale))
  
}

rollapplyr_lm <- function(x, y, width, intercept) {
  
  if (!requireNamespace("zoo", quietly = TRUE)) {
    stop("zoo package required for this function")
  }
  
  if (is.matrix(x) || is.matrix(y) || intercept) {
    
    if (!is.matrix(x)) {
      
      temp_attr <- attributes(x)
      x <- as.matrix(zoo::coredata(x))
      attr(x, "dimnames") <- NULL
      attr(x, "index") <- temp_attr[["index"]]
      attr(x, "class") <- temp_attr[["class"]]
      
    }
    
    if (!is.matrix(y)) {
      
      temp_attr <- attributes(y)
      y <- as.matrix(zoo::coredata(y))
      attr(y, "dimnames") <- NULL
      attr(y, "index") <- temp_attr[["index"]]
      attr(y, "class") <- temp_attr[["class"]]
      
    }
    
    n_rows_xy <- nrow(x)
    n_cols_x <- ncol(x)
    
    if (intercept) {
      n_cols_x <- n_cols_x + 1
    }
    
    result <- list("coefficients" = matrix(as.numeric(NA), n_rows_xy, n_cols_x),
                   "r.squared" = matrix(as.numeric(NA), n_rows_xy, 1),
                   "std.error" = matrix(as.numeric(NA), n_rows_xy, n_cols_x))
    
    if (zoo::is.zoo(x)) {
      
      x_attr <- attributes(x)
      x_attr[["dim"]] <- NULL
      x_attr[["dimnames"]] <- NULL
      
    } else if (zoo::is.zoo(y)) {
      
      x_attr <- attributes(y)
      x_attr[["dim"]] <- NULL
      x_attr[["dimnames"]] <- NULL
      
    }
    
    for (i in 1:n_rows_xy) {
      
      x_subset <- x[max(1, i - width + 1):i, , drop = FALSE]
      y_subset <- y[max(1, i - width + 1):i, , drop = FALSE]
      data <- as.data.frame(cbind(y_subset, x_subset))
      
      if (intercept) {
        
        # "Unparseable 'response'; use is deprecated.  Use as.name(.) or `..`!"
        fit <- lm(reformulate(termlabels = ".", response = as.name(names(data)[1])), data = data)
        
      } else {
        fit <- lm(reformulate(termlabels = ".-1", response = as.name(names(data)[1])), data = data)
      }
      
      summary_fit <- summary(fit)
      summary_fit_coef <- coef(summary_fit)[ , "Estimate"]
      
      if (nrow(coef(summary_fit)) == n_cols_x) {
        
        result[["coefficients"]][i, ] <- summary_fit_coef
        
        # "In summary.lm(fit) : essentially perfect fit: summary may be unreliable"
        if (!(isTRUE(all.equal(as.numeric(rep(summary_fit_coef[1], length(y_subset))), as.numeric(y_subset))) &&
              isTRUE(all.equal(as.numeric(summary_fit_coef[-1]), rep(0, length(summary_fit_coef[-1])))))) {
          
          result[["r.squared"]][i, ] <- summary_fit$r.squared
          result[["std.error"]][i, ] <- coef(summary_fit)[ , "Std. Error"]
          
        } else {
          
          result[["r.squared"]][i, ] <- as.numeric(NA)
          result[["std.error"]][i, ] <- as.numeric(NA)
          
        }
        
      }
      
    }
    
    if (exists("x_attr")) {
      
      attributes(result[["coefficients"]]) <- x_attr
      attributes(result[["r.squared"]]) <- x_attr
      attributes(result[["std.error"]]) <- x_attr
      
    }
    
    attr(result[["coefficients"]], "dim") <- c(n_rows_xy, n_cols_x)
    attr(result[["r.squared"]], "dim") <- c(n_rows_xy, 1)
    attr(result[["std.error"]], "dim") <- c(n_rows_xy, n_cols_x)
    
    x_dimnames <- dimnames(x)
    y_dimnames <- dimnames(y)
    
    attr(result[["coefficients"]], "dimnames") <- dimnames_lm_x(x_dimnames, n_cols_x, intercept)
    if (length(x_dimnames) > 1) {
      attr(result[["r.squared"]], "dimnames") <- list(x_dimnames[[1]], "R-squared")
    } else {
      attr(result[["r.squared"]], "dimnames") <- list(NULL, "R-squared")
    }
    attr(result[["std.error"]], "dimnames") <- dimnames_lm_x(x_dimnames, n_cols_x, intercept)
    
  } else {
    
    n_rows_xy <- length(x)
    n_cols_x <- 1
    
    result <- list("coefficients" = rep(as.numeric(NA), n_rows_xy),
                   "r.squared" = rep(as.numeric(NA), n_rows_xy),
                   "std.error" = rep(as.numeric(NA), n_rows_xy))
    
    if (zoo::is.zoo(x)) {

      x_attr <- attributes(x)
      x_attr[["dim"]] <- NULL
      x_attr[["dimnames"]] <- NULL

    } else if (zoo::is.zoo(y)) {

      x_attr <- attributes(y)
      x_attr[["dim"]] <- NULL
      x_attr[["dimnames"]] <- NULL

    }
    
    for (i in 1:n_rows_xy) {
      
      x_subset <- x[max(1, i - width + 1):i]
      y_subset <- y[max(1, i - width + 1):i]
      data <- as.data.frame(cbind(y_subset, x_subset))
      
      if (intercept) {
        fit <- lm(reformulate(termlabels = ".", response = names(data)[1]), data = data)
      } else {
        fit <- lm(reformulate(termlabels = ".-1", response = names(data)[1]), data = data)
      }
      
      summary_fit <- summary(fit)
      summary_fit_coef <- coef(summary_fit)[ , "Estimate"]
      
      if (nrow(coef(summary_fit)) == n_cols_x) {
        
        result[["coefficients"]][i] <- summary_fit_coef
        
        # "In summary.lm(fit) : essentially perfect fit: summary may be unreliable"
        if (!(isTRUE(all.equal(as.numeric(rep(summary_fit_coef[1], length(y_subset))), as.numeric(y_subset))) &&
              isTRUE(all.equal(as.numeric(summary_fit_coef[-1]), rep(0, length(summary_fit_coef[-1])))))) {
          
          result[["r.squared"]][i] <- summary_fit$r.squared
          result[["std.error"]][i] <- coef(summary_fit)[ , "Std. Error"]
          
        } else {
          
          result[["r.squared"]][i] <- as.numeric(NA)
          result[["std.error"]][i] <- as.numeric(NA)
          
        }
        
      }
      
    }
    
    if (exists("x_attr")) {
      
      attributes(result[["coefficients"]]) <- x_attr
      attributes(result[["r.squared"]]) <- x_attr
      attributes(result[["std.error"]]) <- x_attr
      
    }
    
  }
  
  return(result)
  
}

Try the roll package in your browser

Any scripts or data that you put into this service are public.

roll documentation built on May 29, 2024, 6:02 a.m.