Models for event counts

## Warning: package 'merTools' was built under R version 4.2.3

In this session, we will be analyzing hierarchical data.

We will also familiarize with the ggeffects-package that will allow you to do predictions from various models and estimate their uncertainty.

We will familiarize with:

  • how to estimate hierarchical models using the lme4 package
    • with random intercepts
    • with level-2 variables
    • with cross-level interactions

We will use the data on Members of the European Parliament (MEPs). We are interested in whether MEPs invest less resources in their circumscription when they are more electorally safe (i.e. their party received a high vote share in the last election).

This time, the data includes several groupings; it is a time-series cross-sectional (imbalanced panel) data in which each MEP is observed every 6 months over a 5-year period. Thus, I have observations of individuals (ID) in up to 10 periods (Period). These individuals also have a nationality (Nationality). Thus, I have at least 3 “groupings” that are relevant for different reasons.

  • accounting for individual groupings allows me to adjust the standard error of my predictors.
  • accounting for nationality captures potential confounders at the member state level, including the simpsons paradox
  • accounting for time may allow me to control away potential periodic changes

I begin by loading in my data. I also rename a few of my variables for ease.

load("MEP.rda")
df <- MEP
df$y <- df$ShareOfLocalAssistants

Introduction

Exploring the variables when our data has a hierarchical structure requires a few additional steps. In addition to exploring the distribution of the \(x\) and \(y\) in the pooled (entire) data set, I also explore their distribution within the groups.

Univariate statistics

Dependent variable: Size of local staff

I approximate local investment/constituency work by the number of staffers that MEPs emply. They vary from observation to observation.

I begin by my usual histogram. My variable is continuous, albeit limited at 0 (noone can have a negative number of local staffers). I will ignore that here, and rely on linear modeling anyways. It means I’d have to check my assumptions later.

df %>%
  ggplot +
  geom_histogram(aes(y)) +
  ggtitle("Size of local staff among MEPs")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Second, I check the same distribution within my main grouping, the member states. If staff size vary greatly by nationality, we may want to add random intercepts to the model to account for additional variance. It also tells us if we have within-group variation.

df %>%
  ggplot +
  geom_histogram(aes(y)) +
  ggtitle("Size of local staff among MEPs") +
  facet_wrap(~Nationality)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

I can make the histograms comparable by reporting the proportions rather than raw count.

df %>%
  ggplot +
  geom_histogram(aes(x = y,
                     y = ..ncount..)) +
  ggtitle("Size of local staff among MEPs") +
  facet_wrap(~Nationality)
## Warning: The dot-dot notation (`..ncount..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(ncount)` instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

A relevant question is whether staff size varies within MEPs. Will I be able to leverage within-individual changes to assess whether MEPs respond to different electoral incentives? Or will I have to rely on cross-sectional variation (between-individuals)?

Here, I do this by counting the number of MEPs that have hired or fired at least one staffer (i.e. they experienced a change during the period of study).

tab_change <-
  df %>%
  group_by(ID) %>%
  summarize(change = length(unique(y))>1) %>%
  group_by(change) %>%
  summarize(N = n())
tab_change
I find that I have

MEPs that have experienced at least one within-individual change in staff size. That’s cool!

I can also chech how many changes individuals experience

  df %>%
    arrange(Period) %>%
  group_by(ID) %>%
  mutate(change = lag(y) != y) %>%
  group_by(ID) %>%
  summarize(N = sum(change, na.rm = T)) %>%
    dplyr::select(N) %>%
    ggplot +
    geom_histogram(aes(N)) +
  ggtitle("Alterations in local staff size by MEPs") +
  xlab("Number of hiring/firing decisions")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Here we see that quite a few MEPs make several changes in their staff during the period of study. It means I can leverage changes over time (if I want to) to observe within-indual effects.

Main predictor

My main independent variable is vote share acquired last election. That one too, is metric.

df %>%
  ggplot +
  geom_histogram(aes(VoteShare_LastElection)) +
  ggtitle("Vote share to MEPs' national party last election")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1074 rows containing non-finite values (`stat_bin()`).

A relevant question is whether this one too varies within MEPs; i.e. will I be able to leverage within-individual changes in incentives (vote share) to assess whether they change their behavior (staff size)? Or will I have to rely on cross-sectional variation (between-individuals)?

tab_change <-
  df %>%
  group_by(ID) %>%
  summarize(change = length(unique(VoteShare_LastElection))>1) %>%
  group_by(change) %>%
  summarize(N = n())
tab_change
I find that I have

within-individual changes in vote share. This is getting cooler and cooler!

Bivariate statistics

Let’s explore the bivariate relationship in a pooled data first.

Here, I plot both a smoothed line and a bivariate linear regression.

df %>%
  ggplot +
  geom_smooth(aes(x = VoteShare_LastElection,
                  y = y))+
  geom_smooth(aes(x = VoteShare_LastElection,
                  y = y),
              method = "lm", 
              lty = 2,
              color = "black") +
  ylab("Size of local staff") +
  xlab("Vote-share last election") +
  ggtitle("Bivariate relationship between vote share and local investment")
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## Warning: Removed 1074 rows containing non-finite values (`stat_smooth()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 1074 rows containing non-finite values (`stat_smooth()`).

The loess line indicates first a drop, then an increase. What is going on? The liner model says I should expect more local investment among MEPs that received large support last election; not really what I expected.

Let’s chech how this looks within member states. Member states have different party systems and different ways of aggregating votes to seats, so we may be tapping in to some cross-national variation unrelated to electoral uncertainty as such.

With many groups, it is easier to get a sense of the relationship when I don’t plot the uncertainty. The same goes for linear predictions instead of loess lines. The larger the data, the longer it will take to calculate a loess line as well.

df %>%
  ggplot +
  geom_smooth(aes(x = VoteShare_LastElection,
                  y = y,
                  color = Nationality),
              se = F, method = "lm") +
  ggtitle("Within-state relationship between vote share and local staff size")
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 1074 rows containing non-finite values (`stat_smooth()`).

We get a clear sense that the effect is in fact negative, but was masked by Malta.

Aggregating data

Another way of exploring hierarchical data is to look at the grouped relationship. Let’s create a new data frame with individual-level covariates. Now, we can explore between-individual variation

I already aggregated data to the individual level when I counted the number of changes in staff. dplyr is our friend. Here, I aggregate the dependent variable and include some groupings that dont change over time.

df_id <-
  df %>%
  group_by(Nationality, Female, ID) %>%
  summarize(staff = mean(y),
            vote = mean(VoteShare_LastElection, na.rm = T))
## `summarise()` has grouped output by 'Nationality', 'Female'. You can override
## using the `.groups` argument.

Let’s chech the between- individual variation in local investment and vote share.

df_id %>%
  ggplot +
  geom_smooth(aes(y = staff,
                  x = vote),
              method = "loess", se = T) 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 151 rows containing non-finite values (`stat_smooth()`).

We seem to have a curvy relationship again. Is this masked by member states?

df_id %>%
  ggplot +
  geom_smooth(aes(y = staff,
                  x = vote),
              method = "lm", se = F) +
  facet_wrap(~Nationality)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 151 rows containing non-finite values (`stat_smooth()`).

It’s less clear, but most lines seem either flat or decreasing. I’d expect quite a lot of the variation to happen within individuals then; and not between.

