--- title: "Penalized Factor Analysis: grid search" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Penalized Factor Analysis: grid search} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{css, echo=FALSE} pre { max-height: 600px; overflow-y: auto; } pre[class] { max-height: 600px; } ``` ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE) ``` ## Introduction **Aim**. This vignette shows how to fit a penalized factor analysis model with the scad and mcp penalties using the routines in the `penfa` package. The employed penalty is aimed at encouraging sparsity in the factor loading matrix. Since the scad and mcp penalties cannot be used with the automatic tuning procedure (see Geminiani et al., 2021 for details), the optimal tuning parameter will be found through grid searches. **Data**. For illustration purposes, we use the cross-cultural data set `ccdata` containing the standardized ratings to 12 items concerning organizational citizenship behavior. Employees from different countries were asked to rate their attitudes towards helping other employees and giving suggestions for improved work conditions. The items are thought to measure two latent factors: **help**, defined by the first seven items (`h1` to `h7`), and **voice**, represented by the last five items (`v1` to `v5`). See `?ccdata` for details. This data set is a standardized version of the one in the [`ccpsyc`](https://github.com/Jo-Karl/ccpsyc/) package, and only considers employees from Lebanon and Taiwan (i.e., `"LEB"`, `"TAIW"`). This vignette is meant as a demo of the capabilities of `penfa`; please refer to Fischer et al. (2019) and Fischer and Karl (2019) for a description and analysis of these data. Let us load and inspect `ccdata`. ```{r dataset, R.options = list(width = 100)} library(penfa) data(ccdata) summary(ccdata) ``` ## Model specification Before fitting the model, we need to write a model syntax describing the relationships between the items and the latent factors. To facilitate its formulation, the rules for the syntax specification broadly follow the ones required by [lavaan](https://CRAN.R-project.org/package=lavaan/). The syntax must be enclosed in single quotes `' '`. ```{r syntax} syntax = 'help =~ h1 + h2 + h3 + h4 + h5 + h6 + h7 + 0*v1 + v2 + v3 + v4 + v5 voice =~ 0*h1 + h2 + h3 + h4 + h5 + h6 + h7 + v1 + v2 + v3 + v4 + v5' ``` The factors `help` and `voice` appear on the left-hand side, whereas the observed variables on the left-hand side. Following the rationale in Geminiani et al. (2021), we only specify the minimum number of identification constraints. We are setting the scales of the factors by fixing their factor variances to 1. This can be done in one of two ways: 1) by adding `'help ~~ 1*help'` and `'voice ~~ 1*voice'` to the syntax above; or 2) by setting the argument `std.lv = TRUE` in the fitting function (see below). To avoid rotational freedom, we fix one loading per factor to zero. Parameters can be easily fixed to user-defined values through the pre-multiplication mechanism. By default, unique variances are automatically added to the model, and the factors are allowed to correlate. These specifications can be modified by altering the syntax (see `?penfa` for details on how to write the model syntax). ## Model fitting The core of the package is given by the `penfa` function, a short form for *PENalized Factor Analysis*, that implements the framework discussed in Geminiani et al. (2021). If users decide for either the `scad` or `mcp` penalties, they need to use `strategy = "fixed"` since the automatic procedure is not feasible. ### Scad We start off with the scad. In the function call, we now specify `pen.shrink = "scad"`, and we provide through the `eta` argument the fixed value of the tuning parameter to be employed during optimization (here, for instance, 0.05). The name given to the starting value (here, the factor loading matrix `"lambda"`) reflects the parameter matrix to be penalized. All of its elements are penalized, which means here that the penalization is applied to all factor loadings (except the ones fixed for identification). The scad penalty relies on an additional shape parameter, which is set by default to 3.7 (Fan and Li 2001). This value can be conveniently modified through the `a.scad` argument. See `?penfaOptions` for additional details on the possible options. ```{r scad} scad.fit <- penfa(## factor model model = syntax, data = ccdata, std.lv = TRUE, ## penalization pen.shrink = "scad", eta = list(shrink = c("lambda" = 0.05), diff = c("none" = 0)), ## fixed tuning strategy = "fixed") ``` #### Grid search In order to find the optimal value of the tuning parameter, a grid search needs to be conducted, and the optimal model is the one with the lowest GBIC (Generalized Bayesian Information Criterion). For demo purposes, we use a grid of 51 values evenly spaced between 0 and 0.15. However, for accurate analyses, please consider finer grids (of at least 200 elements) with an upper bound reasonable for the data at hand. ```{r scad_grid_search} # Grid of values for tuning parameter eta.grid <- seq(0, 0.15, length.out = 51) # Return GBIC from a converged and admissible penfa model with fixed tuning penfa.fixedTun <- function(eta, penalty, ...){ fitted <- penfa(model = syntax, data = ccdata, std.lv = TRUE, pen.shrink = penalty, eta = list(shrink = c("lambda" = eta), diff = c("none" = 0)), strategy = "fixed", verbose = FALSE, ...) if(all(fitted@Vcov$solution) & fitted@Vcov$admissibility) return(BIC(fitted)) } # additional penfaOptions can be passed GBIC.scad <- sapply(eta.grid, penfa.fixedTun, penalty = "scad") ``` The optimal tuning parameter is the one generating the penalized model with the lowest `GBIC`. ```{r scad_plot, fig.height=4, fig.width=6, message=FALSE, warning=FALSE} optimtun.scad <- eta.grid[[which.min(GBIC.scad)]] # To plot GBIC across tuning values # p <- plotly::plot_ly(x = eta.grid, y = GBIC.scad, type = 'scatter', mode = 'lines') # plotly::layout(p, xaxis = list(showline = TRUE), # title = list(text = "GBIC values across tuning parameters")) ``` We can ultimately fit the model with the optimal tuning parameter (here, `r optimtun.scad`). The `summary` method details information on the model characteristics, the optimization and penalization procedures as well as the parameter estimates with associated standard errors and confidence intervals. The *Type* column distinguishes between the *fixed* parameters set to specific values for identification, the *free* parameters that have been estimated through ordinary maximum likelihood, and the penalized (*pen*) parameters. The standard errors here have been computed as the square root of the inverse of the penalized Fisher information matrix (Geminiani et al., 2021). The last columns report 95% confidence intervals (CI) for the model parameters. Standard errors and CI of the penalized parameters shrunken to zero are not displayed. A different significance level can be specified through the `level` argument in the `summary` call. ```{r scad_fit} scad.fit <- penfa(## factor model model = syntax, data = ccdata, std.lv = TRUE, ## penalization pen.shrink = "scad", # optimal tuning eta = list(shrink = c("lambda" = optimtun.scad), diff = c("none" = 0)), strategy = "fixed", verbose = FALSE) summary(scad.fit) ``` The penalty matrix can be inspected and plotted as shown in ["plotting-penalty-matrix"](https://egeminiani.github.io/penfa/articles/articles/plotting-penalty-matrix.html). #### Implied moments The implied moments (here, the covariance matrix) can be found via the `fitted` method. ```{r scad_implied} implied <- fitted(scad.fit) implied ``` ### Factor scores Lastly, the factor scores can be calculated via the `penfaPredict` function. ```{r fscores} fscores <- penfaPredict(scad.fit) head(fscores) ``` ### Mcp We can fit a penalized factor model with the mcp penalty in a way similar to the scad. By default the shape parameter of the mcp is set to 3. This value can be conveniently modified through the `a.mcp` argument. See `?penfaOptions` for additional details on the possible options. ```{r mcp_grid_search, fig.height=4, fig.width=6, message=FALSE, warning=FALSE} GBIC.mcp <- sapply(eta.grid, penfa.fixedTun, penalty = "mcp") optimtun.mcp <- eta.grid[[which.min(GBIC.mcp)]] optimtun.mcp ``` The tuning value equal to `r optimtun.mcp` generated the model with the lowest GBIC. ```{r} mcp.fit <- penfa(## factor model model = syntax, data = ccdata, std.lv = TRUE, ## penalization pen.shrink = "mcp", # optimal tuning eta = list(shrink = c("lambda" = optimtun.mcp), diff = c("none" = 0)), strategy = "fixed", verbose = FALSE) summary(mcp.fit) ``` ## Conclusion Implementing the above approach in the multiple-group case would imply carrying out grid searches in three dimensions to find the optimal tuning parameter vector. This is clearly not advisable, as it would introduce further complications and possibly new computational problems and instabilities. For this reason, we suggest users to rely on the automatic tuning parameter selection procedure, whose applicability is demonstrated in the vignettes for single (`vignette("automatic-tuning-selection")`) and multiple-group (["multiple-group-analysis"](https://egeminiani.github.io/penfa/articles/articles/multiple-group-analysis.html)) penalized models. ## R Session ```{r} sessionInfo() ``` ## References * Fan, J., & Li, R. 2001. "Variable Selection via Nonconcave Penalized Likelihood and Its Oracle Properties." Journal of the American Statistical Association 96(456), 1348–60. [https://doi.org/10.1198/016214501753382273](https://doi.org/10.1198/016214501753382273) * Fischer, R., Ferreira, M. C., Van Meurs, N. et al. (2019). "Does Organizational Formalization Facilitate Voice and Helping Organizational Citizenship Behaviors? It Depends on (National) Uncertainty Norms." Journal of International Business Studies, 50(1), 125-134. [https://doi.org/10.1057/s41267-017-0132-6](https://doi.org/10.1057/s41267-017-0132-6) * Fischer, R., & Karl, J. A. (2019). "A Primer to (Cross-Cultural) Multi-Group Invariance Testing Possibilities in R." Frontiers in psychology, 10, 1507. [https://doi.org/10.3389/fpsyg.2019.01507](https://doi.org/10.3389/fpsyg.2019.01507) * Geminiani, E. (2020). "A Penalized Likelihood-Based Framework for Single and Multiple-Group Factor Analysis Models." PhD thesis, University of Bologna. [http://amsdottorato.unibo.it/9355/](http://amsdottorato.unibo.it/9355/) * Geminiani, E., Marra, G., & Moustaki, I. (2021). "Single- and Multiple-Group Penalized Factor Analysis: A Trust-Region Algorithm Approach with Integrated Automatic Multiple Tuning Parameter Selection." Psychometrika, 86(1), 65-95. [https://doi.org/10.1007/s11336-021-09751-8](https://doi.org/10.1007/s11336-021-09751-8)