--- title: "Introduction to psfm" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction to psfm} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Load Package First let us load the sfa package. ```{r setup1} library(sfa) ``` ## Which panel model? An overview of `psfm()` `psfm()` is the panel-data entry point of the sfa package (`sfm()`, `zsfm()`, and `ttsfm()` are its cross-sectional counterparts -- see their own help pages, `?sfm`/`?zsfm`/`?ttsfm`). It fits one of twenty-one `model_name` choices, selected via the `model_name` argument. The most commonly used are listed here; `?psfm` has the full list, including the classical panel estimators (`"CSS"`, `"LS"`, `"KSS"`, `"SSRE"`, `"SSCRE"`) added since this vignette was written: | `model_name` | Individual effect | Estimator | Notes | |---|---|---|---| | `"TRE"` | random | simulated ML (Halton draws) | homoskedastic $\sigma_u$ | | `"TRE_Z"` | random | simulated ML | $\sigma_u$ a function of covariates (1 pipe) | | `"GTRE"` | random, persistent + transient | simulated ML | homoskedastic $\sigma_u$, $\sigma_h$ | | `"GTRE_Z"` | random, persistent + transient | simulated ML | $\sigma_u$, $\sigma_h$ both covariate-driven (2 pipes) | | `"TFE"` | fixed | ML (Chen, Schmidt & Wang 2014) | no separate noise-vs-effect random draw for $r_i$ | | `"FD"` | random, persistent | first-differenced GLS (Wang & Ho 2010) | $\sigma_u$ time-invariant, scaled by covariates (1 pipe) | | `"GTRE_SEQ1"`, `"GTRE_SEQ2"` | random | sequential moment-based | fast, but not maximum likelihood | | `"SSFE"` | fixed | LSDV (Schmidt & Sickles 1984) | deterministic, no optimizer | | `"PL80"` | random, persistent | ML, closed form | Pitt & Lee (1980) | | `"BC92"` | random, time-varying | ML, closed form | Battese & Coelli (1992) | | `"CSS"` | fixed, time-varying | firm-specific quadratic | Cornwell, Schmidt & Sickles (1990) | | `"LS"` | fixed, common time path | rank-one factor, ALS | Lee & Schmidt (1993) | | `"KSS"` | fixed, factor structure | factor model, rank chosen by IC | Kneip, Sickles & Song (2012) | | `"SSRE"`, `"SSCRE"` | random / correlated random | GLS, Mundlak means for `SSCRE` | Schmidt & Sickles (1984) | The last four assume no distribution for inefficiency at all -- they identify it from the firm effects rather than from a one-sided likelihood -- so they carry no `sigma_u` and report slopes only. They are the right comparison when you want to know how much of your efficiency ranking is coming from the half-normal assumption. Every model shares the same `formula`/`data`/`individual` calling convention, so switching between them is usually just a change of `model_name` (and, for the `_Z`/`FD` models, adding the relevant formula pipe -- see "Heteroskedastic inefficiency" below). This vignette walks through the most commonly used model, `GTRE`, in full, then shows the heteroskedastic (`GTRE_Z`) variant and the parametric bootstrap. See `?psfm` and each model's own documentation (e.g. `?data_gen_p`, which lists every DGP column and which model it targets) for the rest. ## Generalized True Random Effects Model We are interested in estimating the Generalized True Random Effects Model (GTRE) model of Filippini and Greene (2016) using simulated maximum likelihood. We will begin by using a simulated data set. We use the `data_gen_p` call to create a simulated data set. The `psfm` call runs the likelihood through three successive optimizer routines, only updating the parameter values if there is an improvement in the likelihood. We have opted for 150 iterations of the "bobyqa" procedure, 10 for "psoptim", and 10 for "optim". `rand` sets the seed (via `set.seed()`) for replication. Other arguments set the true $\sigma$'s and $\beta$'s used to generate the data, with `cons` for $\beta_0$/the constant. `rand.gtre` and `rand.psoptim` seed the simulated-ML Halton draws and the particle-swarm stage respectively; without them a model fit with `PSopt = TRUE` is not reproducible from one session to the next, which is why every fit below sets them. Every fit in this vignette is deliberately sized to build quickly rather than to be a serious estimation exercise: a 70-firm, 6-period panel, only 50 Halton draws (`halton_num`), and a handful of optimizer iterations. Real work should use a genuine panel and the defaults -- `halton_num` alone defaults to `ceiling(sqrt(nrow(data))) + 100`, several times what is used here. Expect the numbers below to move if you rerun them at a realistic size. ### Which estimator: `"fiml"` or `"sml"` `model_name = "GTRE"` names the *model*; `estimator` chooses how it is fitted. Since 1.1.4 the default is `"fiml"`, full-information maximum likelihood through the model's closed-skew-normal representation, which integrates the persistent block exactly and needs no draws at all. The fit below pins `estimator = "sml"` -- simulated maximum likelihood over Halton draws -- so that the `halton_num`/`rand.gtre` arguments have something to do. Which to use: - **`"fiml"`** is the default, is deterministic, and is several times faster. It requires a **balanced** panel; on an unbalanced one `psfm()` warns and falls back to `"sml"`. - **`"sml"`** handles unbalanced panels, and on the evidence of the package's own benchmarks returns better technical-efficiency *scores* (see `horserace/` in the development repository). It is the slower route. One consequence worth knowing: `psfm(model_name = "GTRE")` returns an object whose `$model_name` is `"GTRE_FML"` under the default. That is the name to use when you look it up in tables like the one above. ```{r setup2} data_trial <- data_gen_p(t=6, N=70, rand = 16, sig_u = 0.3, sig_v = 0.1, sig_r = 0.1, sig_h = 0.3, cons = 0.5, beta1 = 0.5, beta2 = 0.5) p.gtre_sml <- psfm(formula = y_gtre ~ x1 + x2, model_name = "GTRE", estimator = "sml", data = data_trial, individual = "name", PSopt = TRUE, optHessian = TRUE, halton_num = 50, rand.gtre = 1, rand.psoptim = 1, maxit.bobyqa = 150, maxit.psoptim= 10, maxit.optim = 10) summary(p.gtre_sml) mean(p.gtre_sml$U) mean(p.gtre_sml$H) ``` ## GTRE Results The results give the model parameter estimates in the classic $\lambda$-$\sigma$ framework as well as the mean efficiency scores. We see that most parameters are estimated well: $\hat\lambda \approx 2.7$ against a true $\lambda = \sigma_u/\sigma_v = 0.3/0.1 = 3$. The optimizer struggles more with $\sigma_r$ specifically -- its estimate (about 0.06) sits below the true value of 0.1, with a $t$-value near 2.2 against 40 and 49 for $\beta_1$ and $\beta_2$. The persistent random effect is the component a short panel speaks least clearly about, and on a panel this small it is the estimate to treat with the most caution: at other seeds it can collapse toward zero altogether. $\sigma_h$ comes out somewhat above its true value (0.40 against 0.3), while $\beta_0$, $\beta_1$, and $\beta_2$ are all close to their true values with large $t$-values. Increasing the number of optimizer iterations (via `maxit.bobyqa`/`maxit.psoptim`/`maxit.optim`) typically improves accuracy at the cost of run time; for `GTRE`/`GTRE_Z`/`TRE`/`TRE_Z`, which use simulated ML over Halton draws, each objective-function evaluation is meaningfully more expensive than in the closed-form cross-sectional models, so this trade-off is worth tuning to your own patience/hardware rather than increasing iterations blindly. We also see that the mean transient technical efficiency (`U`) is around 0.81 and the mean persistent technical efficiency (`H`) is around 0.76. We may be interested in plotting the densities of these efficiency scores: ```{r setup3} plot(density(p.gtre_sml$U),main="Density of Transient TE") plot(density(p.gtre_sml$H),main="Density of Persistent TE") ``` To get the total technical efficiency, we would simply multiply the TE's of U and H in the following way: ```{r setup4} total_te <- rep(p.gtre_sml$H, each=6) * p.gtre_sml$U plot(density(total_te),main="Density of Total TE") ``` NB: Using the constant 6 in rep(p.gtre_sml$H, each=6) only works for a balanced panel with t=6 for each individual. For an unbalanced panel, the each argument in rep would not work. Instead, use the times argument and a vector of length N (number of individuals), with each of the time period lengths, e.g. rep(c(2,5),times=c(5,2)). ## When one of the persistent scales comes back as zero `GTRE` splits the time-invariant part of the error into two pieces: a symmetric firm random effect $r_i$ and a one-sided persistent inefficiency $h_i$. Only the *asymmetry* of $h$ distinguishes them, and in a fair share of samples the likelihood cannot do it and merges the two instead -- one scale goes to zero and the other absorbs its variation. `psfm()` reports this in `$sigh_at_bound` and `$sigr_at_bound`, warns at fit time, and repeats it in `print()`/`summary()`: ```{r setup_bnd} c(sigh_at_bound = p.gtre_sml$sigh_at_bound, sigr_at_bound = p.gtre_sml$sigr_at_bound) ``` **This is usually the correct answer, not a failed fit.** On one simulated replication the boundary solution had a log-likelihood *higher than the true parameter vector* by 3.66 -- a likelihood ratio of about 39. It is the panel counterpart of cross-sectional wrong skewness (Waldman 1982), where $\sigma_u = 0$ is the maximum likelihood estimate rather than an optimizer failure. So `psfm()` reports it instead of preventing it: bounding either scale away from zero would corrupt exactly the samples where the boundary is the answer. It is also common. On simulated panels with both components genuinely present, one of the two collapses in something like a quarter to a third of samples at moderate $N$ and $T$, falling as the panel widens. What to do when it happens: 1. **Read the two persistent scales together, not separately.** Their *sum* is still informative; the division between them is not identified in that sample. 2. **Treat the intercept with the same caution.** The frontier intercept and the mean of persistent inefficiency are the same nuisance, since $E[y] = \beta_0 - E[h] - E[u]$, so a collapse displaces $\beta_0$ by about $E[h] = \sigma_h\sqrt{2/\pi}$. Any *level* read off `coef()` moves with it. Slopes are unaffected. 3. **Do not bootstrap that fit.** `psfm_bootstrap()` resamples from the *fitted* model, so it would draw from a world in which the collapsed component is absent, and the bootstrap is inconsistent on the boundary of the parameter space regardless. It warns rather than proceeding quietly. 4. **Do not simply refit with more draws until the zero goes away.** It is tempting, and it is a selection rule: retrying in one direction only, and stopping when you like the answer, biases the reported estimate away from a boundary that is often correct. If you want to know whether a particular collapse is simulation noise, fit the same data with `estimator = "fiml"`, which uses no draws at all -- if that also lands on the boundary, the boundary is the maximum. 5. **What actually helps is a wider panel**, and how much more firms help than periods depends on the estimator; `?psfm` gives the measured elasticities for both. You can tell in advance whether a panel can support the split at all. The quantity that governs it is *not* $\sigma_h/\sigma_r$ -- raising that only decides which of the two scales dies. It is the persistent variance in the firm means relative to the transient variance still contaminating them: $$ S/N \;=\; \frac{\sigma_r^2 + \sigma_h^2(1-2/\pi)}{\left[\sigma_u^2(1-2/\pi) + \sigma_v^2\right]/T} $$ ```{r setup_sn} sn_gtre <- function(sig_u, sig_v, sig_r, sig_h, T) { k <- 1 - 2/pi (sig_r^2 + sig_h^2 * k) / ((sig_u^2 * k + sig_v^2) / T) } # the fit above, read off its own estimates sn_gtre(sig_u = 0.3, sig_v = 0.1, sig_r = 0.1, sig_h = 0.3, T = 6) ``` Across six simulated designs varying all four scales, this ordered the collapse rate perfectly (Spearman $-1$, both estimators): 54% of fits merged at $S/N = 0.67$, 31% at $1.30$, and 6% at $11.7$. **Below about 1, expect roughly half your fits to merge the two components; it takes something like $S/N > 5$ to get under 10%.** Note that $T$ is in the denominator, so a longer panel raises it directly. A quick check on your own data, before fitting anything: average the OLS residuals within each firm and look at their skew. A production frontier implies it should be negative. When it comes out positive the sample carries no between-firm evidence of persistent inefficiency, and that single statistic predicts which samples lose $\sigma_h$ with an ROC area of 0.955. ## Extracting results: standard R modeling generics `psfm()` (like `sfm()`, `zsfm()`, and `ttsfm()`) returns an object of class `"sfareg"`, which supports the usual R modeling generics rather than requiring you to reach into `$out` by hand: ```{r setup5} coef(p.gtre_sml) # named vector of point estimates vcov(p.gtre_sml) # variance-covariance matrix (from the Hessian) logLik(p.gtre_sml) # log-likelihood, with df/nobs attributes set AIC(p.gtre_sml); BIC(p.gtre_sml) # available "for free" once logLik() works ``` `logLik()` (and therefore `AIC()`/`BIC()`) returns `NA` with a warning for the moment-based/LSDV models (`GTRE_SEQ1`, `GTRE_SEQ2`, `SSFE`), since those aren't fit by maximizing a likelihood in the first place. ## Heteroskedastic inefficiency: the pipe formula and `GTRE_Z` The models above assume homoskedastic inefficiency -- a single $\sigma_u$ (and, for GTRE, $\sigma_h$) shared by every observation. `psfm()`'s `_Z` models instead let $\sigma_u$ (and, for `GTRE_Z`, $\sigma_h$) depend on covariates, via a formula with one or two extra parts separated by `|`: `y ~ x1 + x2 | z_u | z_h`. The first part is the usual frontier equation; the second parameterizes $\sigma_u$; the third (only for `GTRE_Z`) parameterizes $\sigma_h$. `GTRE` and `TRE` (no `_Z` suffix) do not accept any pipes -- if you want covariate-driven inefficiency, use `GTRE_Z`/`TRE_Z` and name them explicitly (an older version of this package let a pipe on `model_name = "GTRE"` silently upgrade to `GTRE_Z`; this package now requires the `_Z` name to be written explicitly instead, so it's always clear from the call itself which model was fit). ```{r setup6} ## data_trial already holds every column data_gen_p() produces, including the ## y_gtre_zz/z_gtre/zp_gtre trio this model needs -- no need to simulate again. p.gtre_z <- psfm(formula = y_gtre_zz ~ x1 + x2 | z_gtre | zp_gtre, model_name = "GTRE_Z", data = data_trial, individual = "name", PSopt = TRUE, optHessian = TRUE, halton_num = 50, rand.gtre = 1, rand.psoptim = 1, maxit.bobyqa = 150, maxit.psoptim= 10, maxit.optim = 10) summary(p.gtre_z) ``` Note the coefficient layout: `sigv`, `sigr`, the frontier ($x$) block, then the $\sigma_u$ block (its own intercept plus `z_gtre`), then the $\sigma_h$ block (its own intercept plus `zp_gtre`) -- each pipe segment gets its own intercept, which is why `(Intercept x)`/`(Intercept u)`/`(Intercept h)` are labeled separately rather than sharing one `(Intercept)` row. This `GTRE_Z`/`TRE_Z` link function is $\sigma = \sqrt{\exp(z'\delta)}$ ($\delta$ parameterizes the *variance*, not the standard deviation) -- a genuine, existing difference from `sfm()`'s `NHN_Z`/`NE_Z`, which use $\sigma = \exp(z'\delta)$ directly. Check which convention applies before interpreting a fitted $z$-coefficient. ## Parametric bootstrap for inference: `psfm_bootstrap()` The Hessian-based standard errors above rely on standard asymptotic MLE theory. `psfm_bootstrap()` offers a parametric-bootstrap alternative: starting from a fitted model, it repeatedly (a) simulates a new response from the fitted parameters (using each model's own assumed data-generating process), (b) re-estimates the same model on the simulated data, and (c) uses the spread of those re-estimates as the standard error. It supports `GTRE_Z`, `TRE_Z`, `GTRE`, `GTRE_FML`, `TRE`, `TFE`, `TFE_WMLE` and `FD`. `PL80` and `BC92` are excluded because they do not expose the `$U`/`$H` efficiency structure the function reads, and the moment-based and LSDV estimators (`GTRE_SEQ1`, `GTRE_SEQ2`, `SSFE`) because they are not maximum likelihood. Note that `psfm(model_name = "GTRE")` returns an object whose `$model_name` is `"GTRE_FML"` under the default `estimator = "fiml"`, which is why both names appear in that list. ```{r setup7} data_trial_tre <- data_gen_p(t=5, N=30, rand=16, sig_u=0.3, sig_v=0.1, sig_r=0.1, sig_h=0.3, cons=0.5, beta1=0.5, beta2=0.5) p.tre <- psfm(formula = y_tre ~ x1 + x2, model_name = "TRE", data = data_trial_tre, individual = "name", halton_num = 50, rand.gtre = 1, maxit.bobyqa = 300) set.seed(1) boot <- psfm_bootstrap(p.tre, numCores = 2, BOOT = 5, # a real analysis should use far more, e.g. 199-999 individual = "name", inefdec = TRUE, maxit.bobyqa = 150, maxit.psoptim = 30) boot$se # bootstrap standard errors, one per parameter in coef(p.tre) boot$model$out # a copy of p.tre$out with bootstrap SEs/t-values written in ``` `BOOT = 5` above is only for a fast-building vignette -- a real analysis should use at least 199, and ideally more, replications. Five is far too few to read the standard errors below as anything but a demonstration of the calling convention. Each replication refits the full model, so runtime scales roughly linearly with `BOOT` (parallelized across `numCores`); note that the refits use the model's default Halton draw count rather than the reduced `halton_num` of the original fit. See `?psfm_bootstrap` for the full per-model data-generating process assumptions, and for `boot_eff`/ `boot_eff_h` (bootstrap draws of the technical-efficiency scores themselves, not just the parameters). ## Practical tips - **`inefdec`**: `TRUE` (default) fits a production-frontier model (inefficiency lowers output, $y = x'\beta + v - u$); `FALSE` fits a cost-frontier model (inefficiency raises cost, $y = x'\beta + v + u$). Make sure this matches how your outcome variable is defined. - **`PSopt`**: enables the middle (`psoptim`, particle-swarm) optimization stage. It can help escape a poor starting point, but its wide, undirected search can occasionally land on a degenerate solution (e.g. a variance component collapsing toward its numerical floor) that the default two-stage (`bobyqa` then `optim`) path would not have found. If a fit's standard errors look implausible (e.g. far larger than the point estimate, or a variance parameter near zero), try refitting with `PSopt = FALSE` before concluding the model doesn't fit the data. - **Pipe counts are enforced per model**: `psfm()` checks the formula's pipe count against what the requested `model_name` actually supports (see the table above) and stops with a clear error if they don't match, rather than silently ignoring extra pipes or extra covariates. - **`individual`**: the panel ID column name (a character string, e.g. `"name"` above), required by every model in this vignette. - **Robust and clustered standard errors** come from the `sandwich` package once the fit is made with `keep_objective = TRUE`: `sandwich::sandwich(fit)` for Huber--White, `sandwich::vcovCL(fit, cluster = d$firm)` for clustered, either of which `lmtest::coeftest()` accepts. `sandwich::vcovHC()` does *not* work and is not meant to -- its finite-sample corrections are built from hat values, which are not defined for a nonlinear likelihood whose parameters include variance components. - **Other things worth knowing about**, each with its own help page: `marginal_effects()` for how a covariate moves inefficiency in the `_Z` models, `efficiency_ci()` for Horrace--Schmidt intervals around the efficiency scores rather than point predictions, and `skewness_test()` for a p-value on cross-sectional wrong skewness. That last one applies to `sfm()` fits; it is *not* valid on between-firm residuals, because the firm mean carries $-\bar u_i$ as well as $-h_i$ and the mean of $T$ half-normals is skewed the same way.