There are three types of missing values

It is always useful to think about the data generating process – especially for missing values – because it tells us how to model them.

Dependent variable: Attendance

I approximate legislative investment by the number of committee meetings an MEP has attended during the semester (period of observation.)

df %>%
  ggplot +
  geom_histogram(aes(y)) +
  ggtitle("Number of committee meetings (y)") +
  xlab("Number of committe meetings")

The variable is a count and naturally bounded at 0 (i.e. none can attend a negative number of meetings). The canoncial model for count data is the poisson moel, of course, but for simplicity, I’ll stick to the linear model here.

For the purposes of this exercise, I am mainly interested in the differences in means between MEPs with an expressed reelection intent and those that don’t have such a clear ambition. This is a comparison of means between categorical variables, so the linear model can be defended. (I could also make the bet that the conditional distribution (and thus my residuals) will be approximately normal once I have specified my full model. If so, I could also run a linear model. )

Independent variable: Reelection intent

The EPRG survey is distributed at the beginning of each legislative term and asks where an MEP sees him/herself in 10 years. One of the possible answers is the “European Parliament”. This is my operationalization of a re-election seeker (FutureInEP==1).

df$future <- ifelse(df$x == 1, "reelection", "exit")
df %>%
  ggplot +
  geom_bar(aes(future)) +
  xlab("Intention to seek reelection") +
  ggtitle("Re-election intent upon election",
          subtitle = "EPRG Survey distributed in 2009 and 2014")

As is clear from the barplot, most MEPs did not reply to the survey. I have a lot of missing data (70 %). To what extent is this a problem?

Moderating variable : electoral system

Our moderating variable is the electoral system in which the MEP would compete if they were to seek reelection. It is a binary variable, and the MEPs are approxiamately equeally distributed accross both systems. That’s good news: we have variation.

df %>%
  mutate(
    ClosedList = if_else(ClosedList == 1, "Party-centered", "Candidate-centered")
  ) %>%
  ggplot +
  geom_bar(aes(ClosedList)) +
  xlab("Electoral system") +
  ggtitle("Electoral system in which MEPs would compete",
          subtitle = "EPRG Survey distributed in 2009 and 2014")

Potential confounder : Sources of missing observations

There are two potential problems with missing data: statistical power (we lose data) and bias. The key question to ask ourselves is whether the assignment mechanism that allocates NAs is systematic. Is it related to x or y, or both? The problem – and opportunities – depend on their structure and our theory about the data generating process (the one that generates NAs).

Ward and Ahlquist are on a mission in chapter 12: We have litte to lose – and a lot to gain – from imputing the missing values. Once we think about the assumptions surrounding most strategies for addressing missing values, imputation is the better option. This only holds as long as we can claim that the data is either MCAR (missing completely at random) or MAR (missing at random). If the data is MNAR (missing not at random), we need to re-theorize and / or look towards selection models (Heckman regression), sensitivity analysis and what not (fun stuff, but beyond this course).

na_response <- mean(is.na(df$FutureInEP)) * 100
na_response

## [1] 70.38596

#Name my variables
vars <- c("y", "x", "ClosedList")

df0 <- 
  #Store results in a new data frame
  df %>%
  #Select them
  dplyr::select(all_of(vars))

library(finalfit)
missing_plot(df0)

library(mice)

miss_pat <- missing_pattern(df0)

miss_pat

##     ClosedList   y    x     
## 395          1   1    1    0
## 923          1   1    0    1
## 27           1   0    1    1
## 80           1   0    0    2
##              0 107 1003 1110

missing_compare(df0,
                dependent = "x",
                explanatory = c("y", "ClosedList"))

#NA in attendance vs. party-centered systems
cor.test(as.numeric(is.na(df$y)), df$ClosedList)

## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(is.na(df$y)) and df$ClosedList
## t = 2.5355, df = 1423, p-value = 0.01134
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.01518664 0.11857821
## sample estimates:
##        cor 
## 0.06706245

#NA in attendance vs. future in the EP (reelection ambition)
cor.test(as.numeric(is.na(df$y)), df$x)

## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(is.na(df$y)) and df$x
## t = 0.63699, df = 420, p-value = 0.5245
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.06458376  0.12615157
## sample estimates:
##        cor 
## 0.03106673

#NA in ambition vs. attendance
cor.test(as.numeric(is.na(df$x)), df$y)

## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(is.na(df$x)) and df$y
## t = -2.5653, df = 1316, p-value = 0.01042
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.12406177 -0.01660529
## sample estimates:
##         cor 
## -0.07053817

Conditional NAs

A quick way to check if missing observations are conditional on other predictors is to regress the NAs on the rest of the predictors.

I begin by recoding my variable into 0 (not missing) and 1 (missing).