If I have cateegorical variables, I can explore the distribution of my numeric outcome (staff) for different categories (men and women respectively).

df_id %>%
  filter(!is.na(Female)) %>%
  ggplot() +
  geom_density(aes(staff,
                   fill = Female %>% as.factor),
               alpha = 0.4)

There is not a lot to report there; there seems to be slightly more men with larger staff.

Over time

In panel data, we may have bumps and drops in the outcome of our dependent variable over time. If I don’t expect a linear trend, I can still calculate group averages.

df %>%
  group_by(DatePeriod) %>%
  summarize(staff = mean(y)) %>%
  ggplot +
  geom_line(aes(x = DatePeriod %>% as.Date,
                y = staff))

We can note a few thins straight away: We had a jump in the Spring term of 2014 (EU election) and a significant dump right after (new parliament), then a drop after 2016 (reform capping MEPs’ local spending). Plotting variables over time always jolts the imagination.

Models

Let’s run some models instead.

Random intercepts

I usually start with a simple pooled analysis with one predictor.

mod.pool <- lm(y ~ VoteShare_LastElection ,
               df)

summary(mod.pool)
## 
## Call:
## lm(formula = y ~ VoteShare_LastElection, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.681 -1.857 -0.513  0.678 40.661 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.276853   0.070836  32.143   <2e-16 ***
## VoteShare_LastElection 0.006895   0.002967   2.324   0.0202 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.874 on 6067 degrees of freedom
##   (1074 observations deleted due to missingness)
## Multiple R-squared:  0.0008894,  Adjusted R-squared:  0.0007247 
## F-statistic: 5.401 on 1 and 6067 DF,  p-value: 0.02016

I find indications that MEPs hold a bigger local staff when they receive a larger vote share. Strange…

I add individual intercepts.

mod.id <- lmer(y ~ VoteShare_LastElection + (1|ID),
               df)

summary(mod.id)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ VoteShare_LastElection + (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 25166.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.3248 -0.3010 -0.0621  0.2829 24.7955 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 6.848    2.617   
##  Residual             2.291    1.514   
## Number of obs: 6069, groups:  ID, 1023
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            2.310758   0.127163  18.172
## VoteShare_LastElection 0.008105   0.004641   1.746
## 
## Correlation of Fixed Effects:
##             (Intr)
## VtShr_LstEl -0.746

I get a stronger effect in the unexpected direction. Yet, I leverage the right kind of variation, don’t I?

Fixed effects to leverage within-group variation

I can always go for a fixed-effect model to leverage only within-individual variation (i.e. over time).

mod.id.fix <- lm(y ~ VoteShare_LastElection + 
                     as.factor(ID),
               df)

stargazer(mod.id.fix,
          type = "html",
          omit = "as.factor")
Dependent variable:
y
VoteShare_LastElection 0.003
(0.006)
Constant 1.546***
(0.489)
Observations 6,069
R2 0.771
Adjusted R2 0.725
Residual Std. Error 1.509 (df = 5045)
F Statistic 16.604*** (df = 1023; 5045)
Note: p<0.1; p<0.05; p<0.01

Nope. If the true within-individual effect is masked, it is by something that varies over time. The EP capped the local spending in 2016, so let’s control for that. the same goes for seniority and age.

mod.id.fix.t <- lm(y ~ VoteShare_LastElection + 
                     Reform2016 +
                     Incumbent +
                     Age +
                     as.factor(ID),
                   df)

stargazer(mod.id.fix, mod.id.fix.t,
          type = "html",
          omit = "as.factor")
Dependent variable:
y
(1) (2)
VoteShare_LastElection 0.003 -0.010
(0.006) (0.006)
Reform2016 -0.438***
(0.069)
Incumbent 0.179*
(0.095)
Age -0.176***
(0.027)
Constant 1.546*** 12.707***
(0.489) (1.709)
Observations 6,069 6,069
R2 0.771 0.781
Adjusted R2 0.725 0.737
Residual Std. Error 1.509 (df = 5045) 1.475 (df = 5042)
F Statistic 16.604*** (df = 1023; 5045) 17.569*** (df = 1026; 5042)
Note: p<0.1; p<0.05; p<0.01

Now, we’re talking! It’s tempting to add the electoral period of 2014 as a control. However, since I assume the lower local investment among MEPs from parties with strong support is due to electoral incentives, I’d be controlling away my own effect.

Several groups of random effects

Let’s leverage both within- and between-individual variation (for good sports) by reverting to random intercepts. However, I’d like to add yet another group variable: nationality. Nationality controls for other potential confounders. For example, both vote share and staff size are reasonably linked to nationality: One because of the party system, the other because of labor cost.

mod.ran <- lmer(y ~ VoteShare_LastElection
                + (1|ID)
                + (1|Nationality),
               df)

summary(mod.ran)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ VoteShare_LastElection + (1 | ID) + (1 | Nationality)
##    Data: df
## 
## REML criterion at convergence: 24698.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.1269 -0.2978 -0.0416  0.2914 24.4319 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  ID          (Intercept) 3.654    1.912   
##  Nationality (Intercept) 3.695    1.922   
##  Residual                2.304    1.518   
## Number of obs: 6069, groups:  ID, 1023; Nationality, 28
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)             2.457426   0.383950    6.40
## VoteShare_LastElection -0.007239   0.004388   -1.65
## 
## Correlation of Fixed Effects:
##             (Intr)
## VtShr_LstEl -0.238

This time, the effect flips because of my group control.

Let’s add the time-varying controls as well.

mod.ran.t <- lmer(y ~ 
                  VoteShare_LastElection
                  + Reform2016
                  + Age
                + (1|ID)
                + (1|Nationality),
               df)

summary(mod.ran.t)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ VoteShare_LastElection + Reform2016 + Age + (1 | ID) + (1 |  
##     Nationality)
##    Data: df
## 
## REML criterion at convergence: 24497.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4733 -0.2973 -0.0445  0.2993 24.7773 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  ID          (Intercept) 3.648    1.910   
##  Nationality (Intercept) 3.734    1.932   
##  Residual                2.212    1.487   
## Number of obs: 6069, groups:  ID, 1023; Nationality, 28
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)             3.935761   0.494324   7.962
## VoteShare_LastElection -0.011815   0.004349  -2.717
## Reform2016             -0.709434   0.053272 -13.317
## Age                    -0.023262   0.005745  -4.049
## 
## Correlation of Fixed Effects:
##             (Intr) VtS_LE Rf2016
## VtShr_LstEl -0.207              
## Reform2016   0.071  0.060       
## Age         -0.625  0.036 -0.166
stargazer(mod.pool,  mod.id, mod.ran, mod.ran.t,
          title = "Random effects (on individual MEPs and their nationality)",
          type = "html")
Random effects (on individual MEPs and their nationality)
Dependent variable:
y
OLS linear
mixed-effects
(1) (2) (3) (4)
VoteShare_LastElection 0.007** 0.008* -0.007* -0.012***
(0.003) (0.005) (0.004) (0.004)
Reform2016 -0.709***
(0.053)
Age -0.023***
(0.006)
Constant 2.277*** 2.311*** 2.457*** 3.936***
(0.071) (0.127) (0.384) (0.494)
Observations 6,069 6,069 6,069 6,069
R2 0.001
Adjusted R2 0.001
Log Likelihood -12,583.270 -12,349.030 -12,248.610
Akaike Inf. Crit. 25,174.550 24,708.070 24,511.230
Bayesian Inf. Crit. 25,201.390 24,741.620 24,558.210
Residual Std. Error 2.874 (df = 6067)
F Statistic 5.401** (df = 1; 6067)
Note: p<0.1; p<0.05; p<0.01

