# R/asym.R In panelr: Regression Models and Utilities for Repeated Measures and Panel Data

#### Documented in asymtidy.asym

```#' @title Estimate asymmetric effects models using first differences
#' @description The function fits the asymmetric effects first difference
#'  model described in Allison (2019) using GLS estimation.
#' @inheritParams wbm
#' @param variance One of `"toeplitz-1"`, `"constrained"`, or `"unconstrained"`.
#'   The toeplitz variance specification estimates a single error variance and
#'   a single lag-1 error correlation with other lags having zero correlation.
#'   The constrained model assumes no autocorrelated errors or heteroskedastic
#'   errors. The unconstrained option allows separate variances for every
#'   period as well as every lag of autocorrelation. This can be very
#'   computationally taxing as periods increase and will be inefficient when
#'   not necessary. See Allison (2019) for more.
#' @param error.type Either "CR2" or "CR1S". See the `clubSandwich` package for
#'  more details.
#' @param ... Ignored.
#' @references
#' Allison, P. D. (2019). Asymmetric fixed-effects models for panel data.
#' *Socius*, *5*, 1-12. https://doi.org/10.1177/2378023119826441
#' @examples
#'
#' \dontrun{
#' data("teen_poverty")
#' # Convert to long format
#' teen <- long_panel(teen_poverty, begin = 1, end = 5)
#' model <- asym(hours ~ lag(pov) + spouse, data = teen)
#' summary(model)
#' }
#'
#' @rdname asym
#' @importFrom stats BIC AIC na.omit
#' @export
asym <- function(formula, data, id = NULL, wave = NULL, use.wave = FALSE,
min.waves = 1,
variance = c("toeplitz-1", "constrained", "unconstrained"),
error.type = c("CR2", "CR1S"), ...) {

if (!requireNamespace("nlme")) need_package("nlme")
if (!requireNamespace("clubSandwich")) need_package("clubSandwich")
if (!requireNamespace("car")) need_package("car")

the_call <- match.call()
the_env <- parent.frame()

formula <- Formula::Formula(formula)

# Send to helper function for data prep
prepped <- diff_data(formula = formula, data = data, id = id,
wave = wave,  min.waves = min.waves,
weights = NULL, use.wave = use.wave,
asym = TRUE)

e <- prepped\$e
pf <- prepped\$pf
data <- e\$data
wave <- prepped\$wave
id <- prepped\$id
dv <- prepped\$dv
weights <- prepped\$weights

# Use helper function to generate formula to pass to lme4
fin_formula <- as.formula(e\$fin_formula)
if (!is.null(pf\$constants)) {
constants <- paste(pf\$constants, collapse = " - ")
up_form <- as.formula(paste(". ~ . -", constants))
fin_formula <- update(fin_formula, up_form)
}

if (variance[1] == "unconstrained") {
cor_form <- as.formula(paste("~ 1 |", id))
var_form <- as.formula(paste("~ 1 |", wave))
corr <- nlme::corSymm(form = cor_form)       # unstructured covariance
weights <- nlme::varIdent(form = var_form)
} else if (variance[1] == "constrained") {
cor_form <- as.formula(paste(" ~ 1 |", id))
corr <- nlme::corARMA(value = -.9999, form = cor_form, p = 0, q = 1,
fixed = TRUE)
weights <- NULL
} else if (variance[1] == "toeplitz-1") {
cor_form <- as.formula(paste(" ~ 1 |", id))
corr <- nlme::corARMA(form = cor_form, p = 0, q = 1)
weights <- NULL
}

gls_mod <- nlme::gls(fin_formula, data = data,
na.action = na.omit,
correlation = corr, weights = weights)

gls_mod\$call\$model <- substitute(fin_formula)

the_vcov <- vcov_CR(gls_mod, cluster = data[[id]], type = error.type[1],
data = data)
coef_table <- as.data.frame(clubSandwich::coef_test(gls_mod, vcov = the_vcov,
test = "naive-t", cluster = data[[id]]))
coef_table <- as.data.frame(coef_table[c("beta","SE","tstat","p_t")])
if ("tstat" %nin% names(coef_table)) { # old version of clubSandwich
names(coef_table) <- c("estimate", "std.error", "p.value")
coef_table["statistic"] <- coef_table\$estimate / coef_table\$std.error
} else {
names(coef_table) <- c("estimate", "std.error", "statistic", "p.value")
coef_table["term"] <- rownames(coef_table)
}
coef_table\$term <- gsub("plus__", "\\+", coef_table\$term)
coef_table\$term <- gsub("minus__", "\\-", coef_table\$term)

if (length(e\$asym_list) > 0) {
asym_diffs <- test_asyms(gls_mod, e\$asym_list, the_vcov)
for (var in names(e\$asym_list)) {
if (stringr::str_detect(var, "^lag_.*_\$")) {
var <- stringr::str_replace(var, "^lag_(.*)_\$", "\\1")
}
rownames(coef_table) <- stringr::str_replace(
rownames(coef_table), paste0("lag_", var, "_"), paste0("lag(", var, ")")
)
asym_diffs\$term <- stringr::str_replace(
asym_diffs\$term, paste0("lag_", var, "_"), paste0("lag(", var, ")")
)
}
} else {
asym_diffs <- NULL
}

mod_info <- list(dv = dv, min.wave = prepped\$minwave,
max.wave = prepped\$maxwave,
num_distinct = prepped\$num_distinct,
AIC = AIC(gls_mod), BIC = BIC(gls_mod),
variance = variance[1],
errors = paste0("Cluster-robust (", error.type[1], ")"))
gls_mod\$mod_info <- mod_info
gls_mod\$coef_table <- coef_table
gls_mod\$asym_diffs <- asym_diffs
colnames(the_vcov) <- rownames(coef_table)
rownames(the_vcov) <- rownames(coef_table)
gls_mod\$vcov <- the_vcov

# Make it a fdm object
class(gls_mod) <- c("asym", "fdm", "gls")
gls_mod
}

#' @export
summary.asym <- function(object, ...) {

dots <- list(...)
if ("digits" %in% names(dots)) {
digits <- dots\$digits
} else {
digits <- getOption("jtools-digits", 2)
}

x <- object

mod_fit <- paste0(bold("MODEL FIT:\n"),
italic("AIC = "), round(x\$mod_info\$AIC, digits),
italic(", BIC = "), round(x\$mod_info\$BIC, digits),  "\n")

if (x\$mod_info\$variance == "toeplitz-1") {
corstruct <- object\$modelStruct\$corStruct
class(corstruct) <- "corARMA"
theta <- coef(corstruct, unconstrained = FALSE)
variance <- paste0(italic("Variance structure: "), x\$mod_info\$variance,
" (theta = ", round(theta, digits), ")")
} else {
variance <- paste0(italic("Variance structure: "), x\$mod_info\$variance)
}

mod_info <- paste0(bold("MODEL INFO:\n"),
italic("Entities: "), x\$mod_info\$num_distinct, "\n",
italic("Time periods: "), paste0(x\$mod_info\$min.wave, "-",
x\$mod_info\$max.wave), "\n",
italic("Dependent variable: "), x\$mod_info\$dv, "\n",
italic("Model type: "),
"Linear asymmetric effects (first differences)\n",
variance
)

mod_info_list <- list(N = x\$mod_info\$num_distinct,
min_wave = x\$mod_info\$min.wave,
max_wave = x\$mod_info\$max.wave,
variance = x\$mod_info\$variance, AIC = x\$mod_info\$AIC,
BIC = x\$mod_info\$BIC)

coef_table <- x\$coef_table
names(coef_table) <- sapply(names(coef_table), function(x) {switch(x,
"estimate" = "Est.",
"std.error" = "S.E.",
"p.value" = "p",
"statistic" = "t val.",
x
)})
rownames(coef_table) <- coef_table\$term
# rownames(coef_table) <- coef_table\$term
coef_table <- coef_table[c("Est.", "S.E.", "t val.", "p")]

if (!is.null(x\$asym_diffs)) {
names(x\$asym_diffs) <- c("variable", "chi^2", "p")
}

out <- list(mod_info = mod_info, coef_table = coef_table, digits = digits,
asym_diffs = x\$asym_diffs, errors = x\$mod_info\$errors,
mod_info_list = mod_info_list)
class(out) <- "summary.asym"
out
}

#' @export
print.summary.asym <- function(x, ...) {
cat(x\$mod_info, "\n\n")

cat(italic("Standard errors:"), x\$errors, "\n")
print(md_table(x\$coef_table, digits = x\$digits, sig.digits = FALSE,
format = getOption("panelr.table.format", "multiline")))

if (!is.null(x\$asym_diffs)) {
cat("\n", bold("Tests of asymmetric effects:\n"), sep = "")
rownames(x\$asym_diffs) <- x\$asym_diffs\$variable
x\$asym_diffs <- x\$asym_diffs %not% "variable"
print(md_table(x\$asym_diffs, digits = x\$digits, sig.digits = FALSE,
format = getOption("panelr.table.format", "multiline")))
}
}

test_asyms <- function(model, asyms, vcov = NULL, escape = TRUE) {
output <- data.frame(term = rep(NA, length(asyms)),
statistic = rep(NA, length(asyms)),
p.value = rep(NA, length(asyms)))
for (var in names(asyms)) {
index <- which(names(asyms) == var)
output[index, "term"] <- var
if (escape) {
test_term <- paste0("plus__", var, " = -minus__", var)
} else {
test_term <- paste0("`+", var, "` = -`-", var, "`")
}
the_test <- as.data.frame(car::lht(model, test_term, vcov. = vcov))
output[index, "statistic"] <- the_test[2, "Chisq"]
output[index, "p.value"] <- the_test[2, "Pr(>Chisq)"]
}
output
}

#' @rdname fdm_tidiers
#' @rawNamespace
#' if (getRversion() >= "3.6.0") {
#'   S3method(generics::tidy, asym)
#' } else {
#'   export(tidy.asym)
#' }

tidy.asym <- function(x, conf.int = FALSE, conf.level = .95, ...) {

if (!requireNamespace("generics")) {
stop_wrap("You must have the generics package to use tidy methods.")
}

params <- x\$coef_table
# params\$term <- rownames(params)

# Getting confidence intervals if requested
if (conf.int == TRUE) {
ints <- as.data.frame(confint(x, level = conf.level))
# Renaming the columns to fit the tidy model
names(ints) <- c("conf.low", "conf.high")
# Renaming the terms to remove the backticks to match the params d.f.
ints\$oterm <- stringr::str_remove_all(rownames(ints), "`")
params\$oterm <- rownames(params)
# Put things together
params <- dplyr::left_join(params, ints, by = "oterm")
}
return(tibble::as_tibble( # Return a tibble
# Only return the relevant columns
params %just% c("term", "estimate", "statistic", "std.error",
"conf.low", "conf.high", "p.value", "group")
))
}
```

## Try the panelr package in your browser

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

panelr documentation built on Aug. 22, 2023, 5:08 p.m.