Nothing
Fixed Effects Individual Slopes using feisr"
In feisr: Estimating Fixed Effects Individual Slope Models
\newcommand{\Exp}{\mathrm{E}}
\newcommand\given[1][]{\:#1\vert\:}
\newcommand{\Cov}{\mathrm{Cov}}
\newcommand{\Var}{\mathrm{Var}}
\newcommand{\rank}{\mathrm{rank}}
\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
The main purpose of the package feisr is the estimation of fixed effects individual slope models and respective test statistics. The fixed effects individual slope (FEIS) estimator is a more general version of the well-known fixed effects estimator (FE), which allows to control for heterogeneous slopes in addition to time-constant heterogeneity [e.g. @Bruderl.2015.387; @Frees.2001.0; @Lemieux.1998; @Polachek.1994.0; @Ruttenauer.2020; @Wooldridge.2010.384]. Formally, the FEIS estimator can be expressed as
$$
\begin{align}
\bm y_{i} =& \bm X_{i}\bm\beta + \bm W_i \bm\alpha_i + \bm \epsilon_{i},
\end{align}
$$
where $\bm y_{i}$ is $T \times 1$, $\bm X_{i}$ is $T \times K$, and $\bm\epsilon_{i}$ is $T \times 1$. $\bm W_i$ is a $T \times J$ matrix of slope variables, and $\bm\alpha_i$ a $J \times 1$ vector of individual-specific slope parameters, for $J$ slope parameters including a constant term. If $\bm W_i$ consists of a constant term only, $\bm W_i = \bm 1$, thus $\bm\alpha_i$ reduces to $\alpha_{i1}$, and the above equation represents the well-known formula of a conventional FE model with individual fixed effects.
As with the conventional FE, FEIS can be estimated using lm() by including $N-1$ individual-specific dummies and interaction terms of each slope variable with the $N-1$ individual-specific dummies ($(N-1) *J$ controls). This is however highly inefficient. As with the conventional FE estimator, we can achieve the same result by running an lm() on pre-transformed data. Therefore, specify the 'residual maker' matrix $\bm M_i = \bm I_T - \bm W_i(\bm W^\intercal_i \bm W_i)^{-1}\bm W^\intercal_i$, and estimate
$$
\begin{align}
y_{it} - \hat{y}{it} =& (\bm x{it} - \hat{\bm x}{it})\bm\beta + \epsilon{it} - \hat{\epsilon}{it}, \
\bm M_i \bm y_i =& \bm M_i \bm X_i\bm\beta + \bm M_i \bm \epsilon{i}, \
\tilde{\bm y}{i} =& \tilde{\bm X}{i}\bm\beta + \tilde{\bm \epsilon}{i},
\end{align}
$$
where $\tilde{\bm y}{i}$, $\tilde{\bm X}{i}$, and $\tilde{\bm \epsilon}{i}$ are the residuals of regressing $\bm y_{i}$, each column-vector of $\bm X_{i}$, and $\bm \epsilon_{i}$ on $\bm W_i$. Intuitively, we (1) estimate the individual-specific predicted values for the dependent variable and each covariate based on an individual intercept and the additional slope variables of $\bm W_i$, (2) 'detrend' the original data by these individual-specific predicted values, and (3) run an OLS model on the residual data. Similarly, we can estimate a correlated random effects (CRE) model [@Chamberlain.1982; @Mundlak.1978.0; @Wooldridge.2010.384] including the individual specific predictions $\hat{\bm X}{i}$ to obtain the FEIS estimator:
$$
\begin{align}
\bm y{i} =& \bm X_{i}\bm\beta + \hat{\bm X}{i}\bm\rho + \bm \epsilon{i}.
\end{align}
$$
We use this transformation of the FEIS estimator to derive a Hausman-like Artificial Regression Test (ART) for hetergeneous slopes in panel or otherwise nested models.
The main functions of the feisr package are:
-
feis(): Estimates fixed effects individual slope estimators by applying linear lm models to 'detrended' data.
-
feistest(): Estimates a regression-based Hausman test for fixed effects individual slope models.
-
bsfeistest(): Estimates a bootstrapped Hausman test for fixed effects individual slope models.
-
slopes(): Extracts the individual slopes from feis objects.
The functions included in the R package feisr are also available in the xtfeis ado for Stata. The package plm [@Croissant.2008.0; @Croissant.2019] provides estimation of related models, as for instance the mean group (MG) or common correlated effects mean groups (CCEMG) estimator via pmg() function, and the estimator of models with variable coefficients via pvcm(). The package fixest also provides estimation of FE models with individual slopes in its function `feols', including higher dimensional FEIS models and instrumental variable approaches.
Estimating feis models
We demonstrate the most important functions of the feisr package by conducting an exemplary replication of the results presented in @Ludwig.2018.0. We therefore use the mwp panel data, containing information on wages and family status of 268 men. This is a random sample drawn from the National Longitudinal Survey of Youth [@NLSY79.2012], and more details on the selection of observations and variable construction can be found in @Ludwig.2018.0. To load the package and data use:
library(feisr)
data("mwp", package = "feisr")
head(mwp)
The data set contains a unique person identifier (id) and survey year indicator (year). Furthermore, we have information about the log hourly wage rate (lnwage), work experience (exp) and its square (expq), family status (marry), enrollment in current education (enrol), years of formal education education (yeduc), age (age), birth cohort (cohort), and a grouped year indicator (yeargr).
To illustrate the usage of the feisr package, we exemplary investigate the 'marriage wage premium': we analyze whether marriage leads to an increase in the hourly wage for men. We use the function feis to estimate fixed effects individual slope models to control for the hypothesis that those men who are more likely to marry or marry earlier, also have a steeper wage growth over time. Similar to the plm function, feis requires to indicate a unique person / group identifier. To include individual-specific slopes, feis uses two-part formulas (expr | slope_expr), where slope_expr gives the expression for modelling the individual slopes.
In the most simple setting, we just regress the wage (lnw) on the dichotomous marriage indicator (marry) and the control variable year. Instead of controlling for a common linear time trend (lnw ~ marry + year), we use year as slope variable in feis() and thus allow each individual in our dataset to have their own linear time trend (lnw ~ marry | year):
wages.feis0 <- feis(lnw ~ marry | year, data = mwp, id = "id")
summary(wages.feis0)
The summary() function prints the model output with the conventional information. The result of this simple model suggests that marriage increases the wage above what we would have expected based on each individual's linear time trend. Note however that the assumption of a linear time trend or a linear effect of other control variables might be spurious. In our example, it is highly unlikely that wages are linearly increasing with time / age. Thus, more complex transformations of the slope variable(s) might be necessary in many applications [see e.g. @Kneip.2009; @Wolfers.2006 for further examples].
An advantage of FEIS is that it allows to use any observed variable to be specified as individual-specific covariate. The formula in feis() also allows 'as is' transformations of original variables using I(). To replicate the analysis of @Ludwig.2018.0, we use work experience (exp) and squared work experience as the slope variables:
wages.feis <- feis(lnw ~ marry + enrol + yeduc + as.factor(yeargr)
| exp + I(exp^2), data = mwp, id = "id")
summary(wages.feis)
coef1 <- unname(round(wages.feis$coefficients[1], 3))
In the example above, we can see that marriage has only a small effect of r coef1, and is not statistically significant. Note that the (first stage) slope regression, 'detrending' the data, always includes an intercept to control for constant person-specific differences. The regression on the detrended data, by default, does not contain an overal intercept (intercept = FALSE). By default, the feis() functions returns conventional standard errors for homoscedastic errors (the scaled variance-covariance matrix is attached to the output objects and can be extracted by vcov()). Using the option robust = TRUE instead returns panel-robust standard errors [e.g. @Arellano.1993; @Berger2017; @Wooldridge.2010.384]:
wages.feis.rob <- feis(lnw ~ marry + enrol + yeduc + as.factor(yeargr)
| exp + I(exp^2), data = mwp, id = "id",
robust = TRUE)
To compare the results with a conventional fixed effects and random effects models, we use the well-known plm [@Croissant.2008.0] package and estimate a within and a random effects model. Note that the feis() function would produce identical results to a plm(..., model = "whitin", effect = "individual") model if specified with an intercept in the slope variables only (without any additional slope variables). However, for this case, we recommend using the well-established plm package instead.
library(plm)
wages.fe <- plm(lnw ~ marry + enrol + yeduc + as.factor(yeargr)
+ exp + I(exp^2), data = mwp, index = c("id", "year"),
model = "within", effect = "individual")
wages.re <- plm(lnw ~ marry + enrol + yeduc + as.factor(yeargr)
+ exp + I(exp^2), data = mwp, index = c("id", "year"),
model = "random", effect = "individual")
The difference between the three models produced above can be visualized using the screenreg() and texreg() functions of the texreg package [@Leifeld.2013]. Since version 1.37.1, the texreg packages comes with an extract method for feis models.
# Fallback for production of vignettes if texreg < 1.37.1 (internal extract.feis)
if(utils::packageVersion("texreg") < "1.37.1"){
library(texreg)
setMethod("extract", signature = className("feis", "feisr"),
definition = feisr:::extract.feis)
}
library(texreg)
screenreg(list(wages.feis, wages.feis.rob, wages.fe, wages.re), digits = 3,
custom.model.names = c("FEIS", "FEIS robust", "FE", "RE"))
coef2 <- unname(round(wages.fe$coefficients[1], 3))
coef3 <- unname(round(wages.re$coefficients[2], 3))
diff <- coef2 - coef1
diff2 <- coef3 - coef1
rel <- round(coef3/coef1, 1)
r2.fe <- unname(round(summary(wages.fe)$r.squared[1], 3))
r2.feis <- unname(round(wages.feis$r2, 3))
The most striking difference here is the effect of marriage on wage. The FE returns a highly significant effect of marriage on the logarithmic wage of r coef2, which exceeds the prediction of the FEIS by r diff. Consequently, according to the FE, we would conclude that there really is a 'marital wage premium' for men. The random effects model produces results, which are relatively similar to the results of the conventional FE model. However, the FEIS revises this conclusion by showing that a substantial part of the effect stems from the fact that men with a steeper wage growth tend to marry at higher rates (those men with a stronger effect of experience on wage also exhibit a stronger effect of experience on marriage). In consequence, we can conclude that the 'marital wage premium' mainly stems from selection on growth rather than from a causal effect of marriage on wage. The models also differ in terms of explained variance: $R^2$ of r r2.fe in the FE compared to r r2.feis in the FEIS model. However, this seems quite plausible, given that we have 'discarded' all the variance in wages, which can be explained by individual work experience in the FEIS model.
Hausman-like specification tests
The example above illustrates how heterogeneous growth curves related to the covariates can drastically influence the conclusions drawn in conventional fixed effects models. For applied research, we thus propose to test for the presence of heterogeneous slopes when using panel models. Therefore, the function feistest() provides a Hausman-like artificial regression test [as discussed in @Ruttenauer.2020], which estimates the Mundlak specification [@Mundlak.1978.0] of the FEIS model and performs a Wald $\chi^2$ test on the coefficients of the individual-specific predicted values using the aod package [@Lesnoff.2012]. The test can be performed on the FEIS model estimated by feis(), and no further model specifications are necessary. By default, the test is jointly performed on all covariates given on the right side of the formula. If required, the test can also be restricted to specified covariates only using the option terms (e.g. feistest(wages.feis, terms = c("marry", "enrol"))). With the option type = "all", the function feistest() provides a test of FEIS against FE, the regression-based version of the well-known Hausman test comparing FE against RE [@Hausman.1978.0; @Wooldridge.2010.384], and a third test comparing the FEIS directly against the RE model. In the example below, we use cluster-robust standard errors for the artificial regression test. A summary method is provided in the package:
ht <- feistest(wages.feis, robust = TRUE, type = "all")
summary(ht)
chi2_feis <- round(unname(ht$wald_feis$result$chi2[1]), 3)
chi2_fe <- round(unname(ht$wald_fe$result$chi2[1]), 3)
p_fe <- round(unname(ht$wald_fe$result$chi2[3]), 3)
chi2_re <- round(unname(ht$wald_re$result$chi2[1]), 3)
p_re <- round(unname(ht$wald_re$result$chi2[3]), 3)
Alternatively, we can also use a bootstrapped Hausman test to compare the FEIS, the FE, and the RE model. The function bsfeistest() performs bootstrapping by pairs cluster resampling from the original dataset with replacement, meaning that we run a feis, a plm within model, and a plm random model for each rep random sample. Note that resampling is perfomed on the cluster ids. Thus, for unbalanced samples, only the number of clusters is identical in each replication while the total number of observations can vary. This method is analogue to the 'bootstrap-se' method described in @Cameron.2008, allowing us to receive an empirical estimate of the variance-covariance matrix ${\bm V}_{\hat{\bm \beta}_1 - \hat{\bm \beta}_0}$ and to perform the original version of the Hausman test [@Hausman.1978.0].
bsht <- bsfeistest(wages.feis, type = "all", rep = 100, seed = 91020104)
summary(bsht)
bschi2_feis <- round(unname(bsht$wald_feis$result$chi2[1]), 3)
bschi2_fe <- round(unname(bsht$wald_fe$result$chi2[1]), 3)
bsp_fe <- round(unname(bsht$wald_fe$result$chi2[3]), 3)
bschi2_re <- round(unname(bsht$wald_re$result$chi2[1]), 3)
bsp_re <- round(unname(bsht$wald_re$result$chi2[3]), 3)
The first test comparing FEIS against FE offers a clear rejection of the null- hypothesis that FE is consistent. A highly significant $\chi^2$ of r chi2_feis (r bschi2_feis in the bootstrapped version) indicates that estimates of the conventional FE are inconsistent because of heterogeneous slopes. Thus, we should use FEIS instead of FE. Interestingly, the second test of FE against RE gives a non-significant $\chi^2$ of r chi2_fe with p=r p_fe (or r bschi2_fe with p=r bsp_fe respectively), indicating that we could use an RE model instead of a conventional FE model. This result is in line with the fact that RE and FE produce relatively similar coefficients in our example.
Thus, when RE is the preferred model, we highly recommend to test the RE against the FE and against the FEIS model. The third test statistic of feistest() and bsfeistest() offers a direct comparison of FEIS against RE: both versions -- artificial regression based and bootstrapped -- show a highly significant $\chi^2$ of r chi2_re with p=r p_re (or r bschi2_re with p=r bsp_re respectively), thereby indicating that we should reject the null-hypothesis of consistent estimates in the RE model. In sum, when ignoring the possibility of heterogeneous slopes, we would lean towards using a conventional RE model, thereby relying on an effect of marriage which equals r rel times the effect obtained in a model accounting for individual-specific slopes.
Extracting individual slopes
Occasionally, estimates of individual slopes might be of interest in themselves. For example, we could compute the average conditional slopes over the sample of $i$ as estimates for $E(\bm \alpha_{i2})$. Similarly, we might compare estimated slopes across treatment groups. Using the function slopes() on an object of class "feis" returns the individual slope estimates for each id:
alpha <- slopes(wages.feis)
head(alpha)
As shown above, we receive a matrix which contains the individual ids in the rownames. For each individual, we get an intercept (the individual fixed effects as in a conventional FE), and a coefficient for each slope variable (in our case exp and squared exp). We can use the rownames to merge the individual slopes to the original data.
colnames(alpha) <- paste0("slp_", colnames(alpha))
alpha.df <- data.frame(alpha, id = rownames(alpha))
mwp <- merge(mwp, alpha.df, by = "id")
Thus, we can also use the individual slope parameters for further transformation or analysis. As an example, we plot the predicted ln wage trends by marriage status using ggplot2 [@Wickham.2016], and also add the average trends by marriage group (ever married vs. never married).
### Individual predicted trend of lnw (based on slopes)
mwp$pred <- mwp$slp_.Intercept. + mwp$exp * mwp$slp_exp + mwp$exp^2 * mwp$slp_I.exp.2.
### Average value by evermarry and exp bins
mwp$exp_gr <- round(mwp$exp, 0)
aggr <- aggregate(mwp$pred, by = list(exp = mwp$exp_gr, evermarry = mwp$evermarry), mean)
### Plot
library(ggplot2)
zp1 <- ggplot(data = mwp, aes(x = exp, y = pred)) +
geom_line(aes(group = id, col = as.factor(marry)),
linetype = "solid", lwd = 0.5, alpha = 0.4) +
geom_line(data = aggr, aes(y = x, linetype = as.factor(evermarry)),
alpha = 1, size = 1.2) +
labs(color = "Married", linetype = "Ever-married") +
guides(colour = guide_legend(override.aes = list(alpha = 1))) +
labs(y = "Predicted log wage growth", x = "Work experience") +
theme_bw()
zp1
This example shows that the data exhibits strong heterogeneity in the predicted wage growth over experience. It also shows, that even conditional on the overall marriage effect -- marry is controlled for in the model from which the slopes are extracted -- ever-married individuals, on average, have a stronger wage growth over the years than never-married men. Furthermore, the average trends already start to diverge when most respondents are still unmarried. This leads to a violation of the parallel trends assumption and biased point estimates in conventional FE models.
For more information on FEIS models, the respective test statistics, and applied examples, see also @Ruttenauer.2020.
References
Try the feisr package in your browser
Any scripts or data that you put into this service are public.
feisr documentation built on April 1, 2022, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
\newcommand{\Exp}{\mathrm{E}} \newcommand\given[1][]{\:#1\vert\:} \newcommand{\Cov}{\mathrm{Cov}} \newcommand{\Var}{\mathrm{Var}} \newcommand{\rank}{\mathrm{rank}} \newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}}
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
The main purpose of the package feisr is the estimation of fixed effects individual slope models and respective test statistics. The fixed effects individual slope (FEIS) estimator is a more general version of the well-known fixed effects estimator (FE), which allows to control for heterogeneous slopes in addition to time-constant heterogeneity [e.g. @Bruderl.2015.387; @Frees.2001.0; @Lemieux.1998; @Polachek.1994.0; @Ruttenauer.2020; @Wooldridge.2010.384]. Formally, the FEIS estimator can be expressed as
$$ \begin{align} \bm y_{i} =& \bm X_{i}\bm\beta + \bm W_i \bm\alpha_i + \bm \epsilon_{i}, \end{align} $$ where $\bm y_{i}$ is $T \times 1$, $\bm X_{i}$ is $T \times K$, and $\bm\epsilon_{i}$ is $T \times 1$. $\bm W_i$ is a $T \times J$ matrix of slope variables, and $\bm\alpha_i$ a $J \times 1$ vector of individual-specific slope parameters, for $J$ slope parameters including a constant term. If $\bm W_i$ consists of a constant term only, $\bm W_i = \bm 1$, thus $\bm\alpha_i$ reduces to $\alpha_{i1}$, and the above equation represents the well-known formula of a conventional FE model with individual fixed effects.
As with the conventional FE, FEIS can be estimated using lm() by including $N-1$ individual-specific dummies and interaction terms of each slope variable with the $N-1$ individual-specific dummies ($(N-1) *J$ controls). This is however highly inefficient. As with the conventional FE estimator, we can achieve the same result by running an lm() on pre-transformed data. Therefore, specify the 'residual maker' matrix $\bm M_i = \bm I_T - \bm W_i(\bm W^\intercal_i \bm W_i)^{-1}\bm W^\intercal_i$, and estimate
$$
\begin{align}
y_{it} - \hat{y}{it} =& (\bm x{it} - \hat{\bm x}{it})\bm\beta + \epsilon{it} - \hat{\epsilon}{it}, \
\bm M_i \bm y_i =& \bm M_i \bm X_i\bm\beta + \bm M_i \bm \epsilon{i}, \
\tilde{\bm y}{i} =& \tilde{\bm X}{i}\bm\beta + \tilde{\bm \epsilon}{i},
\end{align}
$$
where $\tilde{\bm y}{i}$, $\tilde{\bm X}{i}$, and $\tilde{\bm \epsilon}{i}$ are the residuals of regressing $\bm y_{i}$, each column-vector of $\bm X_{i}$, and $\bm \epsilon_{i}$ on $\bm W_i$. Intuitively, we (1) estimate the individual-specific predicted values for the dependent variable and each covariate based on an individual intercept and the additional slope variables of $\bm W_i$, (2) 'detrend' the original data by these individual-specific predicted values, and (3) run an OLS model on the residual data. Similarly, we can estimate a correlated random effects (CRE) model [@Chamberlain.1982; @Mundlak.1978.0; @Wooldridge.2010.384] including the individual specific predictions $\hat{\bm X}{i}$ to obtain the FEIS estimator:
$$
\begin{align}
\bm y{i} =& \bm X_{i}\bm\beta + \hat{\bm X}{i}\bm\rho + \bm \epsilon{i}.
\end{align}
$$
We use this transformation of the FEIS estimator to derive a Hausman-like Artificial Regression Test (ART) for hetergeneous slopes in panel or otherwise nested models.
The main functions of the feisr package are:
-
feis(): Estimates fixed effects individual slope estimators by applying linearlmmodels to 'detrended' data. -
feistest(): Estimates a regression-based Hausman test for fixed effects individual slope models. -
bsfeistest(): Estimates a bootstrapped Hausman test for fixed effects individual slope models. -
slopes(): Extracts the individual slopes fromfeisobjects.
The functions included in the R package feisr are also available in the xtfeis ado for Stata. The package plm [@Croissant.2008.0; @Croissant.2019] provides estimation of related models, as for instance the mean group (MG) or common correlated effects mean groups (CCEMG) estimator via pmg() function, and the estimator of models with variable coefficients via pvcm(). The package fixest also provides estimation of FE models with individual slopes in its function `feols', including higher dimensional FEIS models and instrumental variable approaches.
Estimating feis models
We demonstrate the most important functions of the feisr package by conducting an exemplary replication of the results presented in @Ludwig.2018.0. We therefore use the mwp panel data, containing information on wages and family status of 268 men. This is a random sample drawn from the National Longitudinal Survey of Youth [@NLSY79.2012], and more details on the selection of observations and variable construction can be found in @Ludwig.2018.0. To load the package and data use:
library(feisr) data("mwp", package = "feisr") head(mwp)
The data set contains a unique person identifier (id) and survey year indicator (year). Furthermore, we have information about the log hourly wage rate (lnwage), work experience (exp) and its square (expq), family status (marry), enrollment in current education (enrol), years of formal education education (yeduc), age (age), birth cohort (cohort), and a grouped year indicator (yeargr).
To illustrate the usage of the feisr package, we exemplary investigate the 'marriage wage premium': we analyze whether marriage leads to an increase in the hourly wage for men. We use the function feis to estimate fixed effects individual slope models to control for the hypothesis that those men who are more likely to marry or marry earlier, also have a steeper wage growth over time. Similar to the plm function, feis requires to indicate a unique person / group identifier. To include individual-specific slopes, feis uses two-part formulas (expr | slope_expr), where slope_expr gives the expression for modelling the individual slopes.
In the most simple setting, we just regress the wage (lnw) on the dichotomous marriage indicator (marry) and the control variable year. Instead of controlling for a common linear time trend (lnw ~ marry + year), we use year as slope variable in feis() and thus allow each individual in our dataset to have their own linear time trend (lnw ~ marry | year):
wages.feis0 <- feis(lnw ~ marry | year, data = mwp, id = "id") summary(wages.feis0)
The summary() function prints the model output with the conventional information. The result of this simple model suggests that marriage increases the wage above what we would have expected based on each individual's linear time trend. Note however that the assumption of a linear time trend or a linear effect of other control variables might be spurious. In our example, it is highly unlikely that wages are linearly increasing with time / age. Thus, more complex transformations of the slope variable(s) might be necessary in many applications [see e.g. @Kneip.2009; @Wolfers.2006 for further examples].
An advantage of FEIS is that it allows to use any observed variable to be specified as individual-specific covariate. The formula in feis() also allows 'as is' transformations of original variables using I(). To replicate the analysis of @Ludwig.2018.0, we use work experience (exp) and squared work experience as the slope variables:
wages.feis <- feis(lnw ~ marry + enrol + yeduc + as.factor(yeargr) | exp + I(exp^2), data = mwp, id = "id") summary(wages.feis)
coef1 <- unname(round(wages.feis$coefficients[1], 3))
In the example above, we can see that marriage has only a small effect of r coef1, and is not statistically significant. Note that the (first stage) slope regression, 'detrending' the data, always includes an intercept to control for constant person-specific differences. The regression on the detrended data, by default, does not contain an overal intercept (intercept = FALSE). By default, the feis() functions returns conventional standard errors for homoscedastic errors (the scaled variance-covariance matrix is attached to the output objects and can be extracted by vcov()). Using the option robust = TRUE instead returns panel-robust standard errors [e.g. @Arellano.1993; @Berger2017; @Wooldridge.2010.384]:
wages.feis.rob <- feis(lnw ~ marry + enrol + yeduc + as.factor(yeargr) | exp + I(exp^2), data = mwp, id = "id", robust = TRUE)
To compare the results with a conventional fixed effects and random effects models, we use the well-known plm [@Croissant.2008.0] package and estimate a within and a random effects model. Note that the feis() function would produce identical results to a plm(..., model = "whitin", effect = "individual") model if specified with an intercept in the slope variables only (without any additional slope variables). However, for this case, we recommend using the well-established plm package instead.
library(plm) wages.fe <- plm(lnw ~ marry + enrol + yeduc + as.factor(yeargr) + exp + I(exp^2), data = mwp, index = c("id", "year"), model = "within", effect = "individual") wages.re <- plm(lnw ~ marry + enrol + yeduc + as.factor(yeargr) + exp + I(exp^2), data = mwp, index = c("id", "year"), model = "random", effect = "individual")
The difference between the three models produced above can be visualized using the screenreg() and texreg() functions of the texreg package [@Leifeld.2013]. Since version 1.37.1, the texreg packages comes with an extract method for feis models.
# Fallback for production of vignettes if texreg < 1.37.1 (internal extract.feis) if(utils::packageVersion("texreg") < "1.37.1"){ library(texreg) setMethod("extract", signature = className("feis", "feisr"), definition = feisr:::extract.feis) }
library(texreg) screenreg(list(wages.feis, wages.feis.rob, wages.fe, wages.re), digits = 3, custom.model.names = c("FEIS", "FEIS robust", "FE", "RE"))
coef2 <- unname(round(wages.fe$coefficients[1], 3)) coef3 <- unname(round(wages.re$coefficients[2], 3)) diff <- coef2 - coef1 diff2 <- coef3 - coef1 rel <- round(coef3/coef1, 1) r2.fe <- unname(round(summary(wages.fe)$r.squared[1], 3)) r2.feis <- unname(round(wages.feis$r2, 3))
The most striking difference here is the effect of marriage on wage. The FE returns a highly significant effect of marriage on the logarithmic wage of r coef2, which exceeds the prediction of the FEIS by r diff. Consequently, according to the FE, we would conclude that there really is a 'marital wage premium' for men. The random effects model produces results, which are relatively similar to the results of the conventional FE model. However, the FEIS revises this conclusion by showing that a substantial part of the effect stems from the fact that men with a steeper wage growth tend to marry at higher rates (those men with a stronger effect of experience on wage also exhibit a stronger effect of experience on marriage). In consequence, we can conclude that the 'marital wage premium' mainly stems from selection on growth rather than from a causal effect of marriage on wage. The models also differ in terms of explained variance: $R^2$ of r r2.fe in the FE compared to r r2.feis in the FEIS model. However, this seems quite plausible, given that we have 'discarded' all the variance in wages, which can be explained by individual work experience in the FEIS model.
Hausman-like specification tests
The example above illustrates how heterogeneous growth curves related to the covariates can drastically influence the conclusions drawn in conventional fixed effects models. For applied research, we thus propose to test for the presence of heterogeneous slopes when using panel models. Therefore, the function feistest() provides a Hausman-like artificial regression test [as discussed in @Ruttenauer.2020], which estimates the Mundlak specification [@Mundlak.1978.0] of the FEIS model and performs a Wald $\chi^2$ test on the coefficients of the individual-specific predicted values using the aod package [@Lesnoff.2012]. The test can be performed on the FEIS model estimated by feis(), and no further model specifications are necessary. By default, the test is jointly performed on all covariates given on the right side of the formula. If required, the test can also be restricted to specified covariates only using the option terms (e.g. feistest(wages.feis, terms = c("marry", "enrol"))). With the option type = "all", the function feistest() provides a test of FEIS against FE, the regression-based version of the well-known Hausman test comparing FE against RE [@Hausman.1978.0; @Wooldridge.2010.384], and a third test comparing the FEIS directly against the RE model. In the example below, we use cluster-robust standard errors for the artificial regression test. A summary method is provided in the package:
ht <- feistest(wages.feis, robust = TRUE, type = "all") summary(ht)
chi2_feis <- round(unname(ht$wald_feis$result$chi2[1]), 3) chi2_fe <- round(unname(ht$wald_fe$result$chi2[1]), 3) p_fe <- round(unname(ht$wald_fe$result$chi2[3]), 3) chi2_re <- round(unname(ht$wald_re$result$chi2[1]), 3) p_re <- round(unname(ht$wald_re$result$chi2[3]), 3)
Alternatively, we can also use a bootstrapped Hausman test to compare the FEIS, the FE, and the RE model. The function bsfeistest() performs bootstrapping by pairs cluster resampling from the original dataset with replacement, meaning that we run a feis, a plm within model, and a plm random model for each rep random sample. Note that resampling is perfomed on the cluster ids. Thus, for unbalanced samples, only the number of clusters is identical in each replication while the total number of observations can vary. This method is analogue to the 'bootstrap-se' method described in @Cameron.2008, allowing us to receive an empirical estimate of the variance-covariance matrix ${\bm V}_{\hat{\bm \beta}_1 - \hat{\bm \beta}_0}$ and to perform the original version of the Hausman test [@Hausman.1978.0].
bsht <- bsfeistest(wages.feis, type = "all", rep = 100, seed = 91020104) summary(bsht)
bschi2_feis <- round(unname(bsht$wald_feis$result$chi2[1]), 3) bschi2_fe <- round(unname(bsht$wald_fe$result$chi2[1]), 3) bsp_fe <- round(unname(bsht$wald_fe$result$chi2[3]), 3) bschi2_re <- round(unname(bsht$wald_re$result$chi2[1]), 3) bsp_re <- round(unname(bsht$wald_re$result$chi2[3]), 3)
The first test comparing FEIS against FE offers a clear rejection of the null- hypothesis that FE is consistent. A highly significant $\chi^2$ of r chi2_feis (r bschi2_feis in the bootstrapped version) indicates that estimates of the conventional FE are inconsistent because of heterogeneous slopes. Thus, we should use FEIS instead of FE. Interestingly, the second test of FE against RE gives a non-significant $\chi^2$ of r chi2_fe with p=r p_fe (or r bschi2_fe with p=r bsp_fe respectively), indicating that we could use an RE model instead of a conventional FE model. This result is in line with the fact that RE and FE produce relatively similar coefficients in our example.
Thus, when RE is the preferred model, we highly recommend to test the RE against the FE and against the FEIS model. The third test statistic of feistest() and bsfeistest() offers a direct comparison of FEIS against RE: both versions -- artificial regression based and bootstrapped -- show a highly significant $\chi^2$ of r chi2_re with p=r p_re (or r bschi2_re with p=r bsp_re respectively), thereby indicating that we should reject the null-hypothesis of consistent estimates in the RE model. In sum, when ignoring the possibility of heterogeneous slopes, we would lean towards using a conventional RE model, thereby relying on an effect of marriage which equals r rel times the effect obtained in a model accounting for individual-specific slopes.
Extracting individual slopes
Occasionally, estimates of individual slopes might be of interest in themselves. For example, we could compute the average conditional slopes over the sample of $i$ as estimates for $E(\bm \alpha_{i2})$. Similarly, we might compare estimated slopes across treatment groups. Using the function slopes() on an object of class "feis" returns the individual slope estimates for each id:
alpha <- slopes(wages.feis) head(alpha)
As shown above, we receive a matrix which contains the individual ids in the rownames. For each individual, we get an intercept (the individual fixed effects as in a conventional FE), and a coefficient for each slope variable (in our case exp and squared exp). We can use the rownames to merge the individual slopes to the original data.
colnames(alpha) <- paste0("slp_", colnames(alpha)) alpha.df <- data.frame(alpha, id = rownames(alpha)) mwp <- merge(mwp, alpha.df, by = "id")
Thus, we can also use the individual slope parameters for further transformation or analysis. As an example, we plot the predicted ln wage trends by marriage status using ggplot2 [@Wickham.2016], and also add the average trends by marriage group (ever married vs. never married).
### Individual predicted trend of lnw (based on slopes) mwp$pred <- mwp$slp_.Intercept. + mwp$exp * mwp$slp_exp + mwp$exp^2 * mwp$slp_I.exp.2. ### Average value by evermarry and exp bins mwp$exp_gr <- round(mwp$exp, 0) aggr <- aggregate(mwp$pred, by = list(exp = mwp$exp_gr, evermarry = mwp$evermarry), mean) ### Plot library(ggplot2) zp1 <- ggplot(data = mwp, aes(x = exp, y = pred)) + geom_line(aes(group = id, col = as.factor(marry)), linetype = "solid", lwd = 0.5, alpha = 0.4) + geom_line(data = aggr, aes(y = x, linetype = as.factor(evermarry)), alpha = 1, size = 1.2) + labs(color = "Married", linetype = "Ever-married") + guides(colour = guide_legend(override.aes = list(alpha = 1))) + labs(y = "Predicted log wage growth", x = "Work experience") + theme_bw() zp1
This example shows that the data exhibits strong heterogeneity in the predicted wage growth over experience. It also shows, that even conditional on the overall marriage effect -- marry is controlled for in the model from which the slopes are extracted -- ever-married individuals, on average, have a stronger wage growth over the years than never-married men. Furthermore, the average trends already start to diverge when most respondents are still unmarried. This leads to a violation of the parallel trends assumption and biased point estimates in conventional FE models.
For more information on FEIS models, the respective test statistics, and applied examples, see also @Ruttenauer.2020.
References
Try the feisr package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.