df0 <-
  df0 %>%
  mutate(x_na = if_else(is.na(x), 1, 0))

Now, I can estimate a binomial model to see if the data is balanced.

mod.na <- glm(x_na ~ ClosedList,
              family = "binomial",
              df0)
summary(mod.na)

## 
## Call:
## glm(formula = x_na ~ ClosedList, family = "binomial", data = df0)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.63908    0.08186   7.807 5.85e-15 ***
## ClosedList   0.44045    0.11665   3.776 0.000159 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1731.5  on 1424  degrees of freedom
## Residual deviance: 1717.2  on 1423  degrees of freedom
## AIC: 1721.2
## 
## Number of Fisher Scoring iterations: 4

… they are not. The non-responses to the survey are systematically correlated with the electoral system. Here, I find that the probability of a missing observation on x increases by 55.3407789 % when we go from candidate-centered to party-centered systems.

If we assume the data is MAR rather than MNAR, then conditioning on electoral system in the main model would ideally address the bias induced by the missing observations. We can also use the predictors in the model to impute non-responses.

Ward and Ahlquist (2018) also argue that we can use the outcome variable (y) to impute values on x. Their approach is that we have a (theoretical) relationship between the two, and that this is reflected in their correlation. The correlation is symmetric, so if we can predict y using x, then they argue that we can also use y to predict x.

Before we go that far, let’s consider some low-stakes options.

Removing data (listwise exclusion / complete sample analysis)

My ideal model looks like this:

mod.base <- lm(y ~
                 x * ClosedList,
               df0)
summary(mod.base)

## 
## Call:
## lm(formula = y ~ x * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.839  -7.408  -1.358   6.324  47.642 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   19.4082     0.8460  22.942   <2e-16 ***
## x             -1.0500     1.5119  -0.694    0.488    
## ClosedList    -0.4642     1.2479  -0.372    0.710    
## x:ClosedList   2.9452     2.2375   1.316    0.189    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.26 on 391 degrees of freedom
##   (1030 observations deleted due to missingness)
## Multiple R-squared:  0.00507,    Adjusted R-squared:  -0.002564 
## F-statistic: 0.6641 on 3 and 391 DF,  p-value: 0.5745

From the base-line model, we can see that MEPs that already know that they will seek reelection attend on average 3 committee meetings more than their colleagues in party-centered systems without such ambition. The result is not statistically significant, though.

If I do nothing more at the estimation stage, R will proceed to listwise exclusion of all observations with missing values on any of the covariates included in the model. .

mod.base$residuals %>% length

## [1] 395

Here, we see that the length of my data is 395, while my data frame contains 1425 observations. We lost quite a lot of data on the other variables in the regression.

Re-weighing data (regression weights)

A common problem / feature among survey researchers is that samples are not representative of the population / original data. Researchers fix the symptom, by re-balancing the data. They bank on the assumption that the observations that are complete, are representative enough for their group to “stand in” for non-responses.

If the main problem is just over/under-representation of some groups of observations. If so, we can include weights in the regression to compensate for the overrepresenation.

df %>%
  group_by(Nationality) %>%
  summarize(prop = mean(is.na(x))) %>%
  dplyr::select(Nationality, prop) %>%
  ggplot +
  geom_bar(aes(x = prop,
               y = Nationality),
           stat = "identity") +
  ggtitle("Proportion of non-responses by nationality") +
  xlab("")

Distribution of non-responses on the EPRG survey

From the Figure, we see that there is some variation in response rate across nationalities. If that is all, and we don’t think the same mechanism will affect the values taken on by my \(y\), then we can simply weight the observations according to their probability of occurring in the survey. E.g. we can see that Hungary and Spain are overrepresented, while Denmark has a lower response rate.

I begin by creating a dummy indicator variable for whether the MEP responded to the survey.

df$future_na <- df$FutureInEP %>% is.na %>% as.numeric

I then fit a model to predict the probability of a non-response as a function of nationality.

mod.weight <- glm(future_na ~ Nationality, 
                  family = "binomial",
                  df)
summary(mod.weight)

