The `panelr`

package contributes two categories of things:

- A
`panel_data`

object and some tools to create/manipulate them. - A series of regression modeling functions for panel data.

`panel_data`

framesCheck out the other vignette for a lot of detail on how to take your raw data
and reshape it into a `panel_data`

format. Here's a short version, using some
example data provided by this package.

```
library(panelr)
data("teen_poverty")
teen_poverty
```

```
# A tibble: 1,151 x 28
id pov1 mother1 spouse1 inschool1 hours1 pov2 mother2 spouse2 inschool2
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 22 1 0 0 1 21 0 0 0 1
2 75 0 0 0 1 8 0 0 0 1
3 92 0 0 0 1 30 0 0 0 1
4 96 0 0 0 0 19 1 1 0 0
5 141 0 0 0 1 0 0 0 0 1
6 161 0 0 0 1 0 0 0 0 1
7 220 0 0 0 1 6 0 0 0 1
8 229 0 0 0 1 0 1 0 0 1
9 236 0 0 0 1 0 0 0 0 1
10 240 0 0 0 1 18 1 0 0 1
# … with 1,141 more rows, and 18 more variables: hours2 <dbl>, pov3 <dbl>,
# mother3 <dbl>, spouse3 <dbl>, inschool3 <dbl>, hours3 <dbl>, pov4 <dbl>,
# mother4 <dbl>, spouse4 <dbl>, inschool4 <dbl>, hours4 <dbl>, age <dbl>,
# black <dbl>, pov5 <dbl>, mother5 <dbl>, spouse5 <dbl>, inschool5 <dbl>,
# hours5 <dbl>
```

These data come from a subset of young women surveyed as part of the
National Longitudinal Survey of Youth starting in 1979. The `teen_poverty`

data come in “wide” format, meaning there is one row per respondent and
each of the repeated measures is in a separate column for each wave.

We need to convert this to “long” format, in which you have one row for each
respondent in each wave of the 5-wave survey. We'll use `long_panel()`

for
that.

```
teen <- long_panel(teen_poverty, begin = 1, end = 5, label_location = "end")
teen
```

```
# Panel data: 5,755 x 9
# entities: id [1151]
# wave variable: wave [1, 2, 3, ... (5 waves)]
id wave age black pov mother spouse inschool hours
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 22 1 16 0 1 0 0 1 21
2 22 2 16 0 0 0 0 1 15
3 22 3 16 0 0 0 0 1 3
4 22 4 16 0 0 0 0 1 0
5 22 5 16 0 0 0 0 1 0
6 75 1 17 0 0 0 0 1 8
7 75 2 17 0 0 0 0 1 0
8 75 3 17 0 0 0 0 1 0
9 75 4 17 0 0 0 0 1 4
10 75 5 17 0 1 0 0 1 0
# … with 5,745 more rows
```

Now we have a `panel_data`

object! It is a special version of a `tibble`

,
which is itself a special kind of `data.frame`

. `panel_data`

objects
work very hard to make sure you never accidentally drop the variables that
are the identifiers for each respondent and the indicators for which wave the
row corresponds to. `panel_data`

objects also try to stay in order by ID and
wave.

Note that if your raw data are already in long format, you can use the
`panel_data()`

function to convert them to `panel_data`

format.

```
data("WageData")
wages <- panel_data(WageData, id = id, wave = t)
```

`panel_data`

frames are designed to work with `tidyverse`

packages,
particularly `dplyr`

. When used inside `mutate()`

, functions like `lag()`

work
properly by taking the previous value for the specific respondent. If you ever
need to do something that is easier to do with a “regular” data frame, you can
just use the `unpanel()`

function to convert the `panel_data`

frame back to
normal.

The original motivation to create this package was to automate the process of fitting “within-between” models, sometimes called “between-within” or “hybrid” models (see Allison, 2009; Bell & Jones, 2015). These combine the benefits of what econometricians call “fixed effects” models — robustness to time-invariant confounding chief among them — as well as what they call “random effects” models, which allow the inclusion of time-invariant coefficients. Within-between models include coefficients that are identical to the fixed effects equivalent, but the flexibility to also include the random effects and other time-invariant predictors (this was noticed by Mundlak, 1978). They are fit via multilevel models which allow for some other nice possibilities like inclusion of random slopes and generalized linear model specifications.

