Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
44,155 characters · 17 sections · 18 citation commands
Estimating Semiparametric and Nonparametric Fixed Effects Panel Data Models with mgcv
Keywords: semiparametric regression, nonparametric regression, panel data, fixed effects, generalized additive models, splines
JEL codes: C14, C15, C23, C87
Empirical work in economics frequently uses panel data to account for unobserved unit-level heterogeneity through fixed effects. At the same time, many applied questions call for flexible functional forms. For example, borri_martins_filho_kalatzis_2022 estimate a nonlinear investment curve, while baglan_yoldas_2014 estimate an inflation--growth relationship semiparametrically. Such empirical applications remain relatively rare, however, despite a well-developed theoretical literature on semiparametric and nonparametric panel data models.
A central practical obstacle is the need to handle fixed effects while also choosing basis dimension, series length, bandwidths, or smoothing penalties. Nonparametric methods that are popular in cross-sectional settings, such as kernel or nearest-neighbor estimators, do not naturally accommodate additive fixed effects. Series-based approaches can incorporate fixed effects more directly, but they require the researcher to specify basis functions, choose interaction terms, and select the number of series terms.
This paper provides a practical guide to estimating semiparametric and nonparametric fixed effects panel data models using penalized splines. The approach connects three literatures. The first is the econometric literature on series estimation for semiparametric panel data models with fixed effects baltagi_li_2002,an_et_al_2016. The second is the broader econometric literature on sieve and penalized sieve estimation and inference for nonparametric and semiparametric models ai_chen_2003,chen_2007,chen_pouzo_2009,chen_pouzo_2012,chen_pouzo_2015. The third is the statistical literature on generalized additive models, smoothing parameter selection, and penalized regression splines wood_2003,wood_2011,wood_2017_gam.
Conceptually, the proposed approach is grounded in the econometric theory of series, sieve, and penalized sieve estimation. Computationally, it uses the fast and stable implementations developed in the generalized additive model literature, making flexible estimation practical even in panel settings with high-dimensional fixed effects.
The paper gives concrete implementation guidance in R using the mgcv package wood_mgcv_package. I first describe the class of fully nonparametric, partially linear, and varying-coefficient panel models covered by the approach. I then show how the resulting penalized-spline estimators can be implemented using unit indicators or first differences, how a penalized unit-effect specification can be used as a regularized alternative, and how to construct penalty-adjusted cluster-robust standard errors and confidence bands from the fitted bam object. Monte Carlo simulations evaluate estimation accuracy for the unknown function, empirical size and power for finite-dimensional parameters, and coverage of confidence bands for the nonparametric component.
The proposed approach can handle a variety of fixed effects panel data models. One example is the fully nonparametric model
where \(X_{it}\) may be scalar or vector-valued, \(g(\cdot)\) is unknown, and \(\mu_i\) is a unit fixed effect.
Another example is the partially linear model
where \(X_{it} = (X_{1,it}', X_{2,it}')'\), \(\beta\) is an unknown finite-dimensional parameter, and \(g(\cdot)\) is an unknown function.
A third example is a varying-coefficient model,
where \(\beta_0(\cdot)\) is an unknown intercept function and \(\beta_1(\cdot)\) is an unknown smooth coefficient function, or vector of coefficient functions when \(X_{1,it}\) is vector-valued.
In these additive fixed-effects models, the level of the unknown smooth is not separately identified from the intercept and unit effects. The identified object is a normalized version of the function. For population statements one can write \(g^c(x)=g(x)-E[g(X_{2,it})]\), with the obvious modification in the fully nonparametric model. The tables below use the operational grid analog \[ g_G^c(x_j)=g(x_j)-\frac{1}{J}\sum_{\ell=1}^J g(x_\ell), \qquad j=1,\ldots,J, \] on the evaluation grid \(G=\{x_1,\ldots,x_J\}\). For RMSE calculations the same normalization is applied on the independent test set. For first-differenced estimators, the data identify smooth contrasts \(g(X_{2,it})-g(X_{2,i,t-1})\); reported level curves are reconstructed in the fitted basis and then normalized using the same grid-centering convention.
This section describes the penalized objective that underlies the estimators used in this paper and the role of fast restricted maximum likelihood (fREML) in selecting the penalty weights. A Gaussian identity-link gam or bam fit can be viewed as a penalized least squares problem with the amount of penalization chosen automatically. After each smooth term has been represented by a finite set of basis functions, the model is linear in its unknown coefficients. Let $R$ denote the resulting mgcv model matrix, including parametric regressors, spline basis columns, and any factor indicators, and let $\theta$ collect all coefficients. For a fixed vector of smoothing parameters $\lambda=(\lambda_1,\ldots,\lambda_J)'$, mgcv estimates
The matrices $S_j$ are known positive semidefinite penalty matrices constructed from the chosen smooth bases; each $S_j$ has nonzero entries only for the coefficients belonging to the corresponding smooth term. This is directly analogous to the familiar ridge-type objective \[ \min_{\beta_0,\ldots,\beta_k} \sum_i \left(Y_i - \beta_0 - \sum_{j=1}^{k} X_{ij}\beta_j\right)^2 + \lambda \sum_{j=1}^{k}\beta_j^2. \] Here the intercept is left unpenalized, while larger values of \(\lambda\) shrink the slope coefficients toward zero. In the spline case, however, the row \(R_i\) contains the usual regressors, spline basis functions, fixed-effect indicators, and any other columns generated by the model formula, rather than only the raw covariates \(X_i\).
For common cubic and thin-plate spline bases, the quadratic form $\theta'S_j\theta$ is a finite-dimensional measure of wiggliness, such as an integrated squared derivative. Thus \(\lambda_j\) controls how costly wiggliness is for the \(j\)th smooth. A large $\lambda_j$ pushes the smooth toward the null space of its penalty, such as a low-order polynomial component; with select = TRUE, mgcv adds an additional shrinkage penalty on this null space, so that a smooth can be penalized all the way toward zero. Ordinary parametric regressors and factor fixed effects are not penalized unless an explicit parametric penalty is supplied. By contrast, a term such as s(unit, bs = "re") uses a ridge-type penalty on the unit effects, corresponding to a Gaussian random-effect representation.
The label bs = "re" refers to a computational random-effect representation of a ridge penalty on the unit indicators, not to the classical random-effects assumption that the unit effect is uncorrelated with the regressors. Conditional on the smoothing parameter, the estimator is penalized least squares with a quadratic penalty on the unit effects. Because this specification shrinks the unit effects toward zero, however, the resulting estimator differs numerically from including unpenalized factor fixed effects. When the unit effects are correlated with the regressors and \(T\) is fixed, this shrinkage can transmit part of that correlation into the estimated smooth or finite-dimensional coefficient. I therefore treat s(unit, bs = "re") as a predictive or regularized alternative, not as a substitute for the fixed-effects estimand delivered by unpenalized unit indicators or first differencing.
It is important to distinguish the coefficient objective in equation (ref) from the role of method = "REML" or method = "fREML". Conditional on \(\lambda\), the coefficient estimates solve the penalized least squares problem above; REML or fREML is then used to choose \(\lambda\) and the scale parameter. In the Gaussian case, this can be understood by treating the penalized spline coefficients as Gaussian random effects with precision proportional to the smoothing parameters. For a trial value of \(\lambda\), mgcv solves the penalized coefficient problem and evaluates a restricted marginal likelihood criterion, with the scale parameter profiled or updated. The outer smoothing-parameter step then searches over \(\lambda\). Schematically, up to constants and rank adjustments, the negative REML criterion has the form \[ \mathcal{V}_{\mathrm{REML}}(\lambda,\sigma^2) = \frac{1}{2\sigma^2}Q_{\lambda}(\hat{\theta}_{\lambda}) +\frac{1}{2}\log |R'R+S_{\lambda}| -\frac{1}{2}\log |S_{\lambda}|_{+} +\frac{N-d_0}{2}\log \sigma^2, \] where $N=nT$, $|\cdot|_+$ denotes the product of positive eigenvalues, and $d_0$ is the dimension of the unpenalized component after imposing the relevant identifiability constraints. The data-fit contribution enters through \(Q_{\lambda}(\hat{\theta}_{\lambda})\) and the scale term. The determinant terms adjust the marginal likelihood for the design, penalty, and dimension of the penalized and unpenalized coefficient spaces, thereby discouraging overly flexible smooths when the smoothing parameters are chosen. The bam option fREML is a computationally efficient REML-type implementation for large data sets wood_2011,wood_goude_shaw_2015,wood_li_shaddick_augustin_2017. Thus, for the estimators used below, the coefficient objective is “RSS plus a quadratic wiggliness penalty,” while REML/fREML selects the penalty weights rather than requiring the researcher to choose them by hand.
While REML/fREML is the default choice in mgcv::bam, one could use a more familiar generalized cross-validation (GCV) criterion. REML/fREML does not require a fully parametric likelihood model for the substantive econometric object of interest; instead, it is used as a data-driven rule for selecting the smoothing parameters. Other smoothing selection criteria, including GCV/Cp-type criteria, are also available in mgcv. I use fREML in the simulations because REML-type criteria are often more stable in finite samples and less prone to undersmoothing than GCV-type criteria reiss_ogden_2009, but the same implementation framework can be used with GCV if desired.
Standard mgcv returns a Bayesian or frequentist covariance matrix that assumes independent, homoskedastic observations. For panel data with within-unit dependence, a cluster-robust analogue is needed. This section develops a penalty-adjusted cluster-robust covariance estimator from the fitted bam object.
Let $R$ denote the fitted mgcv linear predictor matrix, $S_{\lambda}$ the penalty matrix, and $\hat{u}$ the residual vector. For Gaussian identity-link models without observation weights, the mgcv working-weight matrix is the identity matrix, so the penalized Hessian matrix is
A penalty-adjusted clustered sandwich then has the form
where $R_{it}$ is the row of $R$ for observation $(i,t)$. $\hat{\Omega}$ is the unit-cluster robust covariance component: it first sums the score contributions $R_{it}'\hat{u}_{it}$ over time within each unit and then forms outer products of the unit-level sums. This has the same unit-level outer-product form as the robust covariance estimator for within-group estimators in arellano_1987, and therefore allows arbitrary heteroskedasticity and serial dependence within units under independence across units. In the present setting, this expression can be viewed as a penalized-spline analog of the score covariance term in series-based semiparametric panel estimators. Setting $S_{\lambda}=0$ gives the corresponding series-style cluster sandwich that ignores the penalty.
For a scalar finite-dimensional parameter \(\beta\), or for one component of a vector \(\beta\), let \(\hat{\sigma}_{\beta}^2\) denote the corresponding diagonal element of \(\widehat{V}\). Under the Gaussian approximation, \[ \frac{\hat{\beta}-\beta}{\hat{\sigma}_{\beta}} \approx N(0,1). \] Thus an approximate \(1-\alpha\) confidence interval for \(\beta\) is \[ \left[ \hat{\beta}-z_{1-\alpha/2}\hat{\sigma}_{\beta},\, \hat{\beta}+z_{1-\alpha/2}\hat{\sigma}_{\beta} \right]. \]
For inference on the nonparametric component, let \(G=\{x_1,\ldots,x_J\}\) be the evaluation grid. Let \(R_G\) denote the grid lpmatrix, and let \(R_G^{(g)}\) be the same matrix after all columns not belonging to the smooth of interest have been set to zero. Because the level of \(g\) is not separately identified from the intercept and unit effects, the target is the grid-centered smooth. Let \[ C_G=I_J-\mathbbm{1}\mathbbm{1}'/J \] be the grid-centering matrix, where \(\mathbbm{1}\) is the \(J\)-vector of ones. Premultiplication by \(C_G\) subtracts the average over the \(J\) grid points. The matrix that maps fitted coefficients into centered smooth values is \[ L_G=C_G R_G^{(g)}. \] Thus \[ \hat{g}_G^c=L_G\hat{\theta}, \qquad \hat{g}_G^c(x_j) = \hat{g}(x_j)-\frac{1}{J}\sum_{\ell=1}^J \hat{g}(x_\ell), \] and the estimated covariance matrix for the centered grid vector is \[ \widehat{\Sigma}_G=L_G\widehat{V}L_G'. \] Pointwise and uniform bands use the same estimated grid covariance matrix but answer different questions. A pointwise \(1-\alpha\) interval at grid point \(x_j\) is \[ \hat{g}_G^c(x_j)\pm z_{1-\alpha/2}\hat{\sigma}_j, \qquad \hat{\sigma}_j^2=(\widehat{\Sigma}_G)_{jj}. \] This interval is calibrated separately at each \(x_j\). For each individual point, \[ P\left( g_G^c(x_j) \in \left[ \hat{g}_G^c(x_j)-z_{1-\alpha/2}\hat{\sigma}_j,\, \hat{g}_G^c(x_j)+z_{1-\alpha/2}\hat{\sigma}_j \right] \right) \approx 1-\alpha. \] This pointwise statement does not imply that all intervals cover the curve simultaneously. The distinction between pointwise component intervals and whole-function coverage is standard in spline and GAM inference; see, for example, marra_wood_2012 on coverage of GAM component intervals.
A uniform band instead controls the probability that the entire reported grid curve is covered simultaneously: \[ P\left( g_G^c(x_j)\in \left[\hat{g}_G^c(x_j)-c_{1-\alpha}\hat{\sigma}_j,\, \hat{g}_G^c(x_j)+c_{1-\alpha}\hat{\sigma}_j\right] \text{ for all } j=1,\ldots,J \right)\approx 1-\alpha. \] The critical value \(c_{1-\alpha}\) is the \(1-\alpha\) quantile of the maximum absolute standardized error over the grid. Each standardized component is marginally standard normal, but the maximum depends on their joint correlation structure. For a spline fit, nearby grid points are typically strongly correlated, so the maximum is neither the maximum of independent standard normals nor a single standard normal. Except in special cases, there is no simple scalar closed-form critical value; simultaneous bands for penalized splines are therefore commonly constructed from the approximating Gaussian distribution krivobokova_kneib_claeskens_2010. In the simulations I approximate this critical value from an auxiliary Gaussian vector \(Z_G\sim N(0,\widehat{\Sigma}_G)\), where \(Z_G\) is simulated rather than observed in the data: \[ c_{1-\alpha} = \operatorname{quantile}_{1-\alpha} \left( \max_{1\leq j\leq J} \left|Z_{G,j}/\hat{\sigma}_j\right| \right). \] This is a simultaneous band on the reported grid. Coverage is evaluated over the 50-point grid used in the simulation.
This section translates the preceding estimation and inference formulas into mgcv code. Suppose that one wants to estimate the model from equation (ref), with scalar $X_1$ and $X_2$. There are several ways to implement the estimator, and the same post-estimation covariance calculations can then be applied to the fitted object.
The main implementation uses factor fixed effects:
An alternative regularized specification uses a penalized unit effect:
This specification is often useful for prediction or as a shrinkage comparison, but it changes the estimator when \(\mu_i\) is correlated with the regressors. The fixed-effects interpretation in the usual panel-data sense comes from the factor-indicator or first-differenced specifications.
Fixed effects can also be removed by first differencing. A constrained smooth of the form \[ g(X_{it}) - g(X_{i,t-1}) \] can be represented in mgcv using linear functional smooths. The implementation constructs the matrix of current and lagged covariates and the corresponding $+1/-1$ weights.
The core first-differenced bam call used in the simulations is:
Here dy and dx1 are first differences, while X_pair stores $(X_{2,it},X_{2,i,t-1})$ and L_pair stores the corresponding $(1,-1)$ weights. The term s(X_pair, by = L_pair) therefore estimates a single smooth $g(\cdot)$ evaluated at both arguments and constrained to enter as $g(X_{2,it})-g(X_{2,i,t-1})$.
By default, s() smooths in mgcv use thin plate regression spline bases, specified by bs = "tp" wood_2003. Researchers can also choose other spline bases, such as cubic regression splines with bs = "cr" or P-splines with bs = "ps".
The same syntax extends naturally to varying-coefficient models. For example, a model in which the coefficient on \(X_{1,it}\) varies smoothly with \(X_{2,it}\), \[ Y_{it} = \beta_0(X_{2,it}) + X_{1,it} \beta_1(X_{2,it}) + \mu_i + \varepsilon_{it}, \] can be estimated using a by-variable smooth:
The simulations below focus on one-dimensional smooths; the varying-coefficient and tensor-product examples are included to show how the same fixed-effects syntax extends to richer specifications.
For a fully nonparametric bivariate component \(g(X_{1,it},X_{2,it})\), one can use a tensor-product smooth:
If the goal is to separate main effects from interactions, ti() provides the corresponding tensor-product interaction:
The basis dimension \(k\) sets the maximum flexibility of the smooth and is chosen by the researcher; conditional on this basis, mgcv chooses the smoothing parameter by fREML. If \(k\) is omitted, mgcv uses default basis dimensions, but these defaults should be checked in applications.
For bam, the default smoothing parameter selection criterion is fREML, corresponding to method = "fREML". Researchers who wish to use the GCV/Cp criterion instead can specify method = "GCV.Cp". The default choices are best viewed as convenient and computationally efficient options, rather than binding or unusual constraints imposed on the estimator.
After fitting any of the fixed-effects, penalized-unit-effects, or first-differenced specifications, standard errors are computed from the fitted bam object rather than from the default mgcv covariance matrix. The fitted linear predictor matrix is obtained with
and the unit-level score sums are formed by summing the row-wise products of the design matrix and residuals within unit:
This command produces $\hat{\Omega}$ from equation (ref).
For the penalty-adjusted covariance estimator, one needs to obtain $B_{\lambda}^{-1}$ from equation (ref). In the Gaussian identity-link case this is obtained from the Bayesian covariance matrix returned by mgcv, rescaled by the estimated residual variance:
The simulations also apply the usual finite-sample cluster correction for a one-way cluster-robust covariance matrix cameron_miller_2015. In this panel setting, each unit is one cluster, so the number of clusters is \(n\). Let \(N\) denote the total number of observations and let \(d\) denote the substantive degrees of freedom of the fitted model. The correction used below is \[ \frac{n}{n-1}\frac{N-1}{N-d}. \] For a balanced panel, \(N=nT\); for an unbalanced panel, \(N\) is simply the total number of observed unit-time cells. For the fixed-effects specification, \(d\) excludes nuisance unit indicators, matching the within-estimator convention rather than inflating the adjustment by counting every unit effect. In code, the correction is:
For pointwise and uniform bands for $g(\cdot)$, the fitted smooth is evaluated on a grid by calling predict(..., type = "lpmatrix") on the grid data. As in Section (ref), the full prediction matrix should not be used unchanged, because it includes the intercept, linear regressors, and unit effects. In code, the construction of \(L_G=C_G R_G^{(g)}\) corresponds to selecting the columns for s(x2), zeroing the remaining columns, and centering the selected grid rows:
Pointwise intervals only require the diagonal of C_grid:
These intervals are appropriate for statements about a fixed grid point, such as the value of the centered smooth at \(x=0.5\). A plotted collection of pointwise intervals is visually useful, but it is not a simultaneous 95% band for the curve.
For a uniform band, first convert the covariance matrix into the correlation matrix of the standardized grid errors, simulate draws from the corresponding Gaussian approximation, and take the 95th percentile of the maximum absolute draw. The simulations below use \(499\) Gaussian draws for this step:
The only difference between the pointwise and uniform bands is the critical value: \(1.96\) for each separate point versus the simulated sup-norm critical value for the whole grid. The number of Gaussian draws controls only the simulation error in this critical value; using more draws, such as \(999\) or \(1{,}999\), gives a slightly more stable cutoff at additional computational cost. In the first-differenced implementation, the grid data evaluate the linear functional smooth with X_pair = x_grid and L_pair = 1; the same zeroing, centering, and pointwise-or-uniform critical-value steps then report the normalized level curve \(g_G^c\).
The simulation section studies both estimation accuracy and inference. For estimation of $g(\cdot)$, the benchmark loss is out-of-sample root mean squared error on an independently generated test set. For inference on $\beta$, the benchmark outcomes are empirical size under the null and empirical power under the alternative. For inference on $g(\cdot)$, the benchmark outcome is coverage of confidence bands.
For inference on $g(\cdot)$, I primarily use a penalty-adjusted cluster sandwich estimator. A series-style cluster sandwich estimator that ignores the penalty is too conservative, so its results are reported in the appendix. One could consider a wild cluster bootstrap that refits the bam model in each bootstrap draw, but it is left for future work.
The simulation designs deliberately correlate the unit effects with the regressors. For each unit, draw a latent index \(\alpha_i\sim \mathrm{Beta}(2.2,2.2)\) and set \(\mu_i=1.1(\alpha_i-\bar{\alpha})+\nu_i\), where \(\nu_i\sim N(0,1.25^2)\) is independent noise. In the one-regressor designs, \(X_{it}\) is a noisy draw around \(\alpha_i\), clipped to lie in $[0, 1]$. In the partially linear designs, \(X_{2,it}\) is generated the same way and \(X_{1,it}\) is correlated with both \(X_{2,it}\) and \(\alpha_i\). Thus the unobserved unit effect is not independent of the covariates.
The tables use a heteroskedastic AR(1) error design as the primary specification. Conditional variances are proportional to \(1+2\exp(0.75 Z_{it})\) and normalized to keep the average error variance fixed across designs. Here \(Z_{it}\) denotes the covariate index used in the variance function. In the estimation-accuracy designs and nonparametric coverage simulations, there is a single nonparametric regressor, so \(Z_{it}=X_{it}\). The partially linear outcome adds \(X_{1,it}\beta\), and for that design \(Z_{it}=X_{1,it}+X_{2,it}\). Within units, errors have AR(1) correlation with \(\rho=0.5\). All reported tables use \(M=1{,}000\) Monte Carlo draws.
I consider three data generating processes (DGPs):
and the outcome is then generated as \[ Y_{it} = g_j(X_{it}) + \mu_i + \varepsilon_{it}. \]
Table (ref) reports estimation accuracy for the three sample-size designs. Each entry reports $100$ times the root mean squared error for the centered function on an independently generated test set with $1{,}000$ test points per draw. The fixed spline estimators use $K=7$ basis terms. This fixed-\(K\) benchmark can perform well when its number of series terms “matches” the smoothness of the unknown function, but is unable to adapt to the smoothness automatically. In the sine design, the fixed-\(K\) spline performs especially well, but it overfits in the other two designs and is dominated by the bam estimators.
The semiparametric inference simulations use the partially linear model \[ Y_{it} = X_{1,it}\beta + g(X_{2,it}) + \mu_i + \varepsilon_{it}. \] The first exercise reports empirical rejection probabilities for tests of $H_0:\beta=1$, both under the null and under a fixed alternative.
Table (ref) reports size results for tests of $H_0:\beta=1$ in the partially linear model. The design uses $g(x)=\sin(2\pi x)$ and the heteroskedastic AR(1) error process described above.
Table (ref) reports power results for the same test when the true value is \(\beta=1.075\).
Tables (ref) and (ref) report penalty-adjusted inference for the unknown function $g$. Table (ref) continues with the partially linear model used for the \(\beta\) results, now focusing on the smooth component. Table (ref) reports the corresponding results for the fully nonparametric model \(Y_{it}=g(X_{it})+\mu_i+\varepsilon_{it}\). In both cases, the target is the centered function evaluated on a grid of 50 points, with $g(x)=\sin(2\pi x)$. The “Avg. pointwise” column averages marginal coverage over the grid points, while “Min. pointwise” reports the worst grid-point coverage. The “Uniform” column reports the fraction of Monte Carlo draws in which the entire grid curve is covered simultaneously by the sup-$t$ band, using \(499\) Gaussian draws to approximate the critical value in each fitted model. The corresponding series-style bands that ignore the penalty are reported in Appendix Table (ref).
Several findings are worth highlighting. For inference on \(\beta\), the three bam methods yield similar finite-sample behavior. Under the null hypothesis, empirical rejection rates tend to lie between about 5% and 6%, close to the nominal 5% level. Inference on \(g(\cdot)\) is more demanding. The bam FE and bam FD estimators tend to have coverage probabilities in the 92% to 94% range, while the penalized-unit-effect specification tends to have lower coverage, sometimes falling below 90%. This pattern is consistent with the shrinkage interpretation of s(unit, bs = "re"): when unit effects are correlated with regressors, shrinking them toward zero can introduce bias that the confidence bands do not account for.
The modest deviations from nominal coverage are consistent with the broader literature on cluster-robust inference. imbens_kolesar_2016 and pustejovsky_tipton_2018 show that conventional heteroskedasticity- and cluster-robust procedures can have finite-sample distortions in various settings, motivating degrees-of-freedom and leverage-based corrections. cameron_miller_2015 and mackinnon_nielsen_webb_2023 emphasize that cluster-robust inference is asymptotic in the number of clusters, so modest over-rejection can remain in finite samples even when the number of clusters is not especially small. The penalized-spline setting adds a further source of finite-sample distortion because the covariance calculation treats the estimated smoothing parameter as fixed, ignoring the additional variability from its selection.
Another possible reason for deviations from nominal coverage is smoothing bias in nonparametric series and spline estimators. Penalized series and spline estimators can have bias that is non-negligible relative to the standard error, and confidence bands built from the sandwich covariance matrix do not correct for this bias. The REML/fREML criterion targets prediction performance rather than coverage-optimal undersmoothing or bias correction, so the resulting bands can be slightly too narrow even when the variance estimate is accurate marra_wood_2012. This effect is most directly relevant to inference on \(g(\cdot)\) rather than to tests of finite-dimensional parameters, which is consistent with the presented simulation evidence.
A wild cluster bootstrap that refits the bam model in each draw could improve finite-sample properties and is left for future work.
This paper provides a practical guide to estimating semiparametric and nonparametric panel data models with fixed effects using the mgcv package in R. The proposed approach connects the penalized sieve and spline literature in theoretical econometrics with REML/fREML implementations of penalized splines in statistics. The fixed-effects estimators use unit indicators or first differencing, while penalized unit effects provide a regularized alternative with a distinct shrinkage interpretation. The same syntax can be applied to partially linear, varying-coefficient, fully nonparametric, and related specifications.
The focus is implementation. I provide mgcv::bam commands for estimation and show how to construct penalty-adjusted cluster-robust standard errors and confidence bands from the fitted model object. Monte Carlo simulations with correlated unit effects, heteroskedasticity, and serially correlated errors show good finite-sample performance in the one-dimensional smooth designs studied here: the estimators adapt reasonably well to the smoothness of the unknown function, tests for finite-dimensional parameters have empirical size close to nominal, and confidence bands for the centered nonparametric component have accurate coverage.