## 
## Call:
## glm(formula = future_na ~ Nationality, family = "binomial", data = df)
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                1.34993    0.42414   3.183  0.00146 **
## NationalityBelgium        -0.58779    0.53353  -1.102  0.27060   
## NationalityBulgaria       -0.94446    0.54676  -1.727  0.08410 . 
## NationalityCroatia        -0.94446    0.77237  -1.223  0.22140   
## NationalityCyprus         -1.16761    0.73930  -1.579  0.11426   
## NationalityCzech Republic -1.19578    0.53204  -2.248  0.02461 * 
## NationalityDenmark        -1.26988    0.58322  -2.177  0.02945 * 
## NationalityEstonia        -0.94446    0.77237  -1.223  0.22140   
## NationalityFinland        -0.71394    0.59146  -1.207  0.22740   
## NationalityFrance         -0.25131    0.47364  -0.531  0.59570   
## NationalityGermany        -0.34011    0.45458  -0.748  0.45435   
## NationalityGreece         -0.21852    0.55871  -0.391  0.69570   
## NationalityHungary         1.24034    0.73364   1.691  0.09090 . 
## NationalityIreland        -1.24457    0.62530  -1.990  0.04655 * 
## NationalityItaly          -0.60271    0.46113  -1.307  0.19120   
## NationalityLatvia         -0.99325    0.65019  -1.528  0.12661   
## NationalityLithuania      -0.30847    0.63670  -0.484  0.62804   
## NationalityLuxembourg     -0.65678    0.74491  -0.882  0.37795   
## NationalityMalta          -1.01345    0.72302  -1.402  0.16100   
## NationalityNetherlands    -0.47446    0.52377  -0.906  0.36502   
## NationalityPoland         -0.50263    0.47698  -1.054  0.29199   
## NationalityPortugal       -0.34662    0.55149  -0.629  0.52966   
## NationalityRomania        -0.85745    0.50312  -1.704  0.08833 . 
## NationalitySlovakia       -1.26988    0.58322  -2.177  0.02945 * 
## NationalitySlovenia       -1.34993    0.65566  -2.059  0.03951 * 
## NationalitySpain          -0.04652    0.48754  -0.095  0.92398   
## NationalitySweden         -0.41562    0.55364  -0.751  0.45283   
## NationalityUnited Kingdom -0.21285    0.46770  -0.455  0.64904   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1731.5  on 1424  degrees of freedom
## Residual deviance: 1681.8  on 1397  degrees of freedom
## AIC: 1737.8
## 
## Number of Fisher Scoring iterations: 5

df$weight <- predict(mod.weight, df, "response")

The weight variable multiplies in practice our dependent variable to increase or decrease according to weight it is given. Now, I weigh up the observations from under-represented groups and down the over-represented groups.

mod.w <- lm(Attendance ~ 
              x * ClosedList,
            weight = 1/weight,
            df)

summary(mod.w)

## 
## Call:
## lm(formula = Attendance ~ x * ClosedList, data = df, weights = 1/weight)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.625  -9.453  -1.720   7.607  64.608 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   19.4783     0.8397  23.197   <2e-16 ***
## x             -0.8878     1.4916  -0.595    0.552    
## ClosedList    -0.4212     1.2888  -0.327    0.744    
## x:ClosedList   2.9464     2.2863   1.289    0.198    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.78 on 391 degrees of freedom
##   (1030 observations deleted due to missingness)
## Multiple R-squared:  0.00509,    Adjusted R-squared:  -0.002543 
## F-statistic: 0.6668 on 3 and 391 DF,  p-value: 0.5728

This is a common way to address some of the information biases we find in surveys such as the European Social Survy, World Values Survey etc.

stargazer(mod.base, mod.w,
          #Esthetics: redefine the dependent variable and column names
          dep.var.caption = "Committee attendance",
          dep.var.labels.include = F,
          model.numbers = F,
          column.labels = c("Listwise excl", "Weights"),
          title = "Listwise exclusion and regression weights.",
          type = "html")

Mean imputation

I can simply replace my missing observations by the mean of the variable. Here, the “mean” is the proportion of MEPs that answer the survey question by the positive (x == 1).

df0 <- 
  df0 %>%
  mutate(x_bar = if_else(is.na(x), mean(x, na.rm = T), x))

df0$x_bar %>% hist

The histogram shows what I just did: I now have my original observations with 0s and 1s, but I also have a big group of “imputed” observations: the variable mean.

I can run my model with this new “imputed” variable. There is no new information in the model; to the contrary, I am “watering down” the correlation with noise. My N is higher, however, meaning that the standard errors will be artificially low. I.e. I’ll be underpredicting while being overly certain about the result.

mod.m <- lm(y ~ 
              x_bar * ClosedList,
            df0)

summary(mod.m)

## 
## Call:
## lm(formula = y ~ x_bar * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.479  -8.022  -1.206   6.380  74.794 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       18.3558     0.6645  27.623   <2e-16 ***
## x_bar             -1.0602     1.5988  -0.663    0.507    
## ClosedList        -0.7361     0.9540  -0.772    0.440    
## x_bar:ClosedList   2.9194     2.3659   1.234    0.217    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.85 on 1314 degrees of freedom
##   (107 observations deleted due to missingness)
## Multiple R-squared:  0.001267,   Adjusted R-squared:  -0.001013 
## F-statistic: 0.5557 on 3 and 1314 DF,  p-value: 0.6444