From here, I'll give a somewhat technical description of these models. If you just want to look at how to estimate them in R, skip ahead to the next mini-section.

Note that fixed effects models can be fit using individual demeaning. That is, you can subtract the entity's own mean for each predictor and the dependent variable and fit a model via OLS that is equivalent to the so-called least squares dummy variable approach (in which dummy variables for every entity ID are included as predictors).

Let's get a bit more technical. We have entities \(i = 1, …, n\) who are measured at times \(t = 1, …, T\). We have as our dependent variable \(y_{it}\), the variable \(y\) for individual \(i\) at time \(t\). We have predictors that vary over time \(x_{it}\), variables that do not vary over time \(z_i\), and variables we did not measure that do not vary over time \(\alpha_i\) as well as random error \(\epsilon_{it}\)

The fixed effects model, then, looks like this:

\[ y_{it} = \mu_t + \beta_1x_{it} + \gamma z_i + \alpha_i + \epsilon_{it} \]

Although \(\alpha_i\) is not observed, it can be estimated by including a dummy variable for each \(i\). The \(\gamma\) is undefined because the \(z_i\) are perfectly collinear with the \(\alpha_i\) dummy variables.

The individual-mean-centered version of the fixed effects models is based on calculating a mean of \(y\) and \(x\) for each \(i\) — so \(\bar{y_i}\) and \(\bar{x_i}\) and subtracting it from each \(y_{it}\) and \(x_{it}\). The model can be expressed like this, including \(\bar{z_i}\) and \(\bar{\alpha_i}\) for demonstration:

\[ y_{it} - \bar{y_i} = \mu_t + \beta_1(x_{it} - \bar{x_i}) + (z_i - \bar{z_i} = 0) + (\alpha_i - \bar{\alpha_i} = 0) + (\epsilon_{it} - \bar{\epsilon_i}) \]

By de-meaning everything, all the time-invariant variables drop out:

\[ y_{it} - \bar{y_i} = \mu_t + \beta_1(x_{it} - \bar{x_i}) + (\epsilon_{it} - \bar{\epsilon_i}) \]

This is often called the “within” estimator. You can take these de-meaned variables and fit an OLS regression and get valid estimates (with some adjustments to the standard errors).

You can also do something slightly different and get the same results with multilevel models. Take this, for example:

\[ y_{it} = \beta_{0i} + \beta_1(x_{it} - \bar{x_i}) + (\epsilon_{it} - \bar{\epsilon_i}) \]

Where \(\beta_{0i}\) is a random intercept estimated for each \(i\). This is equivalent to subtracting \(\bar{y_i}\) in terms of the estimation of \(\beta_1\). But in the multilevel modeling framework, we can include those time-invariant \(z_i\) as well. Conceptually, they are basically being included in a model predicting \(\beta_{0i}\):

\[ \beta_{0i} = \beta_0 + \gamma z_i + u_{0i} \]

Where \(u_{0i}\) is the random error of the model predicting \(\beta_{0i}\).

In fact, we can include the \(\bar{x_i}\) in our multilevel model as well and they are used just like the \(z_i\):

\[ \beta_{0i} = \beta_0 + \beta_2 \bar{x_i} + \gamma z_i + u_{0i} \]

Now we can substitute into the previous multilevel equation and we have our within-between model:

$$
y*{it} = \beta*{0} + \beta*1(x*{it} - \bar{x_i})

- \beta
*2 \bar{x_i} + \gamma z_i + u*{0i} + \epsilon_{it} $$

The \(\beta_1\) has the same interpretation as in the fixed effects model, these are the effects of within-entity deviations of \(x\) on within-entity deviations of \(y\). The \(\beta_2\) is basically predicting the \(\bar{y_i}\), however, so these coefficients are helpful for predicting differences in mean levels across entities. The same is true for the \(z_i\).

A similar model that I call the “contextual” model because this is how it is often interpreted (see, e.g., Raudenbush & Bryk, 2002). Here we do not demean the \(x_i\):

$$
y*{it} = \beta*{0} + \beta*1 x*{it}

- \beta
*2 \bar{x_i} + \gamma z_i + u*{0i} + \epsilon_{it} $$

Believe it or not, the \(\beta_1\) is unchanged in this model; it is the \(\beta_2\)
that changes. The interpretation of \(\beta_2\) becomes a the *difference*
between the within- and between-entities effects. A significant coefficient for
\(\beta_2\) means significant differences between the within- and between-entity
effects. For those who are familiar, this is like a variable-by-variable Hausman
test. Substantively, \(\beta_2\) is often interpreted as a *contextual* effect.