Group-level effects

I can add group-level variables to this. The random intercepts are a way to cluster the residuals. The level-two variables are then regressed on these group means.

In R it looks the same, though. Here, I add gender to the mix. However, its’ effect was already subsumed by the individual-level intercepts, so the effect of vote share is small.

mod.ran.t.c <- lmer(y ~ 
                  VoteShare_LastElection 
                  + Female
                  + Reform2016
                  + Age
                + (1|ID)
                + (1|Nationality),
               df)

summary(mod.ran.t.c)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ VoteShare_LastElection + Female + Reform2016 + Age + (1 |  
##     ID) + (1 | Nationality)
##    Data: df
## 
## REML criterion at convergence: 24487
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4650 -0.2975 -0.0436  0.3013 24.7727 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  ID          (Intercept) 3.638    1.907   
##  Nationality (Intercept) 3.660    1.913   
##  Residual                2.213    1.487   
## Number of obs: 6067, groups:  ID, 1021; Nationality, 28
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)             4.118588   0.498547   8.261
## VoteShare_LastElection -0.011755   0.004350  -2.702
## Female                 -0.281310   0.134398  -2.093
## Reform2016             -0.705550   0.053314 -13.234
## Age                    -0.024684   0.005779  -4.271
## 
## Correlation of Fixed Effects:
##             (Intr) VtS_LE Female Rf2016
## VtShr_LstEl -0.203                     
## Female      -0.162 -0.022              
## Reform2016   0.076  0.060 -0.024       
## Age         -0.634  0.035  0.102 -0.168

The clustering means that the standard errors are also corrected to account for the correlation within the group (i.e. that observations are not i.i.d.). Here, it also means that the uncertainty of gender accounts for the fact that MEPs in general didn’t change gender during the period of study.

Cross-level interactions

We can make interactions accross levels. Here, I check whether men and women respond to vote share in the same way.

mod.ran.t.c.i <- lmer(y ~ 
                  VoteShare_LastElection 
                  * Female
                  + Reform2016
                  + Age
                + (1|ID)
                + (1|Nationality),
               df)

