a number of changes to robust_isatpanel and some smaller changes to t…

…esting, vignette, and making t.pval explicit in isatpanel

NEWS.md

@@ -2,21 +2,21 @@
 
 ## Major changes
 
--   Implemented the argument `uis` in `isatpanel()` to enable user-specified indicator saturation
+- Implemented the argument `uis` in `isatpanel()` to enable user-specified indicator saturation
+- Added a new argument `t.pval` to mirror the setting of a target p-value to the gets::isat() function.
 - Added new function `get_indicators()`
 
-
 ## Minor changes
 
 -   Set the default argument for plotting in `isatpanel` to `TRUE` to enable automatic plotting
 -   Changed the default plot in `plot.isatpanel()` to now be a combination of time series and `plot_grid()` using the `cowplot` package
-
+-   Automated testing for `robust_isatpanel()` using the Koch et al. 2022 results
 
 ## Bug fixes
 
--   Bugfix in robust_isatpanel()
+-   Bugfix in robust_isatpanel(). Before it appeared like HAC S.E. with clusters were possible. This has been fixed. Now also the estimation type is taken from the inital model rather than specified as a separate argument; hence 'effect' has been removed as an argument.  
 -   Bugfixes in plotting to work with uis specifications
--   Improved Plotting when there is just a single identified break (now the colour scale adjusts)
+-   Improved plotting when there is just a single identified break (now the colour scale adjusts)
 
 # getspanel 0.1.4
 

R/isatpanel.R

@@ -28,6 +28,7 @@
 #' @param csis Logical. Use Coefficient Step Indicator Saturation. Constructed by Default is FALSE.
 #' @param cfesis Logical. Use Coefficient-Fixed Effect Indicator Saturation. Default is FALSE.
 #' @param uis Matrix or List. This can be used to include a set of UIS (User Specified Indicators). Must be equal to the sample size (so it is recommended to use this only with datasets without \code{NA} values. Default is \code{NULL}. See the reference by Genaro Sucarrat (2020) below for an explanation of the UIS system.
+#' @param t.pval numeric value between 0 and 1. The significance level used for the two-sided regressor significance t-tests
 #' @param ... Further arguments to [gets::isat()]
 #' @param data The input data.frame object.
 #' @param formula Please specify a formula argument. The dependent variable will be the left-most element, separated by a ~ symbol from the remaining regressors. Note the intercept will always be removed, if effect is not "none" - this means that if any fixed effects are specified, the intercept will always be removed.