From this framework, we can do cross-level interactions, random slopes, generalized linear models, and all kinds of interesting stuff.

In the fixed effects framework, it is generally considered wrong to
operationalize an interaction between two time-varying variables (let's call
them \(w\) and \(x\)) by taking
the product of their individual-demeaned forms. That is, you are **not**
supposed to generate the interaction term \(xw_{it}\) by doing this:

\[ xw_{it} = (x_{it} - \bar{x_{i}}) \times (w_{it} - \bar{w_i}) \]

Instead, the conventional wisdom goes, you should first take the product of the observed variables and subtract the individual-level mean of that product, like so:

\[ xw_{it} = x_{it}w_{it} - \overline{xw}_i \]

Where \(\overline{xw}_i\) can also be expressed as \(\frac{\sum_{t=1}^{T_i}{x_{it}w_{it}}}{T_i}\), the sum of all products for each \(i\) divided by the number of time points for each \(i\), \(T_i\).

Giesselmann and Schmidt-Catran (2020)
show that this conventional method for generating \(xw_{it}\) does not have the
unbiasedness that the individual terms do. I'll leave it to them to explain
why exactly this is, but the solution is to start with the first, wrong
version of \(xw_{it}\), which I'll call \(xw_{it}^*\), and subtract *its* mean too:

\[ xw_{it}^* = (x_{it} - \bar{x_{i}}) \times (w_{it} - \bar{w_i}) \ xw_{it} = xw_{it}^* - \overline{xw_i^*} \]

I call this the “double-demeaning” approach to interactions, in contrast to
the one-time demeaning in the conventional approach. By default, `wbm()`

calculates interactions via the double-demeaning method. You can change this
via the `interaction.style`

argument if you need your results to match other
software.

The workhorse function for within-between models is `wbm()`

, which is built on
top of `lme4`

's `lmerMod()`

and `glmerMod()`

. It is not so hard to understand
how to treat your data to estimate within-between models, but the programming
can be a challenge to those who aren't skilled with R (or whatever else they
might use) and is error-prone in any case.

The main thing to know in order to use `wbm()`

is how the model formula
works, because it's a little different from your typical regression model.
It is split into up to 3 parts, each for a different kind of variable. Each
part is separated by a `|`

. The pattern is like this:

```
dependent ~ time_varying | time_invariant | cross_lev_interactions + (random_slopes | id)
```

So you start with your dependent variable on the left-hand side like normal and
then what comes next are variables that vary over time. You will only get
within-entity estimates for these variables. Next are time-invariant variables;
the between-entity terms for the time-varying variables are added automatically
so no need to try to include them here. Finally, in the third part you can
specify cross-level interactions (i.e., within-entity by
between-entity/time-invariant) as well as additional random effects terms
using the `lme4`

-style syntax. By default, `(1 | id)`

(or whatever the ID
variable is) is added internally for a random intercept so you do not need to
include it yourself.

Let's walk through an example with the `wages`

data we looked at briefly
earlier. We'll predict the logarithm of wages (`lwage`

) using weeks worked
(`wks`

), union membership (`union`

), marital status (`ms`

),
blue (vs. white) collar job status (`occ`

), black race (`blk`

), and
female sex (`fem`

).

```
model <- wbm(lwage ~ wks + union + ms + occ | blk + fem, data = wages)
summary(model)
```

```
MODEL INFO:
Entities: 595
Time periods: 1-7
Dependent variable: lwage
Model type: Linear mixed effects
Specification: within-between
MODEL FIT:
AIC = 2036.78, BIC = 2119.13
Pseudo-R² (fixed effects) = 0.27
Pseudo-R² (total) = 0.69
Entity ICC = 0.57
WITHIN EFFECTS:
----------------------------------------------------
Est. S.E. t val. d.f. p
----------- ------- ------ -------- --------- ------
wks 0.00 0.00 1.06 3566.00 0.29
union 0.06 0.03 2.53 3566.00 0.01
ms -0.08 0.03 -2.57 3566.00 0.01
occ -0.08 0.02 -3.32 3566.00 0.00
----------------------------------------------------
BETWEEN EFFECTS:
----------------------------------------------------------
Est. S.E. t val. d.f. p
------------------ ------- ------ -------- -------- ------
(Intercept) 6.30 0.20 30.85 588.00 0.00
imean(wks) 0.01 0.00 2.25 588.00 0.02
imean(union) 0.15 0.03 4.67 588.00 0.00
imean(ms) 0.17 0.05 3.07 588.00 0.00
imean(occ) -0.41 0.03 -13.31 588.00 0.00
blk -0.15 0.05 -2.81 588.00 0.01
fem -0.32 0.06 -4.96 588.00 0.00
----------------------------------------------------------
p values calculated using Satterthwaite d.f.
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
id (Intercept) 0.2992
Residual 0.2589
------------------------------------
```

