La identificación del efecto causal \(D \rightarrow Y\) mediante técnicas que, como la regresión múltiple o las diversas variedades de matching, controlan por características potencialmente confusoras puede verse limitada debido a la sospecha de confusores omitidos en el análisis. Este es un problema que puede resultar desconcertante, en principio, pues estos confusores omitidos pueden ser o bien variables que se prestan poco a la medición (como la ‘’motivación’’, en estudios educativos) o bien variables que se sospecha intervienen en el proceso generador de datos, pero no se tiene una idea precisa de qué tipo de factor podría ser. O sencillamente por el hecho de que siempre puede existir una variable confusora cuya existencia ignoramos, pero no por eso deja de influir en nuestros resultados.

Así planteada, la dificultad parece insalvable, pues sólo se puede controlar por las variables de existencia conocida y que pueden ser medidas… a menos que se puedan reunir observaciones para cada unidad en distintos momentos en el tiempo. Los llamados datos panel permiten descontar el efecto de las características que distinguen a las unidades entre sí pero se mantienen constantes en el tiempo, sin necesidad de observarlas (es decir, sin necesidad de medirlas e incluirlas de manera explícita en el modelo).

Modelo

Sean \(Y_{it}\) el valor de la variable de respuesta en la unidad \(i\) en el momento \(t\) y \(D_{it}\) la variable de tratamiento, que puede ser numérica o indicadora, donde \(i = 1, \dots, n\) y \(t = 1, \dots, T\). Nótese que el número de observaciones en la base de datos es \(n \times T\).

La generación del valor de \(Y_{it}\) responde a la siguiente ecuación: \[\begin{align} Y_{it} = \alpha + \mathbf{A}'_i \boldsymbol{\gamma} + \lambda_t + \delta D_{it} + \mathbf{X}'_{it} \boldsymbol{\beta} + \varepsilon_{it}, \tag{1} \end{align}\]

donde \(\mathbf{X}_{it}\) es un vector de características observadas que cambian en el tiempo y \(\mathbf{A}_i\) es un vector de características no observadas y que, para cada unidad \(i\), mantienen el mismo valor durante todos los momentos \(t\) (es en ese sentido que son invariantes, si bien se espera que estas características varíen entre las unidades). Finalmente, \(\lambda_t\) representa una tendencia temporal de \(Y\) que es igual para todas las unidades.

El efecto causal de interés es el parámetro \(\delta\).

Con mediciones repetidas en el tiempo, se puede dar cuenta de las características invariantes no observadas en \(\mathbf{A}_i\) de dos maneras. La primera de ellas comienza por expresar la ecuación en los siguientes términos: \[\begin{align} Y_{it} = \alpha_i + \lambda_t + \delta D_{it} + \mathbf{X}'_{it} \boldsymbol{\beta} + \varepsilon_{it}, \tag{2} \end{align}\] donde \[\alpha_i = \alpha + \mathbf{A}'_i \boldsymbol{\gamma}.\] En esta especificación, \(\alpha_i\) representa un intercepto para cada una de las unidades que incluye las características de las unidades que son invariantes en el tiempo.