mod.ind <- lm(y ~ 
              + x_na +
              + x_bar * ClosedList,
              df0)
summary(mod.ind)

## 
## Call:
## lm(formula = y ~ +x_na + +x_bar * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.763  -7.762  -1.436   6.238  75.238 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       19.4714     0.7890  24.679  < 2e-16 ***
## x_na              -1.7045     0.6534  -2.609  0.00919 ** 
## x_bar             -1.0493     1.5953  -0.658  0.51079    
## ClosedList        -0.6018     0.9533  -0.631  0.52794    
## x_bar:ClosedList   2.9426     2.3607   1.246  0.21281    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.82 on 1313 degrees of freedom
##   (107 observations deleted due to missingness)
## Multiple R-squared:  0.006417,   Adjusted R-squared:  0.00339 
## F-statistic:  2.12 on 4 and 1313 DF,  p-value: 0.07615

stargazer(mod.base, mod.w, mod.m, mod.ind, #mod.sam,
          dep.var.labels.include = F,
          dep.var.caption = "Attendance",
          column.labels = c("listwise exclusion", "weights", "mean imp.", "mean imp. + indicator"),
          type = "html")

Information from other covariates

df0 <- 
  df0 %>%
  group_by(ClosedList) %>%
  mutate(x_bar_g = if_else(is.na(x), mean(x, na.rm = T), x))

mod.g <- lm(y ~ 
              x_bar_g * ClosedList,
            df0)

summary(mod.g)

## 
## Call:
## lm(formula = y ~ x_bar_g * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.571  -8.016  -1.198   6.388  74.802 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         18.3871     0.6753  27.227   <2e-16 ***
## x_bar_g             -1.1271     1.5981  -0.705    0.481    
## ClosedList          -0.7753     0.9497  -0.816    0.414    
## x_bar_g:ClosedList   3.0864     2.3652   1.305    0.192    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.85 on 1314 degrees of freedom
##   (107 observations deleted due to missingness)
## Multiple R-squared:  0.001407,   Adjusted R-squared:  -0.0008732 
## F-statistic: 0.617 on 3 and 1314 DF,  p-value: 0.604

Including new covariates

Let’s instead use an entirely new covariate.

In fact, let’s add several predictors to the data set we use. My candidate predictors are the legislator’s Nationality, EPGroup (parliamentary group), EP (legislative period/survey round) and Age when they responded to the survey.

#Center age (for convenience later)
df <- 
  df %>%
  mutate(
    Age_c = Age - 40,
    Age2_c = Age_c^2
  )

##Select variables
vars <- c("y", "x", "ClosedList",
          "EPGroup", 
          "Nationality", 
          "Age", "Age_c", "Age2_c",
          "position", "lrecon")

df0 <- df %>%
  dplyr::select(all_of(vars))

Usually, I’d check the patterns of missing data once again (i.e. are there missing observations on the new covariates?).

df0 <- 
  df0 %>%
  group_by(Nationality) %>%
  mutate(
    x_bar_g2 = if_else(is.na(x), mean(x, na.rm = T), x)) %>%
  #make sure to ungroup (useful for some of the sampling functions)
  ungroup()

mod.g2 <- lm(y ~ 
              x_bar_g2 * ClosedList,
            df0)

summary(mod.g2)

## 
## Call:
## lm(formula = y ~ x_bar_g2 * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.744  -8.027  -1.277   6.405  74.827 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          18.2853     0.6447  28.362   <2e-16 ***
## x_bar_g2             -0.8313     1.5065  -0.552    0.581    
## ClosedList           -0.6906     0.8988  -0.768    0.442    
## x_bar_g2:ClosedList   2.9393     2.2243   1.321    0.187    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.84 on 1314 degrees of freedom
##   (107 observations deleted due to missingness)
## Multiple R-squared:  0.001561,   Adjusted R-squared:  -0.0007184 
## F-statistic: 0.6848 on 3 and 1314 DF,  p-value: 0.5613

stargazer(mod.base,  mod.g, mod.g2,
          column.labels = c("listwise exclusion", "group mean/ClosedList", "group mean/Nationality"),
          type = "html")

mod.imp.age <- glm(x ~ Age_c + Age2_c,
                   family = "binomial",
                   df0)

summary(mod.imp.age)

## 
## Call:
## glm(formula = x ~ Age_c + Age2_c, family = "binomial", data = df0)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.299117   0.166678  -1.795   0.0727 .  
## Age_c        0.054337   0.025004   2.173   0.0298 *  
## Age2_c      -0.004572   0.001119  -4.085  4.4e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 525.96  on 421  degrees of freedom
## Residual deviance: 487.40  on 419  degrees of freedom
##   (1003 observations deleted due to missingness)
## AIC: 493.4
## 
## Number of Fisher Scoring iterations: 5