As you can see, the output distinguishes within- and between-entity effects.
When you see `imean()`

around a variable, that is the between-entity effect
represented as the individual mean.

Here, we see there seems to be a wage penalty for switching from white collar
to blue collar work (`occ`

) and although married people earn more (`imean(ms)`

),
just becoming married (`ms`

) coincides with a drop in earnings. We also see a
boost in earnings from joining a union (`union`

).

Maybe we think the timing of the marriage effect is off and the true effect
occurs the time period after a person becomes married. We can ask for the
lagged effect using `lag()`

.

```
model <- wbm(lwage ~ wks + union + lag(ms) + occ | blk + fem, data = wages)
summary(model)
```

```
MODEL INFO:
Entities: 595
Time periods: 2-7
Dependent variable: lwage
Model type: Linear mixed effects
Specification: within-between
MODEL FIT:
AIC = 1247.06, BIC = 1327.41
Pseudo-R² (fixed effects) = 0.28
Pseudo-R² (total) = 0.74
Entity ICC = 0.64
WITHIN EFFECTS:
------------------------------------------------------
Est. S.E. t val. d.f. p
------------- ------- ------ -------- --------- ------
wks -0.00 0.00 -1.43 2999.92 0.15
union 0.05 0.03 1.82 2991.78 0.07
lag(ms) -0.04 0.03 -1.19 2971.04 0.23
occ -0.06 0.02 -2.56 2992.05 0.01
------------------------------------------------------
BETWEEN EFFECTS:
------------------------------------------------------------
Est. S.E. t val. d.f. p
-------------------- ------- ------ -------- -------- ------
(Intercept) 6.36 0.21 30.51 588.00 0.00
imean(wks) 0.01 0.00 2.25 588.00 0.02
imean(union) 0.15 0.03 4.48 588.00 0.00
imean(lag(ms)) 0.15 0.06 2.80 588.01 0.01
imean(occ) -0.41 0.03 -13.16 588.02 0.00
blk -0.15 0.05 -2.92 587.99 0.00
fem -0.33 0.06 -5.17 588.01 0.00
------------------------------------------------------------
p values calculated using Satterthwaite d.f.
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
id (Intercept) 0.3068
Residual 0.2325
------------------------------------
```

Well that doesn't change the direction of the estimate, but it also moved it sufficiently close to 0 that we can't say much about it one way or another.

Keep in mind that you do not have to stick to linear models. Using the `family`

argument (just like `glm()`

), you can estimate logit (`family = binomial`

),
probit (`family = binomal(link = "probit")`

), poisson (`family = poisson`

), or
other model families and links as needed.

Now maybe we want to include an effect of time since wages tend to go up for
everyone, on average, over time. We can just include the time variable in the
formula or set `use.wave`

to `TRUE`

.

```
model <- wbm(lwage ~ wks + union + ms + occ | blk + fem, data = wages, use.wave = TRUE)
summary(model)
```