Esta formulación permite estimar el efecto de \(D\) sobre \(Y\) tratando a los interceptos \(\alpha_i\) (los “efectos fijos”) como parámetros a estimar con variables indicadoras. El efecto tiempo también puede ser descontado con variables indicadoras para cada momento en que se tomaron las mediciones. La ecuación se estima con mínimos cuadrados mediante: \[\begin{align*} Y_{it} = \alpha + \mathbf{G}'_{i} \boldsymbol{\Gamma} + \mathbf{T}'_{t} \boldsymbol{\lambda} + \delta D_{it} + \mathbf{X}'_{it} \boldsymbol{\beta} + \varepsilon_{it}, \end{align*}\] donde \(\mathbf{G}_i\) es un vector de variables indicadoras para \(n-1\) unidades y \(\mathbf{T}_t\) es un vector de variables indicadoras para \(T-1\) períodos. Nótese que los interceptos por unidad, \(\alpha + \Gamma_i = \alpha_i\), contienen el efecto de las características invariantes.

¿Qué resuelve esto? Si sólo se tienen datos para un momento en el tiempo, la estimación de \(\delta\) en la ecuación exgiría tener mediciones de todas las variables en \(\mathbf{X}\) y \(\mathbf{A}\). Si alguna o todas las variables en \(\mathbf{A}\) no pueden ser observadas porque su medición es difícil o porque desconocemos su existencia, el estimador estaría irremediablemente sesgado.

Al disponer de distintas observaciones en el tiempo, con los coeficientes correspondientes a la indicadoras de cada unidad se “descuenta” el efecto de todas las características invariantes (conocidas y desconocidas) sin necesidad de incluirlas en el modelo (y, por lo tanto, de medirlas).

La otra manera de descontar el efecto de estas características consiste en tomar la media de \(Y\) para la unidad \(i\), \[\begin{align*} \overline{Y}_i = \alpha_i + \bar{\lambda} + \delta \overline{D}_i + \mathbf{\overline{X}}'_{i} \boldsymbol{\beta} + \bar{\varepsilon}_i, \end{align*}\] y sustraerla de \(Y_{it}\) en la ecuación , \[\begin{align} Y_{it} - \overline{Y}_i = \lambda_t - \bar{\lambda} + \delta \left(D_{it} - \overline{D}_i \right) + \left( \mathbf{X}_{it} - \mathbf{\overline{X}}_{i} \right)' \boldsymbol{\beta} + \left(\varepsilon_{it} - \bar{\varepsilon}_i \right). \tag{3} \end{align}\]

Es decir que el efecto de las características invariantes se cancela al sustraer la media de la unidad, por lo que, una vez más, no es necesario medir las variables en \(\mathbf{A}\) para descontar su efecto. La ecuación anterior también se estima con mínimos cuadrados, con una precaución Si el modelo incluye \(K\) regresores, los grados de libertad asumidos en la regresión múltiple son \(nT - K -1\). Sin embargo, se pierde un grado de libertad por cada efecto fijo suprimido (las \(\alpha_i\)), de modo que los grados de libertad son \(nT - K - n = n \left(T-1\right) - K\), igual que en la especificación con variables indicadoras. Hecha esta corrección, los resultados del modelo con interceptos por unidad son idénticos a los del modelo diferencia con respecto a la media.

Supuesto de identificación

Si se asume que todas las características confusoras son invariantes en el tiempo o, de manera más razonable (pero no siempre enteramente convicente), que se ha controlado por las características que cambian en el tiempo para cada unidad, entonces el tratamiento es independiente de los estados potenciales y la identificación por efectos fijos proporciona un estimador insesgado del efecto causal.

Ejemplo

En esta sección se utilizarán los siguientes paquetes:

## Importar de Stata
install.packages("foreign")

## Análisis
install.packages("plm")
install.packages("lmtest")
install.packages("car")

## Gráficas 
install.packages("ggplot2")
install.packages("ggstance")
install.packages("jtools")
install.packages("wesanderson")

Para ilustrar la técnica de efectos fijos, usaremos el estudio de Finkel y Smith (2011) sobre un programa de educación cívica implementado en Kenia entre 2001 y 2002. El estudio valora el impacto del programa sobre los actitudes y valores de los participantes. Como la participación en el programa fue de carácter voluntario, cabe esperar que los participantes compartan características que los hagan sistemáticamente distintos de los no participantes, un caso típico de autoselección.

La información se encuentra en el archivo EF.data.1

library(foreign)
ef <- read.dta("EF.dta")

La variable de resultado es tol_w, un índice de tolerancia formado por el promedio de las respuestas a cuatro preguntas, cada una de ellas medida en una escala de 1 a 4. Valores mayores indican actitudes más tolerantes. En el estudio participaron 2301 sujetos, 1,489 de los cuales participaron en al menos uno de los talleres de educación cívica del mencionado programa. La variable multce_w reporta el número de talleres a los que cada entrevistado asistió, pero el valor máximo, el 4, representa la asistencia a cuatro o más cursos. En el análisis, esta variable es tratada como numérica.

A fin de ilustrar la manera en la que el modelo de efectos fijos hace una diferencia en la estimación, primero calculamos el estimador inocente, que en este caso sería el coeficiente de la regresión del índice de tolerancia sobre el número de talleres. Para ello, usamos sólo el subconjunto de la base de datos con las mediciones tomadas después de concluido el programa.

ols <- lm(tol_w ~ multce_w, data = subset(ef, ef$waveint == 1)) 
summary(ols)
## 
## Call:
## lm(formula = tol_w ~ multce_w, data = subset(ef, ef$waveint == 
##     1))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2579 -0.8928 -0.3928  0.7288  2.2288 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.77115    0.02966  59.711  < 2e-16 ***
## multce_w     0.12169    0.01753   6.943 4.97e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.026 on 2297 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.02056,    Adjusted R-squared:  0.02013 
## F-statistic: 48.21 on 1 and 2297 DF,  p-value: 4.97e-12

El diseño de investigación incluyó la aplicación de dos cuestionarios a cada uno de los participantes en el estudio: uno antes de la realización de los talleres y otro cuando el programa concluyó.2 Esto quiere decir que se dispone de dos observaciones por sujeto. Cada variable es medidas en dos ocasiones, antes del inicio de los talleres y después de la conclusión del programa. La variable waveint asigna un valor de 0 a la primera medición y de 1 a la segunda. La variable id es un número que identifica a los sujetos del estudio.

Para ver la estructura de los datos panel, revisemos las primeras 25 observaciones. Las columnas uno a la cuatro contienen las variables antes comentadas. Se puede ver que a cada sujeto le corresponden dos filas, una para cada medición. Así, por ejemplo, el individuo identificado con el número 10107 obtuvo una puntuación de 1.50 en la primera medición del índice de tolerancia y de 2.25 en la segunda. El número de talleres en la primera medición para todos los sujetos es de cero, y este sujeto en particular había asistido a uno al concluir el programa.

head(ef[, c(1:4)], n = 20)
##       id waveint tol_w multce_w
## 1  10101       0  1.00        0
## 2  10101       1  1.00        3
## 3  10103       0  1.00        0
## 4  10103       1  1.00        3
## 5  10104       0  1.00        0
## 6  10104       1  1.00        1
## 7  10105       0  4.00        0
## 8  10105       1  4.00        1
## 9  10107       0  1.50        0
## 10 10107       1  2.25        1
## 11 10121       0  3.50        0
## 12 10121       1  3.00        3
## 13 10122       0  2.50        0
## 14 10122       1  1.00        4
## 15 10131       0  1.00        0
## 16 10131       1  1.00        0
## 17 10133       0  1.00        0
## 18 10133       1  1.50        1
## 19 10134       0  1.00        0
## 20 10134       1  1.25        1

Ahora vamos a descontar el efecto de la unidad con un modelo de efectos fijos, usando la función plm(). En el vector del argumento index se debe especificar, primero, la variable con la que se identifica a las unidades y, después, la que identifica el momento de la medición. Con el argumento model = "within" se toman los efectos fijos por unidad. La función plm() implementa la especificación de diferencias con respecto a la media de la unidad (ecuación ), con el ajuste apropiado en los grados de libertad.

library("plm")
fixed <- plm(tol_w ~ multce_w,
             data = ef, 
             index = c("id", "waveint"),
             model = "within")
summary(fixed)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = tol_w ~ multce_w, data = ef, model = "within", 
##     index = c("id", "waveint"))
## 
## Unbalanced Panel: n = 2301, T = 1-2, N = 4595
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -1.59046 -0.38631  0.00000  0.38631  1.59046 
## 
## Coefficients:
##          Estimate Std. Error t-value Pr(>|t|)   
## multce_w 0.045229   0.016399  2.7581 0.005861 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2033.6
## Residual Sum of Squares: 2026.9
## R-Squared:      0.0033065
## Adj. R-Squared: -0.99686
## F-statistic: 7.60689 on 1 and 2293 DF, p-value: 0.0058607

Al descontar las características invariantes de los sujetos, el coeficiente se reduce considerablemente. Un ajuste adicional se consigue incorporando la tendencia temporal. Para ello, se añade el argumento effect = "twoways".

fixed.t <- plm(tol_w ~ multce_w,
            data = ef,
            index = c("id", "waveint"),
            model = "within",
            effect = "twoways")
summary(fixed.t)
## Twoways effects Within Model
## 
## Call:
## plm(formula = tol_w ~ multce_w, data = ef, effect = "twoways", 
##     model = "within", index = c("id", "waveint"))
## 
## Unbalanced Panel: n = 2301, T = 1-2, N = 4595
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -1.65763 -0.39451  0.00000  0.39451  1.65763 
## 
## Coefficients:
##          Estimate Std. Error t-value  Pr(>|t|)    
## multce_w 0.131555   0.022593  5.8229 6.594e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2029.9
## Residual Sum of Squares: 2000.3
## R-Squared:      0.014578
## Adj. R-Squared: -0.97514
## F-statistic: 33.9067 on 1 and 2292 DF, p-value: 6.594e-09

Ahora el coeficiente de interés incrementó su valor. Se puede mostrar fácilmente que esto es debido a que la tendencia temporal del índice de tolerancia es decreciente, mientras que el valor del número de talleres es, por fuerza, creciente.

Un último ajuste en el marco del modelo de efectos fijos se consigue con la inclusión de características que varían en el tiempo. Finkel y Smith controlan por exposición a los medios (med_rcdw), interés en la política (med_rcdw), pertenencia a grupos (group_rcdw) y discusión política (polpart1_w). Todas las variables son tratadas por ellos como numéricas. Por este motivo se usa as.numeric() en la última de estas variables, un factor con tres niveles.

fixed.c <- plm(tol_w ~ multce_w + med_rcdw + int_rcdw + group_rcdw +
                 as.numeric(polpart1_w),
            data = ef,
            index = c("id", "waveint"),
            model = "within",
            effect = "twoways")
summary(fixed.c)
## Twoways effects Within Model
## 
## Call:
## plm(formula = tol_w ~ multce_w + med_rcdw + int_rcdw + group_rcdw + 
##     as.numeric(polpart1_w), data = ef, effect = "twoways", model = "within", 
##     index = c("id", "waveint"))
## 
## Unbalanced Panel: n = 2301, T = 1-2, N = 4586
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -1.67590 -0.38497  0.00000  0.38497  1.67590 
## 
## Coefficients:
##                        Estimate Std. Error t-value  Pr(>|t|)    
## multce_w               0.129335   0.022757  5.6832 1.491e-08 ***
## med_rcdw               0.230858   0.094846  2.4340   0.01501 *  
## int_rcdw               0.036416   0.090634  0.4018   0.68788    
## group_rcdw             0.102548   0.115846  0.8852   0.37613    
## as.numeric(polpart1_w) 0.027229   0.031462  0.8654   0.38689    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2026.3
## Residual Sum of Squares: 1988
## R-Squared:      0.018865
## Adj. R-Squared: -0.97389
## F-statistic: 8.76392 on 5 and 2279 DF, p-value: 3.02e-08

Asumiendo que se ha controlado por todas las características relevantes con variación temporal, la estimación del efecto del programa de educación cívica está dada por el coeficiente 0.129.

Por último, se debe dar cuenta del hecho de que en un panel, compuesto de dos o más observaciones de una misma unidad, cabe esperar que exista autocorrelación en los residuos. Es más apropiado el uso de errores estándar robustos para las unidades (en clusters). La función coeftest() genera automáticamente los errores estándar cuando se incluye en su argumento un objeto de plm.

library(lmtest)
fixed.r <- coeftest(fixed.c, vcov = vcovHC, type = "HC1")
fixed.r
## 
## t test of coefficients:
## 
##                        Estimate Std. Error t value  Pr(>|t|)    
## multce_w               0.129335   0.023874  5.4174 6.681e-08 ***
## med_rcdw               0.230858   0.092792  2.4879   0.01292 *  
## int_rcdw               0.036416   0.089381  0.4074   0.68374    
## group_rcdw             0.102548   0.109704  0.9348   0.35000    
## as.numeric(polpart1_w) 0.027229   0.032153  0.8469   0.39717    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dado el objetivo del estudio, lo relevante de los resultados es que la significancia del coeficiente para multce_w se sostiene aún después de corregir los errores estándar.

La conclusión, asumiendo que se ha controlado por las características relevantes que varían individualmente en el tiempo, cada taller adicional produce un incremento de 0.129 unidades en el índice de tolerancia.

Un comentario final

En lo personal, encuentro innecesario el supuesto de linealidad entre la variable de tratamiento y la de resultado. Podríamos indagar sobre el efecto que tiene la asistencia a una cantidad dada de talleres en comparación con no asistir a ninguno. Para ello, sólo tenemos que tratar la variable de tratamiento como variable categórica, generando una variable indicadora para cada uno de sus valores y dejando como categoría de referencia la no asistencia a ningún taller.

fixed.c1 <- plm(tol_w ~ factor(multce_w) + med_rcdw + int_rcdw + group_rcdw +
                 as.numeric(polpart1_w),
               data = ef,
               index = c("id", "waveint"),
               model = "within",
               effect = "twoways")
fixed.r1 <- coeftest(fixed.c1, vcov = vcovHC, type = "HC1")
fixed.r1
## 
## t test of coefficients:
## 
##                        Estimate Std. Error t value  Pr(>|t|)    
## factor(multce_w)1      0.107579   0.063235  1.7012 0.0890337 .  
## factor(multce_w)2      0.361737   0.100660  3.5936 0.0003330 ***
## factor(multce_w)3      0.435529   0.107031  4.0692 4.879e-05 ***
## factor(multce_w)4      0.429891   0.118769  3.6196 0.0003015 ***
## med_rcdw               0.230167   0.092645  2.4844 0.0130486 *  
## int_rcdw               0.033331   0.089854  0.3709 0.7107112    
## group_rcdw             0.107544   0.109799  0.9795 0.3274565    
## as.numeric(polpart1_w) 0.029004   0.032272  0.8987 0.3688905    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cada coeficiente, bajo los supuestos de identificación del modelo, nos habla del cambio medio en el valor del índice de tolerancia, generado, respectivamente, por la asistencia a uno, dos, tres o cuatro o más talleres, en comparación con no haber asistido a ninguno. Para facilitar la interpretación de los resultados, generamos una gráfica de los coeficientes con sus respectivos intervalos de confianza de 95% (línea delgada) y 90% (línea gruesa).

library(ggplot2)
library(ggstance)
library(jtools)
library(wesanderson)
plot_summs(fixed.r1, ci_level = .95, inner_ci_level = .90,
           coefs = c("1 taller" = "factor(multce_w)1", "2 talleres" = "factor(multce_w)2",
                     "3 talleres" = "factor(multce_w)3", "4 o más" = "factor(multce_w)4"),
           colors = wes_palette("Darjeeling2")[2]) + labs(x = "Estimación", y = "")

A un nivel del significancia del 95%, con prueba de dos colas, asistir a un solo taller no parece hacer una diferencia con respecto a no asistir a ninguno.3 El efecto es claramente detectable sólo a partir de la asistencia a dos talleres. Ahora bien, los efectos correspondientes a dos, tres y cuatro o más talleres, si bien significativos, presentan valores muy similares entre sí.

Para evaluar esta impresión de manera más sistemática, realizamos un test de la hipótesis nula de igualdad de estos tres coeficientes. De acuerdo con el resultado, no podemos descartar que el efecto sea el mismo para dos talleres, tres talleres y cuatro o más.

VC <- vcovHC(fixed.c1, type = "HC1") # Matriz de covarianzas para errores estándar robustos
library(car)
linearHypothesis(fixed.c1, c("factor(multce_w)2 = factor(multce_w)3",
                              "factor(multce_w)3 = factor(multce_w)4"),
                              test = "F", vcov = VC)
## Linear hypothesis test
## 
## Hypothesis:
## factor(multce_w)2 - factor(multce_w)3 = 0
## factor(multce_w)3 - factor(multce_w)4 = 0
## 
## Model 1: restricted model
## Model 2: tol_w ~ factor(multce_w) + med_rcdw + int_rcdw + group_rcdw + 
##     as.numeric(polpart1_w)
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df      F Pr(>F)
## 1   2278                 
## 2   2276  2 0.1912  0.826

Para saber más…

Algunas introducciones del modelo de efectos fijos son Allison (2009), Brüderl y Ludwig (2015) y Halaby (2004). La referencia obligada es Wooldridge (2010). En este texto se sigue de cerca a Angrist y Pischke (2009). Croissant y Millo (2008) exponen los usos del paquete plm y en Hanck et al. (2019) hay algunos ejemplos aplicados.

Cómo citar

Los Breviarios Digitales tienen fines pedagógicos y de divulgación del conocimiento. Favor de referir este breviario de la siguiente forma:

Salazar-Elena, Rodrigo (2021): “Efectos fijos”. Breviarios Digitales. FLACSO. Disponible por Internet en https://sites.google.com/flacso.edu.mx/labdem/breviarios-digitales?authuser=0

Referencias

Allison, Paul D. 2009. Fixed effects regression models. Thousand Oaks, CA: Sage.
Angrist, Joshua D., y Jörn-Steffen Pischke. 2009. Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton: Princeton University Press.
Brüderl, Josef, y Volker Ludwig. 2015. «Fixed-effects panel regression». En The SAGE Handbook of Regression Analysis and Causal Inference, editado por Henning Best y Christof Wolf, 7-29. Thousand Oaks: Sage.
Croissant, Yves, y Giovanni Millo. 2008. «Panel Data Econometrics in R: The plm Package». Journal of Statistical Software 27 (2): 1-43. https://www.jstatsoft.org/v027/i02.
Finkel, Steven E., y Amy Erica Smith. 2011. «Civic education, political discussion, and the social transmission of democratic knowledge and values in a new democracy: Kenya 2002». American Journal of Political Science 55 (2): 417-35.
Halaby, Charles N. 2004. «Panel models in sociological research: Theory into Practice». Annual Review of Sociology 30: 507-44.
Hanck, Christoph, Martin Arnold, Alexander Gerber, y Martin Schmelzer. 2019. Introduction to Econometrics with R. Essen: Department of Business Administration; Economics, University of Duisburg-Essen. https://www.econometrics-with-r.org.
Wooldridge, Jeffrey M. 2010. Econometric Analysis of Cross Section and Panel Data. 2.ª ed. Cambridge, MA: The MIT Press.

  1. La base de datos está disponible en el repositorio de los breviarios. Se replica el análisis que es reportado en la quinta columna del cuadro 2 (p. 427).↩︎

  2. El cuestionario posterior a la intervención se realizó en dos períodos distintos. Esa diferencia no es tomada en cuenta para efectos del análisis que aquí se presenta.↩︎

  3. Sin embargo, al mismo nivel de significancia, para la hipótesis alternativa de efecto positivo (una cola), se rechaza la hipótesis de efecto nulo para la asistencia a un sólo taller.↩︎