df0$pred.age <- predict(mod.imp.age, df0, type = "response")
df0 <-
  df0 %>%
  mutate(x.age = if_else(is.na(x), pred.age, x))

mod.age <- lm(y ~ 
             + x.age * ClosedList,
            df0)

summary(mod.age)

## 
## Call:
## lm(formula = y ~ +x.age * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.523  -7.942  -1.247   6.262  75.184 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       18.5225     0.6369  29.082   <2e-16 ***
## x.age             -1.5927     1.4834  -1.074    0.283    
## ClosedList        -0.7840     0.9046  -0.867    0.386    
## x.age:ClosedList   3.0683     2.1582   1.422    0.155    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.84 on 1314 degrees of freedom
##   (107 observations deleted due to missingness)
## Multiple R-squared:  0.001618,   Adjusted R-squared:  -0.0006612 
## F-statistic: 0.7099 on 3 and 1314 DF,  p-value: 0.5461

mod.imp.epgr <- glm(x ~ Age_c + Age2_c
                    +  EPGroup
                    ,
                   family = "binomial",
                   df)
#Check the results
# summary(mod.imp.epgr)

df0$pred.epgr <- predict(mod.imp.epgr, df0, type = "response")
df0 <-
  df0 %>%
  mutate(x_gr = if_else(is.na(x), pred.epgr, x))

mod.gr <- lm(y ~ 
             + x_gr * ClosedList,
            df0)

summary(mod.gr)

## 
## Call:
## lm(formula = y ~ +x_gr * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.793  -7.945  -1.307   6.370  75.369 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      18.3403     0.6147  29.836   <2e-16 ***
## x_gr             -1.0334     1.4144  -0.731   0.4651    
## ClosedList       -0.8833     0.8723  -1.013   0.3114    
## x_gr:ClosedList   3.3691     2.0246   1.664   0.0963 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.84 on 1314 degrees of freedom
##   (107 observations deleted due to missingness)
## Multiple R-squared:  0.002448,   Adjusted R-squared:  0.0001704 
## F-statistic: 1.075 on 3 and 1314 DF,  p-value: 0.3587

library(lme4)
mod.imp.ran <- glmer(x ~ Age_c + Age2_c
                    +  (1|EPGroup),
                   family = "binomial",
                   df0)

summary(mod.imp.ran)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: x ~ Age_c + Age2_c + (1 | EPGroup)
##    Data: df0
## 
##      AIC      BIC   logLik deviance df.resid 
##    485.3    501.5   -238.7    477.3      418 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1683 -0.7179 -0.4391  0.9765  7.4477 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  EPGroup (Intercept) 0.3495   0.5912  
## Number of obs: 422, groups:  EPGroup, 8
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.596654   0.294265  -2.028   0.0426 *  
## Age_c        0.047631   0.026319   1.810   0.0703 .  
## Age2_c      -0.004556   0.001157  -3.938 8.21e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) Age_c 
## Age_c  -0.209       
## Age2_c  0.042 -0.888
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00296411 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

df0$pred.ran <- predict(mod.imp.ran, df0, type = "response", allow.new.levels = T)
df0 <-
  df0 %>%
  mutate(x.ran = if_else(is.na(x), pred.ran, x))

mod.ran <- lm(y ~ 
             + x.ran * ClosedList,
            df0)

summary(mod.ran)

## 
## Call:
## lm(formula = y ~ +x.ran * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.781  -7.989  -1.316   6.376  75.371 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       18.3368     0.6216  29.500   <2e-16 ***
## x.ran             -1.0212     1.4446  -0.707    0.480    
## ClosedList        -0.8685     0.8829  -0.984    0.325    
## x.ran:ClosedList   3.3340     2.0755   1.606    0.108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.84 on 1314 degrees of freedom
##   (107 observations deleted due to missingness)
## Multiple R-squared:  0.002277,   Adjusted R-squared:  -5.978e-07 
## F-statistic: 0.9997 on 3 and 1314 DF,  p-value: 0.3921

stargazer(mod.base,  mod.gr, mod.ran,
          column.labels = c("listwise exclusion", "fixeff", "raneff"),
          type = "html")

Random sampling

Absent any new information, we can of course fill in the gaps using information from the model, but without any assumptions about a trand; that is, we can do random sampling. It has the advantage that it allows us to inflate the standard error, reflecting the uncertainty in the estimation of x.

I can use the sample() function to draw randomly from the observed part of my x.

sam <- sample(
  #Draw values that are not NA
  df0$x[!is.na(df0$x)],
  #Make as many draws as I have data
  size = nrow(df0),
  #For reach draw, replace so that I can draw anew
  replace = T)

#Replace the missing observations by my draws