```
MODEL INFO:
Entities: 595
Time periods: 1-7
Dependent variable: lwage
Model type: Linear mixed effects
Specification: within-between
MODEL FIT:
AIC = -1688.87, BIC = -1600.19
Pseudo-R² (fixed effects) = 0.44
Pseudo-R² (total) = 0.89
Entity ICC = 0.8
WITHIN EFFECTS:
----------------------------------------------------
Est. S.E. t val. d.f. p
----------- ------- ------ -------- --------- ------
wks 0.00 0.00 1.91 3565.00 0.06
union 0.03 0.02 2.29 3565.00 0.02
ms -0.03 0.02 -1.71 3565.00 0.09
occ -0.03 0.01 -1.81 3565.00 0.07
----------------------------------------------------
BETWEEN EFFECTS:
-----------------------------------------------------------
Est. S.E. t val. d.f. p
------------------ ------- ------ -------- --------- ------
(Intercept) 5.91 0.20 28.95 588.64 0.00
imean(wks) 0.01 0.00 2.25 588.00 0.02
imean(union) 0.15 0.03 4.67 588.00 0.00
imean(ms) 0.17 0.05 3.07 588.00 0.00
imean(occ) -0.41 0.03 -13.31 588.00 0.00
blk -0.15 0.05 -2.81 588.00 0.01
fem -0.32 0.06 -4.96 588.00 0.00
t 0.10 0.00 81.29 3565.00 0.00
-----------------------------------------------------------
p values calculated using Satterthwaite d.f.
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
id (Intercept) 0.3094
Residual 0.1533
------------------------------------
```

Including `t`

wipes out some of those previously observed effects. Believe it
or not, we just fit a growth curve model!

Now, we might think people have different trajectories. We can include that as a random slope, which will go in the third part of the formula.

```
model <- wbm(lwage ~ wks + union + ms + occ | blk + fem | (t | id), use.wave = TRUE, data = wages)
summary(model)
```

```
MODEL INFO:
Entities: 595
Time periods: 1-7
Dependent variable: lwage
Model type: Linear mixed effects
Specification: within-between
MODEL FIT:
AIC = -2064.42, BIC = -1963.07
Pseudo-R² (fixed effects) = 0.43
Pseudo-R² (total) = 0.92
Entity ICC = 0.84
WITHIN EFFECTS:
----------------------------------------------------
Est. S.E. t val. d.f. p
----------- ------- ------ -------- --------- ------
wks 0.00 0.00 1.47 3498.40 0.14
union 0.02 0.01 1.61 3561.26 0.11
ms -0.04 0.02 -1.97 3416.40 0.05
occ -0.02 0.01 -1.42 3563.73 0.16
----------------------------------------------------
BETWEEN EFFECTS:
----------------------------------------------------------
Est. S.E. t val. d.f. p
------------------ ------- ------ -------- -------- ------
(Intercept) 5.93 0.20 29.69 588.52 0.00
imean(wks) 0.01 0.00 2.12 587.97 0.03
imean(union) 0.16 0.03 4.90 587.96 0.00
imean(ms) 0.17 0.05 3.25 588.06 0.00
imean(occ) -0.39 0.03 -13.16 588.01 0.00
blk -0.13 0.05 -2.49 588.00 0.01
fem -0.31 0.06 -4.90 588.06 0.00
t 0.10 0.00 54.67 594.66 0.00
----------------------------------------------------------
p values calculated using Satterthwaite d.f.
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
id (Intercept) 0.3057
id t 0.03499
Residual 0.1334
------------------------------------
```

And now we have a latent growth curve model. The general effect on the other
coefficients is more uncertainty and attenuated estimates. It's worth
keeping in mind that it is sometimes wrong to use a growth curve model like this
if you think the variables in your model *cause* the time trend; if you think
wages are going up because more people are moving into white collar work, then
including the growth curve will make it harder for you to see the true effect
of `occ`

.

By default, `wbm()`

does as the name suggests. But if you'd rather have the
contextual model described earlier, in which the means are not subtracted from
the time varying variables, that's an option too.

```
model <- wbm(lwage ~ wks + union + ms + occ | blk + fem, data = wages, model = "contextual")
summary(model)
```

```
MODEL INFO:
Entities: 595
Time periods: 1-7
Dependent variable: lwage
Model type: Linear mixed effects
Specification: contextual
MODEL FIT:
AIC = 2036.78, BIC = 2119.13
Pseudo-R² (fixed effects) = 0.27
Pseudo-R² (total) = 0.69
Entity ICC = 0.57
WITHIN EFFECTS:
----------------------------------------------------
Est. S.E. t val. d.f. p
----------- ------- ------ -------- --------- ------
wks 0.00 0.00 1.06 3566.00 0.29
union 0.06 0.03 2.53 3566.00 0.01
ms -0.08 0.03 -2.57 3566.00 0.01
occ -0.08 0.02 -3.32 3566.00 0.00
----------------------------------------------------
CONTEXTUAL EFFECTS:
-----------------------------------------------------------
Est. S.E. t val. d.f. p
------------------ ------- ------ -------- --------- ------
(Intercept) 6.30 0.20 30.85 588.00 0.00
imean(wks) 0.01 0.00 1.93 660.18 0.05
imean(union) 0.09 0.04 2.15 1411.72 0.03
imean(ms) 0.25 0.06 3.95 1047.31 0.00
imean(occ) -0.33 0.04 -8.55 1401.80 0.00
blk -0.15 0.05 -2.81 588.00 0.01
fem -0.32 0.06 -4.96 588.00 0.00
-----------------------------------------------------------
p values calculated using Satterthwaite d.f.
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
id (Intercept) 0.2992
Residual 0.2589
------------------------------------
```