summary(mod.ran.t.c.i)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ VoteShare_LastElection * Female + Reform2016 + Age + (1 |  
##     ID) + (1 | Nationality)
##    Data: df
## 
## REML criterion at convergence: 24486
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4273 -0.3010 -0.0464  0.3021 24.7573 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  ID          (Intercept) 3.664    1.914   
##  Nationality (Intercept) 3.623    1.903   
##  Residual                2.207    1.486   
## Number of obs: 6067, groups:  ID, 1021; Nationality, 28
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    4.247646   0.499609   8.502
## VoteShare_LastElection        -0.019631   0.005113  -3.840
## Female                        -0.803357   0.222335  -3.613
## Reform2016                    -0.701027   0.053278 -13.158
## Age                           -0.024237   0.005797  -4.181
## VoteShare_LastElection:Female  0.025357   0.008592   2.951
## 
## Correlation of Fixed Effects:
##             (Intr) VtS_LE Female Rf2016 Age   
## VtShr_LstEl -0.217                            
## Female      -0.166  0.406                     
## Reform2016   0.079  0.036 -0.037              
## Age         -0.632  0.014  0.037 -0.168       
## VtShr_LsE:F  0.085 -0.524 -0.795  0.027  0.031

They don’t. Men are more reactive to the party’s support than women.

Random slopes

We can estimate separate slopes for each country. This is a form of interaction effect. However, we specify it differently. It means we will estimate mini-models within each group. It requires a lot of variation (and will “explain” a lot of variation), so I remove some of the more redundant covariates. Controlling for gender when I already have individual-level intercepts is over the top.

mod.slope <- lmer(y ~ 
                  + Reform2016
                  + Age
                + (1|ID)
                + (VoteShare_LastElection |Nationality),
               df)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00274656 (tol = 0.002, component 1)
summary(mod.slope)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ +Reform2016 + Age + (1 | ID) + (VoteShare_LastElection |  
##     Nationality)
##    Data: df
## 
## REML criterion at convergence: 24426.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.1624 -0.3076 -0.0409  0.2968 22.4524 
## 
## Random effects:
##  Groups      Name                   Variance Std.Dev. Corr 
##  ID          (Intercept)            4.02625  2.0066        
##  Nationality (Intercept)            8.24286  2.8710        
##              VoteShare_LastElection 0.01097  0.1047   -0.88
##  Residual                           2.12779  1.4587        
## Number of obs: 6069, groups:  ID, 1023; Nationality, 28
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  3.245224   0.425967   7.618
## Reform2016  -0.723314   0.052859 -13.684
## Age         -0.020990   0.006012  -3.491
## 
## Correlation of Fixed Effects:
##            (Intr) Rf2016
## Reform2016  0.115       
## Age        -0.749 -0.180
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00274656 (tol = 0.002, component 1)

We can’t see the random effects, but they are there. Let’s check them out.

Interpretation

Random intercepts

The random intercepts are stored, but not immediately reported.

We can always take a look at the coefficients using the coefficients() function. It reports the coefficients for all groupings: Their (varying) intercepts, and their slopes regardless of whether they are “fixed” or “vary” by group.

Here, I look at the first coefficients for the 10 first MEPs.

coefficients(mod.slope)$ID[1:10,]

And here I can see the varying slopes i estimated.

coefficients(mod.slope)$Nationality[1:10,]

The intercepts are interpreted s deviations from the (grand) mean. Here, we see that Austria has a higer estimated number of local staff that the adjused average (grand mean).

I can separate out the “fixed” effects (those that are common for all groups) and the “random” effects (the varying slopes and intercepts) using two functions.

fixef(mod.slope)
## (Intercept)  Reform2016         Age 
##  3.24522415 -0.72331434 -0.02099024

And here are the varying ones

ranef(mod.slope)
## $ID
##          (Intercept)
## 840     0.4470575887
## 988    -0.0443138152
## 997    -0.4410573470
## 1023    2.1425511331
## 1059   -0.0313271030
## 1073    2.7250304459
## 1122   -0.2465788689
## 1129    0.6101969169
## 1164    1.0607100738
## 1179   -0.8150482473
## 1183   -0.0808875249
## 1185   -0.5062726168
## 1186    0.5596368872
## 1204    0.2679528549
## 1309    0.0075546391
## 1351   -1.0856090026
## 1394    3.5558456510
## 1403   -1.0588202717
## 1405    0.0725310629
## 1407   -0.8879596175
## 1419    0.4333059485
## 1431   -1.1274524137
## 1473   -0.3743233812
## 1654    2.1587010361
## 1665    1.4042854733
## 1746    0.2711182345
## 1832   -0.9515346855
## 1849   -0.7649496774
## 1854   -0.1811013145
## 1892   -0.7271466490
## 1902    0.1102401567
## 1906    1.4064731688
## 1907   -0.5013601785
## 1909    0.4387397912
## 1911   -0.6396720497
## 1913    0.1584064909
## 1914    0.4878663357
## 1917    0.0956977668
## 1933   -0.2658182813
## 1941    0.6981827255
## 1945    1.5696880780
## 1946    0.6740922608
## 1956   -0.9768078033
## 1965   -0.3078798497
## 1985   -1.3495663021
## 1993   -1.5035926502
## 2002   -1.6291974042
## 2025   -0.9169845734
## 2037   -0.3723345532
## 2054    0.1944949942
## 2073   -0.7157727272
## 2081   -1.3785521279
## 2097    2.3009974342
## 2109    1.7139553575
## 2119    0.4666370592
## 2125    0.6610649235
## 2128    1.5029636432
## 2155    2.0148004602
## 2173   -1.2235197434
## 2187   -1.1377034302
## 2212    2.7118243420
## 2229    0.8576996108
## 2278    0.8477705123
## 2295    0.2439532086
## 2307   -0.5841256148
## 2309   -1.2779334249
## 2315    0.1652977521
## 2319    0.7348133776
## 2327   -0.0840095216
## 2338    0.3916982071
## 4246    1.8950630673
## 4253   -0.8693110987
## 4260   -1.2264411788
## 4262    0.0522539242
## 4265   -0.3592341105
## 4267    0.4907536529
## 4270   -0.1911874754
## 4274   -0.0167167702
## 4289    3.1414838944
## 4293   -0.5608552605
## 4294    0.1115578907
## 4308    2.5407109405
## 4318   -1.1860424508
## 4319   -0.8152653780
## 4326    0.4731641989
## 4328    1.2676901821
## 4334   -0.6925487436
## 4335    1.5477937341
## 4342   -1.1353946741
## 4345    0.3403590814
## 4346   -0.5313477593
## 4347    0.4346942785
## 4348   -0.6337038905
## 4393   -1.8586802003
## 4412   -0.5741401576
## 4423   -0.7781318839
## 4429   -0.0429259200
## 4432   -0.2559653302
## 4436    0.9980634297
## 4440    0.2502791719
## 4462    3.1622568836
## 4482    0.6539019872
## 4507    3.4383215015
## 4513   -1.4718261598
## 4514    0.3055840685
## 4516   -0.9838777790
## 4519   -1.6613223025
## 4521   -1.2455033866
## 4525    0.5650448813
## 4529    0.1205938289
## 4531    0.1776424693
## 4532   -0.1326550541
## 4533    0.8553021709
## 4536    0.2066974646
## 4537   -1.3075845040
## 4538    0.6624011032
## 4540   -1.4302094057
## 4542   -1.0929419234
## 4545   -1.0564926090
## 4546   -1.3486206223
## 4550   -0.0033188317
## 4553   -1.1461973623
## 4554    2.2281884714
## 4555   -2.6632350405
## 4556    0.4612062853
## 4560   -0.0016735338
## 4746    0.0795088602
## 4751   -0.8420979748
## 5537   -0.6601407217
## 5565    0.8860201835
## 5729    0.5335578288
## 5737    0.3242446385
## 5846    0.5914607632
## 5956   -0.8581427979
## 21331  -0.9948710616
## 21817  -1.3704447230
## 21818  -0.4276127809
## 21918   0.2171286613
## 22097   1.0856357490
## 22418  -1.0699237188
## 23413  -0.1436659265
## 23691   0.2999685039
## 23693  -0.7619853349
## 23695   2.2485849670
## 23699  -0.4970585968
## 23704  -0.1751767660
## 23705  -2.2691798131
## 23706   0.5261245895
## 23707  -1.0837280050
## 23712   0.5015198047
## 23768   4.3838202607
## 23781  -3.7524186060
## 23782  -1.0312951150
## 23784   3.7224693900
## 23785   3.7580728634
## 23787   2.2206931533
## 23788  -3.4152812103
## 23791   2.0984468362
## 23792  -2.9163721190
## 23805  -1.3034651030
## 23808  -1.6004656355
## 23840   1.3514792308
## 23866  -1.4707707375
## 23868  -1.7589127654
## 23873   1.1684117398
## 23894   1.0694156612
## 23938  -1.2628721304
## 24030   2.7080112176
## 24505  -0.0472191081
## 24594  -1.3290653774
## 24942  -0.8119997666
## 25704  -0.0861354404
## 25718  -0.3819680278
## 25919  -0.3893518728
## 26829   0.8017119423
## 26837   0.3329820098
## 26851   0.5568292994
## 28112  -0.3622463841
## 28114   0.5746239977
## 28115   0.3124206306
## 28117   1.5730063217
## 28120   0.9695663038
## 28122  -0.3847108430
## 28123   0.8013022604
## 28124   0.6977922874
## 28125  -0.5132265947
## 28126  -0.7490308742
## 28130  -0.8643054584
## 28131  -0.5367417255
## 28132  -1.7103996060
## 28139   1.0799129905
## 28148  -0.3903778940
## 28153  -0.7606678981
## 28154  -0.5201658541
## 28155  -1.7687186826
## 28156   1.5056361520
## 28161   0.4036705296
## 28165   1.0571181118
## 28166  -1.1449864511
## 28167  -0.0520805197
## 28170  -1.4728495530
## 28171  -0.3884039043
## 28174  -0.5627222378
## 28177  -0.2879177319
## 28178   0.4144880627
## 28182  -0.8227508596
## 28192   2.4232636255
## 28193  -1.1442883649
## 28206  -0.1869740372
## 28208   0.0201261433
## 28210  -0.5575823611
## 28214  -1.1926739928
## 28227   1.0095287613
## 28228   1.8428520852
## 28229  -0.3704012460
## 28233   0.3045605984
## 28238  -0.4960746062
## 28240  -1.3107284650
## 28241  -1.7496130517
## 28242  -1.1665176419
## 28243  -1.6627039110
## 28244  -1.5571031380
## 28246  -1.1119555637
## 28247   1.5079199387
## 28248   0.5287735521
## 28251   1.3882245615
## 28252  -0.2938345043
## 28253  -0.9083495082
## 28256  -0.6483149407
## 28257   0.6254893867
## 28266   0.5179007768
## 28269  -2.8307364538
## 28275   2.2103087358
## 28276  -7.7109157017
## 28278  -0.5920957164
## 28279   5.3751201765
## 28280   1.1494260198
## 28284  -2.9355959241
## 28288   0.0598464262
## 28291  -0.1450510172
## 28292  -0.7791353959
## 28295  -0.8196924935
## 28297  -2.4880365840
## 28298   7.3244446462
## 28299   1.2149168823
## 28300   1.1140994591
## 28301   0.2362413969
## 28306   0.2372814229
## 28307  -0.4822731403
## 28308  -0.5182461627
## 28310   1.3508816896
## 28314   1.6697484723
## 28316  -0.8334412771
## 28320   2.9117895614
## 28321   0.0789452447
## 28323   0.4958462796
## 28324   1.4539707890
## 28330   3.6056317849
## 28331  -0.4567510513
## 28336  -0.2489481975
## 28338  -0.6501910759
## 28340  -1.2791026100
## 28341   0.3671189713
## 28342   0.0538313910
## 28346  -1.8630137000
## 28349  -0.7346557947
## 28352  -1.5869544748
## 28353   2.4225647047
## 28354  -0.7053790268
## 28358   6.0811112073
## 28362   2.0622368881
## 28365  -1.5732256827
## 28367  -0.8841017573
## 28370   1.9491966349
## 28372   5.5183199483
## 28377  -1.2182864845
## 28380  -0.4472227621
## 28389  -1.3699674680
## 28390  -0.4473706190
## 28393  -0.2976000890
## 28394  -2.3941892663
## 28397   0.2191367512
## 28398  -1.7315531452
## 28399  -0.5148977867
## 28400  -1.5440458078
## 28404  -1.4470741646
## 28407  -0.2677873664
## 28422   2.3138525963
## 28424   0.5699772959
## 28429  -1.1354768744
## 28463  -0.0291575641
## 28474   0.1513418799
## 28477   0.5808543062
## 28481  -0.7408483516
## 28493  -0.8838498621
## 28497  -1.1821425048
## 28506   1.0553877023
## 28508  -0.7118802570
## 28513   0.7996634148
## 28514  -1.4313231934
## 28572  -1.1986545247
## 28584   6.3944710736
## 28586   1.4050481856
## 28615  -0.9923480518
## 28617   0.0250244466
## 28618  -1.8842342104
## 28619   4.8380382184
## 29074   2.7361657407
## 29579  -0.7634740426
## 30190  -0.7586573952
## 30475  22.1562130432
## 33569  -1.0054408850
## 33570   0.6205380989
## 33775   1.5701267021
## 33982  -0.8912956827
## 33984   0.4635630923
## 33989   0.2866795271
## 33990  -0.4768678906
## 33997  -0.1419022725
## 34231   2.7874082588
## 34232   1.3935694097
## 34234  -1.2966090570
## 34235  -0.7688813204
## 34249  -0.2604135776
## 34250  -1.4778280158
## 34254  -0.0794140970
## 34728   1.1890240362
## 35135   1.0207960597
## 35743  -0.7519459332
## 36281   2.6439200349
## 36392  -0.0089043901
## 36396   0.2000226493
## 37229  -0.5411644467
## 37312  -0.7644304195
## 37324   2.2580808086
## 37524  -0.5229948809
## 38398   0.7542791650
## 38420  -1.8053208614
## 38596   0.2423602030
## 38601   0.4041333022
## 38602   0.3577645423
## 38605   0.4247028744
## 38610  -1.3761993285
## 38613   2.4678355957
## 39317  -1.8518176788
## 39319  -1.8593071426
## 39321  -0.6601250754
## 39711   1.4436225631
## 39713   0.4676627167
## 39714   0.8224253425
## 39717   0.2293146653
## 39721  -0.0006084121
## 39722  -2.3246377918
## 39724   0.0269857334
## 39725  -0.9195963554
## 39726   5.4204160472
## 39916   0.0573727558
## 40224  -1.3516584234
## 40599  -0.6776128997
## 41007  -0.4540660308
## 58758   1.0146633848
## 58789  -0.7738635680
## 72775  -0.7367670518
## 72779  -0.2365363839
## 88882   5.6340700211
## 90110  -0.4194871611
## 93624  -0.3666572785
## 94282   1.0471503761
## 94283   0.9561621878
## 95017  -2.0911638396
## 95074  -0.2743255903
## 95281  -1.5418578679
## 96603   0.0279597234
## 96646  -0.0236307134
## 96650   0.2346421660
## 96651  -0.1008450236
## 96652  -1.5809804405
## 96653  -0.7143887852
## 96654   0.9144868271
## 96655   0.0126522657
## 96656   0.7139478513
## 96658  -1.1202750426
## 96660  -1.8749486995
## 96661   2.2543205053
## 96662  -1.1279049276
## 96663  -0.0138548012
## 96664  -1.0195053627
## 96668  -0.1152293866
## 96670   0.5037952042
## 96671   1.3759906654
## 96672   0.4380405383
## 96673  -0.1699623932
## 96674   0.0835588568
## 96675  -0.6353364327
## 96677  -0.5599602521
## 96678  -0.1746660764
## 96680   1.0121531187
## 96681  -0.3956568660
## 96682   0.0237530925
## 96683  -0.2977967742
## 96684  -0.3292826256
## 96685   0.2479972510
## 96686  -0.4634435064
## 96692  -8.0302844677
## 96693   5.2883896220
## 96694  -2.5631264409
## 96696  -5.8217609566
## 96697  -0.2769457485
## 96698   0.0600233107
## 96701  -0.9590278046
## 96702  -1.5923836698
## 96703   2.5158650804
## 96704  -0.6969064778
## 96705  -0.3810343402
## 96709   0.0495050349
## 96710  -1.0719242917
## 96713  -1.1622421696
## 96714   1.3824457346
## 96715  -1.6617363285
## 96716  -0.1475517325
## 96717  -0.9658163843
## 96718   0.9096136881
## 96719  -2.5192969970
## 96725  -0.3092546979
## 96726  -0.6145655872
## 96728   0.6831415964
## 96729  -0.2741729605
## 96730  -0.7610211574
## 96731  -1.3585353861
## 96733   0.7754279380
## 96734  -0.2408537112
## 96736   1.2721701132
## 96739  -0.4415335793
## 96741  -0.0852312279
## 96746   0.4721936127
## 96747  -0.5557145704
## 96748   0.0805795959
## 96752   0.0288448151
## 96753   3.6205721266
## 96754   0.6874634133
## 96755  -0.0132874497
## 96757  -0.2283196247
## 96758  -1.4873647176
## 96760  -2.0789781023
## 96763  -0.2312724744
## 96765  -0.4374842444
## 96768  -0.6079535254
## 96769  -0.8827062550
## 96770   2.4019603599
## 96771   1.0811898959
## 96773  -3.3941224388
## 96774  -2.3814534083
## 96775   0.2443306243
## 96776  -2.5040459004
## 96777   0.1051204719
## 96778   2.0056724172
## 96779  -3.2389227334
## 96780   0.3996673381
## 96781  -3.9905840280
## 96782  -4.6279768409
## 96783  -2.0758013927
## 96784   3.6470107135
## 96785   2.5740617625
## 96786   6.1757278677
## 96787  -0.9480854511
## 96788  -3.9249844228
## 96789   0.3923382685
## 96790   8.0801555486
## 96791  -2.2808910575
## 96792  -1.5510033532
## 96793  -1.9845078081
## 96794  -3.3138002111
## 96795  -1.0822222338
## 96796   6.6342321982
## 96798  -0.5282539560
## 96799  -1.0138297064
## 96800  -1.4085194309
## 96801   0.6611957337
## 96802   1.9942210294
## 96803  -2.9762308735
## 96805   1.4364136382
## 96806  -1.6713124700
## 96807  -1.0167421680
## 96808  -0.2877874829
## 96809   2.5989002779
## 96810   0.6580207236
## 96811   0.0486544478
## 96812  -0.8629784870
## 96813  -0.4916985503
## 96814  -0.5131273544
## 96815   0.3419494674
## 96816  -1.8037855484
## 96817  -0.7533429382
## 96818   3.0725797288
## 96819  -1.6740400793
## 96820   0.3199697872
## 96823  -0.9968965316
## 96824  -0.1100479537
## 96825  -1.5292358608
## 96831   0.6513778689
## 96832  -1.5124642465
## 96833   2.3531811283
## 96834  -0.3289495540
## 96835  -0.4541277819
## 96836  -0.7698935242
## 96837   0.0700393718
## 96838  -0.2567688874
## 96839   0.2709098668
## 96840   0.8249084736
## 96842   0.4818294017
## 96843   2.2826363006
## 96844  -0.8850079911
## 96845   1.0931001945
## 96846  -1.6816258408
## 96847   0.2524013216
## 96848   0.1976506271
## 96849   0.3661453203
## 96850   0.6082495563
## 96851   0.0876268318
## 96852   2.0171635885
## 96854  -0.4323377532
## 96855  -2.1528829769
## 96857  -0.7017421876
## 96859   0.0838763519
## 96861  -0.9481712251
## 96862   0.6331052872
## 96863   0.6009849270
## 96864  -1.2129126149
## 96865   0.3792614933
## 96866  -0.7304066561
## 96867  -0.7522499966
## 96870   0.7414648470
## 96871   3.1988704498
## 96872  -0.8096921612
## 96873   0.1746335466
## 96874  -1.2526233257
## 96875  -1.7627197272
## 96876  -0.7646864951
## 96877   1.3052015399
## 96878  -2.5579779214
## 96880   1.8931991269
## 96882  -1.1815847779
## 96884   2.9088141999
## 96885  -0.9437468382
## 96886  -0.7483519081
## 96887  -1.0082252706
## 96888  -0.2043430339
## 96889   0.9547591345
## 96890   6.8559261242
## 96891  -1.0089440736
## 96892  -0.6552753924
## 96897   0.6662133102
## 96898   0.3272038667
## 96899   0.9342119647
## 96900  -1.6250863893
## 96901   0.8522313408
## 96902   0.0752869163
## 96903  -0.8817237007
## 96904  -0.5781939814
## 96905  -0.8166881844
## 96906   0.9231538185
## 96907  -0.6955573952
## 96908   0.1913163474
## 96910  -0.7337839288
## 96911  -1.2059915920
## 96912   0.1779814215
## 96914  -0.2849997210
## 96915  -1.3290464059
## 96916   0.5303235589
## 96917  -0.0291717784
## 96918  -0.7172474396
## 96919  -0.6204251502
## 96920  -1.4390252459
## 96922  -1.2296540011
## 96925  -1.7775470367
## 96930  -0.3635356765
## 96931   0.5817288484
## 96932  -0.0567946527
## 96933   0.2070023025
## 96934  -1.2468678083
## 96936  -1.0017516104
## 96937  -0.4947742744
## 96940  -0.6024798603
## 96944  -0.0541659015
## 96945  -0.5436690140
## 96946   0.4481206601
## 96947  -0.5136409529
## 96948  -0.2700808175
## 96949  -0.2100141176
## 96950   0.5269184830
## 96951   2.8606612803
## 96952  -0.4691141530
## 96953   1.4278582746
## 96954  -3.0940286480
## 96955   1.4291660148
## 96956  -1.1828286109
## 96957  -0.1579657664
## 96959   1.8279376649
## 96960   2.2537173909
## 96973   0.0023265404
## 96974  -0.1853364369
## 96975  -0.3292851579
## 96976   0.2205412693
## 96977  -0.5620798078
## 96978   0.1377849955
## 96979  -0.9617380950
## 96986  -0.8742511830
## 96987  -1.0945665581
## 96989  -0.8622367558
## 96990  -0.6603851526
## 96991  -1.1440030049
## 96992  -2.9021880434
## 96993  -0.9236243673
## 96994  -2.7827920721
## 96996   1.4248037181
## 96998  -0.8229064673
## 97004   1.3083363486
## 97007  -1.5620685278
## 97008  -1.7001038325
## 97009  -2.1333489602
## 97011  -1.3228548619
## 97014  -0.8895679609
## 97017   1.0275226816
## 97018  -0.2011581253
## 97019  -0.6646735978
## 97020   0.5450323390
## 97022  -0.6311913060
## 97025   0.4547234774
## 97026  -1.1068208414
## 97050   1.8550737516
## 97058  -0.7473838611
## 97076  -1.3058903313
## 97078  -0.4637037858
## 97086   0.3473265225
## 97125   3.1467652107
## 97130   0.8283071202
## 97131  -1.4420276163
## 97137  -0.9100257806
## 97156  -1.2151859317
## 97193   1.2816391594
## 97196   0.4027956947
## 97197   1.0284997419
## 97198  -2.2391691228
## 97199   3.4123007101
## 97203  -0.9831977794
## 97228  -1.7928066066
## 97229   1.4201435632
## 97230  -1.2009339522
## 97280   0.5020682845
## 97293  -0.7474052080
## 97296  -0.3822633111
## 97308  -0.3834992600
## 97344  -0.5949942509
## 97513  -0.8026477562
## 97968  -0.7140184291
## 98341  -0.7775059483
## 99325  -1.7056689149
## 99416  -1.8489435622
## 99419   7.1870670992
## 99650  -0.2140238024
## 101580 -0.0664364344
## 101954 -4.6037099042
## 102887  1.0858294043
## 102931  1.2542117633
## 103132 -0.6701753124
## 103246 -0.3917693768
## 103381 -1.1314761785
## 103488  6.6947466668
## 105192 -0.7720510981
## 105624 -0.3895658517
## 105849 -0.1486744016
## 106202  1.4048957263
## 107041 -0.7099827571
## 107212 -1.7519906513
## 107385 -0.2502056819
## 107386 -0.7515295508
## 107973  0.2295667668
## 107977 -1.9624953310
## 108080  1.0242252646
## 108329  0.1021181707
## 109337 -0.5830119047
## 109649  0.2087032262
## 110365  0.4792392288
## 110977  0.8195972219
## 110984  2.2500529526
## 110987  4.7349485104
## 111011 -0.8484005587
## 111014 -0.6009218768
## 111017 -0.8788281618
## 111019 -0.3883752790
## 111024 -0.1681599227
## 111026  0.1697353140
## 111027  0.5763146246
## 111068  0.7746852173
## 111126 -1.2512730581
## 111496 -1.9518168952
## 112014  0.7763507151
## 112071 -0.4458855245
## 112611 -0.3258207467
## 112620 -1.6092486087
## 112744 -2.7175202732
## 112748 -1.3236597256
## 112755 -0.4177714139
## 112760  1.6949625096
## 113487 -0.1182820655
## 113597 -0.6876553197
## 113890 -0.3530410157
## 113892 -0.1411314333
## 113959 -0.3981816715
## 114268  0.3288653589
## 114279  0.2884778034
## 114840  2.1242229835
## 114873 -0.6464909880
## 115868 -0.8137768604
## 116663  3.5824068048
## 116816  1.7083130585
## 116823 -1.8009761096
## 117119 -0.7622031444
## 117477  0.3246193306
## 117491 -0.9984645453
## 118144 -1.1801105442
## 118658 -0.1801809745
## 118708  1.5341470370
## 118709 -0.3032814808
## 118710  3.8236485190
## 118858 -2.7442637290
## 118859  0.1540284841
## 118860  3.0533287397
## 118999  3.0738992701
## 119431  2.3687347197
## 119434  1.5606378522
## 119435 -1.1160303825
## 119436  8.1475084343
## 119438 -0.2450534893
## 119439 -2.1711191256
## 119440 -3.0843343330
## 119441 -3.3225086614
## 120478 -0.9893561026
## 122317 -2.0317264589
## 122404 -0.5268094965
## 124287 -0.2826701553
## 124688  0.6746406483
## 124689  0.1186704542
## 124691 -2.3473451356
## 124692  0.3456401097
## 124695  0.2709486663
## 124696 -0.7078181631
## 124697 -0.7527876564
## 124698  0.1033101967
## 124699 -0.8963289885
## 124700  0.4277628256
## 124701 -0.9447511271
## 124702 -0.5090614637
## 124703  0.5107925111
## 124705 -0.0239668103
## 124710  0.0127972418
## 124713  1.6982868016
## 124720 -0.4802194360
## 124725  1.2108685795
## 124726 -0.5876864445
## 124727 -0.1916388834
## 124728 -0.8518395057
## 124729  0.9343100060
## 124730 -0.3199001948
## 124732  1.0430013810
## 124733 -0.8423430665
## 124734 -0.2275701693
## 124735 -1.0989349464
## 124736  0.1717448005
## 124737 -0.3058024133
## 124738 -0.1630073667
## 124739 -0.3755480920
## 124741  0.6167570236
## 124743  0.1294410551
## 124744 -0.0340605311
## 124745  2.1962844246
## 124746 -1.4372419112
## 124747  0.4677477902
## 124748 -0.4380661266
## 124750 -0.3098444589
## 124751 -0.5812536164
## 124753 -0.3411742731
## 124754  2.6922076159
## 124757 -0.6123656914
## 124758 -0.6021718045
## 124759 -3.4088582613
## 124760 -0.1658879764
## 124761 -0.5955563432
## 124763  6.0427799305
## 124764 -0.3797481457
## 124765 -0.0006833202
## 124766 -3.5236503079
## 124767 -0.6489708828
## 124768 -0.7294956012
## 124769 -0.0392396874
## 124770 -0.0730833483
## 124771  0.4783008164
## 124772  1.2081201466
## 124776 -0.4164722049
## 124777 -1.1012640190
## 124778 -2.4600754088
## 124779 -3.1530717418
## 124780 -1.9370385622
## 124781  0.1001671195
## 124784  0.3231571858
## 124787 -2.0547855173
## 124788 -1.0853546566
## 124789 -0.8759928251
## 124791 -1.3410469930
## 124793 -0.3032962409
## 124796 -0.5629361523
## 124797 -2.0773321460
## 124799 -1.0465666033
## 124800 -0.1687764256
## 124801 -1.5218666867
## 124803  0.3383252651
## 124811 -2.1393383539
## 124812 -0.6554723850
## 124813 -1.9783146933
## 124814 -1.6903798110
## 124817  1.5919720697
## 124819 -1.6592576659
## 124821 -0.1132495452
## 124822  0.1023230437
## 124826  0.2489959382
## 124828  4.3692361904
## 124831 -0.3324988561
## 124833 -1.1597335567
## 124835 -0.1616972424
## 124836 -0.2083621243
## 124837  0.2331847637
## 124838 -0.2908024070
## 124839  0.5148940862
## 124840 -1.1725819406
## 124841  1.4866654518
## 124842  3.5093910317
## 124843 -2.8184792196
## 124844 -1.9732695755
## 124846 -1.0251099890
## 124847 -1.1354405878
## 124848 -0.8631923049
## 124849  1.0999331971
## 124850 -1.4031909543
## 124851  5.8655167913
## 124852 -1.5421359646
## 124853  6.4681450140
## 124854  3.8494338016
## 124855 -1.6590819867
## 124856 -0.9853533275
## 124858  0.1814486909
## 124859 -2.0113410795
## 124860  0.1497613502
## 124861 -0.9736987813
## 124864 -0.0629035747
## 124866 -1.6375564204
## 124867 -0.2994041128
## 124868 -1.9066282916
## 124869  2.2966906592
## 124870 -0.2019265652
## 124871  6.9476212315
## 124872  0.4525486134
## 124873  1.9497491039
## 124874  2.2856838098
## 124875 -0.4999763566
## 124877  3.0232429369
## 124878 -0.9136988386
## 124880 -0.9236311611
## 124881 -1.4019661868
## 124882  4.4762181888
## 124883  0.2200995973
## 124884  3.3748120391
## 124887  1.7026689778
## 124888  0.5676525654
## 124889 -1.5710934262
## 124890 -4.0742479851
## 124893 -0.4283226451
## 124894  2.9034138848
## 124895 -1.9309073691
## 124896 -0.9151603535
## 124897 -2.5731541631
## 124898  0.3707747408
## 124899  5.1883400311
## 124900  0.9781493865
## 124902  3.1775397828
## 124903 -2.6845105492
## 124926 -1.5597398616
## 124928  1.8873200265
## 124929  0.7356376632
## 124930 -1.3435818642
## 124932  1.2929730181
## 124934  0.3351080360
## 124935 -0.7847683501
## 124936 -0.8228761672
## 124937 -0.1151032288
## 124939  6.1190081848
## 124940 -0.3818308953
## 124941 -0.8818791962
## 124942  0.2520950307
## 124943 -0.9247766117
## 124944 -1.3024923139
## 124945 -0.5950897952
## 124947 -0.2208612409
## 124948 -0.0961330567
## 124949 -1.1148494769
## 124950  2.9224642821
## 124951  4.7772942731
## 124952 -1.9698676440
## 124953 -2.7888227690
## 124954 -1.2430445530
## 124955  0.6800591022
## 124956  0.7895618021
## 124957  0.4770831099
## 124958  0.6827435749
## 124960 -1.4973792351
## 124961 -0.7050287790
## 124962 -0.0733945182
## 124963 -1.1381920605
## 124964  0.2786590160
## 124965  2.5052110397
## 124966  3.8289425851
## 124967  1.1895783937
## 124968 -0.5019419628
## 124969  0.4521189292
## 124970 -1.7021254199
## 124971  0.2362897847
## 124973  0.0745988373
## 124984 -1.2985716480
## 124986 -0.2253485539
## 124987 -0.0232366838
## 124988 -0.5888092934
## 124989 -0.6171602332
## 124990 -0.1865581882
## 124991 -0.5695307480
## 124992 -0.4559802509
## 124993  0.8484336302
## 124994 -0.4494990245
## 124995 -0.1836273366
## 124996  0.9456108413
## 125002  0.2852878901
## 125004  0.3749405654
## 125005 -0.5800973976
## 125013  2.8708752238
## 125016 -0.0742046578
## 125018 -0.9846181454
## 125019 -0.3391649024
## 125020 -0.1331020138
## 125021 -0.1911874930
## 125022  0.2311108977
## 125023 -0.4425900971
## 125025  0.0575265457
## 125026  1.1320596268
## 125027 -0.0482550037
## 125028 -0.0768158527
## 125029 -0.1816635460
## 125030 -0.1295398808
## 125031  2.4804465815
## 125032 -1.3098392708
## 125035  6.3438758431
## 125037 -0.4353375850
## 125039 -1.0664712854
## 125041 -0.9386092928
## 125043 -0.6413186626
## 125044 -0.7188658976
## 125045 -0.5860061888
## 125046 -1.0402892233
## 125047  2.0885927759
## 125048 -0.8391833672
## 125050  2.3194042502
## 125051 -1.0340014943
## 125052 -0.4928541714
## 125053 -0.9719369875
## 125061  0.5624183201
## 125062  0.2925848955
## 125063  0.7524226946
## 125064  0.1385659350
## 125065 -0.8663627827
## 125067 -0.5204257886
## 125068 -0.1939104596
## 125070 -0.5844204963
## 125071 -0.4344229696
## 125072 -0.6598252181
## 125090  3.3495134601
## 125092  1.4049861000
## 125093  3.4421757186
## 125099  0.0843383589
## 125103 -0.0483708492
## 125104  0.3145884950
## 125105 -0.8755515697
## 125107 -0.0819314206
## 125109 -1.4548063038
## 125110 -0.9175705715
## 125112 -0.4515808580
## 125113 -0.8201080177
## 125128 -1.8287060969
## 125214 -4.8455907125
## 125237 -2.1525489463
## 125670 -0.0629918768
## 125684 -0.5080751087
## 125997  1.9610116250
## 127330  6.9704496716
## 128717 -0.0974378561
## 129073 -0.9670258950
## 129164 -1.0803801089
## 129256  6.3113587054
## 130100 -0.8302579333
## 130833  2.0592278555
## 131051 -1.2531463077
## 131507 -0.4941043224
## 131749 -0.1062493607
## 132366 -0.0703123767
## 132925  2.9761146597
## 133316 -0.2341339909
## 133326  1.8403840588
## 134243 -0.2284126413
## 135490 -0.3350796009
## 135540 -3.6777995586
## 135541 -3.8438756134
## 183022  0.1928020383
## 183338 -0.2271052470
## 183793 -1.2415082795
## 185341 -2.0991330255
## 185619 -0.5053564689
## 185637 -2.0262394761
## 
## $Nationality
##                (Intercept) VoteShare_LastElection
## Austria        -1.45291486           0.0402362300
## Belgium        -1.04815710           0.0027568188
## Bulgaria        1.85648912           0.0023250048
## Croatia        -3.66285342           0.1516031178
## Cyprus          1.45040113          -0.0449145834
## Czech Republic  0.97460725          -0.0366798102
## Denmark        -0.29631871          -0.0350621577
## Estonia         0.07098672          -0.0353425336
## Finland        -1.18159386           0.0209643596
## France         -0.48884079          -0.0040586864
## Germany         0.04240794          -0.0054156075
## Greece          0.26547605          -0.0323747935
## Hungary         0.94336127          -0.0553356247
## Ireland         1.09596316          -0.0401321336
## Italy           1.25581052          -0.0019890206
## Latvia         -0.83072708           0.1087012381
## Lithuania      12.12321547          -0.4524590165
## Luxembourg     -0.25308339          -0.0232906293
## Malta          -0.38699776           0.0573647953
## Netherlands    -1.46108690           0.0008384453
## Poland          3.97034691           0.0125701433
## Portugal       -1.95884039           0.0190642455
## Romania         0.50555933          -0.0145496882
## Slovakia       -0.19411997          -0.0048999490
## Slovenia        0.08805439          -0.0372517840
## Spain           0.60690676          -0.0680279021
## Sweden         -1.47575575          -0.0011136927
## United Kingdom  1.48088227          -0.0239490788
## 
## with conditional variances for "ID" "Nationality"

Remember, the group means/varyin intercepts are assumed to be normally distributed, and so are the residuals they are meant to group. In other words, we have 28 distributions at the member-state level. Together, their means also constitute a normal distribution.

This is a lot of distributions (and assumptions), so there are some controversy over whether we can reliably estimate them. That’s why the lme4 package doesn’t report the random intercepts with any p-values. However, R has other packages.

arm has two twin functions: se.fixef() and se.ranef() that are useful.

merTools allows us to also simulate the random/varying coefficients.

re <- REsim(mod.ran,
            n.sims = 100)
re %>% head

the best way is still to display these intercepts graphically. You can use the simulations to construct your own plots: coefplots or density plots are nice possibiliites.

There are a few “ready-mades” out there too.

  • merToools has a quick and dirty way
plotREsim(re)

It gives us a sense of the variation between the groups. Here, we see that there is still quite some individual-level variation that we haven’t accounted for with other predictors. The groups marked in black have a 95% confidence interval that does not overlap zero. We can also see there are member states that have higher/lower level of local staff than what we might expect.

  • sjPlot reports effects from all kinds of models. It also gives a nice labelling of our random intercepts.

Here, I only ask for the varying intercept for nationality.

library(sjPlot)
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
plot_model(mod.ran,
           type = "re",
           terms = "Nationality")
## [[1]]

## 
## [[2]]

Red are states below average; blue are states above average. The color scheme is in relation to the point estimate, not the confidence interval.

Random slopes

Our random slopes can also be plotted.

We can report the national-level effects of vote shares together with their inercept in sjPlot.

plot_model(mod.slope,
           type = "re",
           terms = "Nationality")
## [[1]]

## 
## [[2]]

Or we can do it our own way:

#simulate
re <- REsim(mod.slope) %>%
  as.data.frame() %>%
  #select national variation; those observations without numbers
  filter(grepl("\\D", groupID ))

tab_slope <-
  re %>%
  #Select my variables
  dplyr::select(groupID, term, mean, sd) %>%
  #Pivot to wide format
  tidyr::pivot_wider(values_from = c(mean, sd),
                     names_from = term) %>%
  dplyr::select(1,2,4,3,5) %>%
  rename("alpha" = 2,
         "se(alpha)" = 3,
         "beta" = 4,
         "se(beta)" = 5) %>%
  as.data.frame()

tab_slope %>%
  arrange(beta) %>%
  ggplot +
  geom_point(aes(y = 1:nrow(tab_slope),
                  x = beta)) +
  geom_errorbar(aes(y = 1:nrow(tab_slope),
                    xmin = beta - `se(beta)`*1.96,
                    xmax = beta + `se(beta)`*1.96)
                ) +
  geom_text(aes(y = 1:nrow(tab_slope),
                x = -0.6,
                label = groupID),
            hjust = 0
                ) +
  geom_vline(aes(xintercept = 0),
             lty = 2) +
  ggtitle("Effect of vote share within member state delegations")

Non-varying slopes

The “fixed”/non-varying effecs of the model can be interpreted and displayes as per usual using either your own simulations or the ggeffects package, for example.

library(ggeffects)

vote_share <- seq(0, 60, 0.1)

dff <- ggpredict(mod.ran,
                 type = "fixed",
                 terms = list(VoteShare_LastElection = vote_share))

dff %>% 
  plot +
  ylim(c(0,4)) +
  ylab("Local staff size")
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

The ggeffects package (as well as the sjPlot) includes also variation from the groups, so the precision of the scenarios we make tend to be surrounded with more uncertainty than if we ran a pooled model.

Regression tables

stargazer() is of course a useful package for model displays. However, sjPlot also offers a table function that provides more model statistics for hierarchical models

library(sjPlot)

tab_model(mod.pool, mod.id, mod.ran,
          file = "sjplot.html")

Literature

Ward, Michael D., and John S. Ahlquist. 2018. Maximum Likelihood for Social Science: Strategies for Analysis. Analytical Methods for Social Research. Cambridge: Cambridge University Press. https://doi.org/10.1017/9781316888544.