df0 <- 
  df0 %>%
  mutate(x.sam = if_else(is.na(x), sam, x))

Note that this random part means that I’ll get a new result every time I run my R-code. Also, it won’t help me make much headway if I have a lot of missing observations and I am actually interested in the effect of this particular variable.

mod.sam <- lm(y ~ 
              x.sam * ClosedList,
            df0)

summary(mod.sam)

## 
## Call:
## lm(formula = y ~ x.sam * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.357  -8.190  -1.190   6.134  74.766 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       17.8656     0.5270  33.901   <2e-16 ***
## x.sam              0.4912     0.9325   0.527    0.598    
## ClosedList         0.3242     0.7198   0.450    0.652    
## x.sam:ClosedList  -0.4469     1.2978  -0.344    0.731    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.85 on 1314 degrees of freedom
##   (107 observations deleted due to missingness)
## Multiple R-squared:  0.0002821,  Adjusted R-squared:  -0.002 
## F-statistic: 0.1236 on 3 and 1314 DF,  p-value: 0.9462

Sampling from a prediction

We can also sample from a prediction if we know the distribution of this prediction. One way is to assume a normal distribution of each prediction.

##Missing values in x are related to the value of other xs
mod2.na <- glm(x ~ 
               + Age_c
               + Age2_c
               + position
               + EPGroup,
               family="binomial",
               df0
)
summary(mod2.na)

## 
## Call:
## glm(formula = x ~ +Age_c + Age2_c + position + EPGroup, family = "binomial", 
##     data = df0)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       -2.366516   1.096786  -2.158   0.0310 * 
## Age_c              0.037068   0.029898   1.240   0.2150   
## Age2_c            -0.004099   0.001283  -3.195   0.0014 **
## position           0.313443   0.169086   1.854   0.0638 . 
## EPGroupECR         0.965138   0.740880   1.303   0.1927   
## EPGroupEFDD       -0.730251   1.285852  -0.568   0.5701   
## EPGroupEPP         0.698014   0.402373   1.735   0.0828 . 
## EPGroupGreens/EFA -1.918548   1.095199  -1.752   0.0798 . 
## EPGroupGUE/NGL     0.433698   0.903595   0.480   0.6312   
## EPGroupNI          1.400315   1.017342   1.376   0.1687   
## EPGroupPES         0.389912   0.418219   0.932   0.3512   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 458.61  on 369  degrees of freedom
## Residual deviance: 398.12  on 359  degrees of freedom
##   (1055 observations deleted due to missingness)
## AIC: 420.12
## 
## Number of Fisher Scoring iterations: 5

I ask for the standard error of the predictions in order to know the spread of the distribution.

preds <- predict(mod2.na, df0, "response", se.fit = T) 

df0$x.glm <- preds$fit
df0$x.se <- preds$se.fit

#Draw from a normal distribution with a mean an a standard deviation
sims <- rnorm(nrow(df0), df0$x.glm, df0$x.se)

#Replace the missing observations by my draws

df0 <- 
  df0 %>%
  mutate(x_sim = if_else(is.na(x), sims, x))

mod.sim <- lm(y ~ 
              x_sim * ClosedList,
            df0)

summary(mod.sim)

## 
## Call:
## lm(formula = y ~ x_sim * ClosedList, data = df0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.319  -7.935  -1.131   6.278  75.158 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       18.0125     0.6196  29.072   <2e-16 ***
## x_sim             -0.2902     1.4190  -0.205    0.838    
## ClosedList        -0.6291     0.8702  -0.723    0.470    
## x_sim:ClosedList   3.2259     2.0178   1.599    0.110    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.84 on 1231 degrees of freedom
##   (190 observations deleted due to missingness)
## Multiple R-squared:  0.003699,   Adjusted R-squared:  0.001271 
## F-statistic: 1.523 on 3 and 1231 DF,  p-value: 0.2067

stargazer(mod.base,  mod.sam, mod.sim,
          column.labels = c("listwise exclusion", "random sample", "simulation"),
          type = "html")

library(sjPlot)
plot_models(mod.base, mod.g, mod.sim)

replace_names <- function(x,
                          vars = c("x", "ClosedList", "x:ClosedList")){
  
  names(x$coefficients)[-1] <- vars
  return(x)
}

mod.m <- replace_names(mod.m)
mod.g <- replace_names(mod.g)
mod.sam <- replace_names(mod.sam)
mod.sim <- replace_names(mod.sim)

plot_models( mod.base,
             mod.m,
             mod.g,
             mod.sam,
             mod.sim,
             m.labels = c("Listwise exclusion",
                          "Mean imputation",
                          "Group-mean imputation",
                          "Random sampling",
                          "Simulation"),
             legend.title = "Model approaches",
             title = "Coeficients from four strategies of NA handling")