Now the individual means have a new interpretation as the difference in effect compared to the within-entity estimates.

If you don't want to use any of the time-invariant variables, you can also just ask for the “within” estimator:

```
model <- wbm(lwage ~ wks + union + ms + occ, data = wages, model = "within")
summary(model)
```

```
MODEL INFO:
Entities: 595
Time periods: 1-7
Dependent variable: lwage
Model type: Linear mixed effects
Specification: within
MODEL FIT:
AIC = 2266.01, BIC = 2310.35
Pseudo-R² (fixed effects) = 0
Pseudo-R² (total) = 0.69
Entity ICC = 0.69
----------------------------------------------------------
Est. S.E. t val. d.f. p
----------------- ------- ------ -------- --------- ------
(Intercept) 6.68 0.02 413.08 594.00 0.00
wks 0.00 0.00 1.06 3566.00 0.29
union 0.06 0.03 2.53 3566.00 0.01
ms -0.08 0.03 -2.57 3566.00 0.01
occ -0.08 0.02 -3.32 3566.00 0.00
----------------------------------------------------------
p values calculated using Satterthwaite d.f.
RANDOM EFFECTS:
------------------------------------
Group Parameter Std. Dev.
---------- ------------- -----------
id (Intercept) 0.3819
Residual 0.2589
------------------------------------
```

This can help declutter your output when you really just don't care about the between-subjects effects.

You don't have to estimate these models using multilevel models and in fact you may get better inferences by avoiding some of the assumptions inherent to multilevel modeling (see McNeish, 2019). You can use the semiparametric generalized estimating equations (GEE) approach to estimation, with the main tradeoff being that you can no longer use random slopes or anything like that. But if you only care about the average effects across all entities, GEE can be a better approach that doesn't require you to be right about the distribution of effects and several other assumptions.

`wbgee()`

builds on `geeglm()`

from the `geepack`

package and works just like
`wbm()`

.

```
model <- wbgee(lwage ~ wks + union + ms + occ | blk + fem, data = wages)
summary(model)
```

```
MODEL INFO:
Entities: 595
Time periods: 1-7
Dependent variable: lwage
Model type: Linear GEE
Variance: ar1 (alpha = 0.79)
Specification: within-between
MODEL FIT:
QIC = 672.19, QICu = 669.6, CIC = 12.29
WITHIN EFFECTS:
------------------------------------------
Est. S.E. z val. p
----------- ------- ------ -------- ------
wks 0.00 0.00 0.07 0.94
union 0.03 0.02 1.32 0.19
ms -0.08 0.03 -2.78 0.01
occ -0.03 0.02 -1.64 0.10
------------------------------------------
BETWEEN EFFECTS:
-------------------------------------------------
Est. S.E. z val. p
------------------ ------- ------ -------- ------
(Intercept) 6.29 0.22 28.48 0.00
imean(wks) 0.01 0.00 2.09 0.04
imean(union) 0.16 0.03 4.74 0.00
imean(ms) 0.17 0.06 3.01 0.00
imean(occ) -0.41 0.03 -12.95 0.00
blk -0.14 0.05 -2.90 0.00
fem -0.30 0.06 -4.89 0.00
-------------------------------------------------
```

This gives us more conservative estimates, in general. Note that by default,
`wbgee()`

uses an AR-1 working error correlation structure in estimation.
This makes sense in general but at times it may make sense to use
“exchangeable” as the argument to `cor.str`

which assumes all within-entity
correlations are equal regardless of time lag. Other options include
“unstructured”, which can be very computationally intensive, and “independence,”
assuming no correlation within entities.

Like `wbm()`

, you can do generalized linear models via the `family`