@@ -104,6 +105,7 @@ isatpanel <- function(
     cfesis_var = colnames(mxreg),
     cfesis_id = NULL,
     uis = NULL,
+    t.pval = 0.001,
 
     plot = TRUE,
     plm_model = "within",
@@ -551,7 +553,7 @@ isatpanel <- function(
 
   options(mc.warning = FALSE)
   #sis <- FALSE # don't allow sis argument - does not make sense in a panel context, only JSIS makes sense
-  ispan <- gets::isat(y, mxreg = mx, iis = iis, sis = FALSE, uis = sispanx, user.estimator = user.estimator, mc = FALSE, ...)
+  ispan <- gets::isat(y, mxreg = mx, iis = iis, sis = FALSE, uis = sispanx, user.estimator = user.estimator, mc = FALSE, t.pval = t.pval, ...)
   #ispan <- isat.short(y, mxreg = mx, iis=iis, sis=FALSE, uis=sispanx, user.estimator = user.estimator, mc=FALSE, ...)
 
   ###############################
@@ -561,13 +563,12 @@ isatpanel <- function(
 
   out$isatpanel.result <- ispan
 
-
   # Create a final data object
   indicators <- out$isatpanel.result$aux$mX
   if (is.null(engine)) {
-    indicators <- indicators[,!colnames(indicators) %in% c(names(estimateddata),if (exists("dummies")) {names(dummies)}else{NULL})]
+    indicators <- indicators[,!colnames(indicators) %in% c(names(estimateddata),if (exists("dummies")) {names(dummies)}else{NULL}), drop = FALSE]
   } else {
-    indicators <- indicators[,!colnames(indicators) %in% names(estimateddata)]
+    indicators <- indicators[,!colnames(indicators) %in% names(estimateddata), drop = FALSE]
   }
 
 

R/robust_isatpanel.R

@@ -1,12 +1,11 @@
 #' Get robust Standard Errors for the isatpanel result
 #'
 #' @param object An isatpanel object
-#' @param robust Logical (TRUE or FALSE). Should the Standard Errors be robustified for Heterogeneity?
-#' @param HAC Should Heteroscedasticity and Autocorrelation Robust Standard Errors be used?
-#' @param lag Number of Lags to be used with HAC in coeftest. Cannot be specified when HAC = FALSE.
-#' @param effect The effects introduced into the plm model, one of "individual", "time", "twoways" (default), or "nested"
-#' @param type Type of Robust procedure e.g. HC0 for White SE or HC3 for Lang
-#' @param cluster 'group' or 'time' or FALSE
+#' @param robust Logical (TRUE or FALSE). Should the Standard Errors be robustified for Heterogeneity? This uses [plm::vcovHC] with the specified type (default is "HC0").
+#' @param HAC Should Heteroscedasticity and Autocorrelation Robust Standard Errors be used? This uses [plm::vcovNW], which uses the Newey-West estimator.
+#' @param lag Maximum Number of Lags to be used with [plm::vcovNW] using the Newey-West estimator. Cannot be specified when HAC = FALSE. Default is \code{NULL}.
+#' @param type Character string. Type of Robust procedure e.g. 'HC0' for White SE or 'HC3' for Lang.
+#' @param cluster Should an object with clustered S.E. be included? Choose between 'group' or 'time' or FALSE. Uses [plm::vcovHC] with the cluster argument.
 #'
 #' @return A list with robust estimates
 #' @export
@@ -36,8 +35,6 @@
 #'   plot = FALSE,
 #'   t.pval = 0.01
 #' )
-#' plot(result)
-#' plot_grid(result)
 #' robust_isatpanel(result)
 #'}
 #'
@@ -48,38 +45,70 @@ robust_isatpanel <- function(object,
                              effect = "twoways",
                              type = "HC0",
                              cluster = "group"){
-  if(!is.logical(cluster)){docluster <- TRUE}
-  if(cluster == TRUE){stop("Please specify the cluster setting i.e. either 'group' or 'time'. TRUE is not allowed. To disable clustering use FALSE. ")}
 
+  # Checks
+  if(cluster == TRUE){stop("Please specify the cluster setting i.e. either 'group' or 'time'. TRUE is not allowed. To disable clustering use FALSE.")}
   if(!is.null(lag)&!HAC){stop("You cannot specify a lag when you select HAC = FALSE.")}
-  out <- list()
 
+
+  # Set-up
   df <- object$finaldata
+  df$id <- as.factor(df$id)
+  df$time <- as.factor(df$time)
+  effect <- object$arguments$effect
   pdata <- plm::pdata.frame(df,index = c("id","time"))
 
+  # parse FE
+  parse_FE <- if (effect == "twoways") {
+    "| id + time"
+  } else if (effect == "time")  {
+    " + time"
+  } else if (effect == "individual") {
+    " + id"
+  } else if (effect == "none") {
+    ""
+  } else {
+    stop("No Fixed Effect Specification was selected. Choose from 'twoways', 'time', 'individual', or 'none'")
+  }
+
   # Let's look at the plm output
   formula <- as.formula(paste0("y ~ ",paste0(names(pdata[,!names(pdata) %in% c("y","id","time")]),collapse = " + ")))
-  plm_object <- plm::plm(formula = formula,data = pdata, effect = effect,model = "within")
+  plm_object <- plm::plm(formula = formula, data = pdata, effect = effect, model = "within")
 
-  out$plm_object <- plm_object
 
-  #summary(plm_object)
-  # This should be exactly the same output as the one from isatpanel
+  out <- list()
+  out$plm_object <- plm_object
 
+  out$robust <- lmtest::coeftest(plm_object, vcov = plm::vcovHC(plm_object, type = type))
 
-  if(robust){
-    if(docluster){
-      out$robust <- lmtest::coeftest(plm_object, vcov=sandwich::vcovHC(plm_object,type=type,cluster=cluster))
-    }else{
-      out$robust <- lmtest::coeftest(plm_object, vcov=sandwich::vcovHC(plm_object,type=type))
-    }
+  if(cluster != FALSE){
+    out$cluster <- lmtest::coeftest(plm_object, vcov = plm::vcovHC(plm_object, type = type, cluster = cluster))
   }
 
   if(HAC){
-    lm_mod <- lm(y ~ as.factor(time) + .-1,df)
-    out$HAC <- lmtest::coeftest(lm_mod, vcov=sandwich::vcovHAC(lm_mod, cluster=cluster, lag = lag))
+    out$HAC <- lmtest::coeftest(plm_object, vcov = plm::vcovNW(plm_object, type = type, maxlag = lag))
   }
 
+
+
+
+
+  #summary(plm_object)
+  # This should be exactly the same output as the one from isatpanel
+
+
+  # if(robust){
+  #   if(docluster){
+  #     out$robust <- lmtest::coeftest(plm_object, vcov=sandwich::vcovHC(plm_object,type=type,cluster=cluster))
+  #   }else{
+  #     out$robust <- lmtest::coeftest(plm_object, vcov=sandwich::vcovHC(plm_object,type=type))
+  #   }
+  # }
+  #
+  # if(HAC){
+  #
+  # }
+
   #
   # # The following also for autocorrelation
   # # Heteroskedasticity- and autocorrelation-consistent (HAC)