stargazer(mod.base, mod.m, mod.g, mod.sam, mod.sim,
          column.labels = c("Base-line", "Mean repl.", "Group repl.", "Random sample", "Imputation"),
          type = "html")

Mice (Multiple Imputation for missing data using Chained Equations)

We can go all in on this approach and simulate x-s, store them in several data frames, run separate models on each of them before we average over the coefficients. This is what mice does.

library(mice)

library(mice)

df0 <- df %>% 
  #Quadric term for age
  mutate(x = as.factor(x),
         ClosedList = as.factor(ClosedList)) %>%
  #Select variables
  dplyr::select(y, x, ClosedList, position,
                EPGroup, Age_c, Age2_c) %>%
  #Recode to factor
  mutate(EPGroup = as.factor(EPGroup))

df0.imp <- mice(df0,
                m = 10,
                # #linear outcome
                # method = c("norm", 
                #            #Binary variables as factors
                #            "logreg", "logreg",
                #            #Metric
                #            "norm",
                #            #Categorical, unordered
                #            "polyreg",
                #            "norm", "norm")
                )

## 
##  iter imp variable
##   1   1  y  x  position
##   1   2  y  x  position
##   1   3  y  x  position
##   1   4  y  x  position
##   1   5  y  x  position
##   1   6  y  x  position
##   1   7  y  x  position
##   1   8  y  x  position
##   1   9  y  x  position
##   1   10  y  x  position
##   2   1  y  x  position
##   2   2  y  x  position
##   2   3  y  x  position
##   2   4  y  x  position
##   2   5  y  x  position
##   2   6  y  x  position
##   2   7  y  x  position
##   2   8  y  x  position
##   2   9  y  x  position
##   2   10  y  x  position
##   3   1  y  x  position
##   3   2  y  x  position
##   3   3  y  x  position
##   3   4  y  x  position
##   3   5  y  x  position
##   3   6  y  x  position
##   3   7  y  x  position
##   3   8  y  x  position
##   3   9  y  x  position
##   3   10  y  x  position
##   4   1  y  x  position
##   4   2  y  x  position
##   4   3  y  x  position
##   4   4  y  x  position
##   4   5  y  x  position
##   4   6  y  x  position
##   4   7  y  x  position
##   4   8  y  x  position
##   4   9  y  x  position
##   4   10  y  x  position
##   5   1  y  x  position
##   5   2  y  x  position
##   5   3  y  x  position
##   5   4  y  x  position
##   5   5  y  x  position
##   5   6  y  x  position
##   5   7  y  x  position
##   5   8  y  x  position
##   5   9  y  x  position
##   5   10  y  x  position

df0.imp %>% names

##  [1] "data"            "imp"             "m"               "where"          
##  [5] "blocks"          "call"            "nmis"            "method"         
##  [9] "predictorMatrix" "visitSequence"   "formulas"        "post"           
## [13] "blots"           "ignore"          "seed"            "iteration"      
## [17] "lastSeedValue"   "chainMean"       "chainVar"        "loggedEvents"   
## [21] "version"         "date"

df0.imp$imp$x %>% head

mod <- lm(y ~ x * ClosedList,
          data = complete(df0.imp, factor = 1))
summary(mod)

## 
## Call:
## lm(formula = y ~ x * ClosedList, data = complete(df0.imp, factor = 1))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.514  -7.981  -1.242   6.486  75.019 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     18.2421     0.5118  35.641   <2e-16 ***
## x1              -0.7513     0.8906  -0.844    0.399    
## ClosedList1     -0.2615     0.6980  -0.375    0.708    
## x1:ClosedList1   1.2847     1.2164   1.056    0.291    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.76 on 1421 degrees of freedom
## Multiple R-squared:  0.0008479,  Adjusted R-squared:  -0.001261 
## F-statistic: 0.402 on 3 and 1421 DF,  p-value: 0.7516

mod <- with(df0.imp, 
            lm(y ~ 
                x
                 * ClosedList)
            )

names(mod)

## [1] "call"     "call1"    "nmis"     "analyses"

res <- pool(mod) %>%
  summary
res

coefs <- list(summary(mod.base)$coefficients[,1],
              res$estimate)

ses <- list(summary(mod.base)$coefficients[,2],
           res$std.error)

stargazer(mod.base, mod$analyses[[1]],
          coef = coefs,
          se = ses,
          title = "Initial results from mice",
          type = "html")

dfNA <- is.na(df0)
#Remove the dependent variable
dfNA[,1] = FALSE

#Recode to categorical variable
df0$x <- as.factor(df0$x)

df1.imp <- mice(df0,
                m = 10,
                #linear outcome
                # method = "norm",
                #Where are the NAs
                where = dfNA,
                printFlag = F
)