argument.
It is for these generalized linear models that GEEs are likely to stand out the
most in terms of added benefit above and beyond the multilevel models, although
this is not a well-tested question to my knowledge.

Sometimes, theory may suggest that increases in a variable have a different effect than decreases in a variable. For instance, getting married and getting divorced are probably not equivalent (in the sense that one is the exact opposite of the other) in their effects on other outcomes. Allison (2019) described a method for estimating models with asymmetric effects based on first differences.

First, you take first differences:

\[ y_{it} - y_{it-1} = (\mu_t - \mu_{t-1}) + \beta(x_{it} - x_{it -1}) + (\epsilon_{it} - \epsilon_{it-1}) \]

We need a slightly different model for asymmetric effects in which we decompose the differences into positive and negative variables.

Our asymmetric effects model will be:

\[ y_{it} - y_{it-1} = (\mu_t - \mu_{t-1}) + \beta^+x_{it}^+ + \beta^-x_{it}^- + (\epsilon_{it} - \epsilon_{it-1}) \]

Where

\[ x_{it}^+ = x_{it} - x_{it -1} \text{ if } (x_{it} - x_{it -1}) > 0, \text{otherwise } 0 \ x_{it}^- = -(x_{it} - x_{it -1}) \text{ if } (x_{it} - x_{it -1}) < 0, \text{otherwise } 0 \]

In other words, if the difference is positive, it becomes part of the \(x_{it}^+\) and if it is negative, it is multiplied by -1 to be made positive and is made part of the \(x_{it}^-\) variable. If the effects are symmetric, \(\beta^+ = -\beta^-\).

After fitting the model via GLS, we can then do a test of the contrasts of the \(\beta^+\) and \(\beta^-\) coefficients as a formal way to assess the presence of asymmetric effects.

Here's how it works with the `panelr`

function, `asym()`

.

```
model <- asym(lwage ~ ms + occ + union + wks, data = wages)
summary(model)
```

```
MODEL INFO:
Entities: 595
Time periods: 2-7
Dependent variable: lwage
Model type: Linear asymmetric effects (first differences)
Variance structure: toeplitz-1 (theta = -0.44)
Standard errors: Cluster-robust (CR2)
------------------------------------------------
Est. S.E. t val. p
----------------- ------- ------ -------- ------
(Intercept) 0.10 0.00 41.12 0.00
+ms -0.04 0.02 -1.95 0.05
-ms 0.04 0.04 1.23 0.22
+occ -0.02 0.02 -0.97 0.33
-occ 0.03 0.02 1.17 0.24
+union 0.01 0.02 0.64 0.52
-union -0.03 0.03 -1.11 0.27
+wks 0.00 0.00 0.48 0.63
-wks -0.00 0.00 -0.35 0.72
------------------------------------------------
Tests of asymmetric effects:
--------------------------
chi^2 p
----------- ------- ------
ms 0.01 0.92
occ 0.19 0.66
union 0.44 0.51
wks 0.00 1.00
--------------------------
```

As you can see, in a model comparable to our within-between model from earlier, the effects seem quite symmetric.

Let's look at the `teen`

data from earlier, where `spouse`

indicates whether
the respondent is living with a spouse, `inschool`

indicates whether the
respondent is enrolled in school, and `hours`

is the hours worked in the week
of the survey.

```
summary(asym(hours ~ spouse + inschool, data = teen))
```

```
MODEL INFO:
Entities: 1151
Time periods: 2-5
Dependent variable: hours
Model type: Linear asymmetric effects (first differences)
Variance structure: toeplitz-1 (theta = -0.54)
Standard errors: Cluster-robust (CR2)
------------------------------------------------
Est. S.E. t val. p
----------------- ------- ------ -------- ------
(Intercept) 1.16 0.15 7.74 0.00
+spouse -4.71 1.15 -4.09 0.00
-spouse -0.61 2.20 -0.28 0.78
+inschool -5.65 1.29 -4.38 0.00
-inschool 7.66 0.69 11.09 0.00
------------------------------------------------
Tests of asymmetric effects:
-----------------------------
chi^2 p
-------------- ------- ------
spouse 6.21 0.01
inschool 2.14 0.14
-----------------------------
```

Here we see an asymmetric effect of marriage: gaining a spouse corresponds with
fewer hours worked, but there's no effect on work hours when a spouse is lost.
You can see in the lower table that this difference in coefficients is
associated with a fairly low *p* value. There is only weak evidence of an
asymmetric effect for entering/leaving school.

