Introduction to psfm

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.

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)
#> Warning in commonArgs(par, fn, control, environment()): maxfun < 10 *
#> length(par)^2 is not recommended.
summary(p.gtre_sml)
#> --- SFA Regression Model Summary ---
#> Formula: y_gtre ~ x1 + x2 
#> Total time: 7.608213 
#> Model Output:
#>                   par     st_err     t-val
#> lambda      2.8394395 0.63905659  4.443174
#> sigma       0.3045838 0.01828817 16.654688
#> sigr        0.0887480 0.03505840  2.531434
#> sigh        0.3707308 0.03483123 10.643631
#> (Intercept) 0.5060797 0.05436173  9.309485
#> x1          0.5009620 0.01028763 48.695554
#> x2          0.4979087 0.01265119 39.356655
#> log likelihood: 16.61189 
#> convergence: 1 -- ITERATION LIMIT REACHED; this is not a converged optimum
#>   optimizer message: NEW_X
#>   a non-zero code does not by itself mean the fit failed --
#>   run sfa_diagnostics() on this fit to see the gradient and Hessian.
mean(p.gtre_sml$U)
#> [1] 0.8079776
mean(p.gtre_sml$H)
#> [1] 0.7692999

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:

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:

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():

c(sigh_at_bound = p.gtre_sml$sigh_at_bound,
  sigr_at_bound = p.gtre_sml$sigr_at_bound)
#> sigh_at_bound sigr_at_bound 
#>         FALSE         FALSE

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:

Read the two persistent scales together, not separately. Their sum is still informative; the division between them is not identified in that sample.
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.
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.
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.
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} \]

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)
#> [1] 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:

coef(p.gtre_sml)                 # named vector of point estimates
#>      lambda       sigma        sigr        sigh (Intercept)          x1 
#>   2.8394395   0.3045838   0.0887480   0.3707308   0.5060797   0.5009620 
#>          x2 
#>   0.4979087
vcov(p.gtre_sml)                 # variance-covariance matrix (from the Hessian)
#>                    lambda         sigma          sigr          sigh
#> lambda       0.4083933287  8.880296e-03  2.925880e-03 -2.355887e-04
#> sigma        0.0088802959  3.344573e-04  1.995673e-05 -2.589498e-05
#> sigr         0.0029258803  1.995673e-05  1.229091e-03 -3.072758e-04
#> sigh        -0.0002355887 -2.589498e-05 -3.072758e-04  1.213215e-03
#> (Intercept)  0.0119139231  3.056385e-04 -3.421815e-04  7.434362e-04
#> x1          -0.0006520449 -1.444698e-05 -9.895864e-06  1.175415e-05
#> x2          -0.0007249687 -1.292746e-05 -6.790614e-05 -1.816977e-06
#>               (Intercept)            x1            x2
#> lambda       0.0119139231 -6.520449e-04 -7.249687e-04
#> sigma        0.0003056385 -1.444698e-05 -1.292746e-05
#> sigr        -0.0003421815 -9.895864e-06 -6.790614e-05
#> sigh         0.0007434362  1.175415e-05 -1.816977e-06
#> (Intercept)  0.0029551982 -1.200822e-04 -4.701230e-04
#> x1          -0.0001200822  1.058354e-04 -7.193147e-06
#> x2          -0.0004701230 -7.193147e-06  1.600527e-04
logLik(p.gtre_sml)               # log-likelihood, with df/nobs attributes set
#> 'log Lik.' 16.61189 (df=7)
AIC(p.gtre_sml); BIC(p.gtre_sml) # available "for free" once logLik() works
#> [1] -19.22379
#> [1] 9.057993

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).

## 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)
#> Warning in commonArgs(par, fn, control, environment()): maxfun < 10 *
#> length(par)^2 is not recommended.
summary(p.gtre_z)
#> --- SFA Regression Model Summary ---
#> Formula: y_gtre_zz ~ x1 + x2 | z_gtre | zp_gtre 
#> Total time: 5.810348 
#> Model Output:
#>                         par     st_err      t-val
#> sigv           1.076345e-01 0.04108009  2.6201128
#> sigr           1.762764e-01 0.00000000        Inf
#> (Intercept x)  3.772770e-01 0.10568547  3.5698092
#> x1             4.554984e-01 0.04351017 10.4687800
#> x2             5.565651e-01 0.05456775 10.1995245
#> (Intercept u)  1.328111e-01 0.38404781  0.3458192
#> z_gtre         7.836690e-01 0.24436487  3.2069627
#> (Intercept h) -1.162703e-05 0.00000000       -Inf
#> zp_gtre        5.307625e-01 0.30017536  1.7681748
#> log likelihood: -667.6627 
#> convergence: 1 -- ITERATION LIMIT REACHED; this is not a converged optimum
#>   optimizer message: NEW_X
#>   a non-zero code does not by itself mean the fit failed --
#>   run sfa_diagnostics() on this fit to see the gradient and Hessian.

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.

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)
#> Warning in commonArgs(par, fn, control, environment()): maxfun < 10 *
#> length(par)^2 is not recommended.

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)
#>      lambda       sigma        sigr (Intercept)          x1          x2 
#>  0.34463500  0.01846778  0.01160290  0.05821697  0.01346153  0.01111306
boot$model$out             # a copy of p.tre$out with bootstrap SEs/t-values written in
#>                   par     st_err     t-val
#> lambda      3.1035669 0.34463500  9.005374
#> sigma       0.2624955 0.01846778 14.213702
#> sigr        0.1098259 0.01160290  9.465387
#> (Intercept) 0.4779295 0.05821697  8.209453
#> x1          0.5271614 0.01346153 39.160589
#> x2          0.4917344 0.01111306 44.248335

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.

Introduction to psfm

Load Package

Which panel model? An overview of psfm()