diff --git a/R/simple_slopes.R b/R/simple_slopes.R index 22f9bc8..0b3c5aa 100644 --- a/R/simple_slopes.R +++ b/R/simple_slopes.R @@ -475,139 +475,140 @@ sim_slopes <- function(model, pred, modx, mod2 = NULL, modx.values = NULL, call[[1]] <- survey::svyglm newmod <- eval(call) } else { - # Creating the model + newmod <- j_update(model, data = dt) } - # Getting SEs, robust or otherwise - if (robust != FALSE && is.null(v.cov)) { - # For J-N + # Handling vcov correctly + if (!is.null(vcov)) { + # Use the externally provided vcov matrix + covmat <- vcov + } else if (robust != FALSE && is.null(v.cov)) { covmat <- get_robust_se(newmod, robust, cluster, dt)$vcov } else if (is.null(v.cov)) { - # For J-N + # Default variance-covariance matrix covmat <- vcov(newmod) } else { vcovargs <- v.cov.args vcovargs[[which(sapply(vcovargs, function(x) length(x[[1]]) == 1 && - x[[1]] == "model"))]] <- - newmod + x[[1]] == "model"))]] <- newmod covmat <- do.call(v.cov, vcovargs) } - - # if (robust == FALSE & is.null()) {covmat <- NULL} - + + # Ensure covmat is a matrix (fixes "closure not subsettable" error) + if (is.function(covmat)) { + covmat <- covmat(newmod) + } + + # Johnson-Neyman interval calculations if (johnson_neyman == TRUE) { args <- list(newmod, pred = substitute(pred), modx = substitute(modx), vmat = covmat, plot = jnplot, alpha = jnalpha, digits = digits) - if (exists("jn_args")) {args <- as.list(c(args, jn_args))} + if (exists("jn_args")) { args <- as.list(c(args, jn_args)) } jn <- do.call("johnson_neyman", args) } else { - jn <- NULL - } - + if (j != 0) { - jns[[j]] <- jn - } - - # Looping so any amount of moderator values can be used - for (i in seq_along(modxvals2)) { - - dt <- d - - # Create new design to modify - if (is_survey == TRUE) {designt <- design} - - # The moderator value-adjusted variable - if (is.numeric(dt[[modx]])) { - dt[[modx]] <- dt[[modx]] - modxvals2[i] - } else { - dt[[modx]] <- factor(dt[[modx]], ordered = FALSE) - dt[[modx]] <- relevel(dt[[modx]], ref = as.character(modxvals2[i])) - dt[[modx]] <- stats::C(dt[[modx]], "contr.treatment") + jns[[j]] <- jn } - - if (!is.null(mod2)) { - - # The moderator value-adjusted variable - if (is.numeric(dt[[mod2]])) { - dt[[mod2]] <- dt[[mod2]] - mod2vals2[j] + + # Looping over moderator values + for (i in seq_along(modxvals2)) { + dt <- d + + # Adjust moderator variable + if (is.numeric(dt[[modx]])) { + dt[[modx]] <- dt[[modx]] - modxvals2[i] } else { - dt[[mod2]] <- factor(dt[[mod2]], ordered = FALSE) - dt[[mod2]] <- relevel(dt[[mod2]], ref = as.character(mod2vals2[j])) - dt[[mod2]] <- stats::C(dt[[mod2]], "contr.treatment") + dt[[modx]] <- factor(dt[[modx]], ordered = FALSE) + dt[[modx]] <- relevel(dt[[modx]], ref = as.character(modxvals2[i])) + dt[[modx]] <- stats::C(dt[[modx]], "contr.treatment") } - - - } - - # Creating the model - if (is_survey == TRUE) { - # Update design - designt$variables <- dt - + + # Handle second moderator + if (!is.null(mod2)) { + if (is.numeric(dt[[mod2]])) { + dt[[mod2]] <- dt[[mod2]] - mod2vals2[j] + } else { + dt[[mod2]] <- factor(dt[[mod2]], ordered = FALSE) + dt[[mod2]] <- relevel(dt[[mod2]], ref = as.character(mod2vals2[j])) + dt[[mod2]] <- stats::C(dt[[mod2]], "contr.treatment") + } + } + # Update model - ## Have to do all this to avoid adding survey to dependencies - call <- getCall(model) - call$design <- designt - call[[1]] <- survey::svyglm - newmod <- eval(call) - } else { - newmod <- j_update(model, data = dt) - } - - # Getting SEs, robust or otherwise - if (robust != FALSE) { - if (!is.null(v.cov)) { - vcovargs <- v.cov.args - vcovargs[[which(sapply(vcovargs, function(x) length(x[[1]]) == 1 && - x[[1]] == "model"))]] <- - newmod - covmat <- do.call(v.cov, vcovargs) + if (is_survey == TRUE) { + designt$variables <- dt + call <- getCall(model) + call$design <- designt + call[[1]] <- survey::svyglm + newmod <- eval(call) } else { - covmat <- NULL + newmod <- j_update(model, data = dt) } - # Use summ to get the coefficients - if (has_summ) { - sum <- summ(newmod, robust = robust, model.fit = FALSE, - confint = TRUE, ci.width = ci.width, vifs = FALSE, - cluster = cluster, which.cols = which.cols, pvals = pvals, - vcov = covmat, ...) + + # Getting SEs, robust or otherwise + if (!is.null(vcov)) { + covmat <- vcov # Ensure external vcov is always used + } else if (robust != FALSE) { + if (!is.null(v.cov)) { + vcovargs <- v.cov.args + vcovargs[[which(sapply(vcovargs, function(x) length(x[[1]]) == 1 && + x[[1]] == "model"))]] <- newmod + covmat <- do.call(v.cov, vcovargs) + } else { + covmat <- get_robust_se(newmod, robust, cluster, dt)$vcov + } } else { - sum <- generics::tidy(newmod, conf.int = TRUE, conf.level = ci.width) + if (is.null(v.cov)) { + covmat <- vcov(newmod) + } else { + vcovargs <- v.cov.args + vcovargs[[which(sapply(vcovargs, function(x) length(x[[1]]) == 1 && + x[[1]] == "model"))]] <- newmod + covmat <- do.call(v.cov, vcovargs) + } } - } else { - if (is.null(v.cov)) { - # For J-N - covmat <- vcov(newmod) - } else { - vcovargs <- v.cov.args - vcovargs[[which(sapply(vcovargs, function(x) length(x[[1]]) == 1 && - x[[1]] == "model"))]] <- - newmod - covmat <- do.call(v.cov, vcovargs) + + # Ensure covmat is a matrix (fixing potential closure error) + if (is.function(covmat)) { + covmat <- covmat(newmod) + } + + # Use summ to get the coefficients + if (has_summ) { + if (!is.null(vcov)) { + # When vcov is provided, do NOT pass vcov argument inside summ() + sum <- summ(newmod, robust = robust, model.fit = FALSE, + confint = TRUE, ci.width = ci.width, vifs = FALSE, + cluster = cluster, which.cols = which.cols, pvals = pvals, + ...) # No vcov here! + } else { + # When vcov is NOT provided, use computed covariance matrix + sum <- summ(newmod, robust = robust, model.fit = FALSE, + confint = TRUE, ci.width = ci.width, vifs = FALSE, + cluster = cluster, which.cols = which.cols, pvals = pvals, + vcov = covmat, ...) + } } + + # Extract coefficients if (has_summ) { - sum <- summ(newmod, model.fit = FALSE, confint = TRUE, - ci.width = ci.width, vifs = FALSE, - which.cols = which.cols, pvals = pvals, vcov = covmat, ...) + summat <- sum$coeftable + if (pvals == FALSE) { summat <- summat[, colnames(summat) %nin% "p"] } + slopep <- summat[pred_names, , drop = FALSE] + intp <- summat["(Intercept)", ] } else { - sum <- generics::tidy(newmod, conf.int = TRUE, conf.level = ci.width) + summat <- sum + if (pvals == FALSE) { summat <- summat %not% "p.value" } + slopep <- summat[summat$term %in% pred_names, , drop = FALSE] + intp <- summat[summat$term == "(Intercept)", ] } - } - if (has_summ) { - summat <- sum$coeftable - if (pvals == FALSE) {summat <- summat[,colnames(summat) %nin% "p"]} - slopep <- summat[pred_names, , drop = FALSE] - intp <- summat["(Intercept)", ] - } else { - summat <- sum - if (pvals == FALSE) {summat <- summat %not% "p.value"} - slopep <- summat[summat$term %in% pred_names, , drop = FALSE] - intp <- summat[summat$term == "(Intercept)", ] - } + + # Have to account for variable amount of rows needed due to factor # predictors diff --git a/test_sim_slopes.R b/test_sim_slopes.R new file mode 100644 index 0000000..2fc767e --- /dev/null +++ b/test_sim_slopes.R @@ -0,0 +1,29 @@ + + +library(interactions) +library(sandwich) +library(lmtest) +library(jtools) + + +model <- lm(mpg ~ hp * wt, data = mtcars) +# Assuming 'cyl' is the clustering variable + +clustered_se <- vcovCL(model, cluster = ~cyl) + +#sim_slopes(model, pred = "hp", modx = "wt", vcov = clustered_se) + + +sim_slopes(model, pred = "hp", modx = "wt", jnplot = TRUE,vcov = clustered_se) + +clustered_se <- vcovCL(model, cluster = ~cyl) +sim_slopes(model, pred = "hp", modx = "wt", jnplot = TRUE, vcov = clustered_se) + +getAnywhere(sim_slopes) +conflicted::conflict_scout() + + + +getwd() +setwd("/home/onyxia/interactions") +getwd()