The downside to the first differences method is that it does not generalize to non-continuous dependent variables — you can't run a logit model with a differenced binary outcome. Allison (2019) showed that you can do a modified form for such situations.

Instead of including the \(x_{it}^+\) and \(x_{it}^-\) as predictors, you instead create new variables \(z_{it}^+\) and \(z_{it}^-\) that are the cumulative sum of all differences prior to time \(t\).

\[ z_{it}^+ = \sum_{s = 1}^{t}{x_{is}^+} \ z_{it}^- = \sum_{s = 1}^{t}{x_{is}^-} \ \]

Note that at \(t = 1\), both are set to 0. I'll leave the details as to *why*
this works to the manuscript, but he shows that we're left with the following
equation:

\[ y_{it} = \mu_t + \beta^+ z_{it}^+ + \beta^-z_{it}^- + \alpha_i + \epsilon_{it} \]

So we can treat this like a fixed effects model in which we just need to address the \(\alpha_i\). For situations like this that call for a conditional logit, as Allison used in his paper, another option is the GEE with logit link.

Let's try with the `teen`

data, which also appears in Allison (2019). Here our
outcome variable is `pov`

, poverty, and there's a new predictor, `mother`

, an
indicator for whether the respondent has ever had any children.

```
model <- asym_gee(pov ~ mother + spouse + inschool + hours, data = teen, family = binomial(link = "logit"),
use.wave = TRUE, wave.factor = TRUE)
```

```
mother does not decrease over time so -mother is not included in the model.
```

```
Unordered factor wave variable was converted to ordered. You should check
that the order is correct.
```

```
summary(model)
```

```
MODEL INFO:
Entities: 1151
Time periods: 2-5
Dependent variable: pov
Model family: binomial, Link: logit
Variance: ar1 (alpha = 0.33)
Specification: Asymmetric effects (via GEE)
MODEL FIT:
QIC = 5898.64, QICu = 5897.59, CIC = 11.52
------------------------------------------------
Est. S.E. z val. p
----------------- ------- ------ -------- ------
(Intercept) -0.30 0.06 -5.19 0.00
+mother 0.72 0.11 6.63 0.00
+spouse -0.70 0.14 -5.14 0.00
-spouse 0.43 0.25 1.71 0.09
+inschool -0.02 0.16 -0.15 0.88
-inschool -0.01 0.09 -0.09 0.93
+hours -0.02 0.00 -8.17 0.00
-hours 0.01 0.00 1.61 0.11
wave.L 0.11 0.07 1.65 0.10
wave.Q -0.00 0.05 -0.02 0.98
wave.C -0.05 0.05 -0.97 0.33
------------------------------------------------
Tests of asymmetric effects:
-----------------------------
chi^2 p
-------------- ------- ------
spouse 1.10 0.29
inschool 0.04 0.85
hours 25.57 0.00
-----------------------------
```

The results are broadly similar in terms of coefficient estimates to those
obtained by Allison. Unlike Allison, we do not have good evidence of an
asymmetric effect in the case of `spouse`

but we do have one in the case of
`hours`

. Note that `mother`

never goes down so the negative version of this
variable is dropped from the model with a message. To match Allison, I also
used `use.wave`

to include the wave variable and `wave.factor`

to make it
a factor variable.

Allison, P. D. (2009). Fixed effects regression models. Thousand Oaks, CA: SAGE Publications. https://doi.org/10.4135/9781412993869.d33

Allison, P. D. (2019). Asymmetric fixed-effects models for panel data.
*Socius*, *5*, 1–12. https://doi.org/10.1177/2378023119826441

Bell, A., & Jones, K. (2015). Explaining fixed effects: Random effects
modeling of time-series cross-sectional and panel data.
*Political Science Research and Methods*, *3*, 133–153.
https://doi.org/10.1017/psrm.2014.7

Giesselmann, M., & Schmidt-Catran, A. W. (2020). Interactions in fixed effects
regression models. *Sociological Methods & Research*, 1–28.
https://doi.org/10.1177/0049124120914934

McNeish, D. (2019). Effect partitioning in cross-sectionally clustered data
without multilevel models. *Multivariate Behavioral Research*, Advance online
publication. https://doi.org/10.1080/00273171.2019.1602504