From b6cf95e998a323dd60f1a7476558f793b8403d17 Mon Sep 17 00:00:00 2001 From: Shu Fai Cheung Date: Tue, 29 Oct 2024 17:49:19 +0800 Subject: [PATCH 1/8] WIP --- rebuild_vignettes.R | 1 + vignettes/betaselectr_lm.Rmd.original | 446 ++++++++++++++++++++++++++ 2 files changed, 447 insertions(+) create mode 100644 vignettes/betaselectr_lm.Rmd.original diff --git a/rebuild_vignettes.R b/rebuild_vignettes.R index c97d3b6..6dd48eb 100644 --- a/rebuild_vignettes.R +++ b/rebuild_vignettes.R @@ -4,6 +4,7 @@ base_dir <- getwd() setwd("vignettes/") knitr::knit("betaselectr_lav.Rmd.original", output = "betaselectr_lav.Rmd") +knitr::knit("betaselectr_lm.Rmd.original", output = "betaselectr_lm.Rmd") setwd(base_dir) diff --git a/vignettes/betaselectr_lm.Rmd.original b/vignettes/betaselectr_lm.Rmd.original new file mode 100644 index 0000000..d759930 --- /dev/null +++ b/vignettes/betaselectr_lm.Rmd.original @@ -0,0 +1,446 @@ +--- +title: "Beta-Select Demonstration: Regression by `lm()`" +date: "`r Sys.Date()`" +output: + rmarkdown::html_vignette: + number_sections: true +vignette: > + %\VignetteIndexEntry{Beta-Select Demonstration: Regression by `lm()`} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.path = "" +) +``` + +```{r, echo = FALSE} +format_str <- function(x, digits = 3) { + formatC(x, digits = digits, format = "f") + } +``` + +# Introduction + +This article demonstrates how to use +`lm_betaselect()` from the package +`betaselectr` +to standardize +selected variables in a model fitted +by `lm()` and forming confidence +intervals for the parameters. + +# Data and Model + +The sample dataset from the package +`betaselectr` will be used for in this +demonstration: + +```{r} +library(betaselectr) +head(data_test_mod_cat) +``` + +This is the regression model, fitted by +`lm()`: + +```{r} +lm_out <- lm(dv ~ iv * mod + cov1 + cat1, + data = data_test_mod_cat) +``` + +The model has a moderator, `mod`, posited +to moderate the effect from `iv` to +`med`. The product term is `iv:mod`. +The variable `cat1` is a categorical variable +with three groups: `gp1`, `gp2`, `gp3`. + +These are the results: + +```{r} +summary(lm_out) +``` + +# Problems With Standardization + +One common type of standardized coefficients, +called "betas" in some programs, is +computed by standardizing *all* terms +in the model. + +First, all variables in the model, +including the product term and dummy +variables, are computed: + +```{r} +data_test_mod_cat_z <- data_test_mod_cat +data_test_mod_cat_z$iv_x_mod <- data_test_mod_cat_z$iv * + data_test_mod_cat_z$mod +data_test_mod_cat_z$cat_gp2 <- as.numeric(data_test_mod_cat_z$cat1 == "gp2") +data_test_mod_cat_z$cat_gp3 <- as.numeric(data_test_mod_cat_z$cat1 == "gp3") +head(data_test_mod_cat_z) +``` + +All the variables are then standardized: + +```{r} +data_test_mod_cat_z <- data.frame(scale(data_test_mod_cat_z[, -5])) +head(data_test_mod_cat_z) +``` + +The regression model is then fitted to the +standardized variables: + +```{r} +lm_std_common <- lm(dv ~ iv + mod + cov1 + cat_gp2 + cat_gp3 + iv_x_mod, + data = data_test_mod_cat_z) +``` + +The "betas" commonly reported are the +coefficients in this model: + +```{r} +lm_std_common_summary <- summary(lm_std_common) +printCoefmat(lm_std_common_summary$coefficients, + digits = 5, + zap.ind = 1, + P.values = TRUE, + has.Pvalue = TRUE, + signif.stars = TRUE) +``` + +```{r} +library(lm.beta) +lm_beta <- lm.beta(lm_out) +lm_beta_summary <- summary(lm_beta) +printCoefmat(lm_beta_summary$coefficients, + digits = 4, + zap.ind = 1, + P.values = TRUE, + has.Pvalue = TRUE, + signif.stars = TRUE) +``` + +```{r, results = FALSE} +lm_beta_x <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat, + do_boot = TRUE, + bootstrap = 5000, + iseed = 4567) +``` + +```{r} +summary(lm_beta_x) +``` + +However, for this model, there are several +problems: + +- The product term, `iv:mod`, is also + standardized. This is inappropriate. + One simple but underused solution is + standardized the variables *before* + forming the product term (Friedrich, 1982). + +- The confidence intervals are formed using + the delta-method, which has been found + to be inferior to methods such as + nonparametric percentile bootstrap + confidence interval for the standardized + solution (Falk, 2018). Although there + are situations in which the delta-method + confidence and the nonparametric + percentile bootstrap confidences can be + similar (e.g., sample size is large + and the sample estimates are not extreme), + it is still safe to at least try both + methods and compare the results. + +- There are cases in which some variabLes + are measured by meaningful units and + do not need to be standardized. for + example, if `cov1` is age measured by + year, then age is more + meaningful than the "standardized age". + +- In path analysis, categorical variables + are usually represented by dummy variables, + each of them having only two possible + values (0 or 1). It is not meaningful + to standardize the dummy variables. + +# Beta-Select `lav_betaselect()` + +The function `lav_betaselect()` can be used +to solve this problem by: + +- Standardizing variables before product + terms are formed. + +- Standardizing only variables for which + standardization can facilitate + interpretation. + +- Forming confidence intervals that take + into account selected standardization. + +We call the coefficients computed by +this kind of standardization *beta*s-Select +($\beta{s}_{Select}$, $\beta_{Select}$ +in singular form), +to differentiate them from coefficients +computed by standardizing all variables, +including product terms. + +## Estimates Only + +Suppose we only need to +solve the first problem, with the product +term computed after `iv` and `mod` +are standardized: + +```{r, results = FALSE} +fit_beta <- lav_betaselect(fit) +``` + +```{r, echo = FALSE} +tmp <- capture.output(print(fit_beta)) +``` + +```{r, eval = FALSE} +fit_beta +``` + +This is the output +if printed using the +default options: + +```{r, echo = FALSE} +cat(tmp[-c(25:51)], sep = "\n") +``` + +```{r, echo = FALSE} +b_std <- standardizedSolution(fit)$est.std[3] +b_select <- fit_beta$est[3] +``` + +Compared to the solution with the product +term standardized, the coefficient of +`iv:mod` changed substantially from +`r format_str(b_std)` to +`r format_str(b_select)`. As shown by +Cheung et al. (2022), the coefficient +of *standardized* product term (`iv:mod`) +can be severely biased estimate of the +properly standardized product term +(the product of standardized `iv` and +standardized `mod`). + +The footnote will also indicate +variables that are standardized, +and remarked that product terms +are formed *after* standardization. + +## Estimates and Bootstrap Confidence Interval + +Suppose we want to address +both the first and the second problems, +with the product +term computed after `iv` and `mod` +standardized and bootstrap confidence +interval used, we can call `lav_betaselect()` +again, with additional arguments +set: + +```{r} +fit_beta <- lav_betaselect(fit, + std_se = "bootstrap", + bootstrap = 5000, + iseed = 2345, + parallel = "snow", + ncpus = 20) +``` + +These are the additional arguments: + +- `std_se`: The method to compute the + standard errors as well as confidence + intervals. Set to `"bootstrap"` for + nonparametric bootstrapping. + +- `iseed`: The seed for the random number + generator used for bootstrapping. Set + this to an integer to make the results + reproducible. + +- `parallel`: The method to be used for + parallel processing. It will be passed + to `lavaan::bootstrapLavaan()`. Supported + values are `"none"`, `"snow"`, and + `"multicore"`. + +- `ncpus`: The number of CPU cores to + use if `parallel` processing is not + `"none"`. Default is + `parallel::detectCores(logical = FALSE) - 1`, + or the number of *physical* cores + minus one. + +```{r, echo = FALSE} +tmp <- capture.output(print(fit_beta)) +``` + +This is the output if +printed with default +options: + +```{r, eval = FALSE} +fit_beta +``` + +```{r, echo = FALSE} +cat(tmp[c(1:27, 55:66)], sep = "\n") +``` + +In this dataset, with 200 cases, +the delta-method confidence +intervals are close to the +bootstrap confidence intervals, +except obviously for the +product term because the coefficient +of the product term has substantially +different values in the two +solutions. + +## Estimates and Bootstrap Confidence Intervals, With Only Selected Variables Standardized + +Suppose we want to address also the +the third issue, and standardize only +some of the variables. This can be +done using either `to_standardize` +or `not_to_standardize`. + +- Use `to_standardize` when +the variables to standardize +is much fewer than the variables +not to standardize. + +- Use `not_to_standardize` +when the variables to standardize +is much more than the +variables not to standardize. + +For example, suppose we only +need to standardize `dv` and +`iv`, `cov1`, and `cov2`, +this is the call to do +this, setting +`to_standardize` to `c("iv", "dv", "cov1", "cov2")`: + +```{r, results = FALSE} +fit_beta_select_1 <- lav_betaselect(fit, + std_se = "bootstrap", + to_standardize = c("iv", "dv", "cov1", "cov2"), + bootstrap = 5000, + iseed = 2345, + parallel = "snow", + ncpus = 20) +``` + +If we want to standardize all +variables except for `dv` +and `mod`, we can use +this call, and set +`not_to_standardize` to `c("mod", "dv")`: + +```{r, results = FALSE} +fit_beta_select_2 <- lav_betaselect(fit, + std_se = "bootstrap", + not_to_standardize = c("mod", "dv"), + bootstrap = 5000, + iseed = 2345, + parallel = "snow", + ncpus = 20) +``` + +The results of these calls are identical, +and only those of the first version are +printed: + +```{r, eval = FALSE} +fit_beta_select_2 +``` + +```{r, echo = FALSE} +tmp <- capture.output(print(fit_beta_select_2)) +cat(tmp[c(2:27, 55:66)], sep = "\n") +``` + +The footnotes confirmed that, by +specifying that `dv` and `mod` are not +standardized, all the other four variables +are standardized: `iv`, `med`, `cov1`, and `cov2`. +Therefore, in this case, it is more +convenient to use `not_to_standardize`. + +For *beta*s-*select*, researchers need +to state which variables +are standardized and which are not. +This can be done in table notes, +or in a column of the parameter estimate +tables. The output can of `lav_betaselect()` +can be printed with `show_Bs.by` set +to `TRUE` to demonstrate the second +approach: + +```{r, eval = FALSE} +print(fit_beta_select_2, + show_Bs.by = TRUE) +``` + +```{r, echo = FALSE} +tmp <- capture.output(print(fit_beta_select_2, + show_Bs.by = TRUE)) +cat(tmp[c(15:27, 55:70)], sep = "\n") +``` + +## Categorical Variables + +When calling `lav_betaselect()`, +variables with only two values in +the dataset are assumed to be categorical +and will not be standardized by default. +This can be overriden by setting +`skip_categorical_x` to `FALSE`, though +not recommended. + +# Conclusion + +In structural equation modeling, there +are situations in which standardizing +all variables is not appropriate, or +when standardization needs to be done +before forming product terms. We are +not aware of tools that can do appropriate +standardization *and* form confidence +intervals that takes into account the +selective Standardization. By promoting +the use of *beta*s-*select* using +`lav_betaselect()`, we hope to make it +easier for researchers to do appropriate +Standardization in when reporting SEM +results. + +# References + +Cheung, S. F., Cheung, S.-H., Lau, E. Y. Y., Hui, C. H., & Vong, W. N. (2022). Improving an old way to measure moderation effect in standardized units. *Health Psychology, 41*(7), 502--505. https://doi.org/10.1037/hea0001188 + +Falk, C. F. (2018). Are robust standard errors the best approach for interval estimation with nonnormal data in structural equation modeling? *Structural Equation Modeling: A Multidisciplinary Journal, 25*(2), 244--266. https://doi.org/10.1080/10705511.2017.1367254 + +Friedrich, R. J. (1982). In defense of multiplicative terms in multiple regression equations. *American Journal of Political Science, 26*(4), 797--833. https://doi.org/10.2307/2110973 + + From 39546934d3bfca9c829bfe5cc2a89ba7a9e1aaac Mon Sep 17 00:00:00 2001 From: Shu Fai Cheung Date: Tue, 29 Oct 2024 19:12:25 +0800 Subject: [PATCH 2/8] 0.0.1.12: Fixed a bug in printing the bootstrap p-values Tests, checks, and build_site() passed. --- DESCRIPTION | 2 +- NEWS.md | 8 ++++++-- R/lm_betaselect_methods.R | 2 ++ README.md | 2 +- 4 files changed, 10 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index e39afd6..8d5aa1f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: betaselectr Title: Selective Standardization in Structural Equation Models -Version: 0.0.1.11 +Version: 0.0.1.12 Authors@R: c(person(given = "Shu Fai", family = "Cheung", diff --git a/NEWS.md b/NEWS.md index cb6009f..d24230d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# betaselectr 0.0.1.11 +# betaselectr 0.0.1.12 - Added `lm_betaselect()` and related methods and helper functions. @@ -63,4 +63,8 @@ - Updated `print.summary.lm_betaselect()` to print confidence intervals by - default, if available. (0.0.1.11) \ No newline at end of file + default, if available. (0.0.1.11) + +- Fixed a bug in printing the bootstrap + *p*-values for `summary.lm_betaselect()`. + (0.0.1.12) \ No newline at end of file diff --git a/R/lm_betaselect_methods.R b/R/lm_betaselect_methods.R index b60bc0e..d9a1c1e 100644 --- a/R/lm_betaselect_methods.R +++ b/R/lm_betaselect_methods.R @@ -949,6 +949,7 @@ summary.lm_betaselect <- function(object, colnames(out$coefficients)[i] <- "z value" if (boot_pvalue_type == "asymmetric") { boot_est_list <- split(boot_est, rownames(boot_est)) + boot_est_list <- boot_est_list[rownames(boot_est)] boot_pvalues <- sapply(boot_est_list, est2p, h0 = 0) @@ -1462,6 +1463,7 @@ summary.glm_betaselect <- function(object, colnames(out$coefficients)[i] <- "z value" if (boot_pvalue_type == "asymmetric") { boot_est_list <- split(boot_est, rownames(boot_est)) + boot_est_list <- boot_est_list[names(boot_est_list)] boot_pvalues <- sapply(boot_est_list, est2p, h0 = 0) diff --git a/README.md b/README.md index f12c209..39afb8e 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ Not ready for use. # betaselectr: Do selective standardization in structural equation models and regression models -(Version 0.0.1.11, updated on 2024-10-29, [release history](https://sfcheung.github.io/betaselectr/news/index.html)) +(Version 0.0.1.12, updated on 2024-10-29, [release history](https://sfcheung.github.io/betaselectr/news/index.html)) It computes Beta_Select, standardization in structural equation models with only From 50830d6488e39e5274c1a5b944755e5544f7a590 Mon Sep 17 00:00:00 2001 From: Shu Fai Cheung Date: Tue, 29 Oct 2024 20:10:18 +0800 Subject: [PATCH 3/8] Add a demo dataset for vignettes --- R/data_test.R | 31 +++++++++++++++++++++++++- data-raw/data_test_03.R | 43 ++++++++++++++++++++++++++++++++++++ data/data_test_mod_cat2.rda | Bin 0 -> 16396 bytes man/data_test_mod_cat2.Rd | 38 +++++++++++++++++++++++++++++++ 4 files changed, 111 insertions(+), 1 deletion(-) create mode 100644 data-raw/data_test_03.R create mode 100644 data/data_test_mod_cat2.rda create mode 100644 man/data_test_mod_cat2.Rd diff --git a/R/data_test.R b/R/data_test.R index c237f5d..6374cf3 100644 --- a/R/data_test.R +++ b/R/data_test.R @@ -56,4 +56,33 @@ #' summary(lm_out) #' #' -"data_test_mod_cat" \ No newline at end of file +"data_test_mod_cat" + +#' @title Test Dataset with Moderator +#' and Categorical Variables (Version 2) +#' +#' @description This dataset has one +#' predictor, one moderator, one +#' control variable, one dependent +#' variable, and a categorical variable. +#' +#' Similar to `data_test_mod_cat` but +#' generated from another population. +#' +#' @format A data frame with 500 rows +#' and five variables: +#' \describe{ +#' \item{dv}{Dependent variable, continuous} +#' \item{iv}{Independent variable, continuous} +#' \item{mod}{Moderator, continuous} +#' \item{cov1}{Control variable, continuous} +#' \item{cat1}{String variable with these values: "gp1", "gp2", and "gp3"} +#' } +#' +#' @examples +#' +#' lm_out <- lm(dv ~ iv * mod + cov1 + cat1, data_test_mod_cat) +#' summary(lm_out) +#' +#' +"data_test_mod_cat2" \ No newline at end of file diff --git a/data-raw/data_test_03.R b/data-raw/data_test_03.R new file mode 100644 index 0000000..00c338e --- /dev/null +++ b/data-raw/data_test_03.R @@ -0,0 +1,43 @@ +# Adapted from stdmod + +set.seed(41543) +n <- 500 +x <- rnorm(n) +w <- .4 * x + sqrt(1 - .4^2) * rnorm(n) +v1 <- rnorm(n) +cat1 <- sample(c("gp1", "gp2", "gp3"), n, replace = TRUE, prob = c(.2, .3, .5)) +y <- .2 * x + .3 * w + .4 * x * w + .12 * v1 + + sapply(cat1, + switch, + gp1 = 0, + gp2 = .4, + gp3 = .6) + + rnorm(n) +dat <- data.frame(dv = y, + iv = x, + mod = w, + cov1 = v1, + cat1, + stringsAsFactors = FALSE) +out <- lm(dv ~ iv * mod + cov1 + cat1, + dat) +summary(out) +apply(dat[, 1:4], 2, sd) +colMeans(dat[, 1:4]) +b <- 1 / c(5, 2, 4, 2) +a <- c(20, 15, 50, 20) * -b +dat0 <- scale(dat[, 1:4], + center = a, + scale = b) +dat0 <- round(dat0, 2) +apply(dat0, 2, sd) +colMeans(dat0) +apply(dat0, 2, range) +dat1 <- data.frame(dat0, cat = dat$cat1) +head(dat1) +out <- lm(dv ~ iv * mod + cov1 + cat1, + dat1) +summary(out) +summary(lm.beta::lm.beta(out)) +data_test_mod_cat2 <- dat +usethis::use_data(data_test_mod_cat2, overwrite = TRUE) diff --git a/data/data_test_mod_cat2.rda b/data/data_test_mod_cat2.rda new file mode 100644 index 0000000000000000000000000000000000000000..e6fa00fd8a6dd15f0f4453118d9a67cb64c85032 GIT binary patch literal 16396 zcmV+nK=Z#sT4*^jL0KkKS-(I*VE_{{fB*mg|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0 z|NsC0|Nr1A-v9*F=LsFTJ}3 z000000-rzt15G$|$TS+50%nBM(qI5g0$|V;KTHB(00J-xfHasX=)i*lU`!KCYMvr` zXkY_DFlYq8$u^n;Kxhe#B@U{e?G*8*hE-4U7!iom(V}@YO)yO`(3u!1S>w-Q&ZDa_9jeeQ`(OqrkIrdH8#>>PilIMnw}zUGM-b(rkSEH(9=q7spSWWwNvm4Jq)Lm#-`ejR3dsr(E=tA$)=eBJprMq=$kY8JY;tK#*o+03d__41s`P%*-Ic5rQK$ z7$Y(=21X3Q5C{MO03b$SAi;t%002e^85tNd002M$KoI~CW@dm!MhMIig9Hr*2+WKb zk&&4g0KfndnV295%*Yu71`PlpfHE*(5t)$yg9HJYAQ%K;24rA@03Zm=j1id`nhX#i z0RuBLAPR`hC^IuN1Vm&IV8NN00w4^*nS%h(W?;<#jDSIaW&r?XN)k;Vnj(D!i)-;A z2tp8pM}!e15poPJgR|=W*jY*$Vnjy!UE0on^Zp9y_u~4!h|gxH%<$~@-&8K0rdMJI zLpws9X8~#u`CHaHLiXMW6vJ$S;_PcG$^;pptL;3o@bCB7RepvK!EGrj)(Tb3hEi1s z91!DPItxL+v-@Z456DIGX$f~+T zvTJv7Y+mMGxzF+nxl^*KZMr+;Ila^G(l32cBG2!FHEmvdoonnig{G4DzJg*7J zOdp1aVW9tS6791d%e>F5i%<&ndv?qbfQd>2$f_+#S9T|TC&vINcj|u)#vH953z%#P zO`e*KNAGjB%0DDQPM%7hfh0bgN*G$SDtS%C_F40Q zb(bng4))^f!E6lDi-7b8n)&QqDLlVbb~g>nV5%vpW*!A=r0)pxXSgPI-eNX>bX&S` zi%r_~>{_6b3#Lt6dSSvRw^V^)3qm)hO~loxi+0-k)z8GjnKt#_l}OFIJv+?OE-lIt zE{(Qy=$Ynu!g+$6^!}^dXca&K0t5&UAapm@DOQd3&XqVCi}Bsa8O3}1oF1`yyAVNZ z_KOa_^-PcpWoCUE)M^y~rmtCWy9}Me^v2*1?JWVhBDQIMZPOI3e6YyQw|C5w&`UAZ z2u?KBpnDET7L~}}<^ItPjmZ9`gWE*Y%aM)X@Tn2k>{E~4w6&c3z-dWJT2Q0Kfylgf zDDAl}tBx61w++bPmpV-?SA;M;8wl7J$V?gt+db{hB25c|11*{I{J6dYgL-s77`*WN z1cdq=+FV0UGg#HiM7fG(0CIMsD)A6J^Z~wZ&70NZWX9%0j7wb5A~#AI9BRKc9PuHk zstX0*s-zywc^KVkn6Ngmlejr~1@#xF3H^4INI>yZ06;iH_q%8fXv$?wOD=x!@Tu(Y zkbV2MdH&IfF)LI6>k-dqJ@!^bMES)LEnrQE0rSVGs9ta?!#$5(etQ*uEpYmN;}yEsq&?b2BaPziMjLdpHT-lP^anWUMqYq~3RJ&Yt-MpHkm8;jKGi4zMc z!xhP#j5>GMh^1@NcJBx!H?j7e)+`)jshEQgwez@e zNTdGDYenr0%`KYGN>hw;j?+Z#D`DPVXlOdPaKr=jc*b{Yg_{(vgHf(1T+)+2i zi#)VYu2hUi)UOx)8ILZaJB-&?5FSpP3`4y_V<1;}?yoPcl&9DoyMP+KF+@IPoFSG282}QBU;QYb7(z@SH;l zVpsGfSnkuP2h_L|-kPh-X>wO+3(^d5!ddX4`)Jlqx2V6F?u<0t{e&B1sb5xG8B$Ho z@HX?hf7K3m5wNPDX&|*p$fV%+Owm$w)stNSw;-SIn-!JM175P={RFg3>i2Z9)q`0Rm}q3ZOXL5HOfB9@p%gSp&k1-CSk+uX zrkTD$AzPX0*XJ z@GG4jV7<2RgD6NYW9OZPElz+d{fwx*|>Ma1pj#aKfZBjh!XH*LQ= zOw7zhlxuh+m05;t^PIC6|D-1GbA5n*fmPGuE#WU7aygYC(j zS&RhaAIl$}th`s#6B2Lm>?2n%;b{4X@ta~7$Sk|!>{%zkSn~eT6?ie;Z@fB?v!2^G zu;cOf)0zfjeR5=OAnecaiN+w`y34Mh>qxrWzbypn{!N{>IW8}s0SF?K?(2S&QV_3q zmB{C7(~!D#bnx)Y=syzrUDFDeFw;XMb(*xDL0uaU1UCa=OsZD?AAODy24Ge1Uz*<9 zaR1czJpZ%KLg?Vt({z4(p-n~sQAwELSp`gF2;uM`6KOvx$G}bJr)u)CC~&YK^;meH-;Gmi=wVE zUy=8TRZA;y3P9tETSh{|deWQ2aB>Wm{fA=TC9awds}!(WWHY-6i`t8t+8Gdq1#j$D z%30wTfBVR6pz&5W33nzRl;_=u<9J0Ohx`ja@NN^KPIt^``}x1c+yPz$7z!J5()5b3 zo&A3v3gpC_3=Ple^Y)um8ik&LFV`S$!-9I(rQ<;zCfCaKo9UKm*>O;V5-$13!e4tw zYx@6RN&B}lhKp{Da8~y*km^Drc|#GF4->dNP`4s<5pMyv=`hW9ipDP@2y#KnLL&q zazCE5#U40Q?w9!b(*7@c*Xeah3<@(U^!&@H`0jSuQt!eVld)Gt9w#?CZN9eaNb9WL zEP9B*{!jAlDK-&cR{bCV1i%0jMV?pE)$2FXirvq|9l3NQ+Y1g`PILyuHn|u%2dAe- zmRL)nF?(3tI?ij8lS4y54FCZ8+dEDnESI06r5ESOqci);ZAb%HNU9IFg?GKWl#!Ex zoqF0h{bl_%MUJ;BnpO3f=MU#-D-aL>0D%Gr#}Yy^`pdKtG#nSTDch@O{0$uvq+;e* z6Yytkkp(0$(E7r0S*YZ`8xlWC&aB9_Nwb_nwk& zO_G3QZvbv~z6l9FFJ@Uc3JU9>HZiLHfll4pqYS2L)e}4b2tv)s)k8|B!tUTKa#vgx z0<3oC&E}FtdRlK@H=T>}OgYUL!;f241RBL6qo37@@bs!;JqQLbrsug_i0#Fr{nw31A1kjZs`Nyj zZ+@vHod{@b7Z5p;|1WB{=6N|cF&wPTLqN{y39<_3A5>gdwHC|G+^n-y$Rdi4aja|{ zQZcVMlVY0Et%r4SHO}L(4}KUuW2}TJb0jcoqrCGoN^|;q0fDnJ{wJn3MglIPMQm#N zEw(Xf-(Hx)van-j4;+%cc|szIG<=XOG!Juz3rlxNSg>L@5QGF}&0bywM57SCv`^f` zi`%X`C~*}pK{R4>(^Tu1lOzdf3gu26(p<&X%oqjO`mH4e?$ zLQmfAgpf*8%|E`7VIUC=S@X2Ev7$3-k1S^Dwx2tPHJW@UlGGoG&B$tvlNQRFUX2=u z0Wp**OTIkN{>fdXg(Hw339&DpZra5_u3kX0UFWUKWc)yaV(zVR^hUr}P zjqEA-wzrz#k_8!1$_wFBnE)2fQ&hn9gV`*w2>f%#4o_QUOTGDt%(v>AWyg(%(qFs# zsrKgjnkOQ4acvKi?lRmuw4$TyMfjo3lkTivOSvPEdBi!z0!8Hpd{i}GCTD#oGJqPb z6IF13!X)fCHvT~xq2%I_o7|o5c%!`m8JB}%&fi{o4APqId+g19|Jf2GAHq$M$XCA)~kkvF5)Ljb!y zM>^tqG^9_m>8gjbUW=bHBQsfc$K6Zl7m^umLNJh?!>sw>Stb~}17|3+TT^0`(NvG@ zL$H5iXBIC^S7i#w+GuOb?44;Atkv6$N4on zc%zSFx3ilW--yIg3c%5ShZ_er4!2p&B{#(L?Ms!+cxnprO zTAJP6JZbNZxMI-k%Ty?W4vuz9WeAJ6(wgh73FJa0Jqw@aXmP#D?s*2uuViOGJwJo0 zE}#MelQeJE(F2g8zTQD5aiyLbye|^IFYt7B582y30@78aN6zkz!mr04(J8~3O|ViW z!23nEYH_idQi+DzhFd>xktnxvC5o*I?o$)By))J=c_*K?{MvQBkH@myrS03t?c{)K zbWLMtf`u@Rs=ylXkh6kXg-Gcw&8)z*>(Mqpe>|fx&!9IF?dgW=Ckb2jlIqc^5X2N|p;EceiqAI`gyL#|HcnhU%y)y3;g;3MrwE;ifJ z{;6jTM;K*&S*u5y!po0iT<;pVo3!WX@bzO{KJ;5DKtddRlH9pZUhPQEOajr-)=!oB zGb^AZ_}IfIRG)}RNm=3jsq|XZ!OQeY8R%DNiGekWe`bfz?!ZX(Y8tmJv(26bd|cvbFJA`LS`O}1La=#)9~9A2*&(#Dnu}t z?-QQZRhx63=kt7fv$zh_AAV68o6%i?gs&M@?Ycodj+)aKj4PF#Qbm8S=!?K97dFyW`F3kU8OWG(pf*J%984V5VWQgY`?s!12(te;dk&$7mbPo>{~dj*2{;dwpBQScP?CwMJ4q+tZX^87PH z<6^nAuI{|T679AYo!75)erBr-`76QBH&?cnOB12zE*Hylk!drhw5mf9*bPFZ_4}{z zq+PfF90rYmQ{|Yo_VPxQE5F9RG&f#`BXXTpdnFIE#S8nR#B?d%9)?ww?18+Fez9Ed z5v+Qi1!eoja7$&s@GUVyRS2xh1j#4+Sdu^pzoc0X;?G~pCfx=*+WrA>UNiL6@wQ)| z1O>AYBHxyTzwItR9mmYZsRFzpP)#=&=dWe1%mTc&NOZA@hUxFP2ePdvp$ z#D~`zeyQ$3vonrcF^B>Phv2AuB=gc8{9|W3Cw=x2>usT(+m~3g{~*S?{e(_LbEd_1 zT}{(L0VOAGe9f+f47ZrGloa?*jhW5!=irB;p7HJCOfTove%W3g-3XG7_Er23gbnSt zgW!c9nmTL&*Zn_OOl^o@v&$WNKbtP4;_xJo-VFt}_v{0#HIsnk!l2C2mI~MhOSTY0e)byr{FBH}rzWcfK;=p!JA>@s-!Wvq1&qJn>3(oa0}7Ql^n&R>SAQ>z*8=FfEXvc}=Ph2y&mdNq z4_QvP$Ub^*U}M(pzn+2d5kIY02&4074$dUH>>GZUFTb;1&p)v*VQ}#1xS?S!$YuqP zh98?{)%d)ayeVc`07w0EF!%$*sm^apH{o#8z|~KWpnHbFRkNwC*pznoBFUWA@=h_b zs4>AYXkr5OIZ4l|-NtG|G3TbW!ucpMQ=2iK29S!#R(#BwNx<%ZLK7zGfwwV~yruoz zPc7eH4`uSct?Mul zcfP}MP9?Eay|6Lt7n+>DV^*DPWsyQp?|2=vycAjV}FS_~HP_z_%G+XJotqJY?bc_Cx)a;m%G>i!!;JIJr4bXp*-p-`1LN z$a-^d>(-A9CU@JFYjL?KwXDmzy}S;o-%9+>E?s(>(>ug71!2z^B*5m&hO{K=-5TAy zI-LZtrUR{*@dZ3LS*{cneFBdNJ(xJOb2_zs3W%_ix%SzBzL1I82}+rCR%JBGb+ZVX zMu=vO4*31b2ZY~*Zw56m67%cK0WiB(I|L64Z?oewiR{Jbb`dM%(e~vi;0CNsG5l5d z72y}iar@1x~hx&xr>6kmwiA-eGN(ygfYAo134)jGUXHwCxN}p!S5* zUOPLqQh>_aH+5=9AOaZ?{pbIBb85D9JjFCT#5by^ec6Y&vnO`KmXON}Ql4@4xa=UW z-Xij5;(&D_eVs0B*Y7 z=y!$kdfX6~QR@|kQdsmQz!|!p5&)@^Zr`8_DkG97f{exTl}dPH5VKks=# z=|VoznbW4L6%uuhkKb*Wu2jgcZA#*UAYU>v`!>gp1V4MoI(U zr?-OpHsBxYCe*X!<_50Af*|V8n@h*alkTMRiumNLV5OMHs%BqF9*ApM1k@QPHt(5M{m3pcXs2@rg@@ueS z^uC!lH}FsKjM0IMHXrJJN{T~r%v!(Zppm`Sb$OtH?nFIm>3uZ3B*M64tEUg-{(FBH zXEbNAElesB$%i`$PY|ynzjZL$7wz28Yh+C)vblOBP~G)UmnUn8gD8y$8RmZeKRtie z$TDfzPgC%bPlgR&LK|j)qW9r|t7`UJt+k_*LfdL{CEr4mOlx@%aY=^M_`230VVKR& zoO)Lxz_QNI4fBf)(9B3pII~5G#5$H$O`I+qLOAep zHL~H`pD}DLBsLbl%YIK^%{~Uw{-3m}g=bHznS834)+Fb$4)B-)HKX}rEiEDLGz|V( zAxEpoYh6lT{1x+j^uJ9U7^#zCi?{lDU;q;UCMG5(yKBwwC)frN(K{t@*LQqv=C;XY z1sMXcn0uA^4*O%Tv!6tn9UP+UWj7_FZk5Tu%X;_pR*Anc#rVk_%0|>trca=p9$v8^ z!DY)4j!wVzQp)8Z65Q0OO(%^PG?#Jqp0#JBh#cq#v zcd*fO(INY>XaK*Jbn49Nd4xDuRSbkD@^}2>#M`TTUn79hr8xeIuWy|RMp6Zb>S}+0 zY?dW;Cqw3Uoz_C<@McExwT(g6pS^)AL3u-=kX)J~toTS+1#ai{EYQO~Z@|3?r8JCx zhhYc45*Ux#ZN0=7i)bGE%)07GsJ(hxnT3&0X%Nb*F-tVQX$OWloH!nBmMt-()LBXD zA8? z!u{D50X(7}aKyD=2~D&Tkh1>1j-5hiK*-2#RTSQSV}|Wyp^X1@Vo!s8lG9ZwC<+Jh zOx1{g2CxJ@XK*}KQ~a;gWK%dy;=?rG1_6=EKhxqA zKw__X@QRd047ZrBm`(bLE92jH9U(|a4AX8t;^_xL>P0#^G>rd(lAI-5=q#{0qkmK3 zJX!M^{n2Wtn--A18&`GAPLkR|JPNjG@Y@EW;aia!zsZ){0U@z`3-}s%$o<0tetg)+ zxy&)ItR52Xf2F$YNM^0Vu<+PK)+<-kev%EN)gGQji|4`b@eqvF(miIRO%~O-m zuwn#GqDKc?m)DC^#(vSGKzj4xvN?B5;vaPl&~H4fdM{4O2))eZ&S*ih(RUiTYZ6~W z%EAljgHopVDjYkX^2c-?o&^fo3~ft+hL2Exm8C|h`rR`PPG9kZh(CG$3}9AeTLiWl z__$9Q^D#Xkec(kwU)`9e4IyJrSrK?bKCReLw5wVrNY_^GUv_X-sE^IN_*78G>Gc>W zk?0%ZD0yAL27}tUV)ay_Zi*C-s9?tDI0)<@+^1**u*XFw(j7cpGD zo-9ew^h!%*MjHIuq&l*zvx|-NAyIy)b01CP^z)8G3&E4qe z7?2sJooh_ep#Mh{5Vnkfr$^;jGs(-d?mnI}i*Mn2bRUA*WH8>HQ}2aEb!>jI zkEDy1np7_IFNGIF`=ZT(KSWu-$4~x^On>8E-*V|Ih~S>rznZi0%8(aft8w~$ zE^>_&8sK=W-W2w|t#`e~oxN)_NI^I_jWy$g5fw)h{d)b@US&UOq0$)_kPnWEkEL9L z-4Jj1c+EvGP)M3FZt;*(JD8_XN7~U=s$e{Scxd2@C(1FAx2c@!G16q2#NEf-MxQud zVK@D-LxyB`C+zi6s{Gt)7AC*#yHyMTJ%2j6gb~=RC!hJKD1O`BJl3u664N5#TrWMX z%<$fs@DVdd-81?YbX zcP4?BoOc@w=(9KF_n(N+81QsTPXebDU#mst>2K~^a6_nCfgxn;hIJ5QEDZ>93{c(k z39r!7^cV0At}bbXY*ZIvZSnJaySf!msRC7s)J39M;KnW}q%ZR#*A{J+SBh@;X#T!g z#;3EV#^T*_JO$6iD`gF1f10F0g&|xZqBiB-t^b}FQnk8C9z5$b3kyPZ{+`zBW(Ffa zRQu3nmc5swYwrpz93NSl{u7~Ae8N^$W_D0O{@lof+9M*51S6j8A%+RzpvyH(^@dtH zR!C!K$_1u^?$l{#!pP4L^fQRH$KI6t3v>KtRG>f`QqlW^WKuV-uno}5Q7{)|T#?JJ z!`A)B)PT7kl;1j3@9`|6B> zjkeE_e62%;4|?cv!Ybi(VTe>Y)aA_d|I#o7oOY89%Fly~caA86LvT(U#rBt;eALrp zzwzchD!k+DWv1uo<$++i%~pH-YB*k5*6Ymd4seD~c_FYSTJ8&Oy;3>nu_p}}P%Hgq z84n$%oc433O6M$Plz>s51#is1qD@|W^7eUd`L)P+lHUPBA7nZ0jD#eh{Iek>d+;Oc zWXo0z(PXF@-D;itEc}n`x)Q(_hrHAbMTIm)Hv&0JZT<}yMswuiv!Q*zL`xb#h|yZ) zQssADbi5mVsVxFnj#0v(kFWa;}-cU29s*Wh{H(**P6XY?K<_XKmW#@BuJM^f_UhJl zFBFsdIcSc(`}Sm_ zwWYs|i@DX^ik4PR`d{KcP zxkz-$?s?s>SSSx`>haSy{KVq=52sPNclryKTe>-JU!wXosc>vQRin3f@=)1ChJMUe zyN~K^Zi?*R5wY{4i+muAGEakt119E@Gn(u673Z{tL)68?w;ELFjATTgHh(XdyzmV? zpFV=;Cl9;5w@uJ(DL_a31C%bS7{3F#T9?tGJn?1yEG$kliQXuht78kZ4{+t0PLWOb zg>SX0VBAskf6RNQlxeQcxCU~*e2g8k9I6uK+Kk+5rZ-=};uwpnD-BahO0R7Zk6ISipETX^O92nNxaGs|J{+8q(aT8sWBYmD z68L1^N1)5lyi@l!7UyFk;J`9};U(_80wZQma>;}Amz;950&FLSt5PH@}5Rhtnl$E*=LyupyHH$rC3v7k#ov8 ziSO5ntW@w$o(*&WpnOG`H-IiU;_Ashv&ayyj8;Sx6Kum{FL_-|Qrv%Gx|1s0q@?uP zG3yFdl)s{#p}QTxa+uZMpYAyYmrCb?X|z^|a8QsKO)0N_2amfzbP#RyumkA~uKgX< zWe1M%6R1^%tK`Vh-{>M7AdZ98XsV=*Ryv8dUAEU}u&-A4s(AH}&_X#8jK1k^Wvw6Y$ zMJ-ycYduGquVv9nuk2-l)RLy^(onk8gD2S<-qx4HC@d5YLQYlOEI>)>`14qthXEHf+%9<$Sa`37dhR? zj$S+&%H&c&_RpvCBh$n^7S1dk*qKa?rcWkfWsgUg=#o>Gf%r0I+VnUg2II_EU6@SZ z>>MXhCq497QQ+5T?c_RHcsu|B0P;IZd7&upjwnjR@r6-7oR=ttv)Uj06F7?ORg#$d z*-txJxqY<&#krVVPlRN9POezq{Y+V;V%C|F4+6bO(C^dM8bhal2-a&K)}lyvqd~n# z$EtNUnU*S8=!FIf>sxByxq&5QTyH9PuEZ~d0b=Us)7(x(k!jE+8bN%(cQ{6=tnM*; zGj%K%ouIUfC?+=B#AJV_>n`M3)VCtJgvL*(oxSycX%D62w1(nK^IDHuDYFZe)VL>3 z?gg@^8fscVRGPZ(bE}>|JQK{fhanSX{V*w?1xl<0hw+6=`40v_*?riqX2(a|_~ME&0( z|Krt&=td~<;4tz;(pe8=m)dOUU`M{P3&%t3E32Lsbjx!Rd(=n87iFTNH&0N|COyl$ z_s`_i)8fflg3aDy;#bK}*>6ymGyng+A>Ly)V6IXiedP~tWeR_fiA#jAC|(&pyi?9? zj*E0LYqwAMd-i6UU27p*Gl}nf>7;3pR|Xnc_w4FEhJ4m6{o zpNcx#JRKhRt+zN=(15!04=ws}cl=X-fmp920W&Wur)lme#WG@fzWkE++T<@-->$SOV z4t=6X9dli80?ce~``TK`2SY599$(t;u1`>gMF*OJg6jEPPaGru8B+MRX>0L~QuSIw4+lYC|sm?3%b#?x0yWs8YqS7DRM=bJv=@Dw@67Ed8_X z9(uq|9|i9&J3`HnGtc2X zm{Z+%yI`AlVhz=OmU%m=s9K$x*avbaje0n8O!>m}ydIfyAH`YHOB)9tCdyxq9!aD5 zX_Rq(y#6Bd0)#bmhI4mV)cM-me88gA)oZy*`~)ozQ@tLroTnR-%N4Xzs8V$0iKBo` zsv}+1O$_8Zm$k1`{QsN^CnG52hu(#%IH;3W+sd<*VTbK3)kn_imB`UC$5$fLK^7|m9r-k4)5x3+~wpDAO zu}jb_TmGT0*fJ_lnoy;>W+SqiIE>U7-DNU~bal>y4&P)70B{&>UgYXpqG29@Au{OQ z=a0(Kknk}mdgIs?@DvI=owyCb|A=h(3XiU@ z>p%dhRLV!jdinw|)CeFCi0uxVHwR+Ld!&Viw@`pz>GRx7rcE(^Zu>|=6gPazOs)nl zPE5ABVlZHrF}6bB#^i3uFI3^n!~BcrQx3+&y!d@pcd&)%XzKvxf0)Sg*Hw+q9{m&3 zhGB2M{{hGyww72GVVQCr{B{R4G8JlO;1bN-x-J@^7>duBNZv{WRODwXzt|!Sot--L z9sE9iUF)N5NLbFtSMk0L%x^LeMBa`+Z1-_uOVWGH3)(L9hAvsE_324~Wv!P%FHDu? zdYg8GibfBslb>6oWN2)_gYd}w3UMaXmgq4y|@67iQ9W+{V4#<)(@oy-c_l4;=O{74e)m2r8$dngiW-{m5Nxk_|?B??Ozmlrs3rb?H;G&jtr-XoD z^;dLCF}h2&m^5IDgCR^g(*KR?&g=9!`*e@uiVi*E{vnQ728Eq3SpA40tVnf-PKYWP z=IvCp3~sM(r@`x|si2d`JMAhD5Hn0v7Hxw5onx<9kDFpL(AcV5j3b>ugKM~)|48_@ zwfLeGoq)DVYP*VW=njUX5XFHy>G>?T#$DvDoHts!OSNl@`K+N>b;8$omr!2UJ^H(< z!0O3{zyB*YC$hxgGc_g7Ij}}aI$1>8os+ZnWgmI|ex<9r+?%uu+0&rc{f=!+!N~^UnqnHyxK&p}y-zcXvE75U&<2CRyt{LC3X4#yQmA(?*gBgod9~=k zs?z%vG^_~`n>Z))8eBTwYA$x zUb4^iS2@_Y{Va<&*B$&@>JykVD49Ma1S#6+DF%(hkL_%9B~7&GH*#b zu1LvBi}xzw%ZmCb4DU5xSJwf>%**s{waE81w`(^}6ZD>92Kso$<&+f%eWY_ep7Zi` zHJ!Kc=K03Rrk*UDm)nCEL zS5^OcIk<>tk+H$2pyeCE^o%}(-`LjoB)aznjC897LhfuAChL_-iie+{nf~(MuA6$_ zH1m1lJ26S<<8zzY?2HFz^}D%pT2nHDwKv92_ZiZ_;G6V!-Ih!BVu~&&=Q6`3UCwEv z)KJvJXA1S(Ptwr5*RMC7g^YG@Lzl&+#`^u&eww%Sn)~~2ZWU`kNJ=fSJ4HaTt>pN# zblAwK`VY20Eq`=Q{D?y?T8cR@w`#IVih389OK-eqk_xqUom~3{T_9b+oy=m@zzkX9 zI+d<0=kYD~OYa;x-;IRAk?Ox6=qQ;*ZwY?!Un+}=i+JrhTFe%N!%*)%u7RKc+s%A}*e6ltkO$y&ItipW+g?3WFSuiRv- zEC-My!jb}|BB$Owt@4TILFPPZATm#284if$gC#+g)K>wpLE@z3`CqfFf zN|uNF=S$|7mmrxe5DJ)QJ@+_+oaberl z_^QgpZhK+xa^$(p^TA2m=Jkx(g7$K$)vKP9`i%#E_^C>-KqtsLVygjSo_}I^e-`sa zIwjo&J@sg(d9$m!MmQ_zwqm7|!3OhjfZ;j83CRm*$EQkycx^!3i+(+tLweVU6r6N! zF4OS-+3&P=7i)^PV7@;M$H85HbB_Xxl9Tz;Dovx8p9Sd+a&mVBmG*$~Gqk>DPgdoS z2!R3stVe+sU`io@G%aUM(HT04mrD{j59wa|s^~Ju zcm>aAogJ?3XxUf*piY1!p@A2SX109i$1M`$sTPZq7uC(@=?o!5%70AuUr}*;$#vQ$ z{?WS$llLBmM97$+`Tg2%uN&En!ne=3>xARAyX%FaHUX05KnKfi>7jomW{3DEJi(7tJ6>Xa?D!M#|~3C32o z-D_T3D>lv{itM40R|cmW+}q~(@E<;Vv04WM@i}*68}N7BNA09mte$SIc&vHO?{I40 zgkT*|tdc>RyR}Y3F}fHcK!$}4P5uMkqFX*B00Yi#;LL(Upc~2GWZ1zpTP1p_{njqa zcH08d)qOulbd8GUw-iK0}66+WZok zgDzC8x>$$=b-$IDF``~e{Whp-aEkP?<&rz!EolyvSJAp&gdhMEmhSe!!pU|MY)Jz` zVHOATKrnn#7zq?As2wR4{ZeZputwM-u3WdJ;}$#)Bgu{6pZee$ldG;lTOsNk@xcYv zymo1=$%IHjVU(=Bo*g7s4>HLh&pj)z<}S<@Ur@^Ulu7=)e@8fG7A}TY6qq0^W6GVN zPMZvSTyd}n43OFLx+yN|%Mz3vDL=~oAKe1k$kI%cyUs-^8)hSdK>{qVN|chlFnip> zZi>UC`~t`l2u&WLCH38OorM(swW5ht|Ap)7mc@c(cn)xV+kXPdL;x-TIe<^wbNe<% zPprOhu9l2DIQC0-!!l7e2l>4A)Q7I}{Fu2!tiLJL-PDeU-QZzG`G{zPB@hE;pG@{q zl&lBZ)P1x4+FcncpFY&$P+n&lErX-(!kID!@Lv*RkAa0sg<_q$_R^R~!U?W3bHlG} z5`c|B$T?9^jXkY!z?2cQJ!C4b0T_=SB`CSfV_T>YAr?N=o~Zqv{r@(#vTCPKY?2Q^ z<*A7z@VNa1^x$K?jS~b4&&Iv`(yaEOl#C_{@;=^XJ(fxE5teR!FC4p0l3S&KZqZ9f z8b**Kj3|^-BM_qauDJi;W(*Z3GbWgUl@3bl6>_O4n7G=Ykxh;5an&2hnsNA8>~6%M z26hgWgJTXSWv|V(3>0EaZMf-8JF)J;+8Ad=f^4bN&%Y`-BI?1Ghk>0=;8zzE6KAod zBwmU0O-tE0ot4biU)}!?RaQx6LTJ%jjn&$h(SxPjyxYvV8Ss>6j(_WE7YAjq(XH@P4!!?#^+bruydO@ z<=XpEOfX#2GAtRMhU?1zi_o=~lMdi)_%`aT-#3+U^bi#;3TxSDXO&@oJe5-q2PX7s1&N;TTx>7ajm(q}8O6+3sLq@F)-NjYU2QOG z?6xBkm~$x1wX(DBpFx+o4Lg%kNeC^V5|JQWllIE9%dG=e%H~JIvP?mgqokA+9vr{D z^64a^d(%J{E^RJ565Aj-F8=w?7N<)U)wf`^?@p`-ElG^+)vQ-0)-h3TWeg6TAMm14 zF%^G5j85OJ$CbIWjh^TJdpCOxMWp%;mLgTLOZs3ETyem5eU0%=MNm@Q%ZC=;agq5$ znEyk`Q7}+2oWzSY@)tHM`0AU&8v4tkXF+pzX)tm8-U2Y%mdZ0)l>mqVQ(St-^L?r- a&KE5^LGb Date: Tue, 29 Oct 2024 20:18:22 +0800 Subject: [PATCH 4/8] Update data_test_mod_cat2 --- data-raw/data_test_03.R | 6 +++--- data/data_test_mod_cat2.rda | Bin 16396 -> 10145 bytes 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/data-raw/data_test_03.R b/data-raw/data_test_03.R index 00c338e..b09be03 100644 --- a/data-raw/data_test_03.R +++ b/data-raw/data_test_03.R @@ -1,18 +1,18 @@ # Adapted from stdmod set.seed(41543) -n <- 500 +n <- 300 x <- rnorm(n) w <- .4 * x + sqrt(1 - .4^2) * rnorm(n) v1 <- rnorm(n) cat1 <- sample(c("gp1", "gp2", "gp3"), n, replace = TRUE, prob = c(.2, .3, .5)) -y <- .2 * x + .3 * w + .4 * x * w + .12 * v1 + +y <- .2 * x + .3 * w + .35 * x * w + .12 * v1 + sapply(cat1, switch, gp1 = 0, gp2 = .4, gp3 = .6) + - rnorm(n) + rnorm(n, 0, 1.4) dat <- data.frame(dv = y, iv = x, mod = w, diff --git a/data/data_test_mod_cat2.rda b/data/data_test_mod_cat2.rda index e6fa00fd8a6dd15f0f4453118d9a67cb64c85032..e49cec3854f38dd1229b9071f10f4fa39d0f52bd 100644 GIT binary patch literal 10145 zcmV;SCtlb>T4*^jL0KkKS^uR$)&L8h|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0 z|NsC0|Nr0_e*g~8JjZXn*12=q*}GddS(T!@x4gU0ZQ<_S=XQ5@dULgSdLG{0yf=H@ zyW4xb-rL^STJ~RixbM4God5s;0z7q3(NXh6=Kul!Oqg4*7 ze}OVI0i?k&)Wtm~rYDe8{Z9mCGisTd6xye>#XnT`ngCP6X@Hdcr>Cf48W}c}*w9Qg zGBTc;8YiP8AZes#rWHR!H4KaeJx|nWsj(-fiXBanqe<#%nrTnaCa0&V>S?tefi{sm zr?oXarh1P-PYP|6%6ecLG@g?gJt)IS$>M3IgF`8agvxC*LSTlPn;J|hreIAclhH;p zWb`IP!$u(4P{lPoLU|_AV9OrZjp}%_)pcC#mFU)Q436OqoX1 zFoq#d2AEAXPbSh0CzDgkc+oX81wT~so|9&Vr14D2iKe0Ic%#$OU{k_RPg4_7u{AcD zhkzzcM~zQZ{RWyu@-mE`qr#r2r>Q+pCQNBDC!^C;8lI4Z6GX%rHBU#V$+Vhe(+TMq z(s~*vs(NX()YS0OW`Z>Io{DT~HlehJQ`)B54H`_CO*GR@GEM3+8x+*_Jk-gWPbrhp z85=~xdMT&0o|6fwwv$HH2&jVs0EAQ&pr9QhDZ~PE@%MYUOb{IVY;(}+5G>#z7jmGb z1p#;Prq&b&0H`RNF;D?Y0v|Nuoi7qdLWGn+fT*aV0#u4Z0U)9RNGO2-BpbO1vJj|{ z2|`Fj1QZBK1Qi4*5E2mrAVm^U0wM%b5~8Fel1T`FASx&k5dswiKv6`J3XuX7B@|GA zNhlB@Ac8?8l0qRw0DuVyh@udpiiuJZ5FiACNdywbBk`TIwA&cl6kSKjarv{;srwto;Y z@I{OCWsbvs{s`^WqtP`l?kv0&#~|^!GYh|rZW754^_ANv8)(5ci_nwhNCVQ@wTM-GmC*; z$MYKUGw<%>cd40CQ7{C_n_1&9_k`;956jhP@ibK1+6%iLZQ5BmgvN%#ihU3%Bok<(Mx|ZS)UP(AUi6Uu7v5H`2|8{75UlJNgE@<-(a|tk47Fc7qp%XIA*fE@wwetKa(YWEc0tL@QN9fijFgqm{T@3oT)+$<@LQS!A)MiO){5= z$#0o^XJv+#Vxz4_;+nw~_XArzN_z04kkHwcFghHpx)}HYtq>loi|UUQPLC|E6>`e$ z2iJ<|m1;Y<`-CxA|2#El1{cJuhEH}lzSnl)3j@mc**3VwlXl6u(us^81&9`mvrXr# z&{rApnjKHjp~K@{lwo4oG5jR}qN(#P^W~}Xn{vmd43G-*TIZw(b&N zhB4KHyu!NOa-^4q#GBO3W^49_&PdX_+5s~a`L-G?*ltd?X3p|(?#ze*<6KoON6?|L zo5yNxeInC&ZlV?vNJa{bRC!3xuho~jO^@D)sn&;>2{TXI8wHZF%dYvyTFmx0^xoa4AU*H!zkDna!7mL;n!I&dm22C%)Ao$`+jg~oN$3bI zw%8B==RoOW{}q;9Bmu`$XptFaTiux2@8aB$AOHY@06a6+2evXNKC$|0GRiV}+VdAq z*XC@k9tZ$}003aUO}|$-TY43&WsJu6co~{yuHlU^RBh}pj1Gz$(=)*7E?|kLLl)Cxv@5IkbWJi$V2`9@IXqLNlC)Ot9o%o;0GQSeHCwGRfo}4{$(qIeG zmqu1vYa3;UwM}?l9^k!T7W3%R!!Y*w_-cyA$uI%r;$p!}#u7bd*MijEJL2q{HulSo z|DbdJt2!sEIP(poC_IWjmia}^Kgs`+V%j}l#zi637K;!_=4z-YsS3!ZXCuk}*-odT zL5<&zi2dgcnob=7XAci$yp}L~rX@=c>b}{Oz=TzNH|*{{2M(X}bV60)a8!1OLTB%? znk8w3-luJgh_^eg&K(Vj8mKv2QEtWi`95_g4Ko+k>wP=8b_MDU#rm8twV(x40LQ)%&wPx-J0 zuHFWO%lBk4?&>Hx)T?Ve5x1HYQ{@~_JKtOU^3BN9^#h9Mhr9^^_gcd5`(ejJ+m$%3 z13ZQ5xHR&;%K8s89Sj=}YA+5Q2Mku7P`NDdRR{(<8|M9PEnr83>~{8be`dRIBhzJE zqNL(QJWR$jASY{>I!d9)D^d5;atUmfzpEvLKVB_3$GrT6_2HbKCR$1s1Cl4di7aa_ zk$cd=_FIXuOy}GGJmEKPvK9>T$2MwtulLG5X(Oj%u<+aY`{3ye<360Mbqk(NuXJqA zc=`YQc2|w26wE-`lo4kV(OYb@2^@T9>#;GkrT$-X#BprdTXAcJA=Qs0ctiBue+qtQ9k>FO-ck^0(`I@i*^?LYK5A>|s_!(tm2~W@hKILN zJ~F$&srV-v>=MFK(aiOC&u1RATVCkID0AIpEIE!G)UYeO%?*D}VrZ*iFS1;QoZ~S3 zBHXnx_hTI2hNHTb?&>j3Wxa}v3d=;h0iCQMC+h@DB|H&SRPfzbJ|zZ4E_uRPUo3Fd zQ!$P2_YzGob|@O+vYXd66W8Qe1ZL8TDZbIg^bf9lXh5;ipj6FMD~S!k?+18CT+Y!d z!g8rI(cJq1V}V7K1U(0UyCr{>!aZiN-iXr>srgQ?J>2P%<{Yvb0m70 ziWr#H5^$Fm46TK(f7B-5J1We$W!9?3sF7Lmq? zDb#!zo}*(++q9z(H`O_jfXTzr6;b&zU*ZSfE!?YH3sZm$+3YD*S!NsHadEQQV>e|4KC6VGc`V18Pv8iIXn<*;lC1CX~YAN60tMIMp|gPb6k=<2p-2y6E) z$N8HjbW(diLWW23l4|R%gR8t&{wJS=n1A*0up5o6{^UqOB2Kp=VU{R?)EDHu+N{ZT z3Q#%A%%JF06#in{OU8$s$Iz2@6lV2y^i~V5^?TK9(k6?p6W2V$@#N@A?}D3pymr-c;!S@3%L=AmPX zBCKVGA*p%OC1N$NhqzCnhB!APbSJm?TTwx%AFC_h# zY1phv*pZxWA%@c_D2epD>nGd%rufQZ;YhE`4`0V>jLvOhI`FdZPMxp3Gc~(^>neT! z$AllkOB5-i$$EU0`*^OCMQ)m=dPzX4+ijNfzjtv+D^)m}dlGDOF!-NIx;Uz=45@x4 zZoVbX6J^BXqs|LTT^o`@3QwAfesCx&k-1oxm_pUv0hiv*vwxhfwaB63RO59R(|^W)L!T^qh&H%!c;eKp~X? z-~&_EeI~H})%%N>&t?D7rq$YBgOX0^GkyhbK?*jL&>dma)2R>Xb;NxQ4tp(=^KnKO zo#dtI-LlknXqhwF&D$Aw;gDg#fG+4=jxtO znxjD&d!s$Zw{VJ6^Jj2%>bLk#i;K(v06_o{0NIzA_#Ie{iwxR%PjbYhv}4QdG^plc ztbBR{m}jl)kP!QgW(dcWg4VIC^EyXBGS^_ubj>bv1bwb+z&aT;?AvUevH#J>mYSL= zy?3ScZf;4lKzQUibs6@uyg~wQ{P-<@#7|{_RZE6nqrC}JAb;5RXBRDZY~8Ij9qSJS z=Bu#e;e}hSz82~3YJ&BR9x_cjpbtD}>b$##Z)5f7xVhvp&E#>Opo~W2VrH7D zB5q3&mwmJI=EH*wAB3D3sGG99_Knis$WBYL8zu18-fn{AYAPKho+Te8<)%f z5HmZN4xMzCcE{+NU@u9eXhuIJ{A7-r!&B5Wi^w0nvmrp+(+6Lrk(UK~F@x6)nNka} zbJ|v9c*o_Cp*KYH>6^C&omjatk?ppRzXu-j?B2L(JDWyT=kYv0p_i#?pxIszm(`sy z@S0zopaiT?i!15yGb-$ z(j0Hy6EnvtJaucm;d3WF5$;{|PcNr?dM9pIpiMR7?bCP*)*xo($q*3xgO$mps^O>qhg)B@Ut?cvT4PH?P(wrW&MXF@lv6u6$ zjTbS^UcwOVBnHU#LVzOy_tD>3_|glgF@Fx!JkWmDmw$TQY#yflNqNVKEDH~$e>6z} z4_OV(n8kpm>@1`x?y5X`;+z*TER&&%$@a5c+CJ^Km$gCiHK|Kx>8gWi)j*c!^TgNn zgq@Ufp`g7i$S;p`+MY8m3L%Qp?+gpAF{M%g~!c`wR0!xW0-=4#Pk|m2X zyFw;T{~>e>Xj0trG@oVHC&bI9Jelt=PDvD%Ly|3UG6b;wZ(i>9<2Fs(Qih$Jfa|ypL3LKA5@u?@Mmk{evNb~(Bq{$U4!xJH&v_!!m%$MZuCNnQ2}=D zR7>vnh>372t(Dks)-&qZo7PDk#G-i?qaJJxXbPf2H={=(GFil>hDw_6n{wg)Ce+KG zDEeq&{VgKuGXT;#Q5*%bC8JuSm!k9nE6P(d88r94!-{lCP%?6%+^>9U2>Gj~mOYb8=YM^{lU97y1l!3jGz1`rJ zGZYuW-(K?mqIH3?qlmmJjoBe537_S^n$a_9WvSsIOY}Aes8r{;Se0)owP#`#j1W+H z{9eIxsNILOoUD%fL*pWH%x8@$3E?RDI=o~N?2Xr$w5_8ykxbGNtBsMw#6Q&4yvtUZ z7nqe@pg(&8NbF=Pe-UnEO98TJ;90xc?bO*mU0C|jkG=WMi3>0{;y(MbGu3Z!4)Fd9 z6Mv~|ba~eVo*gG*UOQJsDlPgAEVW50xGA>%ybIIzqynKm=jY!%9A+h z@wiQ#(Mqo>#CS~wIevv>x!M;zPw9pdaY-Tb4L0|qiJL5+RKWqF!$VO;hYkTPjpRB|hHu5Ty&E2RjIk5VM>uodtcx zQ_ql@#vi)2wj8|33C@BGC_lA}llmzxP~)B<>OS#~q3A5Ba~%g**Xt@ANP+Bi>@-kA zp@OuqUb*Fin%PRUc+fw~<{5H%;gt%}n+(_;h#-{$d%jF-v*)N0(2Ks_^e!mL&C&b~ zIx_T#an{Vp*fTII`HyNa+0`xdU22wkr>s)?Bb*jy($s)-%>1UbsabB9f76AH)S}IO zp7klaHl)^+RjJ zRZS#UH3582x2J(0lr@#+kYy10U-zYf7~XO+;=4++;(ioX|2@jwjkbfdW!8gL8FvrI zqS(f))7VoxYzF;dq;8aZ32heu4%3dlQ2tmcGCQkG zW12R!!u~8?b!Otc^kgBy1|849uiRo;g|o?3)YX$TDwmep#L6iy8}M~PF0baCQfHhr zybOxpR);|SyB=Po4x_xSWSLKdHJl5so<;5fUj_buGy1`ViYiAGHveC$oT5S|Ca?|; zjjA*8cHAH44r6P$l`jXpj@aYJj=*0Ofqke1TNLh&Rp(ga@)$B&{oJF~(frPFL+)~P zZxQOXQ%YyuGCa3)(cd|PUo%W1GLUTfU~hLnP11u1`<;_Odmk1B9s%-J2#Zujct^GZ z7Yn$KH_Gnj=xImo9pUd6XJ?Z#AXqWAMiTEP#39h;m=0;C+BZFt!{eXd&pkMXM&qp>R)le7?SLsT zMSv`H+d7l?kU3yDLX`bKKmfRRwKK-AX9aO-i6-uB zhQ-(r0c{2t}FeiLW{8~P?*ap!LKl#BBE;OX7=AZ z#?B)5OY&M!PFwuf)4i&(7mSOQz}>>*qvsnKAyH`aUW!xQT_^K=^7g933j7tB%b1&A zuAMJL5m(tw6PdgjSAf!hJkTW$GR%k|D4|e=fXl;7VC1JJ7q8w`J^B`ALF6F57Sqmi zVz!L9N}+AW{v9hULk~dSA1*Fd^Gon3VsJKpY&5p)Z9qsvzi7$WyU77&#I}1`WK7F? ziPlYXzq|7DtA8$R8mcdp(3rw@ue6N09T}Ro`Nu7*Gatrf|BESR$?%D14$DNWw0T0x z+{UeIqhIj>M1x+vFXQ`a51|-z@(zKWGwgOfFUAuA7ydqTn>jY$x3|ow_hLZ!K9(3m zA}*;(YsGqbgocHh1DiR$?d*;nKi;W0cXqi{Ib*|xwR?4JjXB<5EefhedxJBF;PGq4 z|IGNmf~AF0`z!7mj`F$z8pXmMVnl}=gTWN;z8{eQ-n0rDLA$Ig3K z@}Rf5VRso-#lirMYm=;{Cj{1lLVknN$0O05L8%1^5#2@b-Zzth7=pK2K85&-L0`iD z?heITeUKDbfDtbh5z22`!G*b|dMl*v{Fl zfP5|{TE>&{KQ$Zf=!VF#YfegeVvm+ZXJSA_>K4d*Ccg zlkK~xbboxx_HJCV_Ku2aY3ut-Uj6!GX2Yw4*|Zl2cL=3l-lrA3h$@`U_#l19{bFda z@9new8%G&IZVPZ!G~d((0je}ZzzlD{z+=AWQakf&yo7e( znY>|{4-#u}lw~#vy5Ps;iM9eObHI$*?8n{k?70fsO8%vF7sB*_N~O3dLLrP0i42F_ zW}ibKU4Mm6%83ia;z|W|)grrrLY|HVcec*1fC%%cXF%0?l+-T}*RJTFxRuF!QCrj? z`sNfWMv0FFcl%ZkrH8C-NJ-x0C@}y=@B*u3yzzue2T(xND{U{rq;s`_sT_Oa^XY`$ zGDH_;ACMz4Gj}|r$o_&_uLMd;l^g$j9EIiD1Vr=Hu^c-Z71O`=Zi%7j&e|Tu&DKuH zQQLbSuV3Tzc+of*cX^DY`(@neH#)oe!Wo{(Ox&P!f<-SwEL1=9l4G3O%33^XSS52d zUoNR*QceJ$kKK>A?TZGPK=w;22se22bV6JmZpRSJ+F%am-PBKYIp3rY`dl?q1bPFZ z@B8deHRlQser1!{9uzPVrP7qf2@sKf^_1_;#tt@QYk$JfO-!qhE$XGm)1#{RIrgLU zwOsMV9jUrz^{V(rWyhC7D?=Rzb6MD$mwIU;rQ(R0EPlJurdB1HiT4jAmFQFw?W|rI z(mta9Y{%J`ucs5SR3T{@J>ax+m2D7Y<%vjg))NQ18#{y&2*8pN_PIDPlm)tSNbs&; zSaOy_6R*2IrR{Npd3xdUYgRu#a%SqE>3Xy4xyZexA)FpZ%lWqxr}?+B-o#$xQ|`Zq z279abA^_;2aWc|XsP}uQ_ZGj!=T5_gjyZ{NGP~R%MYy}|FT<8F3_;85gG6ca3$^^X zwYt_q0Ch6ja;(-r2av0Q08hh$QVnYjI1_y1Ex17ur19u4g3GpAT$2m5UhExl1b>AD zP;`r{!sf9{6Mbl!*SO0(ET6*8Z{?%Vn>x)qO^FSif0y#v=NP z?{Qg7>y)L#r9?QS%S<3-?z{M2_DY92YL0MJ)_1VNat;xBPBJSTM)rL;z2ryo1@3r_vAK zTP<;dwlWu>xaxjfEnROCFKrAa)_~iPAk4#VxtLx1fzg{a3(JtX--&hFTrXy$gyt?o zv6^7KoHgHNWMujGK3}0@MW+0oku3uG+T!eZ6cfG~rs3JTLfBaxdRP?ByOHyMhsia)dzvS}o4PgaX2n%|jeaXLb zb!Qpnn2sr+;#*tmH0eOsqW2V>qiYT zHk)iKyLRzO2=%;v{tF?$pULBY)K&%8#DkE=Ep*5=+NNf%4Op4W za+#UsL}VF1l7GIrv`w6YQo%Hx`zR{YfaM7q0@U@~_M#>okY#Q1b+345xL)HnkK15G zoF=Jbbe?j?fA^m0W_Vt+={u4K<3(1wyLH>)<2B(*o;79$zTTB;3IyC^T70&!`ZhlxMs0;sL3gI55?CqyDB#@dv|c`8|#}Zu`OJgzMkh7 z(Oj{RG~QA_?39@Awi4I&iwVt2-w2e7^Ca$C|5fYT$VBcP<}iYKkM|v4P0wqV_Uy8= zK#MxUfMaggJ?9FI*zXIcF93TxXH}DA5*8{eZCMsK+;+Kioqn%L>BkfI(#ko_#d_ z(@O0Y56$;B_HNIk07`xHutwy!i>sGdVJgwhqS#?a>Ke?AN2f_>_T;ATZ~L!#C~9#Q zmU1-FKgZOLIrvgBbvu8Oo58HvydJnxbU9oW=g;ZJx5qMmdiM=zi^jW6>U(R85=pUM z<$K=KD2Ugb#Ala1INfoNxnQmFEt$vum{v<_r9K8*<0$*xdv`P_91)$$^0CbwHVSI{ zD)~n9{Jl-n20ZRyfS(F@hDt||fM9osk?TNJut@uH}%aiuld>dE+Wr zO<@Q$mR!GX_p0^Yt}lt2INLr4FH?}fI{G{Nt#KVRS^T8&-urKaIU029uv}Fe2MpVA zPJc6WBfc(58LbB8_k!-x)p1vbT0$$J_}a;2Lu{E#sTzaq9Ptj(8;SIbIkn?-v9tNo zEzt#NK_d`Dgi!N*9z}%H@9D06^spv)>6g9aqbf+&{g9f&9(@A4+?p!`(?*TiLT(w< z8j_|QMyFWFcHgedUSTq3pv&ED{Aa72PkxbW1f4EQMU)0az};;wn})y{G~C0n{cTeG z$9jvFGaHq*p7@mqi8Y(+FhC&6?DzNVBj;P4&+KG6IwwQ(i}CEM5|UqRQwXylKyk$f zLk_>Cib)L5io(g64kGZDw^mQGf(7HF&vA!9=Lfs!X|7E%sb zal1c zXo0+gs9I+dQE8#o$lkm)HEMl0DQ@!&}9A)5^JEj*KL8-_SmY|N*00^7KIqL P|NLFa6yZWc{*?zF=LsFTJ}3 z000000-rzt15G$|$TS+50%nBM(qI5g0$|V;KTHB(00J-xfHasX=)i*lU`!KCYMvr` zXkY_DFlYq8$u^n;Kxhe#B@U{e?G*8*hE-4U7!iom(V}@YO)yO`(3u!1S>w-Q&ZDa_9jeeQ`(OqrkIrdH8#>>PilIMnw}zUGM-b(rkSEH(9=q7spSWWwNvm4Jq)Lm#-`ejR3dsr(E=tA$)=eBJprMq=$kY8JY;tK#*o+03d__41s`P%*-Ic5rQK$ z7$Y(=21X3Q5C{MO03b$SAi;t%002e^85tNd002M$KoI~CW@dm!MhMIig9Hr*2+WKb zk&&4g0KfndnV295%*Yu71`PlpfHE*(5t)$yg9HJYAQ%K;24rA@03Zm=j1id`nhX#i z0RuBLAPR`hC^IuN1Vm&IV8NN00w4^*nS%h(W?;<#jDSIaW&r?XN)k;Vnj(D!i)-;A z2tp8pM}!e15poPJgR|=W*jY*$Vnjy!UE0on^Zp9y_u~4!h|gxH%<$~@-&8K0rdMJI zLpws9X8~#u`CHaHLiXMW6vJ$S;_PcG$^;pptL;3o@bCB7RepvK!EGrj)(Tb3hEi1s z91!DPItxL+v-@Z456DIGX$f~+T zvTJv7Y+mMGxzF+nxl^*KZMr+;Ila^G(l32cBG2!FHEmvdoonnig{G4DzJg*7J zOdp1aVW9tS6791d%e>F5i%<&ndv?qbfQd>2$f_+#S9T|TC&vINcj|u)#vH953z%#P zO`e*KNAGjB%0DDQPM%7hfh0bgN*G$SDtS%C_F40Q zb(bng4))^f!E6lDi-7b8n)&QqDLlVbb~g>nV5%vpW*!A=r0)pxXSgPI-eNX>bX&S` zi%r_~>{_6b3#Lt6dSSvRw^V^)3qm)hO~loxi+0-k)z8GjnKt#_l}OFIJv+?OE-lIt zE{(Qy=$Ynu!g+$6^!}^dXca&K0t5&UAapm@DOQd3&XqVCi}Bsa8O3}1oF1`yyAVNZ z_KOa_^-PcpWoCUE)M^y~rmtCWy9}Me^v2*1?JWVhBDQIMZPOI3e6YyQw|C5w&`UAZ z2u?KBpnDET7L~}}<^ItPjmZ9`gWE*Y%aM)X@Tn2k>{E~4w6&c3z-dWJT2Q0Kfylgf zDDAl}tBx61w++bPmpV-?SA;M;8wl7J$V?gt+db{hB25c|11*{I{J6dYgL-s77`*WN z1cdq=+FV0UGg#HiM7fG(0CIMsD)A6J^Z~wZ&70NZWX9%0j7wb5A~#AI9BRKc9PuHk zstX0*s-zywc^KVkn6Ngmlejr~1@#xF3H^4INI>yZ06;iH_q%8fXv$?wOD=x!@Tu(Y zkbV2MdH&IfF)LI6>k-dqJ@!^bMES)LEnrQE0rSVGs9ta?!#$5(etQ*uEpYmN;}yEsq&?b2BaPziMjLdpHT-lP^anWUMqYq~3RJ&Yt-MpHkm8;jKGi4zMc z!xhP#j5>GMh^1@NcJBx!H?j7e)+`)jshEQgwez@e zNTdGDYenr0%`KYGN>hw;j?+Z#D`DPVXlOdPaKr=jc*b{Yg_{(vgHf(1T+)+2i zi#)VYu2hUi)UOx)8ILZaJB-&?5FSpP3`4y_V<1;}?yoPcl&9DoyMP+KF+@IPoFSG282}QBU;QYb7(z@SH;l zVpsGfSnkuP2h_L|-kPh-X>wO+3(^d5!ddX4`)Jlqx2V6F?u<0t{e&B1sb5xG8B$Ho z@HX?hf7K3m5wNPDX&|*p$fV%+Owm$w)stNSw;-SIn-!JM175P={RFg3>i2Z9)q`0Rm}q3ZOXL5HOfB9@p%gSp&k1-CSk+uX zrkTD$AzPX0*XJ z@GG4jV7<2RgD6NYW9OZPElz+d{fwx*|>Ma1pj#aKfZBjh!XH*LQ= zOw7zhlxuh+m05;t^PIC6|D-1GbA5n*fmPGuE#WU7aygYC(j zS&RhaAIl$}th`s#6B2Lm>?2n%;b{4X@ta~7$Sk|!>{%zkSn~eT6?ie;Z@fB?v!2^G zu;cOf)0zfjeR5=OAnecaiN+w`y34Mh>qxrWzbypn{!N{>IW8}s0SF?K?(2S&QV_3q zmB{C7(~!D#bnx)Y=syzrUDFDeFw;XMb(*xDL0uaU1UCa=OsZD?AAODy24Ge1Uz*<9 zaR1czJpZ%KLg?Vt({z4(p-n~sQAwELSp`gF2;uM`6KOvx$G}bJr)u)CC~&YK^;meH-;Gmi=wVE zUy=8TRZA;y3P9tETSh{|deWQ2aB>Wm{fA=TC9awds}!(WWHY-6i`t8t+8Gdq1#j$D z%30wTfBVR6pz&5W33nzRl;_=u<9J0Ohx`ja@NN^KPIt^``}x1c+yPz$7z!J5()5b3 zo&A3v3gpC_3=Ple^Y)um8ik&LFV`S$!-9I(rQ<;zCfCaKo9UKm*>O;V5-$13!e4tw zYx@6RN&B}lhKp{Da8~y*km^Drc|#GF4->dNP`4s<5pMyv=`hW9ipDP@2y#KnLL&q zazCE5#U40Q?w9!b(*7@c*Xeah3<@(U^!&@H`0jSuQt!eVld)Gt9w#?CZN9eaNb9WL zEP9B*{!jAlDK-&cR{bCV1i%0jMV?pE)$2FXirvq|9l3NQ+Y1g`PILyuHn|u%2dAe- zmRL)nF?(3tI?ij8lS4y54FCZ8+dEDnESI06r5ESOqci);ZAb%HNU9IFg?GKWl#!Ex zoqF0h{bl_%MUJ;BnpO3f=MU#-D-aL>0D%Gr#}Yy^`pdKtG#nSTDch@O{0$uvq+;e* z6Yytkkp(0$(E7r0S*YZ`8xlWC&aB9_Nwb_nwk& zO_G3QZvbv~z6l9FFJ@Uc3JU9>HZiLHfll4pqYS2L)e}4b2tv)s)k8|B!tUTKa#vgx z0<3oC&E}FtdRlK@H=T>}OgYUL!;f241RBL6qo37@@bs!;JqQLbrsug_i0#Fr{nw31A1kjZs`Nyj zZ+@vHod{@b7Z5p;|1WB{=6N|cF&wPTLqN{y39<_3A5>gdwHC|G+^n-y$Rdi4aja|{ zQZcVMlVY0Et%r4SHO}L(4}KUuW2}TJb0jcoqrCGoN^|;q0fDnJ{wJn3MglIPMQm#N zEw(Xf-(Hx)van-j4;+%cc|szIG<=XOG!Juz3rlxNSg>L@5QGF}&0bywM57SCv`^f` zi`%X`C~*}pK{R4>(^Tu1lOzdf3gu26(p<&X%oqjO`mH4e?$ zLQmfAgpf*8%|E`7VIUC=S@X2Ev7$3-k1S^Dwx2tPHJW@UlGGoG&B$tvlNQRFUX2=u z0Wp**OTIkN{>fdXg(Hw339&DpZra5_u3kX0UFWUKWc)yaV(zVR^hUr}P zjqEA-wzrz#k_8!1$_wFBnE)2fQ&hn9gV`*w2>f%#4o_QUOTGDt%(v>AWyg(%(qFs# zsrKgjnkOQ4acvKi?lRmuw4$TyMfjo3lkTivOSvPEdBi!z0!8Hpd{i}GCTD#oGJqPb z6IF13!X)fCHvT~xq2%I_o7|o5c%!`m8JB}%&fi{o4APqId+g19|Jf2GAHq$M$XCA)~kkvF5)Ljb!y zM>^tqG^9_m>8gjbUW=bHBQsfc$K6Zl7m^umLNJh?!>sw>Stb~}17|3+TT^0`(NvG@ zL$H5iXBIC^S7i#w+GuOb?44;Atkv6$N4on zc%zSFx3ilW--yIg3c%5ShZ_er4!2p&B{#(L?Ms!+cxnprO zTAJP6JZbNZxMI-k%Ty?W4vuz9WeAJ6(wgh73FJa0Jqw@aXmP#D?s*2uuViOGJwJo0 zE}#MelQeJE(F2g8zTQD5aiyLbye|^IFYt7B582y30@78aN6zkz!mr04(J8~3O|ViW z!23nEYH_idQi+DzhFd>xktnxvC5o*I?o$)By))J=c_*K?{MvQBkH@myrS03t?c{)K zbWLMtf`u@Rs=ylXkh6kXg-Gcw&8)z*>(Mqpe>|fx&!9IF?dgW=Ckb2jlIqc^5X2N|p;EceiqAI`gyL#|HcnhU%y)y3;g;3MrwE;ifJ z{;6jTM;K*&S*u5y!po0iT<;pVo3!WX@bzO{KJ;5DKtddRlH9pZUhPQEOajr-)=!oB zGb^AZ_}IfIRG)}RNm=3jsq|XZ!OQeY8R%DNiGekWe`bfz?!ZX(Y8tmJv(26bd|cvbFJA`LS`O}1La=#)9~9A2*&(#Dnu}t z?-QQZRhx63=kt7fv$zh_AAV68o6%i?gs&M@?Ycodj+)aKj4PF#Qbm8S=!?K97dFyW`F3kU8OWG(pf*J%984V5VWQgY`?s!12(te;dk&$7mbPo>{~dj*2{;dwpBQScP?CwMJ4q+tZX^87PH z<6^nAuI{|T679AYo!75)erBr-`76QBH&?cnOB12zE*Hylk!drhw5mf9*bPFZ_4}{z zq+PfF90rYmQ{|Yo_VPxQE5F9RG&f#`BXXTpdnFIE#S8nR#B?d%9)?ww?18+Fez9Ed z5v+Qi1!eoja7$&s@GUVyRS2xh1j#4+Sdu^pzoc0X;?G~pCfx=*+WrA>UNiL6@wQ)| z1O>AYBHxyTzwItR9mmYZsRFzpP)#=&=dWe1%mTc&NOZA@hUxFP2ePdvp$ z#D~`zeyQ$3vonrcF^B>Phv2AuB=gc8{9|W3Cw=x2>usT(+m~3g{~*S?{e(_LbEd_1 zT}{(L0VOAGe9f+f47ZrGloa?*jhW5!=irB;p7HJCOfTove%W3g-3XG7_Er23gbnSt zgW!c9nmTL&*Zn_OOl^o@v&$WNKbtP4;_xJo-VFt}_v{0#HIsnk!l2C2mI~MhOSTY0e)byr{FBH}rzWcfK;=p!JA>@s-!Wvq1&qJn>3(oa0}7Ql^n&R>SAQ>z*8=FfEXvc}=Ph2y&mdNq z4_QvP$Ub^*U}M(pzn+2d5kIY02&4074$dUH>>GZUFTb;1&p)v*VQ}#1xS?S!$YuqP zh98?{)%d)ayeVc`07w0EF!%$*sm^apH{o#8z|~KWpnHbFRkNwC*pznoBFUWA@=h_b zs4>AYXkr5OIZ4l|-NtG|G3TbW!ucpMQ=2iK29S!#R(#BwNx<%ZLK7zGfwwV~yruoz zPc7eH4`uSct?Mul zcfP}MP9?Eay|6Lt7n+>DV^*DPWsyQp?|2=vycAjV}FS_~HP_z_%G+XJotqJY?bc_Cx)a;m%G>i!!;JIJr4bXp*-p-`1LN z$a-^d>(-A9CU@JFYjL?KwXDmzy}S;o-%9+>E?s(>(>ug71!2z^B*5m&hO{K=-5TAy zI-LZtrUR{*@dZ3LS*{cneFBdNJ(xJOb2_zs3W%_ix%SzBzL1I82}+rCR%JBGb+ZVX zMu=vO4*31b2ZY~*Zw56m67%cK0WiB(I|L64Z?oewiR{Jbb`dM%(e~vi;0CNsG5l5d z72y}iar@1x~hx&xr>6kmwiA-eGN(ygfYAo134)jGUXHwCxN}p!S5* zUOPLqQh>_aH+5=9AOaZ?{pbIBb85D9JjFCT#5by^ec6Y&vnO`KmXON}Ql4@4xa=UW z-Xij5;(&D_eVs0B*Y7 z=y!$kdfX6~QR@|kQdsmQz!|!p5&)@^Zr`8_DkG97f{exTl}dPH5VKks=# z=|VoznbW4L6%uuhkKb*Wu2jgcZA#*UAYU>v`!>gp1V4MoI(U zr?-OpHsBxYCe*X!<_50Af*|V8n@h*alkTMRiumNLV5OMHs%BqF9*ApM1k@QPHt(5M{m3pcXs2@rg@@ueS z^uC!lH}FsKjM0IMHXrJJN{T~r%v!(Zppm`Sb$OtH?nFIm>3uZ3B*M64tEUg-{(FBH zXEbNAElesB$%i`$PY|ynzjZL$7wz28Yh+C)vblOBP~G)UmnUn8gD8y$8RmZeKRtie z$TDfzPgC%bPlgR&LK|j)qW9r|t7`UJt+k_*LfdL{CEr4mOlx@%aY=^M_`230VVKR& zoO)Lxz_QNI4fBf)(9B3pII~5G#5$H$O`I+qLOAep zHL~H`pD}DLBsLbl%YIK^%{~Uw{-3m}g=bHznS834)+Fb$4)B-)HKX}rEiEDLGz|V( zAxEpoYh6lT{1x+j^uJ9U7^#zCi?{lDU;q;UCMG5(yKBwwC)frN(K{t@*LQqv=C;XY z1sMXcn0uA^4*O%Tv!6tn9UP+UWj7_FZk5Tu%X;_pR*Anc#rVk_%0|>trca=p9$v8^ z!DY)4j!wVzQp)8Z65Q0OO(%^PG?#Jqp0#JBh#cq#v zcd*fO(INY>XaK*Jbn49Nd4xDuRSbkD@^}2>#M`TTUn79hr8xeIuWy|RMp6Zb>S}+0 zY?dW;Cqw3Uoz_C<@McExwT(g6pS^)AL3u-=kX)J~toTS+1#ai{EYQO~Z@|3?r8JCx zhhYc45*Ux#ZN0=7i)bGE%)07GsJ(hxnT3&0X%Nb*F-tVQX$OWloH!nBmMt-()LBXD zA8? z!u{D50X(7}aKyD=2~D&Tkh1>1j-5hiK*-2#RTSQSV}|Wyp^X1@Vo!s8lG9ZwC<+Jh zOx1{g2CxJ@XK*}KQ~a;gWK%dy;=?rG1_6=EKhxqA zKw__X@QRd047ZrBm`(bLE92jH9U(|a4AX8t;^_xL>P0#^G>rd(lAI-5=q#{0qkmK3 zJX!M^{n2Wtn--A18&`GAPLkR|JPNjG@Y@EW;aia!zsZ){0U@z`3-}s%$o<0tetg)+ zxy&)ItR52Xf2F$YNM^0Vu<+PK)+<-kev%EN)gGQji|4`b@eqvF(miIRO%~O-m zuwn#GqDKc?m)DC^#(vSGKzj4xvN?B5;vaPl&~H4fdM{4O2))eZ&S*ih(RUiTYZ6~W z%EAljgHopVDjYkX^2c-?o&^fo3~ft+hL2Exm8C|h`rR`PPG9kZh(CG$3}9AeTLiWl z__$9Q^D#Xkec(kwU)`9e4IyJrSrK?bKCReLw5wVrNY_^GUv_X-sE^IN_*78G>Gc>W zk?0%ZD0yAL27}tUV)ay_Zi*C-s9?tDI0)<@+^1**u*XFw(j7cpGD zo-9ew^h!%*MjHIuq&l*zvx|-NAyIy)b01CP^z)8G3&E4qe z7?2sJooh_ep#Mh{5Vnkfr$^;jGs(-d?mnI}i*Mn2bRUA*WH8>HQ}2aEb!>jI zkEDy1np7_IFNGIF`=ZT(KSWu-$4~x^On>8E-*V|Ih~S>rznZi0%8(aft8w~$ zE^>_&8sK=W-W2w|t#`e~oxN)_NI^I_jWy$g5fw)h{d)b@US&UOq0$)_kPnWEkEL9L z-4Jj1c+EvGP)M3FZt;*(JD8_XN7~U=s$e{Scxd2@C(1FAx2c@!G16q2#NEf-MxQud zVK@D-LxyB`C+zi6s{Gt)7AC*#yHyMTJ%2j6gb~=RC!hJKD1O`BJl3u664N5#TrWMX z%<$fs@DVdd-81?YbX zcP4?BoOc@w=(9KF_n(N+81QsTPXebDU#mst>2K~^a6_nCfgxn;hIJ5QEDZ>93{c(k z39r!7^cV0At}bbXY*ZIvZSnJaySf!msRC7s)J39M;KnW}q%ZR#*A{J+SBh@;X#T!g z#;3EV#^T*_JO$6iD`gF1f10F0g&|xZqBiB-t^b}FQnk8C9z5$b3kyPZ{+`zBW(Ffa zRQu3nmc5swYwrpz93NSl{u7~Ae8N^$W_D0O{@lof+9M*51S6j8A%+RzpvyH(^@dtH zR!C!K$_1u^?$l{#!pP4L^fQRH$KI6t3v>KtRG>f`QqlW^WKuV-uno}5Q7{)|T#?JJ z!`A)B)PT7kl;1j3@9`|6B> zjkeE_e62%;4|?cv!Ybi(VTe>Y)aA_d|I#o7oOY89%Fly~caA86LvT(U#rBt;eALrp zzwzchD!k+DWv1uo<$++i%~pH-YB*k5*6Ymd4seD~c_FYSTJ8&Oy;3>nu_p}}P%Hgq z84n$%oc433O6M$Plz>s51#is1qD@|W^7eUd`L)P+lHUPBA7nZ0jD#eh{Iek>d+;Oc zWXo0z(PXF@-D;itEc}n`x)Q(_hrHAbMTIm)Hv&0JZT<}yMswuiv!Q*zL`xb#h|yZ) zQssADbi5mVsVxFnj#0v(kFWa;}-cU29s*Wh{H(**P6XY?K<_XKmW#@BuJM^f_UhJl zFBFsdIcSc(`}Sm_ zwWYs|i@DX^ik4PR`d{KcP zxkz-$?s?s>SSSx`>haSy{KVq=52sPNclryKTe>-JU!wXosc>vQRin3f@=)1ChJMUe zyN~K^Zi?*R5wY{4i+muAGEakt119E@Gn(u673Z{tL)68?w;ELFjATTgHh(XdyzmV? zpFV=;Cl9;5w@uJ(DL_a31C%bS7{3F#T9?tGJn?1yEG$kliQXuht78kZ4{+t0PLWOb zg>SX0VBAskf6RNQlxeQcxCU~*e2g8k9I6uK+Kk+5rZ-=};uwpnD-BahO0R7Zk6ISipETX^O92nNxaGs|J{+8q(aT8sWBYmD z68L1^N1)5lyi@l!7UyFk;J`9};U(_80wZQma>;}Amz;950&FLSt5PH@}5Rhtnl$E*=LyupyHH$rC3v7k#ov8 ziSO5ntW@w$o(*&WpnOG`H-IiU;_Ashv&ayyj8;Sx6Kum{FL_-|Qrv%Gx|1s0q@?uP zG3yFdl)s{#p}QTxa+uZMpYAyYmrCb?X|z^|a8QsKO)0N_2amfzbP#RyumkA~uKgX< zWe1M%6R1^%tK`Vh-{>M7AdZ98XsV=*Ryv8dUAEU}u&-A4s(AH}&_X#8jK1k^Wvw6Y$ zMJ-ycYduGquVv9nuk2-l)RLy^(onk8gD2S<-qx4HC@d5YLQYlOEI>)>`14qthXEHf+%9<$Sa`37dhR? zj$S+&%H&c&_RpvCBh$n^7S1dk*qKa?rcWkfWsgUg=#o>Gf%r0I+VnUg2II_EU6@SZ z>>MXhCq497QQ+5T?c_RHcsu|B0P;IZd7&upjwnjR@r6-7oR=ttv)Uj06F7?ORg#$d z*-txJxqY<&#krVVPlRN9POezq{Y+V;V%C|F4+6bO(C^dM8bhal2-a&K)}lyvqd~n# z$EtNUnU*S8=!FIf>sxByxq&5QTyH9PuEZ~d0b=Us)7(x(k!jE+8bN%(cQ{6=tnM*; zGj%K%ouIUfC?+=B#AJV_>n`M3)VCtJgvL*(oxSycX%D62w1(nK^IDHuDYFZe)VL>3 z?gg@^8fscVRGPZ(bE}>|JQK{fhanSX{V*w?1xl<0hw+6=`40v_*?riqX2(a|_~ME&0( z|Krt&=td~<;4tz;(pe8=m)dOUU`M{P3&%t3E32Lsbjx!Rd(=n87iFTNH&0N|COyl$ z_s`_i)8fflg3aDy;#bK}*>6ymGyng+A>Ly)V6IXiedP~tWeR_fiA#jAC|(&pyi?9? zj*E0LYqwAMd-i6UU27p*Gl}nf>7;3pR|Xnc_w4FEhJ4m6{o zpNcx#JRKhRt+zN=(15!04=ws}cl=X-fmp920W&Wur)lme#WG@fzWkE++T<@-->$SOV z4t=6X9dli80?ce~``TK`2SY599$(t;u1`>gMF*OJg6jEPPaGru8B+MRX>0L~QuSIw4+lYC|sm?3%b#?x0yWs8YqS7DRM=bJv=@Dw@67Ed8_X z9(uq|9|i9&J3`HnGtc2X zm{Z+%yI`AlVhz=OmU%m=s9K$x*avbaje0n8O!>m}ydIfyAH`YHOB)9tCdyxq9!aD5 zX_Rq(y#6Bd0)#bmhI4mV)cM-me88gA)oZy*`~)ozQ@tLroTnR-%N4Xzs8V$0iKBo` zsv}+1O$_8Zm$k1`{QsN^CnG52hu(#%IH;3W+sd<*VTbK3)kn_imB`UC$5$fLK^7|m9r-k4)5x3+~wpDAO zu}jb_TmGT0*fJ_lnoy;>W+SqiIE>U7-DNU~bal>y4&P)70B{&>UgYXpqG29@Au{OQ z=a0(Kknk}mdgIs?@DvI=owyCb|A=h(3XiU@ z>p%dhRLV!jdinw|)CeFCi0uxVHwR+Ld!&Viw@`pz>GRx7rcE(^Zu>|=6gPazOs)nl zPE5ABVlZHrF}6bB#^i3uFI3^n!~BcrQx3+&y!d@pcd&)%XzKvxf0)Sg*Hw+q9{m&3 zhGB2M{{hGyww72GVVQCr{B{R4G8JlO;1bN-x-J@^7>duBNZv{WRODwXzt|!Sot--L z9sE9iUF)N5NLbFtSMk0L%x^LeMBa`+Z1-_uOVWGH3)(L9hAvsE_324~Wv!P%FHDu? zdYg8GibfBslb>6oWN2)_gYd}w3UMaXmgq4y|@67iQ9W+{V4#<)(@oy-c_l4;=O{74e)m2r8$dngiW-{m5Nxk_|?B??Ozmlrs3rb?H;G&jtr-XoD z^;dLCF}h2&m^5IDgCR^g(*KR?&g=9!`*e@uiVi*E{vnQ728Eq3SpA40tVnf-PKYWP z=IvCp3~sM(r@`x|si2d`JMAhD5Hn0v7Hxw5onx<9kDFpL(AcV5j3b>ugKM~)|48_@ zwfLeGoq)DVYP*VW=njUX5XFHy>G>?T#$DvDoHts!OSNl@`K+N>b;8$omr!2UJ^H(< z!0O3{zyB*YC$hxgGc_g7Ij}}aI$1>8os+ZnWgmI|ex<9r+?%uu+0&rc{f=!+!N~^UnqnHyxK&p}y-zcXvE75U&<2CRyt{LC3X4#yQmA(?*gBgod9~=k zs?z%vG^_~`n>Z))8eBTwYA$x zUb4^iS2@_Y{Va<&*B$&@>JykVD49Ma1S#6+DF%(hkL_%9B~7&GH*#b zu1LvBi}xzw%ZmCb4DU5xSJwf>%**s{waE81w`(^}6ZD>92Kso$<&+f%eWY_ep7Zi` zHJ!Kc=K03Rrk*UDm)nCEL zS5^OcIk<>tk+H$2pyeCE^o%}(-`LjoB)aznjC897LhfuAChL_-iie+{nf~(MuA6$_ zH1m1lJ26S<<8zzY?2HFz^}D%pT2nHDwKv92_ZiZ_;G6V!-Ih!BVu~&&=Q6`3UCwEv z)KJvJXA1S(Ptwr5*RMC7g^YG@Lzl&+#`^u&eww%Sn)~~2ZWU`kNJ=fSJ4HaTt>pN# zblAwK`VY20Eq`=Q{D?y?T8cR@w`#IVih389OK-eqk_xqUom~3{T_9b+oy=m@zzkX9 zI+d<0=kYD~OYa;x-;IRAk?Ox6=qQ;*ZwY?!Un+}=i+JrhTFe%N!%*)%u7RKc+s%A}*e6ltkO$y&ItipW+g?3WFSuiRv- zEC-My!jb}|BB$Owt@4TILFPPZATm#284if$gC#+g)K>wpLE@z3`CqfFf zN|uNF=S$|7mmrxe5DJ)QJ@+_+oaberl z_^QgpZhK+xa^$(p^TA2m=Jkx(g7$K$)vKP9`i%#E_^C>-KqtsLVygjSo_}I^e-`sa zIwjo&J@sg(d9$m!MmQ_zwqm7|!3OhjfZ;j83CRm*$EQkycx^!3i+(+tLweVU6r6N! zF4OS-+3&P=7i)^PV7@;M$H85HbB_Xxl9Tz;Dovx8p9Sd+a&mVBmG*$~Gqk>DPgdoS z2!R3stVe+sU`io@G%aUM(HT04mrD{j59wa|s^~Ju zcm>aAogJ?3XxUf*piY1!p@A2SX109i$1M`$sTPZq7uC(@=?o!5%70AuUr}*;$#vQ$ z{?WS$llLBmM97$+`Tg2%uN&En!ne=3>xARAyX%FaHUX05KnKfi>7jomW{3DEJi(7tJ6>Xa?D!M#|~3C32o z-D_T3D>lv{itM40R|cmW+}q~(@E<;Vv04WM@i}*68}N7BNA09mte$SIc&vHO?{I40 zgkT*|tdc>RyR}Y3F}fHcK!$}4P5uMkqFX*B00Yi#;LL(Upc~2GWZ1zpTP1p_{njqa zcH08d)qOulbd8GUw-iK0}66+WZok zgDzC8x>$$=b-$IDF``~e{Whp-aEkP?<&rz!EolyvSJAp&gdhMEmhSe!!pU|MY)Jz` zVHOATKrnn#7zq?As2wR4{ZeZputwM-u3WdJ;}$#)Bgu{6pZee$ldG;lTOsNk@xcYv zymo1=$%IHjVU(=Bo*g7s4>HLh&pj)z<}S<@Ur@^Ulu7=)e@8fG7A}TY6qq0^W6GVN zPMZvSTyd}n43OFLx+yN|%Mz3vDL=~oAKe1k$kI%cyUs-^8)hSdK>{qVN|chlFnip> zZi>UC`~t`l2u&WLCH38OorM(swW5ht|Ap)7mc@c(cn)xV+kXPdL;x-TIe<^wbNe<% zPprOhu9l2DIQC0-!!l7e2l>4A)Q7I}{Fu2!tiLJL-PDeU-QZzG`G{zPB@hE;pG@{q zl&lBZ)P1x4+FcncpFY&$P+n&lErX-(!kID!@Lv*RkAa0sg<_q$_R^R~!U?W3bHlG} z5`c|B$T?9^jXkY!z?2cQJ!C4b0T_=SB`CSfV_T>YAr?N=o~Zqv{r@(#vTCPKY?2Q^ z<*A7z@VNa1^x$K?jS~b4&&Iv`(yaEOl#C_{@;=^XJ(fxE5teR!FC4p0l3S&KZqZ9f z8b**Kj3|^-BM_qauDJi;W(*Z3GbWgUl@3bl6>_O4n7G=Ykxh;5an&2hnsNA8>~6%M z26hgWgJTXSWv|V(3>0EaZMf-8JF)J;+8Ad=f^4bN&%Y`-BI?1Ghk>0=;8zzE6KAod zBwmU0O-tE0ot4biU)}!?RaQx6LTJ%jjn&$h(SxPjyxYvV8Ss>6j(_WE7YAjq(XH@P4!!?#^+bruydO@ z<=XpEOfX#2GAtRMhU?1zi_o=~lMdi)_%`aT-#3+U^bi#;3TxSDXO&@oJe5-q2PX7s1&N;TTx>7ajm(q}8O6+3sLq@F)-NjYU2QOG z?6xBkm~$x1wX(DBpFx+o4Lg%kNeC^V5|JQWllIE9%dG=e%H~JIvP?mgqokA+9vr{D z^64a^d(%J{E^RJ565Aj-F8=w?7N<)U)wf`^?@p`-ElG^+)vQ-0)-h3TWeg6TAMm14 zF%^G5j85OJ$CbIWjh^TJdpCOxMWp%;mLgTLOZs3ETyem5eU0%=MNm@Q%ZC=;agq5$ znEyk`Q7}+2oWzSY@)tHM`0AU&8v4tkXF+pzX)tm8-U2Y%mdZ0)l>mqVQ(St-^L?r- a&KE5^LGb Date: Tue, 29 Oct 2024 20:22:53 +0800 Subject: [PATCH 5/8] Update demo data --- data-raw/data_test_03.R | 2 +- data/data_test_mod_cat2.rda | Bin 10145 -> 3265 bytes 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/data-raw/data_test_03.R b/data-raw/data_test_03.R index b09be03..9f7dfa7 100644 --- a/data-raw/data_test_03.R +++ b/data-raw/data_test_03.R @@ -39,5 +39,5 @@ out <- lm(dv ~ iv * mod + cov1 + cat1, dat1) summary(out) summary(lm.beta::lm.beta(out)) -data_test_mod_cat2 <- dat +data_test_mod_cat2 <- dat1 usethis::use_data(data_test_mod_cat2, overwrite = TRUE) diff --git a/data/data_test_mod_cat2.rda b/data/data_test_mod_cat2.rda index e49cec3854f38dd1229b9071f10f4fa39d0f52bd..a1b5c71e4cb01e68ea11c680b6b5eb033a18df62 100644 GIT binary patch literal 3265 zcmZ9Gc|6mPAICo!rYJXae9bmO$T7>(vCRx)&bd>}(Ogle@b%4*%*~kP`m)X4;fN&q zZW_7L+((gb3E}Jejjz=A+wrT%i+ei?ft}8+xroVo3qOYzkU)=F7=!iSq1>*M`VMN=qpOGN*!dR zoq|-VS}>)`U5d_epQz#tp$LkAWCI`r{KsLJ??Vj$BvPB|zXC`xoLikE9Dr-nOr+`@ zS%8AIZVns(^8f;YWTMSc0P!Kn08n4BfTt@80Ahrc?K2pMg9mva zGa1#XJ?Wf`G*P5<4zh~lOh_Y%8~BvxuH~+sYiJNpJGGFXz(*VOhq}&1e$e;v zs~zxk(pxD@SMjx(G@@&ZKV|c3p_?eL&K;|bvszd7@BypoP-Q+{TQYDY83Gj}2#c_` zA#Ih%Nu7f+FQSlGF{<_xc@QeaI*bdhmKu~QkP;I?r6YMbZPPj$s@w^PLQQ7KS_1%y zJ1$1*yd`zt62PNae4ZdhEi6?Q3Md4HNkK_n3L3JqeB=ZONJl>|LAoEz@x^LGDOLjT zG8pU@(1a3z2Xnh90SuP~1+_3<1t5e{VZ|__hjXWbRXQ{OXahN(4gvvAjrlyRFOp$iEU~a!bl>MO43Mun*&E)4t~ z>lSij*8(diTm2>@wfAv!#OOMy`32v56HGw9U4bJqgqSwZNqiAqH*afN0-ULQMoI5- zS5{<#jaaW`xFS};9&6LIn`hCXlBoIiSf94RgPlKTn1(Zp;U|cCrNltPabBNB%<&kk zxpdJIx_*6?TUQMNj5_{MYfp~{qBz#_v(dnQQdC!|J7#cpLxD%Zdva5@8|> zXSs#df_sTuR?0>J?+xW%C2Ad{za28;yHwv4j>M85R{E9_Kq3%UsgNpYG?B9bh0}8#Dd&nK%%Wa26fS(G+KbjF*8WVysSGV~56Ir&1a=l6^a%acMA-tR z2{}>Kf_$ju8;MhcDD#bBH`gaKn)&_DhF$~(FFd$2)Q_`QiTiNDAlNBgmk&!eQcaXW zg_bKEjJO$|+Fhd^;hpN;-ZNx2^GhheOmc{u1W3tN37$FBRxhuAS{&Jj6ssdMh@=(E z*dqfM8O%@Jt_C=JaVeu4gV@j+b|?Ne5J#$Rh3r?~J)WShvovgiUNekK{FT)q3nHjy zMcCXK@K#<8oUcd~@Go}@g|{gsa+jU0yC341=uYTA)|Tb!CG3ZR#;b32RDQE1W=-TG z@1M&QYq@Lh`)&p3V(cn58D>qPtj0BU*BEx#Pi3ALU^ZrJ9wo;@^jydCCT7U4@rceO zt;BoCEWaxj(U)CMcF7gL7UpJ2mgK9H{XKSZ*HdL!upYaLuHrQ9xnPa8DgBTOOmp_j zt=qi1dNogl6c52!#7wqQto>n*X4)(kQ6gfnR3|?Vf%3y<=-wWnni5Q|vF$S~axm;# zArP`e)l>)j$K1bb7%oQ{kED?^^r~z87DkJPPc?M4%)^^CG-if1i@WiR!idu$>`C)9 zc#e2}+^J)dg03jNm(L3Q?XGeqcpPb_PWIN7g7zLq*LFQg^va~kc2#1b+2NK+;i=kPsU9&cZDP`Tnfms`jzPGMLTm*1cTe48D90P_ zfx7cFU&lNq=;GxAO0_Ha*W;CZ!yMUgnw(`o4I_x0R5kBZC_hC%F-TE;K%6p#9X3#U z#?i9jJ}8$Ho$v%2Lx&bHK zKi+4&QrQ!Blvs4kBCxh=DgtTRkXY$dS0;Dv=+>%ib^%Oc7sOall!Fcwr;dQlhP5B# zlUaYj(s!5rmU@%=KYj>hKK83<_DH9B*rN)NG$g_&5D*;+VK)uS0UnN`kP(zdG$B5& z?kw=KP3XTmimtyip~1-l8|mWoidj)ep&eVE0GSxgpN%rdWnt8dhIsAlnU_2$_$KOK zm6e-Vea^VY3uC8JYQcz4ASH5^mPUO^+5boN4?c*~5M0NQQqCXSx%AE)$5oILQ|~+# zKBpg?S(mT|7KK5;YaeB5sBXEN3^?7o)R_F_G~BnRdZ#DQkz9ZH*y|mCv6zuWv3qvj zD*q;v+-`DUF!6tXmay$gSChT4@D|!H>z^7>Q0Ji$XDxBjy`8$+}Wz%xf1NF+@u?sW@(- zD070mG49(pu-qw9^+X^j-@AQoEG9EWT!*^z=;CvxmFJa1*&ps1QfS#-)2|yfz1iJw zufxx;uRSFNc_sNjYPJEulUmX#w1d<3>;-J9w#(XJb4<{+Mj1bfq;i~xgIXpO zFL~gb&68%#EZ$&vwOhmX<(vg0TI)@jXnb&T z>`U`J%Nuzw^jb{g+22U}>zMrg&&)5mY_O|Lj(z(=-WIo*j&NFYyAjT`nkHO2Vpaa< zYjYk(1dLGyMj?eP@eUrh*sUxTPSs+nA!RNt8)OO&9%A9qfXsuAR=4{7uDMXZ~XGaF`(Pc3U3Dfk@vD2|C|# zm2`Pq*G9#nJJDlZ-mb6Nx18f_VKDXXtFuFveOHW0@%)Zes8{4O}8`usio*rwxC4KxodazSEKdeaFyTX&6&j%VSog0_&u zrApy_P32zjbH|OYb6rJyj#kDuEhYPf$CpQ8^>K#$kPJ5wn1hY!zC(hcOa5@E$}T~X zmzmkqdeJi)8{7H@PNU`0qJu!aW4)&}Kbr1tKW7*6aiHt9t(G|#p)vU5Agor987$OF zInsCfSxx78UbgDHH>8rv`wKF;t*xy&LElkbP*g z9h@P=!|3??BQn+MK4M2F2Q5fs2NM9OE>e4D<0WHfx{@XrnI?mKp%;sYt6>e*Y1r}ZAX}x&-IT-zGXjzRp{A;t=}F4{s(aSz#0Gm literal 10145 zcmV;SCtlb>T4*^jL0KkKS^uR$)&L8h|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0 z|NsC0|Nr0_e*g~8JjZXn*12=q*}GddS(T!@x4gU0ZQ<_S=XQ5@dULgSdLG{0yf=H@ zyW4xb-rL^STJ~RixbM4God5s;0z7q3(NXh6=Kul!Oqg4*7 ze}OVI0i?k&)Wtm~rYDe8{Z9mCGisTd6xye>#XnT`ngCP6X@Hdcr>Cf48W}c}*w9Qg zGBTc;8YiP8AZes#rWHR!H4KaeJx|nWsj(-fiXBanqe<#%nrTnaCa0&V>S?tefi{sm zr?oXarh1P-PYP|6%6ecLG@g?gJt)IS$>M3IgF`8agvxC*LSTlPn;J|hreIAclhH;p zWb`IP!$u(4P{lPoLU|_AV9OrZjp}%_)pcC#mFU)Q436OqoX1 zFoq#d2AEAXPbSh0CzDgkc+oX81wT~so|9&Vr14D2iKe0Ic%#$OU{k_RPg4_7u{AcD zhkzzcM~zQZ{RWyu@-mE`qr#r2r>Q+pCQNBDC!^C;8lI4Z6GX%rHBU#V$+Vhe(+TMq z(s~*vs(NX()YS0OW`Z>Io{DT~HlehJQ`)B54H`_CO*GR@GEM3+8x+*_Jk-gWPbrhp z85=~xdMT&0o|6fwwv$HH2&jVs0EAQ&pr9QhDZ~PE@%MYUOb{IVY;(}+5G>#z7jmGb z1p#;Prq&b&0H`RNF;D?Y0v|Nuoi7qdLWGn+fT*aV0#u4Z0U)9RNGO2-BpbO1vJj|{ z2|`Fj1QZBK1Qi4*5E2mrAVm^U0wM%b5~8Fel1T`FASx&k5dswiKv6`J3XuX7B@|GA zNhlB@Ac8?8l0qRw0DuVyh@udpiiuJZ5FiACNdywbBk`TIwA&cl6kSKjarv{;srwto;Y z@I{OCWsbvs{s`^WqtP`l?kv0&#~|^!GYh|rZW754^_ANv8)(5ci_nwhNCVQ@wTM-GmC*; z$MYKUGw<%>cd40CQ7{C_n_1&9_k`;956jhP@ibK1+6%iLZQ5BmgvN%#ihU3%Bok<(Mx|ZS)UP(AUi6Uu7v5H`2|8{75UlJNgE@<-(a|tk47Fc7qp%XIA*fE@wwetKa(YWEc0tL@QN9fijFgqm{T@3oT)+$<@LQS!A)MiO){5= z$#0o^XJv+#Vxz4_;+nw~_XArzN_z04kkHwcFghHpx)}HYtq>loi|UUQPLC|E6>`e$ z2iJ<|m1;Y<`-CxA|2#El1{cJuhEH}lzSnl)3j@mc**3VwlXl6u(us^81&9`mvrXr# z&{rApnjKHjp~K@{lwo4oG5jR}qN(#P^W~}Xn{vmd43G-*TIZw(b&N zhB4KHyu!NOa-^4q#GBO3W^49_&PdX_+5s~a`L-G?*ltd?X3p|(?#ze*<6KoON6?|L zo5yNxeInC&ZlV?vNJa{bRC!3xuho~jO^@D)sn&;>2{TXI8wHZF%dYvyTFmx0^xoa4AU*H!zkDna!7mL;n!I&dm22C%)Ao$`+jg~oN$3bI zw%8B==RoOW{}q;9Bmu`$XptFaTiux2@8aB$AOHY@06a6+2evXNKC$|0GRiV}+VdAq z*XC@k9tZ$}003aUO}|$-TY43&WsJu6co~{yuHlU^RBh}pj1Gz$(=)*7E?|kLLl)Cxv@5IkbWJi$V2`9@IXqLNlC)Ot9o%o;0GQSeHCwGRfo}4{$(qIeG zmqu1vYa3;UwM}?l9^k!T7W3%R!!Y*w_-cyA$uI%r;$p!}#u7bd*MijEJL2q{HulSo z|DbdJt2!sEIP(poC_IWjmia}^Kgs`+V%j}l#zi637K;!_=4z-YsS3!ZXCuk}*-odT zL5<&zi2dgcnob=7XAci$yp}L~rX@=c>b}{Oz=TzNH|*{{2M(X}bV60)a8!1OLTB%? znk8w3-luJgh_^eg&K(Vj8mKv2QEtWi`95_g4Ko+k>wP=8b_MDU#rm8twV(x40LQ)%&wPx-J0 zuHFWO%lBk4?&>Hx)T?Ve5x1HYQ{@~_JKtOU^3BN9^#h9Mhr9^^_gcd5`(ejJ+m$%3 z13ZQ5xHR&;%K8s89Sj=}YA+5Q2Mku7P`NDdRR{(<8|M9PEnr83>~{8be`dRIBhzJE zqNL(QJWR$jASY{>I!d9)D^d5;atUmfzpEvLKVB_3$GrT6_2HbKCR$1s1Cl4di7aa_ zk$cd=_FIXuOy}GGJmEKPvK9>T$2MwtulLG5X(Oj%u<+aY`{3ye<360Mbqk(NuXJqA zc=`YQc2|w26wE-`lo4kV(OYb@2^@T9>#;GkrT$-X#BprdTXAcJA=Qs0ctiBue+qtQ9k>FO-ck^0(`I@i*^?LYK5A>|s_!(tm2~W@hKILN zJ~F$&srV-v>=MFK(aiOC&u1RATVCkID0AIpEIE!G)UYeO%?*D}VrZ*iFS1;QoZ~S3 zBHXnx_hTI2hNHTb?&>j3Wxa}v3d=;h0iCQMC+h@DB|H&SRPfzbJ|zZ4E_uRPUo3Fd zQ!$P2_YzGob|@O+vYXd66W8Qe1ZL8TDZbIg^bf9lXh5;ipj6FMD~S!k?+18CT+Y!d z!g8rI(cJq1V}V7K1U(0UyCr{>!aZiN-iXr>srgQ?J>2P%<{Yvb0m70 ziWr#H5^$Fm46TK(f7B-5J1We$W!9?3sF7Lmq? zDb#!zo}*(++q9z(H`O_jfXTzr6;b&zU*ZSfE!?YH3sZm$+3YD*S!NsHadEQQV>e|4KC6VGc`V18Pv8iIXn<*;lC1CX~YAN60tMIMp|gPb6k=<2p-2y6E) z$N8HjbW(diLWW23l4|R%gR8t&{wJS=n1A*0up5o6{^UqOB2Kp=VU{R?)EDHu+N{ZT z3Q#%A%%JF06#in{OU8$s$Iz2@6lV2y^i~V5^?TK9(k6?p6W2V$@#N@A?}D3pymr-c;!S@3%L=AmPX zBCKVGA*p%OC1N$NhqzCnhB!APbSJm?TTwx%AFC_h# zY1phv*pZxWA%@c_D2epD>nGd%rufQZ;YhE`4`0V>jLvOhI`FdZPMxp3Gc~(^>neT! z$AllkOB5-i$$EU0`*^OCMQ)m=dPzX4+ijNfzjtv+D^)m}dlGDOF!-NIx;Uz=45@x4 zZoVbX6J^BXqs|LTT^o`@3QwAfesCx&k-1oxm_pUv0hiv*vwxhfwaB63RO59R(|^W)L!T^qh&H%!c;eKp~X? z-~&_EeI~H})%%N>&t?D7rq$YBgOX0^GkyhbK?*jL&>dma)2R>Xb;NxQ4tp(=^KnKO zo#dtI-LlknXqhwF&D$Aw;gDg#fG+4=jxtO znxjD&d!s$Zw{VJ6^Jj2%>bLk#i;K(v06_o{0NIzA_#Ie{iwxR%PjbYhv}4QdG^plc ztbBR{m}jl)kP!QgW(dcWg4VIC^EyXBGS^_ubj>bv1bwb+z&aT;?AvUevH#J>mYSL= zy?3ScZf;4lKzQUibs6@uyg~wQ{P-<@#7|{_RZE6nqrC}JAb;5RXBRDZY~8Ij9qSJS z=Bu#e;e}hSz82~3YJ&BR9x_cjpbtD}>b$##Z)5f7xVhvp&E#>Opo~W2VrH7D zB5q3&mwmJI=EH*wAB3D3sGG99_Knis$WBYL8zu18-fn{AYAPKho+Te8<)%f z5HmZN4xMzCcE{+NU@u9eXhuIJ{A7-r!&B5Wi^w0nvmrp+(+6Lrk(UK~F@x6)nNka} zbJ|v9c*o_Cp*KYH>6^C&omjatk?ppRzXu-j?B2L(JDWyT=kYv0p_i#?pxIszm(`sy z@S0zopaiT?i!15yGb-$ z(j0Hy6EnvtJaucm;d3WF5$;{|PcNr?dM9pIpiMR7?bCP*)*xo($q*3xgO$mps^O>qhg)B@Ut?cvT4PH?P(wrW&MXF@lv6u6$ zjTbS^UcwOVBnHU#LVzOy_tD>3_|glgF@Fx!JkWmDmw$TQY#yflNqNVKEDH~$e>6z} z4_OV(n8kpm>@1`x?y5X`;+z*TER&&%$@a5c+CJ^Km$gCiHK|Kx>8gWi)j*c!^TgNn zgq@Ufp`g7i$S;p`+MY8m3L%Qp?+gpAF{M%g~!c`wR0!xW0-=4#Pk|m2X zyFw;T{~>e>Xj0trG@oVHC&bI9Jelt=PDvD%Ly|3UG6b;wZ(i>9<2Fs(Qih$Jfa|ypL3LKA5@u?@Mmk{evNb~(Bq{$U4!xJH&v_!!m%$MZuCNnQ2}=D zR7>vnh>372t(Dks)-&qZo7PDk#G-i?qaJJxXbPf2H={=(GFil>hDw_6n{wg)Ce+KG zDEeq&{VgKuGXT;#Q5*%bC8JuSm!k9nE6P(d88r94!-{lCP%?6%+^>9U2>Gj~mOYb8=YM^{lU97y1l!3jGz1`rJ zGZYuW-(K?mqIH3?qlmmJjoBe537_S^n$a_9WvSsIOY}Aes8r{;Se0)owP#`#j1W+H z{9eIxsNILOoUD%fL*pWH%x8@$3E?RDI=o~N?2Xr$w5_8ykxbGNtBsMw#6Q&4yvtUZ z7nqe@pg(&8NbF=Pe-UnEO98TJ;90xc?bO*mU0C|jkG=WMi3>0{;y(MbGu3Z!4)Fd9 z6Mv~|ba~eVo*gG*UOQJsDlPgAEVW50xGA>%ybIIzqynKm=jY!%9A+h z@wiQ#(Mqo>#CS~wIevv>x!M;zPw9pdaY-Tb4L0|qiJL5+RKWqF!$VO;hYkTPjpRB|hHu5Ty&E2RjIk5VM>uodtcx zQ_ql@#vi)2wj8|33C@BGC_lA}llmzxP~)B<>OS#~q3A5Ba~%g**Xt@ANP+Bi>@-kA zp@OuqUb*Fin%PRUc+fw~<{5H%;gt%}n+(_;h#-{$d%jF-v*)N0(2Ks_^e!mL&C&b~ zIx_T#an{Vp*fTII`HyNa+0`xdU22wkr>s)?Bb*jy($s)-%>1UbsabB9f76AH)S}IO zp7klaHl)^+RjJ zRZS#UH3582x2J(0lr@#+kYy10U-zYf7~XO+;=4++;(ioX|2@jwjkbfdW!8gL8FvrI zqS(f))7VoxYzF;dq;8aZ32heu4%3dlQ2tmcGCQkG zW12R!!u~8?b!Otc^kgBy1|849uiRo;g|o?3)YX$TDwmep#L6iy8}M~PF0baCQfHhr zybOxpR);|SyB=Po4x_xSWSLKdHJl5so<;5fUj_buGy1`ViYiAGHveC$oT5S|Ca?|; zjjA*8cHAH44r6P$l`jXpj@aYJj=*0Ofqke1TNLh&Rp(ga@)$B&{oJF~(frPFL+)~P zZxQOXQ%YyuGCa3)(cd|PUo%W1GLUTfU~hLnP11u1`<;_Odmk1B9s%-J2#Zujct^GZ z7Yn$KH_Gnj=xImo9pUd6XJ?Z#AXqWAMiTEP#39h;m=0;C+BZFt!{eXd&pkMXM&qp>R)le7?SLsT zMSv`H+d7l?kU3yDLX`bKKmfRRwKK-AX9aO-i6-uB zhQ-(r0c{2t}FeiLW{8~P?*ap!LKl#BBE;OX7=AZ z#?B)5OY&M!PFwuf)4i&(7mSOQz}>>*qvsnKAyH`aUW!xQT_^K=^7g933j7tB%b1&A zuAMJL5m(tw6PdgjSAf!hJkTW$GR%k|D4|e=fXl;7VC1JJ7q8w`J^B`ALF6F57Sqmi zVz!L9N}+AW{v9hULk~dSA1*Fd^Gon3VsJKpY&5p)Z9qsvzi7$WyU77&#I}1`WK7F? ziPlYXzq|7DtA8$R8mcdp(3rw@ue6N09T}Ro`Nu7*Gatrf|BESR$?%D14$DNWw0T0x z+{UeIqhIj>M1x+vFXQ`a51|-z@(zKWGwgOfFUAuA7ydqTn>jY$x3|ow_hLZ!K9(3m zA}*;(YsGqbgocHh1DiR$?d*;nKi;W0cXqi{Ib*|xwR?4JjXB<5EefhedxJBF;PGq4 z|IGNmf~AF0`z!7mj`F$z8pXmMVnl}=gTWN;z8{eQ-n0rDLA$Ig3K z@}Rf5VRso-#lirMYm=;{Cj{1lLVknN$0O05L8%1^5#2@b-Zzth7=pK2K85&-L0`iD z?heITeUKDbfDtbh5z22`!G*b|dMl*v{Fl zfP5|{TE>&{KQ$Zf=!VF#YfegeVvm+ZXJSA_>K4d*Ccg zlkK~xbboxx_HJCV_Ku2aY3ut-Uj6!GX2Yw4*|Zl2cL=3l-lrA3h$@`U_#l19{bFda z@9new8%G&IZVPZ!G~d((0je}ZzzlD{z+=AWQakf&yo7e( znY>|{4-#u}lw~#vy5Ps;iM9eObHI$*?8n{k?70fsO8%vF7sB*_N~O3dLLrP0i42F_ zW}ibKU4Mm6%83ia;z|W|)grrrLY|HVcec*1fC%%cXF%0?l+-T}*RJTFxRuF!QCrj? z`sNfWMv0FFcl%ZkrH8C-NJ-x0C@}y=@B*u3yzzue2T(xND{U{rq;s`_sT_Oa^XY`$ zGDH_;ACMz4Gj}|r$o_&_uLMd;l^g$j9EIiD1Vr=Hu^c-Z71O`=Zi%7j&e|Tu&DKuH zQQLbSuV3Tzc+of*cX^DY`(@neH#)oe!Wo{(Ox&P!f<-SwEL1=9l4G3O%33^XSS52d zUoNR*QceJ$kKK>A?TZGPK=w;22se22bV6JmZpRSJ+F%am-PBKYIp3rY`dl?q1bPFZ z@B8deHRlQser1!{9uzPVrP7qf2@sKf^_1_;#tt@QYk$JfO-!qhE$XGm)1#{RIrgLU zwOsMV9jUrz^{V(rWyhC7D?=Rzb6MD$mwIU;rQ(R0EPlJurdB1HiT4jAmFQFw?W|rI z(mta9Y{%J`ucs5SR3T{@J>ax+m2D7Y<%vjg))NQ18#{y&2*8pN_PIDPlm)tSNbs&; zSaOy_6R*2IrR{Npd3xdUYgRu#a%SqE>3Xy4xyZexA)FpZ%lWqxr}?+B-o#$xQ|`Zq z279abA^_;2aWc|XsP}uQ_ZGj!=T5_gjyZ{NGP~R%MYy}|FT<8F3_;85gG6ca3$^^X zwYt_q0Ch6ja;(-r2av0Q08hh$QVnYjI1_y1Ex17ur19u4g3GpAT$2m5UhExl1b>AD zP;`r{!sf9{6Mbl!*SO0(ET6*8Z{?%Vn>x)qO^FSif0y#v=NP z?{Qg7>y)L#r9?QS%S<3-?z{M2_DY92YL0MJ)_1VNat;xBPBJSTM)rL;z2ryo1@3r_vAK zTP<;dwlWu>xaxjfEnROCFKrAa)_~iPAk4#VxtLx1fzg{a3(JtX--&hFTrXy$gyt?o zv6^7KoHgHNWMujGK3}0@MW+0oku3uG+T!eZ6cfG~rs3JTLfBaxdRP?ByOHyMhsia)dzvS}o4PgaX2n%|jeaXLb zb!Qpnn2sr+;#*tmH0eOsqW2V>qiYT zHk)iKyLRzO2=%;v{tF?$pULBY)K&%8#DkE=Ep*5=+NNf%4Op4W za+#UsL}VF1l7GIrv`w6YQo%Hx`zR{YfaM7q0@U@~_M#>okY#Q1b+345xL)HnkK15G zoF=Jbbe?j?fA^m0W_Vt+={u4K<3(1wyLH>)<2B(*o;79$zTTB;3IyC^T70&!`ZhlxMs0;sL3gI55?CqyDB#@dv|c`8|#}Zu`OJgzMkh7 z(Oj{RG~QA_?39@Awi4I&iwVt2-w2e7^Ca$C|5fYT$VBcP<}iYKkM|v4P0wqV_Uy8= zK#MxUfMaggJ?9FI*zXIcF93TxXH}DA5*8{eZCMsK+;+Kioqn%L>BkfI(#ko_#d_ z(@O0Y56$;B_HNIk07`xHutwy!i>sGdVJgwhqS#?a>Ke?AN2f_>_T;ATZ~L!#C~9#Q zmU1-FKgZOLIrvgBbvu8Oo58HvydJnxbU9oW=g;ZJx5qMmdiM=zi^jW6>U(R85=pUM z<$K=KD2Ugb#Ala1INfoNxnQmFEt$vum{v<_r9K8*<0$*xdv`P_91)$$^0CbwHVSI{ zD)~n9{Jl-n20ZRyfS(F@hDt||fM9osk?TNJut@uH}%aiuld>dE+Wr zO<@Q$mR!GX_p0^Yt}lt2INLr4FH?}fI{G{Nt#KVRS^T8&-urKaIU029uv}Fe2MpVA zPJc6WBfc(58LbB8_k!-x)p1vbT0$$J_}a;2Lu{E#sTzaq9Ptj(8;SIbIkn?-v9tNo zEzt#NK_d`Dgi!N*9z}%H@9D06^spv)>6g9aqbf+&{g9f&9(@A4+?p!`(?*TiLT(w< z8j_|QMyFWFcHgedUSTq3pv&ED{Aa72PkxbW1f4EQMU)0az};;wn})y{G~C0n{cTeG z$9jvFGaHq*p7@mqi8Y(+FhC&6?DzNVBj;P4&+KG6IwwQ(i}CEM5|UqRQwXylKyk$f zLk_>Cib)L5io(g64kGZDw^mQGf(7HF&vA!9=Lfs!X|7E%sb zal1c zXo0+gs9I+dQE8#o$lkm)HEMl0DQ@!&}9A)5^JEj*KL8-_SmY|N*00^7KIqL P|NLFa6yZWc{*?z Date: Tue, 29 Oct 2024 20:26:04 +0800 Subject: [PATCH 6/8] Update demo data --- data-raw/data_test_03.R | 2 +- data/data_test_mod_cat2.rda | Bin 3265 -> 3250 bytes 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/data-raw/data_test_03.R b/data-raw/data_test_03.R index 9f7dfa7..0ccec56 100644 --- a/data-raw/data_test_03.R +++ b/data-raw/data_test_03.R @@ -33,7 +33,7 @@ dat0 <- round(dat0, 2) apply(dat0, 2, sd) colMeans(dat0) apply(dat0, 2, range) -dat1 <- data.frame(dat0, cat = dat$cat1) +dat1 <- data.frame(dat0, cat1 = dat$cat1) head(dat1) out <- lm(dv ~ iv * mod + cov1 + cat1, dat1) diff --git a/data/data_test_mod_cat2.rda b/data/data_test_mod_cat2.rda index a1b5c71e4cb01e68ea11c680b6b5eb033a18df62..c99f2305e42e2b9e3cbff993e8264ac5f7c279d7 100644 GIT binary patch literal 3250 zcmZ8bdpy$%AO4xq$h6IULz-C|=D5sAskzT(?u0n*Y@worZtv#S+-l~M=5kyoM9~m( zj^uue$t9(nWFaLXZ>7$w_w)YsKA+F`dH#96&-0l3(2nStd1KwYQ^+-wz=4#1_VDu; z?a9xbUvha~xbv^Uk2{xMf82?fe?PMzrTSIu>MhSHkp%#dzvW!((fNY+tT)Odp2k)e zA%^q}9^=VC=ryL)n(FKnQ8E+|#WDbqQeY1vt9v9$rU6hf52sy#CZ1g7P+$#MTYv&& z9B%6Z*i@>bH2`M-nVH-G3u~-IqCTxlB6%>4;Q$5z@vL7*1oUZpaWr6D9{@s#L=hf? zI1coThX6DrAX1j2M(jo8eks6gVcojWsuj<+p6x9l@AV@0s;PSlAXRy(jB#-?QlECP z=71ah4)T(R`Tsn`Uxxn&$dUsdY$i@pMbRY;Cl_9!n{_5vmHpq;J@4LvUmvtz1nMoR zAXP+-#YiwL#MBcv6JtEmM=I_=Uvq70wa>g;so zXaxl)f)ue_G6Qy<*d|HiDN+dlo2{qP4Zds+k}c5cnwWF9tz5xP5#&)IrGWf&(sb z0lO{nv~FmuMyXhj*uu!M$wq?+yEb(bJ69?-WqI|&V>r3sZQrGm6|MG*-)kDhVd6~MA_Xm>K$L} z(KvZfU6p@N*z(?~&*|7h!FG7RC?xtrw6mFkxMeFRygsH3?+h@+8ipq{Z#8ek`)rmI zSZTfiz8U#vH9hK>3iM#l{do&kN7gZn0VG8zpp5!^b;PI?lJ1=_v5dcVZ{UsZ-&hAZ zNyyb>g?WgW{KvMVOu*`q1uT^_6re7Od}f*+NN8gPba_Q0j$9_yD1;(?PhDKaTqrz- zVY+p=_HK*p_pYPJ*HX-%>T&}p&RLD6b}+U)0uEAS-*i5~Oj7*@&0G$}K}S z#5pLR^4ywjUA9FFGeBu>h;&#XINfA-d(%g`7`pwjwzFnf@C;=vck$0BY#J{I`&b>RgFZ`9Nf#8#4xm^48GUh-1MOWk2eQvQ{5*{xx=>=g7x_sxi;#Y!Vs7;XCi-{hcO>Pc~oVNOwvy zsp3bKt}bh3y|NyRYs!>P7*t*H)s4~ zgLBo5DAsZI7`AVT*5GB|)@Ez9J|QUQS5OHG<-nR2$Fn^F(828XYG$CDp*OZ{Y}?)n zVAsO(ZVzNTv`vNDKev-TU%_wAll8MDl=Q!<8%5{OXPtVP>!bq?HRh`aF_(gSg`Eb4 zkET*SH<~Day%BzAwbLDD&+=-t_Fl?dZau1{i?8X+OJ0zRR**0(@AVZ}EcXX&kK1aH z?OJP+TdXO%KC$JVa%Q+tXLbARO=qe;@^=hbDUQ3(x?91Qpp}74_TryX^Ei}_!+Cyr zfN8pUT1T5iS?%9v*mONqCmZwD=jm-4e`H}`Q6zZ&fc@hQ9h9AsLC8s0l@*cuxF;HE zLElv>wbx~2BZoxB5{M&5%PYjO&VOH8u^D7>D@Jna>3l=IL_qekNQPdP$LR<*6SMaV18 z(I}R@KIZpWdhXDB*N-iOy$0#z_Pd8USteHCc4$G*j#=0yW#>wcw!`Nr9QTImTq)~} zTAr`AHsgsTTGG|EJ_a&j`Est?*dr}tH0O}}tZhqxXO4Om6IGO$LQjAj^TW^CtAmlH zmX3H z{q3dIVUk_!%sX4!&)gXG5{_E&^PlT8U zW{xZXk7qM_MC|f=r5TWZ!WFU%w7Ma zjePhM`$1(vEjKnZR}hGv_*UZSu!PT33qyz-)ot`AmQd19yBDcJSRv}6-BB6lxw#Dv zS?S7c(Ypq34626@{MIigx_#(*{L6A=9bsn7C>8hD%d_pFCo9otvF}%~HuVsq+S8Wq zB#DEygh|1b4`#RT4LMq z6q<6Z_(|aV1+=OXdda8(+5gHlWJ`~<_Hu2$0htL?wx3t+^x0KZ;e_4USzt8FTh5&Tnv|%=t+_Zx^Fxv$8^n zKv}}dQRJNPW|il*cm$<*tV(i9JFYhKjbohA>z;510|ZMl0F;vS4S<$^%jQ+`O4HD` zk?B!ipKCe__*$#7qqly)nCbzAxE#B1&+|5G z{=D8rCCwOu!^)U8G0CUL@R}MqL-Z#7EQGn?f3+?s(>una_-ODG?bK^fA&~P#>ci;*F5rw zXc(P)FJT7{+pn$9mO%L6Z}w4_^z;p5a(@#>;>j#lTUDgT^re<5cGi!dRzI4xL0@X; zH2?i6F!W_bQTBQn=!9I$L3p%lmw4VlIot@)SkZ{cJPjxRYXef=|WF$xp@3UQP-$u>+J;3Wdz0l z1y(UCkrGy!ALb_>Lbp{D0TdJLT&$$>q#Ll-nV%P{)88e|HO?mZ);w!kfx`mXQ&%*< z4skkq0Z-tcg#Z8m literal 3265 zcmZ9Gc|6mPAICo!rYJXae9bmO$T7>(vCRx)&bd>}(Ogle@b%4*%*~kP`m)X4;fN&q zZW_7L+((gb3E}Jejjz=A+wrT%i+ei?ft}8+xroVo3qOYzkU)=F7=!iSq1>*M`VMN=qpOGN*!dR zoq|-VS}>)`U5d_epQz#tp$LkAWCI`r{KsLJ??Vj$BvPB|zXC`xoLikE9Dr-nOr+`@ zS%8AIZVns(^8f;YWTMSc0P!Kn08n4BfTt@80Ahrc?K2pMg9mva zGa1#XJ?Wf`G*P5<4zh~lOh_Y%8~BvxuH~+sYiJNpJGGFXz(*VOhq}&1e$e;v zs~zxk(pxD@SMjx(G@@&ZKV|c3p_?eL&K;|bvszd7@BypoP-Q+{TQYDY83Gj}2#c_` zA#Ih%Nu7f+FQSlGF{<_xc@QeaI*bdhmKu~QkP;I?r6YMbZPPj$s@w^PLQQ7KS_1%y zJ1$1*yd`zt62PNae4ZdhEi6?Q3Md4HNkK_n3L3JqeB=ZONJl>|LAoEz@x^LGDOLjT zG8pU@(1a3z2Xnh90SuP~1+_3<1t5e{VZ|__hjXWbRXQ{OXahN(4gvvAjrlyRFOp$iEU~a!bl>MO43Mun*&E)4t~ z>lSij*8(diTm2>@wfAv!#OOMy`32v56HGw9U4bJqgqSwZNqiAqH*afN0-ULQMoI5- zS5{<#jaaW`xFS};9&6LIn`hCXlBoIiSf94RgPlKTn1(Zp;U|cCrNltPabBNB%<&kk zxpdJIx_*6?TUQMNj5_{MYfp~{qBz#_v(dnQQdC!|J7#cpLxD%Zdva5@8|> zXSs#df_sTuR?0>J?+xW%C2Ad{za28;yHwv4j>M85R{E9_Kq3%UsgNpYG?B9bh0}8#Dd&nK%%Wa26fS(G+KbjF*8WVysSGV~56Ir&1a=l6^a%acMA-tR z2{}>Kf_$ju8;MhcDD#bBH`gaKn)&_DhF$~(FFd$2)Q_`QiTiNDAlNBgmk&!eQcaXW zg_bKEjJO$|+Fhd^;hpN;-ZNx2^GhheOmc{u1W3tN37$FBRxhuAS{&Jj6ssdMh@=(E z*dqfM8O%@Jt_C=JaVeu4gV@j+b|?Ne5J#$Rh3r?~J)WShvovgiUNekK{FT)q3nHjy zMcCXK@K#<8oUcd~@Go}@g|{gsa+jU0yC341=uYTA)|Tb!CG3ZR#;b32RDQE1W=-TG z@1M&QYq@Lh`)&p3V(cn58D>qPtj0BU*BEx#Pi3ALU^ZrJ9wo;@^jydCCT7U4@rceO zt;BoCEWaxj(U)CMcF7gL7UpJ2mgK9H{XKSZ*HdL!upYaLuHrQ9xnPa8DgBTOOmp_j zt=qi1dNogl6c52!#7wqQto>n*X4)(kQ6gfnR3|?Vf%3y<=-wWnni5Q|vF$S~axm;# zArP`e)l>)j$K1bb7%oQ{kED?^^r~z87DkJPPc?M4%)^^CG-if1i@WiR!idu$>`C)9 zc#e2}+^J)dg03jNm(L3Q?XGeqcpPb_PWIN7g7zLq*LFQg^va~kc2#1b+2NK+;i=kPsU9&cZDP`Tnfms`jzPGMLTm*1cTe48D90P_ zfx7cFU&lNq=;GxAO0_Ha*W;CZ!yMUgnw(`o4I_x0R5kBZC_hC%F-TE;K%6p#9X3#U z#?i9jJ}8$Ho$v%2Lx&bHK zKi+4&QrQ!Blvs4kBCxh=DgtTRkXY$dS0;Dv=+>%ib^%Oc7sOall!Fcwr;dQlhP5B# zlUaYj(s!5rmU@%=KYj>hKK83<_DH9B*rN)NG$g_&5D*;+VK)uS0UnN`kP(zdG$B5& z?kw=KP3XTmimtyip~1-l8|mWoidj)ep&eVE0GSxgpN%rdWnt8dhIsAlnU_2$_$KOK zm6e-Vea^VY3uC8JYQcz4ASH5^mPUO^+5boN4?c*~5M0NQQqCXSx%AE)$5oILQ|~+# zKBpg?S(mT|7KK5;YaeB5sBXEN3^?7o)R_F_G~BnRdZ#DQkz9ZH*y|mCv6zuWv3qvj zD*q;v+-`DUF!6tXmay$gSChT4@D|!H>z^7>Q0Ji$XDxBjy`8$+}Wz%xf1NF+@u?sW@(- zD070mG49(pu-qw9^+X^j-@AQoEG9EWT!*^z=;CvxmFJa1*&ps1QfS#-)2|yfz1iJw zufxx;uRSFNc_sNjYPJEulUmX#w1d<3>;-J9w#(XJb4<{+Mj1bfq;i~xgIXpO zFL~gb&68%#EZ$&vwOhmX<(vg0TI)@jXnb&T z>`U`J%Nuzw^jb{g+22U}>zMrg&&)5mY_O|Lj(z(=-WIo*j&NFYyAjT`nkHO2Vpaa< zYjYk(1dLGyMj?eP@eUrh*sUxTPSs+nA!RNt8)OO&9%A9qfXsuAR=4{7uDMXZ~XGaF`(Pc3U3Dfk@vD2|C|# zm2`Pq*G9#nJJDlZ-mb6Nx18f_VKDXXtFuFveOHW0@%)Zes8{4O}8`usio*rwxC4KxodazSEKdeaFyTX&6&j%VSog0_&u zrApy_P32zjbH|OYb6rJyj#kDuEhYPf$CpQ8^>K#$kPJ5wn1hY!zC(hcOa5@E$}T~X zmzmkqdeJi)8{7H@PNU`0qJu!aW4)&}Kbr1tKW7*6aiHt9t(G|#p)vU5Agor987$OF zInsCfSxx78UbgDHH>8rv`wKF;t*xy&LElkbP*g z9h@P=!|3??BQn+MK4M2F2Q5fs2NM9OE>e4D<0WHfx{@XrnI?mKp%;sYt6>e*Y1r}ZAX}x&-IT-zGXjzRp{A;t=}F4{s(aSz#0Gm From c28001f8663831f8249243ec68291f820490cfa4 Mon Sep 17 00:00:00 2001 From: Shu Fai Cheung Date: Tue, 29 Oct 2024 21:21:24 +0800 Subject: [PATCH 7/8] Add references --- vignettes/apa.csl | 1916 ++++++++++++++++++++++++++++++++++++++ vignettes/references.bib | 116 +++ 2 files changed, 2032 insertions(+) create mode 100644 vignettes/apa.csl create mode 100644 vignettes/references.bib diff --git a/vignettes/apa.csl b/vignettes/apa.csl new file mode 100644 index 0000000..aa4f384 --- /dev/null +++ b/vignettes/apa.csl @@ -0,0 +1,1916 @@ + + diff --git a/vignettes/references.bib b/vignettes/references.bib new file mode 100644 index 0000000..2cf988f --- /dev/null +++ b/vignettes/references.bib @@ -0,0 +1,116 @@ + +@article{falk_are_2018, + title = {Are robust standard errors the best approach for interval estimation with nonnormal data in structural equation modeling?}, + volume = {25}, + issn = {1070-5511}, + url = {https://doi.org/10.1080/10705511.2017.1367254 https://www.tandfonline.com/doi/full/10.1080/10705511.2017.1367254}, + doi = {10.1080/10705511.2017.1367254}, + abstract = {When the multivariate normality assumption is violated in structural equation modeling, a leading remedy involves estimation via normal theory maximum likelihood with robust corrections to standard errors. We propose that this approach might not be best for forming confidence intervals for quantities with sampling distributions that are slow to approach normality, or for functions of model parameters. We implement and study a robust analog to likelihood-based confidence intervals based on inverting the robust chi-square difference test of Satorra (2000). We compare robust standard errors and the robust likelihood-based approach versus resampling methods in confirmatory factor analysis (Studies 1 \& 2) and mediation analysis models (Study 3) for both single parameters and functions of model parameters, and under a variety of nonnormal data generation conditions. The percentile bootstrap emerged as the method with the best calibrated coverage rates and should be preferred if resampling is possible, followed ...}, + number = {2}, + journal = {Structural Equation Modeling: A Multidisciplinary Journal}, + author = {Falk, Carl F.}, + month = mar, + year = {2018}, + keywords = {SEM (Structural Equation Modeling), Criticisms, LBCI (Likelihood-Based Confidence Interval), Nonnormal Variables, Robust, Key References}, + pages = {244--266}, +} + +@article{friedrich_defense_1982, + title = {In defense of multiplicative terms in multiple regression equations}, + volume = {26}, + issn = {00925853}, + url = {http://www.jstor.org/stable/2110973?origin=crossref}, + doi = {10.2307/2110973}, + number = {4}, + journal = {American Journal of Political Science}, + author = {Friedrich, Robert J.}, + month = nov, + year = {1982}, + keywords = {Regression, Moderation (Interaction), Standardized Solution}, + pages = {797--833}, +} + +@article{yuan_biases_2011, + title = {Biases and standard errors of standardized regression coefficients}, + volume = {76}, + issn = {0033-3123}, + url = {https://doi.org/10.1007/s11336-011-9224-6}, + doi = {10.1007/s11336-011-9224-6}, + number = {4}, + journal = {Psychometrika}, + author = {Yuan, Ke-Hai and Chan, Wai}, + month = aug, + year = {2011}, + keywords = {Confidence Intervals, Regression, Nonnormal Variables, Sampling Distribution, Standard Error, Standardized Solution, Fixed Predictors}, + pages = {670--690}, +} + +@article{dudgeon_improvements_2017, + title = {Some improvements in confidence intervals for standardized regression coefficients}, + volume = {82}, + issn = {0033-3123, 1860-0980}, + url = {http://link.springer.com/10.1007/s11336-017-9563-z}, + doi = {10.1007/s11336-017-9563-z}, + language = {en}, + number = {4}, + urldate = {2020-04-18}, + journal = {Psychometrika}, + author = {Dudgeon, Paul}, + month = dec, + year = {2017}, + keywords = {Confidence Intervals, Regression, Simulation Studies, Nonnormal Variables, Standardized Solution}, + pages = {928--951}, +} + +@article{cheung_improving_2022, + title = {Improving an old way to measure moderation effect in standardized units}, + volume = {41}, + issn = {1930-7810, 0278-6133}, + url = {https://doi.org/10.1037/hea0001188}, + doi = {10.1037/hea0001188}, + language = {en}, + number = {7}, + urldate = {2022-08-02}, + journal = {Health Psychology}, + author = {Cheung, Shu Fai and Cheung, Sing-Hang and Lau, Esther Yuet Ying and Hui, C. Harry and Vong, Weng Ngai}, + month = jul, + year = {2022}, + keywords = {SEM (Structural Equation Modeling), Regression, Moderation (Interaction), FTBC, sfcheung, sfcheung (First or Correspondence), Effect Sizes, Standardized Solution}, + pages = {502--505}, +} + +@article{jones_computing_2013, + title = {Computing confidence intervals for standardized regression coefficients}, + volume = {18}, + issn = {1939-1463, 1082-989X}, + url = {http://doi.apa.org/getdoi.cfm?doi=10.1037/a0033269}, + doi = {10.1037/a0033269}, + abstract = {With fixed predictors, the standard method (Cohen, Cohen, West, \& Aiken, 2003, p. 86; Harris, 2001, p. 80; Hays, 1994, p. 709) for computing confidence intervals (CIs) for standardized regression coefficients fails to account for the sampling variability of the criterion standard deviation. With random predictors, this method also fails to account for the sampling variability of the predictor standard deviations. Nevertheless, under some conditions the standard method will produce CIs with accurate coverage rates. To delineate these conditions, we used a Monte Carlo simulation to compute empirical CI coverage rates in samples drawn from 36 populations with a wide range of data characteristics. We also computed the empirical CI coverage rates for 4 alternative methods that have been discussed in the literature: noncentrality interval estimation, the delta method, the percentile bootstrap, and the bias-corrected and accelerated bootstrap. Our results showed that for many data-parameter configurations—for example, sample size, predictor correlations, coefficient of determination (R2), orientation of ␤ with respect to the eigenvectors of the predictor correlation matrix, RX—the standard method produced coverage rates that were close to their expected values. However, when population R2 was large and when ␤ approached the last eigenvector of RX, then the standard method coverage rates were frequently below the nominal rate (sometimes by a considerable amount). In these conditions, the delta method and the 2 bootstrap procedures were consistently accurate. Results using noncentrality interval estimation were inconsistent. In light of these findings, we recommend that researchers use the delta method to evaluate the sampling variability of standardized regression coefficients.}, + language = {en}, + number = {4}, + urldate = {2020-04-18}, + journal = {Psychological Methods}, + author = {Jones, Jeff A. and Waller, Niels G.}, + year = {2013}, + keywords = {Bootstrapping, Confidence Intervals, Regression, Simulation Studies, Standardized Solution, Fixed Predictors}, + pages = {435--453}, +} + +@article{jones_normal-theory_2015, + title = {The normal-theory and asymptotic distribution-free ({ADF}) covariance matrix of standardized regression coefficients: {Theoretical} extensions and finite sample behavior}, + volume = {80}, + copyright = {http://www.springer.com/tdm}, + issn = {0033-3123, 1860-0980}, + shorttitle = {The {Normal}-{Theory} and {Asymptotic} {Distribution}-{Free} ({ADF}) {Covariance} {Matrix} of {Standardized} {Regression} {Coefficients}}, + url = {http://link.springer.com/10.1007/s11336-013-9380-y}, + doi = {10.1007/s11336-013-9380-y}, + language = {en}, + number = {2}, + urldate = {2024-06-03}, + journal = {Psychometrika}, + author = {Jones, Jeff A. and Waller, Niels G.}, + month = jun, + year = {2015}, + keywords = {Nonnormal Variables, Sampling Distribution, Standard Error, Standardized Solution, ADF (Asymptotically Distribution Free)}, + pages = {365--378}, +} From 8c119677e27fc6cf184ccefe298e26006a85061b Mon Sep 17 00:00:00 2001 From: Shu Fai Cheung Date: Tue, 29 Oct 2024 21:22:56 +0800 Subject: [PATCH 8/8] Draft the vignette for lm_betaselect() --- vignettes/betaselectr_lm.Rmd | 619 ++++++++++++++++++++++++++ vignettes/betaselectr_lm.Rmd.original | 361 +++++++-------- 2 files changed, 788 insertions(+), 192 deletions(-) create mode 100644 vignettes/betaselectr_lm.Rmd diff --git a/vignettes/betaselectr_lm.Rmd b/vignettes/betaselectr_lm.Rmd new file mode 100644 index 0000000..7ca2e58 --- /dev/null +++ b/vignettes/betaselectr_lm.Rmd @@ -0,0 +1,619 @@ +--- +title: "Beta-Select Demonstration: Regression by `lm()`" +date: "2024-10-29" +output: + rmarkdown::html_vignette: + number_sections: true +vignette: > + %\VignetteIndexEntry{Beta-Select Demonstration: Regression by `lm()`} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +bibliography: "references.bib" +csl: apa.csl +--- + + + + + +# Introduction + +This article demonstrates how to use +`lm_betaselect()` from the package +`betaselectr` +to standardize +selected variables in a model fitted +by `lm()` and forming confidence +intervals for the parameters. + +# Data and Model + +The sample dataset from the package +`betaselectr` will be used for in this +demonstration: + + +``` r +library(betaselectr) +head(data_test_mod_cat2) +#> dv iv mod cov1 cat1 +#> 1 15.53 13.95 50.75 25.33 gp2 +#> 2 17.69 15.07 49.67 20.96 gp1 +#> 3 28.56 14.43 53.42 19.22 gp3 +#> 4 25.00 11.22 42.55 20.18 gp2 +#> 5 19.33 14.93 52.12 22.82 gp2 +#> 6 20.62 10.22 39.36 18.41 gp1 +``` + +This is the regression model, fitted by +`lm()`: + + +``` r +lm_out <- lm(dv ~ iv * mod + cov1 + cat1, + data = data_test_mod_cat2) +``` + +The model has a moderator, `mod`, posited +to moderate the effect from `iv` to +`med`. The product term is `iv:mod`. +The variable `cat1` is a categorical variable +with three groups: `gp1`, `gp2`, `gp3`. + +These are the results: + + +``` r +summary(lm_out) +#> +#> Call: +#> lm(formula = dv ~ iv * mod + cov1 + cat1, data = data_test_mod_cat2) +#> +#> Residuals: +#> Min 1Q Median 3Q Max +#> -17.0892 -4.6312 0.0057 5.0186 18.7053 +#> +#> Coefficients: +#> Estimate Std. Error t value Pr(>|t|) +#> (Intercept) 90.87211 34.04462 2.669 0.00803 ** +#> iv -6.06932 2.33701 -2.597 0.00988 ** +#> mod -1.61636 0.68840 -2.348 0.01954 * +#> cov1 0.09885 0.19433 0.509 0.61136 +#> cat1gp2 1.71248 1.15064 1.488 0.13775 +#> cat1gp3 2.47838 1.10562 2.242 0.02574 * +#> iv:mod 0.13230 0.04656 2.841 0.00481 ** +#> --- +#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +#> +#> Residual standard error: 6.743 on 293 degrees of freedom +#> Multiple R-squared: 0.1149, Adjusted R-squared: 0.09676 +#> F-statistic: 6.338 on 6 and 293 DF, p-value: 2.759e-06 +``` + +# Problems With Standardization + +One common type of standardized coefficients, +called "betas" in some programs, is +computed by standardizing *all* terms +in the model. + +First, all variables in the model, +including the product term and dummy +variables, are computed: + + +``` r +data_test_mod_cat2_z <- data_test_mod_cat2 +data_test_mod_cat2_z$iv_x_mod <- data_test_mod_cat2_z$iv * + data_test_mod_cat2_z$mod +data_test_mod_cat2_z$cat_gp2 <- as.numeric(data_test_mod_cat2_z$cat1 == "gp2") +data_test_mod_cat2_z$cat_gp3 <- as.numeric(data_test_mod_cat2_z$cat1 == "gp3") +head(data_test_mod_cat2_z) +#> dv iv mod cov1 cat1 iv_x_mod cat_gp2 cat_gp3 +#> 1 15.53 13.95 50.75 25.33 gp2 707.9625 1 0 +#> 2 17.69 15.07 49.67 20.96 gp1 748.5269 0 0 +#> 3 28.56 14.43 53.42 19.22 gp3 770.8506 0 1 +#> 4 25.00 11.22 42.55 20.18 gp2 477.4110 1 0 +#> 5 19.33 14.93 52.12 22.82 gp2 778.1516 1 0 +#> 6 20.62 10.22 39.36 18.41 gp1 402.2592 0 0 +``` + +All the variables are then standardized: + + +``` r +data_test_mod_cat2_z <- data.frame(scale(data_test_mod_cat2_z[, -5])) +head(data_test_mod_cat2_z) +#> dv iv mod cov1 iv_x_mod cat_gp2 +#> 1 -0.9458226 -0.44874323 0.23147783 2.553777460 -0.24816181 1.331109 +#> 2 -0.6414005 0.13926755 -0.03378874 0.390649940 0.06187143 -0.748749 +#> 3 0.8905756 -0.19673861 0.88727574 -0.470641109 0.23249121 -0.748749 +#> 4 0.3888429 -1.88201951 -1.78258317 0.004553953 -2.01026427 1.331109 +#> 5 -0.4102652 0.06576621 0.56797339 1.311340372 0.28829267 1.331109 +#> 6 -0.2284575 -2.40702913 -2.56610202 -0.871586942 -2.58464861 -0.748749 +#> cat_gp3 +#> 1 -0.9401258 +#> 2 -0.9401258 +#> 3 1.0601418 +#> 4 -0.9401258 +#> 5 -0.9401258 +#> 6 -0.9401258 +``` + +The regression model is then fitted to the +standardized variables: + + +``` r +lm_std_common <- lm(dv ~ iv + mod + cov1 + cat_gp2 + cat_gp3 + iv_x_mod, + data = data_test_mod_cat2_z) +``` + +The "betas" commonly reported are the +coefficients in this model: + + +``` r +lm_std_common_summary <- summary(lm_std_common) +printCoefmat(lm_std_common_summary$coefficients, + digits = 5, + zap.ind = 1, + P.values = TRUE, + has.Pvalue = TRUE, + signif.stars = TRUE) +#> Estimate Std. Error t value Pr(>|t|) +#> (Intercept) 0.000000 0.054871 0.0000 1.000000 +#> iv -1.629280 0.627360 -2.5970 0.009878 ** +#> mod -0.927480 0.395006 -2.3480 0.019539 * +#> cov1 0.028140 0.055329 0.5087 0.611359 +#> cat_gp2 0.116040 0.077970 1.4883 0.137753 +#> cat_gp3 0.174620 0.077901 2.2416 0.025735 * +#> iv_x_mod 2.439510 0.858601 2.8413 0.004809 ** +#> --- +#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +``` + +However, for this model, there are several +problems: + +- The product term, `iv:mod`, is also + standardized. This is inappropriate. + One simple but underused solution is + standardizing the variables *before* + forming the product term [@friedrich_defense_1982]. + +- The confidence intervals are formed using + the ordinary least squares (OLS), which does not + take into account the sampling variation + of the standardizers (the sample standard + deviations used in standardization) and + so the standard errors may be + biased [@yuan_biases_2011]. + Although there + are situations in which the OLS + confidence and the nonparametric + percentile bootstrap confidences can be + similar (e.g., sample size is large + and the population values are not extreme), + it is recommended to use bootstrap + confidence intervals when computation + cost is low [@jones_computing_2013]. + +- There are cases in which some variabLes + are measured by meaningful units and + do not need to be standardized. for + example, if `cov1` is age measured by + year, then age is more + meaningful than "standardized age". + +- In regression models, categorical variables + are usually represented by dummy variables, + each of them having only two possible + values (0 or 1). It is not meaningful + to standardize the dummy variables. + +# Beta-Select by `lm_betaselect()` + +The function `lm_betaselect()` can be used +to solve these problems by: + +- standardizing variables before product + terms are formed, + +- standardizing only variables for which + standardization can facilitate + interpretation, and + +- forming bootstrap confidence intervals + that take + into account selected standardization. + +We call the coefficients computed by +this kind of standardization *beta*s-Select +($\beta{s}_{Select}$, $\beta_{Select}$ +in singular form), +to differentiate them from coefficients +computed by standardizing all variables, +including product terms. + +## Estimates Only + +Suppose we only need to +solve the first problem, standardizing all +numeric variables, +with the product +term computed after `iv`, `mod`, and `dv` +are standardized. + + +``` r +lm_beta_select <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat2, + do_boot = FALSE) +``` + +The function `lm_beta_iv_mod()` can be +used as `lm()`, with applicable arguments +such as the model formula and `data` passed +to `lm()`. + +By default, *all* numeric variables will +be standardized before fitting the models. +Terms such as product terms are created +*after* standardization. + +Moreover, categorical variables (factors and +string variables) will not be standardized. + +Bootstrapping is done by default. In this +illustration, `do_boot = FALSE` is added +to disabled it because we only want to +address the first problem. We will do bootstrapping when +addressing the issue with confidence intervals. + +The `summary()` method can be used +ont the output of `lm_betaselect()`: + + +``` r +summary(lm_beta_select) +#> Warning in summary.lm_betaselect(lm_beta_select): Bootstrap estimates not +#> available; 'se_method' changed to 'ls'. +#> Call to lm_betaselect(): +#> betaselectr::lm_betaselect(formula = dv ~ iv * mod + cov1 + cat1, +#> data = data_test_mod_cat2, do_boot = FALSE) +#> +#> Variable(s) standardized: dv, iv, mod, cov1 +#> +#> Call: +#> stats::lm(formula = dv ~ iv * mod + cov1 + cat1, data = betaselectr::std_data(data = data_test_mod_cat2, +#> to_standardize = c("dv", "iv", "mod", "cov1"))) +#> +#> Residuals: +#> Min 1Q Median 3Q Max +#> -2.4085 -0.6527 0.0008 0.7073 2.6363 +#> +#> Coefficients: +#> Estimate CI.Lower CI.Upper Std. Error t value Pr(>|t|) +#> (Intercept) -0.308 -0.576 -0.040 0.136 -2.259 0.02461 * +#> iv 0.140 0.021 0.258 0.060 2.324 0.02082 * +#> mod 0.196 0.078 0.315 0.060 3.264 0.00123 ** +#> cov1 0.028 -0.081 0.137 0.055 0.509 0.61136 +#> cat1gp2 0.241 -0.078 0.561 0.162 1.488 0.13775 +#> cat1gp3 0.349 0.043 0.656 0.156 2.242 0.02574 * +#> iv:mod 0.145 0.044 0.245 0.051 2.841 0.00481 ** +#> --- +#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +#> +#> Residual standard error: 0.9504 on 293 degrees of freedom +#> +#> R-squared : 0.115 +#> Adjusted R-squared : 0.097 +#> ANOVA test of R-squared : F(6, 293) = 6.338, p < 0.001 +#> +#> Note: +#> - Results *after* standardization are reported. +#> - Standard errors are least-squares standard errors. +#> - T values are computed by 'Estimate / Std. Error'. +#> - P-values are usual t-test p-values. +#> - Least squares standard errors, t values, p-values, and confidence +#> intervals (if reported) should not be used for coefficients involved +#> in standardization. +#> - Least squares 95.0% confidence interval reported. +``` + +A warning is always issued if standardization +is done but bootstrapping not requested. + + + + + +Compared to the solution with the product +term standardized, the coefficient of +`iv:mod` changed substantially from +2.440 to +0.145. As shown by +@cheung_improving_2022, the coefficient +of *standardized* product term (`iv:mod`) +can be a severely biased estimate of the +properly standardized product term +(the product of standardized `iv` and +standardized `mod`). + +## Estimates and Bootstrap Confidence Interval + +Suppose we want to address +both the first and the second problems, +with the product +term computed after `iv`, `mod`, and `dv` are +standardized and bootstrap confidence +interval used, we can call `lm_betaselect()` +again, with additional arguments +set: + + +``` r +lm_beta_select_boot <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat2, + bootstrap = 5000, + iseed = 4567) +``` + +These are the additional arguments: + +- `bootstrap`: The number of bootstrap + samples to draw. Default is 100. It should + be set to 5000 or even 10000. + +- `iseed`: The seed for the random number + generator used for bootstrapping. Set + this to an integer to make the results + reproducible. + + + + +This is the output of `summary()` + + +``` r +summary(lm_beta_select_boot) +#> Call to lm_betaselect(): +#> betaselectr::lm_betaselect(formula = dv ~ iv * mod + cov1 + cat1, +#> data = data_test_mod_cat2, bootstrap = 5000, iseed = 4567) +#> +#> Variable(s) standardized: dv, iv, mod, cov1 +#> +#> Call: +#> stats::lm(formula = dv ~ iv * mod + cov1 + cat1, data = betaselectr::std_data(data = data_test_mod_cat2, +#> to_standardize = c("dv", "iv", "mod", "cov1"))) +#> +#> Residuals: +#> Min 1Q Median 3Q Max +#> -2.4085 -0.6527 0.0008 0.7073 2.6363 +#> +#> Coefficients: +#> Estimate CI.Lower CI.Upper Std. Error z value Pr(Boot) +#> (Intercept) -0.308 -0.536 -0.080 0.117 -2.636 0.0056 ** +#> iv 0.140 0.009 0.268 0.066 2.109 0.0384 * +#> mod 0.196 0.075 0.317 0.061 3.208 0.0016 ** +#> cov1 0.028 -0.075 0.131 0.052 0.537 0.5700 +#> cat1gp2 0.241 -0.067 0.540 0.155 1.560 0.1276 +#> cat1gp3 0.349 0.064 0.631 0.146 2.394 0.0152 * +#> iv:mod 0.145 0.059 0.228 0.043 3.356 0.0020 ** +#> --- +#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +#> +#> Residual standard error: 0.9504 on 293 degrees of freedom +#> +#> R-squared : 0.115 +#> Adjusted R-squared : 0.097 +#> ANOVA test of R-squared : F(6, 293) = 6.338, p < 0.001 +#> +#> Note: +#> - Results *after* standardization are reported. +#> - Nonparametric bootstrapping conducted. +#> - The number of bootstrap samples is 5000. +#> - Standard errors are bootstrap standard errors. +#> - Z values are computed by 'Estimate / Std. Error'. +#> - The bootstrap p-values are asymmetric p-values by Asparouhov and +#> Muthén (2021). +#> - Percentile bootstrap 95.0% confidence interval reported. +``` + +By default, 95% percentile bootstrap +confidence intervals are printed +(`CI.Lower` and `CI.Upper`). The *p*-values +(`Pr(Boot)`) are asymmetric bootstrap +*p*-values. + +## Estimates and Bootstrap Confidence Intervals, With Only Selected Variables Standardized + +Suppose we want to address also the +the third issue, and standardize only +some of the variables. This can be +done using either `to_standardize` +or `not_to_standardize`. + +- Use `to_standardize` when +the variables to standardize +is much fewer than the variables +not to standardize. + +- Use `not_to_standardize` +when the variables to standardize +is much more than the +variables not to standardize. + +For example, suppose we only +need to standardize `dv` and +`iv`, +this is the call to do +this, setting +`to_standardize` to `c("iv", "dv")`: + + +``` r +lm_beta_select_boot_1 <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat2, + to_standardize = c("dv", "iv"), + bootstrap = 5000, + iseed = 4567) +``` + +If we want to standardize all +variables except for `mod` +and `cov1`, we can use +this call, and set +`not_to_standardize` to `c("mod", "cov1")`: + + +``` r +lm_beta_select_boot_2 <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat2, + not_to_standardize = c("mod", "cov1"), + bootstrap = 5000, + iseed = 4567) +``` + +The results of these calls are identical, +and only those of the first version are +printed: + + +``` r +summary(lm_beta_select_boot_1) +#> Call to lm_betaselect(): +#> betaselectr::lm_betaselect(formula = dv ~ iv * mod + cov1 + cat1, +#> data = data_test_mod_cat2, to_standardize = c("dv", "iv"), +#> bootstrap = 5000, iseed = 4567) +#> +#> Variable(s) standardized: dv, iv +#> +#> Call: +#> stats::lm(formula = dv ~ iv * mod + cov1 + cat1, data = betaselectr::std_data(data = data_test_mod_cat2, +#> to_standardize = c("dv", "iv"))) +#> +#> Residuals: +#> Min 1Q Median 3Q Max +#> -2.4085 -0.6527 0.0008 0.7073 2.6363 +#> +#> Coefficients: +#> Estimate CI.Lower CI.Upper Std. Error z value Pr(Boot) +#> (Intercept) -2.991 -4.769 -1.196 0.899 -3.326 0.0024 ** +#> iv -1.629 -2.667 -0.573 0.539 -3.021 0.0036 ** +#> mod 0.048 0.019 0.078 0.015 3.199 0.0016 ** +#> cov1 0.014 -0.037 0.066 0.026 0.533 0.5700 +#> cat1gp2 0.241 -0.067 0.540 0.155 1.560 0.1276 +#> cat1gp3 0.349 0.064 0.631 0.146 2.394 0.0152 * +#> iv:mod 0.036 0.015 0.056 0.011 3.366 0.0020 ** +#> --- +#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +#> +#> Residual standard error: 0.9504 on 293 degrees of freedom +#> +#> R-squared : 0.115 +#> Adjusted R-squared : 0.097 +#> ANOVA test of R-squared : F(6, 293) = 6.338, p < 0.001 +#> +#> Note: +#> - Results *after* standardization are reported. +#> - Nonparametric bootstrapping conducted. +#> - The number of bootstrap samples is 5000. +#> - Standard errors are bootstrap standard errors. +#> - Z values are computed by 'Estimate / Std. Error'. +#> - The bootstrap p-values are asymmetric p-values by Asparouhov and +#> Muthén (2021). +#> - Percentile bootstrap 95.0% confidence interval reported. +``` + +For *beta*s-*select*, researchers need +to state which variables +are standardized and which are not. +This can be done in table notes. + +## Categorical Variables + +When calling `lm_betaselect()`, +categorical variables (factors and +string variables) will not be standardized +by default. +This can be overriden by setting +`skip_categorical_x` to `FALSE`, though +not recommended. + +In the example above, the coefficients +of the two dummy variables when both +the dummy variables and the outcome +variables are standardized are +0.116 and +0.175: + + +``` r +printCoefmat(lm_std_common_summary$coefficients[5:6, ], + digits = 5, + zap.ind = 1, + P.values = TRUE, + has.Pvalue = TRUE, + signif.stars = TRUE) +#> Estimate Std. Error t value Pr(>|t|) +#> cat_gp2 0.116041 0.077970 1.4883 0.13775 +#> cat_gp3 0.174623 0.077901 2.2416 0.02574 * +#> --- +#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +``` + +These two values are not interpretable +because it does not make sense to talk +about a one-SD change in the dummy variables. + +The *beta*s-*Select* of the dummy variables, +with only the outcome variable standardized, +are +0.241 and +0.349. + + +``` r +lm_beta_select_boot_summary <- summary(lm_beta_select_boot) +printCoefmat(lm_beta_select_boot_summary$coefficients[5:6, ], + digits = 5, + zap.ind = 1, + P.values = TRUE, + has.Pvalue = TRUE, + signif.stars = TRUE) +#> Estimate CI.Lower CI.Upper Std. Error z value Pr(Boot) +#> cat1gp2 0.241350 -0.067312 0.540477 0.154668 1.5604 0.1276 +#> cat1gp3 0.349290 0.064136 0.630539 0.145927 2.3936 0.0152 * +#> --- +#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +``` + +They +represent the difference between `gp2` +and `gp3` from the reference group, +`gp1`, on the *standardized* outcome variable. +That is, their meanings are the same before +and after standardization. The only difference +is in the unit of the outcome variable. + +# Conclusion + +In regression analysis, there +are situations in which standardizing +all variables is not appropriate, or +when standardization needs to be done +before forming product terms. We are +not aware of tools that can do appropriate +standardization *and* form confidence +intervals that takes into account the +selective standardization. By promoting +the use of *beta*s-*select* using +`lm_betaselect()`, we hope to make it +easier for researchers to do appropriate +standardization when reporting regression +results. + +# References diff --git a/vignettes/betaselectr_lm.Rmd.original b/vignettes/betaselectr_lm.Rmd.original index d759930..c390c48 100644 --- a/vignettes/betaselectr_lm.Rmd.original +++ b/vignettes/betaselectr_lm.Rmd.original @@ -8,6 +8,8 @@ vignette: > %\VignetteIndexEntry{Beta-Select Demonstration: Regression by `lm()`} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} +bibliography: "references.bib" +csl: apa.csl --- ```{r, include = FALSE} @@ -42,7 +44,7 @@ demonstration: ```{r} library(betaselectr) -head(data_test_mod_cat) +head(data_test_mod_cat2) ``` This is the regression model, fitted by @@ -50,7 +52,7 @@ This is the regression model, fitted by ```{r} lm_out <- lm(dv ~ iv * mod + cov1 + cat1, - data = data_test_mod_cat) + data = data_test_mod_cat2) ``` The model has a moderator, `mod`, posited @@ -77,19 +79,19 @@ including the product term and dummy variables, are computed: ```{r} -data_test_mod_cat_z <- data_test_mod_cat -data_test_mod_cat_z$iv_x_mod <- data_test_mod_cat_z$iv * - data_test_mod_cat_z$mod -data_test_mod_cat_z$cat_gp2 <- as.numeric(data_test_mod_cat_z$cat1 == "gp2") -data_test_mod_cat_z$cat_gp3 <- as.numeric(data_test_mod_cat_z$cat1 == "gp3") -head(data_test_mod_cat_z) +data_test_mod_cat2_z <- data_test_mod_cat2 +data_test_mod_cat2_z$iv_x_mod <- data_test_mod_cat2_z$iv * + data_test_mod_cat2_z$mod +data_test_mod_cat2_z$cat_gp2 <- as.numeric(data_test_mod_cat2_z$cat1 == "gp2") +data_test_mod_cat2_z$cat_gp3 <- as.numeric(data_test_mod_cat2_z$cat1 == "gp3") +head(data_test_mod_cat2_z) ``` All the variables are then standardized: ```{r} -data_test_mod_cat_z <- data.frame(scale(data_test_mod_cat_z[, -5])) -head(data_test_mod_cat_z) +data_test_mod_cat2_z <- data.frame(scale(data_test_mod_cat2_z[, -5])) +head(data_test_mod_cat2_z) ``` The regression model is then fitted to the @@ -97,7 +99,7 @@ standardized variables: ```{r} lm_std_common <- lm(dv ~ iv + mod + cov1 + cat_gp2 + cat_gp3 + iv_x_mod, - data = data_test_mod_cat_z) + data = data_test_mod_cat2_z) ``` The "betas" commonly reported are the @@ -113,79 +115,59 @@ printCoefmat(lm_std_common_summary$coefficients, signif.stars = TRUE) ``` -```{r} -library(lm.beta) -lm_beta <- lm.beta(lm_out) -lm_beta_summary <- summary(lm_beta) -printCoefmat(lm_beta_summary$coefficients, - digits = 4, - zap.ind = 1, - P.values = TRUE, - has.Pvalue = TRUE, - signif.stars = TRUE) -``` - -```{r, results = FALSE} -lm_beta_x <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, - data = data_test_mod_cat, - do_boot = TRUE, - bootstrap = 5000, - iseed = 4567) -``` - -```{r} -summary(lm_beta_x) -``` - However, for this model, there are several problems: - The product term, `iv:mod`, is also standardized. This is inappropriate. One simple but underused solution is - standardized the variables *before* - forming the product term (Friedrich, 1982). + standardizing the variables *before* + forming the product term [@friedrich_defense_1982]. - The confidence intervals are formed using - the delta-method, which has been found - to be inferior to methods such as - nonparametric percentile bootstrap - confidence interval for the standardized - solution (Falk, 2018). Although there - are situations in which the delta-method + the ordinary least squares (OLS), which does not + take into account the sampling variation + of the standardizers (the sample standard + deviations used in standardization) and + so the standard errors may be + biased [@yuan_biases_2011]. + Although there + are situations in which the OLS confidence and the nonparametric percentile bootstrap confidences can be similar (e.g., sample size is large - and the sample estimates are not extreme), - it is still safe to at least try both - methods and compare the results. + and the population values are not extreme), + it is recommended to use bootstrap + confidence intervals when computation + cost is low [@jones_computing_2013]. - There are cases in which some variabLes are measured by meaningful units and do not need to be standardized. for example, if `cov1` is age measured by year, then age is more - meaningful than the "standardized age". + meaningful than "standardized age". -- In path analysis, categorical variables +- In regression models, categorical variables are usually represented by dummy variables, each of them having only two possible values (0 or 1). It is not meaningful to standardize the dummy variables. -# Beta-Select `lav_betaselect()` +# Beta-Select by `lm_betaselect()` -The function `lav_betaselect()` can be used -to solve this problem by: +The function `lm_betaselect()` can be used +to solve these problems by: -- Standardizing variables before product - terms are formed. +- standardizing variables before product + terms are formed, -- Standardizing only variables for which +- standardizing only variables for which standardization can facilitate - interpretation. + interpretation, and -- Forming confidence intervals that take +- forming bootstrap confidence intervals + that take into account selected standardization. We call the coefficients computed by @@ -199,122 +181,113 @@ including product terms. ## Estimates Only Suppose we only need to -solve the first problem, with the product -term computed after `iv` and `mod` -are standardized: +solve the first problem, standardizing all +numeric variables, +with the product +term computed after `iv`, `mod`, and `dv` +are standardized. ```{r, results = FALSE} -fit_beta <- lav_betaselect(fit) +lm_beta_select <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat2, + do_boot = FALSE) ``` -```{r, echo = FALSE} -tmp <- capture.output(print(fit_beta)) -``` +The function `lm_beta_iv_mod()` can be +used as `lm()`, with applicable arguments +such as the model formula and `data` passed +to `lm()`. -```{r, eval = FALSE} -fit_beta +By default, *all* numeric variables will +be standardized before fitting the models. +Terms such as product terms are created +*after* standardization. + +Moreover, categorical variables (factors and +string variables) will not be standardized. + +Bootstrapping is done by default. In this +illustration, `do_boot = FALSE` is added +to disabled it because we only want to +address the first problem. We will do bootstrapping when +addressing the issue with confidence intervals. + +The `summary()` method can be used +ont the output of `lm_betaselect()`: + +```{r} +summary(lm_beta_select) ``` -This is the output -if printed using the -default options: +A warning is always issued if standardization +is done but bootstrapping not requested. ```{r, echo = FALSE} -cat(tmp[-c(25:51)], sep = "\n") +tmp <- capture.output(suppressWarnings(print(summary(lm_beta_select)))) ``` ```{r, echo = FALSE} -b_std <- standardizedSolution(fit)$est.std[3] -b_select <- fit_beta$est[3] +b_select <- coef(lm_beta_select) +b_std <- coef(lm_std_common) ``` Compared to the solution with the product term standardized, the coefficient of `iv:mod` changed substantially from -`r format_str(b_std)` to -`r format_str(b_select)`. As shown by -Cheung et al. (2022), the coefficient +`r format_str(b_std["iv_x_mod"])` to +`r format_str(b_select["iv:mod"])`. As shown by +@cheung_improving_2022, the coefficient of *standardized* product term (`iv:mod`) -can be severely biased estimate of the +can be a severely biased estimate of the properly standardized product term (the product of standardized `iv` and standardized `mod`). -The footnote will also indicate -variables that are standardized, -and remarked that product terms -are formed *after* standardization. - ## Estimates and Bootstrap Confidence Interval Suppose we want to address both the first and the second problems, with the product -term computed after `iv` and `mod` +term computed after `iv`, `mod`, and `dv` are standardized and bootstrap confidence -interval used, we can call `lav_betaselect()` +interval used, we can call `lm_betaselect()` again, with additional arguments set: -```{r} -fit_beta <- lav_betaselect(fit, - std_se = "bootstrap", - bootstrap = 5000, - iseed = 2345, - parallel = "snow", - ncpus = 20) +```{r results = FALSE} +lm_beta_select_boot <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat2, + bootstrap = 5000, + iseed = 4567) ``` These are the additional arguments: -- `std_se`: The method to compute the - standard errors as well as confidence - intervals. Set to `"bootstrap"` for - nonparametric bootstrapping. +- `bootstrap`: The number of bootstrap + samples to draw. Default is 100. It should + be set to 5000 or even 10000. - `iseed`: The seed for the random number generator used for bootstrapping. Set this to an integer to make the results reproducible. -- `parallel`: The method to be used for - parallel processing. It will be passed - to `lavaan::bootstrapLavaan()`. Supported - values are `"none"`, `"snow"`, and - `"multicore"`. - -- `ncpus`: The number of CPU cores to - use if `parallel` processing is not - `"none"`. Default is - `parallel::detectCores(logical = FALSE) - 1`, - or the number of *physical* cores - minus one. ```{r, echo = FALSE} -tmp <- capture.output(print(fit_beta)) +tmp <- capture.output(suppressWarnings(print(summary(lm_beta_select_boot)))) ``` -This is the output if -printed with default -options: +This is the output of `summary()` -```{r, eval = FALSE} -fit_beta -``` - -```{r, echo = FALSE} -cat(tmp[c(1:27, 55:66)], sep = "\n") +```{r} +summary(lm_beta_select_boot) ``` -In this dataset, with 200 cases, -the delta-method confidence -intervals are close to the -bootstrap confidence intervals, -except obviously for the -product term because the coefficient -of the product term has substantially -different values in the two -solutions. +By default, 95% percentile bootstrap +confidence intervals are printed +(`CI.Lower` and `CI.Upper`). The *p*-values +(`Pr(Boot)`) are asymmetric bootstrap +*p*-values. ## Estimates and Bootstrap Confidence Intervals, With Only Selected Variables Standardized @@ -336,91 +309,103 @@ variables not to standardize. For example, suppose we only need to standardize `dv` and -`iv`, `cov1`, and `cov2`, +`iv`, this is the call to do this, setting -`to_standardize` to `c("iv", "dv", "cov1", "cov2")`: - -```{r, results = FALSE} -fit_beta_select_1 <- lav_betaselect(fit, - std_se = "bootstrap", - to_standardize = c("iv", "dv", "cov1", "cov2"), - bootstrap = 5000, - iseed = 2345, - parallel = "snow", - ncpus = 20) +`to_standardize` to `c("iv", "dv")`: + +```{r results = FALSE} +lm_beta_select_boot_1 <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat2, + to_standardize = c("dv", "iv"), + bootstrap = 5000, + iseed = 4567) ``` If we want to standardize all -variables except for `dv` -and `mod`, we can use +variables except for `mod` +and `cov1`, we can use this call, and set -`not_to_standardize` to `c("mod", "dv")`: +`not_to_standardize` to `c("mod", "cov1")`: ```{r, results = FALSE} -fit_beta_select_2 <- lav_betaselect(fit, - std_se = "bootstrap", - not_to_standardize = c("mod", "dv"), - bootstrap = 5000, - iseed = 2345, - parallel = "snow", - ncpus = 20) +lm_beta_select_boot_2 <- lm_betaselect(dv ~ iv*mod + cov1 + cat1, + data = data_test_mod_cat2, + not_to_standardize = c("mod", "cov1"), + bootstrap = 5000, + iseed = 4567) ``` The results of these calls are identical, and only those of the first version are printed: -```{r, eval = FALSE} -fit_beta_select_2 -``` - -```{r, echo = FALSE} -tmp <- capture.output(print(fit_beta_select_2)) -cat(tmp[c(2:27, 55:66)], sep = "\n") +```{r} +summary(lm_beta_select_boot_1) ``` -The footnotes confirmed that, by -specifying that `dv` and `mod` are not -standardized, all the other four variables -are standardized: `iv`, `med`, `cov1`, and `cov2`. -Therefore, in this case, it is more -convenient to use `not_to_standardize`. - For *beta*s-*select*, researchers need to state which variables are standardized and which are not. -This can be done in table notes, -or in a column of the parameter estimate -tables. The output can of `lav_betaselect()` -can be printed with `show_Bs.by` set -to `TRUE` to demonstrate the second -approach: - -```{r, eval = FALSE} -print(fit_beta_select_2, - show_Bs.by = TRUE) -``` - -```{r, echo = FALSE} -tmp <- capture.output(print(fit_beta_select_2, - show_Bs.by = TRUE)) -cat(tmp[c(15:27, 55:70)], sep = "\n") -``` +This can be done in table notes. ## Categorical Variables -When calling `lav_betaselect()`, -variables with only two values in -the dataset are assumed to be categorical -and will not be standardized by default. +When calling `lm_betaselect()`, +categorical variables (factors and +string variables) will not be standardized +by default. This can be overriden by setting `skip_categorical_x` to `FALSE`, though not recommended. +In the example above, the coefficients +of the two dummy variables when both +the dummy variables and the outcome +variables are standardized are +`r format_str(b_std["cat_gp2"])` and +`r format_str(b_std["cat_gp3"])`: + +```{r} +printCoefmat(lm_std_common_summary$coefficients[5:6, ], + digits = 5, + zap.ind = 1, + P.values = TRUE, + has.Pvalue = TRUE, + signif.stars = TRUE) +``` + +These two values are not interpretable +because it does not make sense to talk +about a one-SD change in the dummy variables. + +The *beta*s-*Select* of the dummy variables, +with only the outcome variable standardized, +are +`r format_str(b_select["cat1gp2"])` and +`r format_str(b_select["cat1gp3"])`. + +```{r} +lm_beta_select_boot_summary <- summary(lm_beta_select_boot) +printCoefmat(lm_beta_select_boot_summary$coefficients[5:6, ], + digits = 5, + zap.ind = 1, + P.values = TRUE, + has.Pvalue = TRUE, + signif.stars = TRUE) +``` + +They +represent the difference between `gp2` +and `gp3` from the reference group, +`gp1`, on the *standardized* outcome variable. +That is, their meanings are the same before +and after standardization. The only difference +is in the unit of the outcome variable. + # Conclusion -In structural equation modeling, there +In regression analysis, there are situations in which standardizing all variables is not appropriate, or when standardization needs to be done @@ -428,19 +413,11 @@ before forming product terms. We are not aware of tools that can do appropriate standardization *and* form confidence intervals that takes into account the -selective Standardization. By promoting +selective standardization. By promoting the use of *beta*s-*select* using -`lav_betaselect()`, we hope to make it +`lm_betaselect()`, we hope to make it easier for researchers to do appropriate -Standardization in when reporting SEM +standardization when reporting regression results. # References - -Cheung, S. F., Cheung, S.-H., Lau, E. Y. Y., Hui, C. H., & Vong, W. N. (2022). Improving an old way to measure moderation effect in standardized units. *Health Psychology, 41*(7), 502--505. https://doi.org/10.1037/hea0001188 - -Falk, C. F. (2018). Are robust standard errors the best approach for interval estimation with nonnormal data in structural equation modeling? *Structural Equation Modeling: A Multidisciplinary Journal, 25*(2), 244--266. https://doi.org/10.1080/10705511.2017.1367254 - -Friedrich, R. J. (1982). In defense of multiplicative terms in multiple regression equations. *American Journal of Political Science, 26*(4), 797--833. https://doi.org/10.2307/2110973 - -