mod1 <- with(df1.imp, 
             lm(y ~ 
                  x
                * ClosedList)
)
res1 <- pool(mod1)

coefs <- list(summary(mod.base)$coefficients[,1],
              res$estimate,
              res1$estimate)

ses <- list(summary(mod.base)$coefficients[,2],
           res$std.error,
           res1$std.error)

stargazer(mod.base, mod$analyses[[1]],mod1$analyses[[1]],
          coef = coefs,
          se = ses,
          title = "Initial results from mice",
          type = "html")

#All NAs
dfNA <- is.na(df0)
#Only two first variables are missing/to be imputed
dfNA[, 3:5] = FALSE

#Two blocks
blocks <- make.blocks(c( "y", "x"))

#Matrix with predictors
predictorMatrix = make.predictorMatrix(df0,
                                       blocks = blocks)
predictorMatrix

##   y x ClosedList position EPGroup Age_c Age2_c
## y 0 1          1        1       1     1      1
## x 1 0          1        1       1     1      1

df2.imp <- mice(df0,
                m = 10,
                #linear outcome
                method = c("norm", "logreg"),
                #Blocks of variables
                blocks = blocks,
                #Data frame specifying where imputations should be made
                where = dfNA,
                #Matrix flagging the predictors
                predictorMatrix = predictorMatrix,
                printFlag = F
)


mod2 <- with(df2.imp, 
            lm(y ~ 
                  x * ClosedList
)
            )
res2 <- pool(mod2) %>%
  summary
res2

library(ggeffects)

#Create an empty object
predictions = NULL

#Index from 1 to 10 (number of datasets/models)
  i = 1
  
  for(i in 1:10){  
    m <- lm(y ~ 
             x * ClosedList,
              data = complete(df2.imp, action = i))
    
  predictions[[i]] <-  ggpredict(m, 
                                 terms = list("x" = c(0,1),
                                              "ClosedList" = c(0,1))
  )
  }

effs <-  pool_predictions(predictions)

effs %>%
  plot +
  ylab("Number of meetings attended") +
  xlab("Election-seeker") +
  ggtitle("Effect of ambition on committee attendance")

In-sample predictions

You can do your usual cross-correlations between observed and predicted values for the imputed variable.

When using mice, there are a few additional functions. For continuous variables, we can plot the density (distribution) of the 10 simulated data (in red) against the observed density (blue).

densityplot(df2.imp, ~ y)

The xyplot() plots the imputed variable agains whatever covariate you have used to impute; just to check if the two overlap appropriately.

xyplot(df2.imp, y ~ EPGroup )

df0$x_miss <- df0$x
id <- sample(row.names(df0 %>% 
                         filter(!is.na(x))), 
             100)
df0$x_miss[row.names(df0) %in% id ] <- NA

#All NAs
dfNA <- is.na(df0)
#Only two first variables are missing/to be imputed
dfNA[, 2:6] = FALSE
dfNA %>% head

##          y     x ClosedList position EPGroup Age_c Age2_c x_miss
## [1,] FALSE FALSE      FALSE    FALSE   FALSE FALSE  FALSE   TRUE
## [2,] FALSE FALSE      FALSE    FALSE   FALSE FALSE  FALSE   TRUE
## [3,] FALSE FALSE      FALSE    FALSE   FALSE FALSE  FALSE   TRUE
## [4,]  TRUE FALSE      FALSE    FALSE   FALSE FALSE  FALSE   TRUE
## [5,] FALSE FALSE      FALSE    FALSE   FALSE FALSE  FALSE  FALSE
## [6,] FALSE FALSE      FALSE    FALSE   FALSE FALSE  FALSE   TRUE

#Two blocks
blocks <- make.blocks(c( "y", "x_miss"))

#Matrix with predictors
predictorMatrix = make.predictorMatrix(df0,
                                       blocks = blocks)

#Remove x as a predictor
predictorMatrix[,"x"] = 0
predictorMatrix

##        y x ClosedList position EPGroup Age_c Age2_c x_miss
## y      0 0          1        1       1     1      1      1
## x_miss 1 0          1        1       1     1      1      0

#Recode to factors
df0 <- 
  df0 %>%
  mutate(EPGroup = as.factor(EPGroup),
         x_miss = as.factor(x_miss))

df2.amp <- mice(df0,
                m = 10,
                #linear outcome
                method = c("norm", "logreg"),
                #Blocks of variables
                blocks = blocks,
                #Data frame specifying where imputations should be made
                where = dfNA,
                #Matrix flagging the predictors
                predictorMatrix = predictorMatrix,
                printFlag = F
)

test <- mice::complete(df2.amp, factor = 3, include = T)

table(test$x, test$x_miss)

##    
##       0   1
##   0 543   6
##   1   2 253