Research Methods in Political Science II
Day 1-1: Linear Regression: A Review
Yuki Yanai
October 22, 2016
Preparation
First, let’s load the packages we are going to use.
library('MASS')
library('dplyr')
library('ggplot2')
library('rethinking')For illustration, we use hr96-09.csv, the dataset provided by Asano and Yanai (2013). Download the dataset and put it in the “data” folder in your project. Then, read the dataset by read.csv().
HR <- read.csv('data/hr96-09.csv', na = '.')Check the variables in the data set.
glimpse(HR)## Observations: 5,614
## Variables: 16
## $ year (int) 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996...
## $ ku (fctr) aichi, aichi, aichi, aichi, aichi, aichi, aichi, ai...
## $ kun (int) 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3...
## $ party (int) 1000, 800, 1001, 305, 1014, 1038, 1, 1000, 800, 1001...
## $ name (fctr) KAWAMURA, TAKASHI, IMAEDA, NORIO, SATO, TAISUKE, IW...
## $ age (int) 47, 72, 53, 43, 51, 51, 45, 51, 71, 30, 31, 44, 61, ...
## $ status (int) 2, 3, 2, 1, 1, 1, 1, 2, 3, 1, 1, 1, 1, 1, 1, 1, 2, 1...
## $ nocand (int) 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7...
## $ wl (int) 1, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0, 0, 0, 0, 1, 0, 2...
## $ rank (int) 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3...
## $ previous (int) 2, 3, 2, 0, 0, 0, 0, 2, 1, 1, 0, 0, 0, 0, 0, 1, 3, 1...
## $ vote (int) 66876, 42969, 33503, 22209, 616, 566, 312, 56101, 44...
## $ voteshare (dbl) 40.0, 25.7, 20.1, 13.3, 0.4, 0.3, 0.2, 32.9, 26.4, 2...
## $ eligible (int) 346774, 346774, 346774, 346774, 346774, 346774, 3467...
## $ turnout (dbl) 49.22, 49.22, 49.22, 49.22, 49.22, 49.22, 49.22, 51....
## $ exp (int) 9828097, 9311555, 9231284, 2177203, NA, NA, NA, 1294...Let’s add labels to some variables.
party_names <- c("independent", "JCP", "LDP", "CGP", "oki",
"tai", "saki", "NFP", "DPJ", "SDP",
"LF", "NJSP", "DRF", "kobe", "nii",
"sei", "JNFP", "bunka", "green", "LP",
"RC", "muk", "CP", "NCP", "ND",
"son", "sek", "NP", "NNP", "NPJ",
"NPD", "minna", "R", "H", "others")
HR <- HR %>%
mutate(wl = factor(wl, levels = 0:2, labels = c("lost", "won", "zombie")),
status = factor(status, levels = 1:3,
labels = c("challenger", "incumbent", "ex")),
party = factor(party, labels = party_names))Now, let’s create a dummy variable indicating a cadidate’s prior experience of MP and a variable measuring the electoral expnditure in million yen, extract the observations of the 2009 election, and save them as HR09.
HR09 <- HR %>%
mutate(experience = as.numeric(status == "incumbent" | status == "ex"),
expm = exp / 10^6) %>%
filter(year == 2009)Let’s save this labelled dataset in “data” folder for future use.
save(HR, file = 'data/HR-labelled.Rda')To load this dataset, run load('data/HR-labelled.Rda') (Do not run HR <- load('data/HR-labelled.Rda'))
Take a look at the data set.
glimpse(HR09)## Observations: 1,139
## Variables: 18
## $ year (int) 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 200...
## $ ku (fctr) aichi, aichi, aichi, aichi, aichi, aichi, aichi, a...
## $ kun (int) 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, ...
## $ party (fctr) DPJ, LDP, JCP, SDP, H, DPJ, LDP, JCP, H, DPJ, LDP,...
## $ name (fctr) SATO, YUKO, SHINODA, YOSUKE, KIMURA, EMI, HIRAYAMA...
## $ age (int) 46, 36, 59, 61, 42, 43, 47, 53, 67, 51, 52, 36, 32,...
## $ status (fctr) challenger, incumbent, challenger, challenger, cha...
## $ nocand (int) 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, ...
## $ wl (fctr) won, lost, lost, lost, lost, won, lost, lost, lost...
## $ rank (int) 1, 2, 3, 4, 5, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, ...
## $ previous (int) 1, 1, 0, 0, 1, 5, 0, 0, 0, 5, 1, 0, 0, 4, 1, 2, 0, ...
## $ vote (int) 122348, 78691, 14485, 6082, 3352, 162237, 58225, 18...
## $ voteshare (dbl) 54.4, 35.0, 6.4, 2.7, 1.5, 66.6, 23.9, 7.8, 1.7, 62...
## $ eligible (int) 369526, 369526, 369526, 369526, 369526, 378272, 378...
## $ turnout (dbl) 61.96, 61.96, 61.96, 61.96, 61.96, 65.59, 65.59, 65...
## $ exp (int) 3559575, 9617887, 1802968, 1642095, 1772906, 932959...
## $ experience (dbl) 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, ...
## $ expm (dbl) 3.559575, 9.617887, 1.802968, 1.642095, 1.772906, 9...Linear Regression
We need a statistical model to analyze data and make an inference. As components of a linear regression model, we need:
- outcome variable,
- likelihood function—how the outcome is generated,
- predictor (explanatory) variable(s),
- the relationship between the likelihood and the predictor(s), and
- prior probability (distribution) of the parmeters.
For instance, when we regress the outcome \(y\) on a sigle predictor \(x\), we might construct the following model.
\[y_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i\] \[\alpha \sim \mbox{Normal}(0, 100)\] \[\beta \sim \mbox{Normal}(0, 10)\] \[\sigma \sim \mbox{Cauchy}^{+}(0, 1)\]
We can fit this model to data by rethinking::map() function in R.
Simple Regression
Binary Predictor
Let’s consider a model that explains the vote share (the outcome variable) \(v\) by the exprience of MP (the explanatory variable) \(x\). The experience is a binary (dummy) variable indicting if a candidate has the prior exprience of MP. The variable takes the value of 1 if a candidate is the incumbents or an ex-MP and 0 otherwise.
We can connstruct a statistical model as follows.
\[v_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i\] \[\alpha \sim \mbox{Normal}(0, 100)\] \[\beta \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]
In this model, we assume that the outcome follows the normal distribution whose center is a (linear) function of the explanatory variable. Here, we estimate three paramters \(\alpha\), \(\beta\), and \(\sigma\). In R, we can get MAP (maximum a posteriori) estimates by rethinking::map().
model1 <- alist(
voteshare ~ dnorm(mu, sigma),
mu <- alpha + beta*experience,
alpha ~ dnorm(0, 100),
beta ~ dnorm(0, 100),
sigma ~ dunif(0, 50)
)
fit1 <- map(model1, data = HR09)Now the estimation results have been saved in fit1. We can view the summary of the results by rethinking::precis().
precis(fit1)## Mean StdDev 5.5% 94.5%
## alpha 13.98 0.62 12.99 14.98
## beta 30.86 0.98 29.29 32.43
## sigma 16.24 0.34 15.70 16.79In the output, the column Mean shows the mean of the estimates for each parameter. As shown, we got \(\hat{\alpha} =\) 13.98, \(\hat{\beta}=\) 30.86, and \(\hat{\sigma}=\) 16.24. Thus, the linear relationship between the outcome and the explanatory variable can be summarized as \[\hat{v}_i = 13.98 + 30.86 \cdot x_i.\]
Let’s visualize the result. To display the unceratainty of our inference, we will first sample the partermeter values.
post <- mvrnorm(10000, mu = coef(fit1), Sigma = vcov(fit1)) %>%
as.data.frame()
head(post)## alpha beta sigma
## 1 13.88953 29.24265 15.69308
## 2 13.06442 32.57933 16.15492
## 3 13.92751 32.28971 15.89614
## 4 13.55726 30.70425 16.43296
## 5 14.51539 31.12089 16.35108
## 6 14.11258 29.98933 16.06091Next, we will calculate \(\mu\) using these samples. \(x\) is either 0 or 1.
mu <- sapply(c(0, 1), function(x) post$alpha + post$beta*x)Then, compute the mean of \(\mu\) and the 89% HPDI (highest posterior density interval) corresponding to each value of \(x\) .
mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI)Visualize the result!
df <- data.frame(experience = c(0, 1), mu_mean,
lb = mu_HPDI[1, ], ub = mu_HPDI[2, ])
plt_m1 <- ggplot(HR09, aes(x = experience)) +
geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
geom_jitter(aes(y = voteshare), color = 'royalblue',
position = position_jitter(width = 0.05)) +
geom_line(data = df, aes(y = mu_mean)) +
scale_x_continuous(breaks = c(0, 1))
print(plt_m1)
The intercept of the line 13.98 is the average vote share (the predicted vote share) of an unexprienced candidate. This is evident by plugging 0 into \(x\) of the predicted value formula above. The average (predicted) vote share of an experienced candidate can be obtained by plugging 1 into \(x\): \(13.98 + 30.86 \cdot 1 = 44.84\).
Let’s calculate the average vote share by experience in R and verify that these equal the figures above.
HR09 %>%
filter(experience == 0) %>%
with(mean(voteshare)) %>%
round(2)## [1] 13.98HR09 %>%
filter(experience == 1) %>%
with(mean(voteshare)) %>%
round(2)## [1] 44.84As these numbers show, the predicted values are the mean of the outcome variable given a spcific value of the explanatory variable.
Continours Predictor
Similarly, let’s think about a model that epxlains the vote share \(v\) by the campaign spending \(m\). This model should look like: \[v_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \gamma m_i\] \[\alpha \sim \mbox{Normal}(0, 100)\] \[\gamma \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]
We estimate three paramters: \(\alpha\), \(\gamma\), and \(\sigma\). We fit this model as follows.
model2 <- alist(
voteshare ~ dnorm(mu, sigma),
mu <- alpha + gamma*exp,
alpha ~ dnorm(0, 100),
gamma ~ dnorm(0, 100),
sigma ~ dunif(0, 50)
)
fit2 <- map(model2, data = HR09)## Error in map(model2, data = HR09): initial value in 'vmmin' is not finite
## The start values for the parameters were invalid. This could be caused by missing values (NA) in the data or by start values outside the parameter constraints. If there are no NA values in the data, try using explicit start values.We got an error! When R gives you an error, don’t panic. You should take a breath and calm down. Then, read the error message. In the message above, R tells you that ``[t]his is could be caused by missing values (NA) in the data.’’ So, let’s examinie if there is an NA.
summary(HR09)## year ku kun party
## Min. :2009 tokyo :110 Min. : 1.000 H :292
## 1st Qu.:2009 osaka : 80 1st Qu.: 2.000 LDP :291
## Median :2009 kanagawa: 77 Median : 4.000 DPJ :271
## Mean :2009 aichi : 54 Mean : 5.658 JCP :152
## 3rd Qu.:2009 saitama : 54 3rd Qu.: 8.000 independent: 70
## Max. :2009 chiba : 51 Max. :25.000 SDP : 31
## (Other) :713 (Other) : 32
## name age status nocand
## NISHINO, AKIRA : 2 Min. :25.00 challenger:683 Min. :2.000
## SATO, MASAYUKI : 2 1st Qu.:41.00 incumbent :394 1st Qu.:3.000
## TAKAGI, YOSHIAKI: 2 Median :50.00 ex : 62 Median :4.000
## ABE, MASATO : 1 Mean :50.08 Mean :4.005
## ABE, SHINZO : 1 3rd Qu.:59.00 3rd Qu.:4.000
## ABE, TAKAYUKI : 1 Max. :85.00 Max. :9.000
## (Other) :1130
## wl rank previous vote
## lost :742 Min. :1.000 Min. : 0.000 Min. : 177
## won :300 1st Qu.:1.000 1st Qu.: 0.000 1st Qu.: 5992
## zombie: 97 Median :2.000 Median : 0.000 Median : 62034
## Mean :2.496 Mean : 1.745 Mean : 61940
## 3rd Qu.:3.000 3rd Qu.: 3.000 3rd Qu.:107292
## Max. :9.000 Max. :16.000 Max. :201461
##
## voteshare eligible turnout exp
## Min. : 0.10 Min. :211750 Min. :61.72 Min. : 10024
## 1st Qu.: 2.40 1st Qu.:298265 1st Qu.:66.13 1st Qu.: 1794542
## Median :30.00 Median :352167 Median :68.62 Median : 4809437
## Mean :26.34 Mean :349973 Mean :69.34 Mean : 6118181
## 3rd Qu.:47.30 3rd Qu.:405333 3rd Qu.:72.24 3rd Qu.: 9109114
## Max. :95.30 Max. :487837 Max. :79.86 Max. :25354069
## NA's :15
## experience expm
## Min. :0.0000 Min. : 0.01002
## 1st Qu.:0.0000 1st Qu.: 1.79454
## Median :0.0000 Median : 4.80944
## Mean :0.4004 Mean : 6.11818
## 3rd Qu.:1.0000 3rd Qu.: 9.10911
## Max. :1.0000 Max. :25.35407
## NA's :15It turned out theat there are 15 NAs in the variable exp, which is the predictor of our model. We have to deal with them before we move forward.
Here, let’s plug the mean value of the variable into NAs—this might not the ideal, but one way to overcome the error we got.
mis <- is.na(HR09$exp)
HR09_imp <- HR09
HR09_imp$exp <- ifelse(mis, mean(HR09$exp, na.rm = TRUE), HR09$exp)Let’s fit the model with the imputed data set.
fit2 <- map(model2, data = HR09_imp)
precis(fit2)## Mean StdDev 5.5% 94.5%
## alpha 7.54 0.76 6.33 8.76
## gamma 0.00 0.00 0.00 0.00
## sigma 16.13 0.34 15.59 16.67We got the estimate of zero for gamma. Therefore, we did not find any effect of money of vote share. Well, not exactly.
Let’s display more digits.
precis(fit2, digits = 8)## Mean StdDev 5.5% 94.5%
## alpha 7.54468298 0.7587308 6.33208464 8.75728131
## gamma 0.00000307 0.0000001 0.00000292 0.00000323
## sigma 16.13001504 0.3379543 15.58989887 16.67013121Now we see the tiny but positive effect of gamma. The reason the effect is so small is that \(exp\) is measured by 1 yen and the gamma shows the effect of 1 yean. Naturally, one yen makes a tiny difference at best. When the effect is so small compared to the estimates of other estimates, we should consider rescaling the variable. Here, we will create a new variable exp_m, which measures the campaign spending in one million yen and fit the same model again.
HR09_imp <- HR09_imp %>%
mutate(exp_m = exp / 10^6)
model3 <- alist(
voteshare ~ dnorm(mu, sigma),
mu <- alpha + gamma*exp_m,
alpha ~ dnorm(0, 100),
gamma ~ dnorm(0, 100),
sigma ~ dunif(0, 50)
)
fit3 <- map(model3, data = HR09_imp)
precis(fit3)## Mean StdDev 5.5% 94.5%
## alpha 7.55 0.76 6.33 8.76
## gamma 3.07 0.10 2.92 3.23
## sigma 16.13 0.34 15.59 16.67We got the effect seemingly 1,000,000 times larger than before. Of course, the sunstantive meaning of the result does not differ between two results. However, we can understand the new results more easily. We expect the vote share increases by 3 points as the spending increases by 1 million yen.
Then what is the meaning of 7.55, the estimate of alpha? This is the expected vote share when the spending is 0. Howerver, there are no candidates who did not spend at all. So it’s really hard to understand what 7.55 means.
To get a better estimate, let’s center (de-mean) the predictor and the fit the model again.
HR09_imp <- HR09_imp %>%
mutate(exp_m_c = exp_m - mean(exp_m))
## Let's save this data set for the future
save(HR09_imp, file = 'data/HR09-imputed.Rda')
model4 <- alist(
voteshare ~ dnorm(mu, sigma),
mu <- alpha + gamma*exp_m_c,
alpha ~ dnorm(0, 100),
gamma ~ dnorm(0, 100),
sigma ~ dunif(0, 50)
)
fit4 <- map(model4, data = HR09_imp)
precis(fit4)## Mean StdDev 5.5% 94.5%
## alpha 26.34 0.48 25.57 27.10
## gamma 3.07 0.10 2.92 3.23
## sigma 16.13 0.34 15.59 16.67The estimates for \(\gamma\) and \(\sigma\) are exactly same as the previoud model, but the value of \(\alpha\) estimate has changed. This value now show the expected vote share of the candidate who spent the average amount of money. That is, we expect that a candidate who is an average person in terms of campaign spending receives 26.34% of votes in the district.
Let’s visualize this result (we use fit3 instead of fit4 becaseu uncentered model is easier to understand when we use visual displays). To add the uncertainty of our inference, simulate parameter values first.
post <- mvrnorm(10000, mu = coef(fit3), Sigma = vcov(fit3)) %>%
as.data.frame()
head(post)## alpha gamma sigma
## 1 8.114647 2.909621 15.88761
## 2 7.118386 3.164406 16.19241
## 3 7.741858 3.086369 16.23153
## 4 7.151619 3.213531 16.60039
## 5 7.619360 3.117216 16.55685
## 6 9.034673 2.887549 15.94148Using these values, calculate \(\mu\). We use the range of the observed data.
m <- seq(min(HR09_imp$exp_m_c), max(HR09_imp$exp_m), length.out = 100)
mu <- sapply(m, function(x) post$alpha + post$gamma*x)Compute the mean of \(mu\) and 89% HPDI corresponding to each m.
mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI)Finally, visualize!
df <- data.frame(exp_m = m, mu_mean,
lb = mu_HPDI[1, ], ub = mu_HPDI[2, ])
plt_m3 <- ggplot(HR09_imp, aes(x = exp_m)) +
geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
geom_jitter(aes(y = voteshare), color = 'royalblue',
position = position_jitter(width = 0.05)) +
geom_line(data = df, aes(y = mu_mean)) +
xlab('expenditure (million yen)')
print(plt_m3)
Multiple Regression
According to the results above, it seems that both the experience and the money affect the vote share. Then, why don’t we include both predictors in a single model? Let’s run multiple regression, which has more than one predictor variables in a model.
Esimation
First, we need a statistical model.
\[v_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i + \gamma m_i\] \[\alpha \sim \mbox{Normal}(0, 100)\] \[\beta \sim \mbox{Normal}(0, 100)\] \[\gamma \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]
Now, the linear model (second line) contains two explanatory variable on the rihgt-hand side. We will estimate four parameters: \(\alpha\), \(\beta\), \(\gamma\), and \(\sigma\).
model5 <- alist(
voteshare ~ dnorm(mu, sigma),
mu <- alpha + beta*experience + gamma*exp_m_c,
alpha ~ dnorm(0, 100),
beta ~ dnorm(0, 100),
gamma ~ dnorm(0, 100),
sigma ~ dunif(0, 50)
)
fit5 <- map(model5, data = HR09_imp)
precis(fit5)## Mean StdDev 5.5% 94.5%
## alpha 19.12 0.66 18.07 20.17
## beta 18.02 1.22 16.07 19.98
## gamma 1.86 0.12 1.66 2.05
## sigma 14.78 0.31 14.29 15.28We don’t like tables, so let’s show graphically the result.
# plot(precis(fit5)) should work but doesn't on my computer
precis.plot(precis(fit5))
Let’s interpret the reuslt. First, the estimate of alpha (the intercept) is 19.13. This is the expected vote share when both experience and exp_m_c are zero. In other words, a new candidate spending the average amount of money is expected to gain 19.13% of votes in the district. However, we don’t know if such a candidate really exist.
Next, think about the beta estimate 18.02. This is the difference in vote share between the candidates with and without the prior experience of MP, other thing euqal. That is, when we compare the candidates who spent the exactly same amount of money, the experienced candidates are expected to receive vote shares 18 points higher than the unexperienced candidates. There is no assurance that we can find a pair of experienced and unexperienced candidates whose spendng are the same.
Lastly, the gamma estimate 1.86. This valus tells us how much the vote share increases by one unit of money, keeping the other variables at the same level. That is, comparing the cadidates with same experience, 1 million yen is expected to increase the vote share by 1.86 points.
In sum, the paramter estimates in multiple regression presents the effect of an explanatory variable on the outcome, other variables equal.
Visualizing Multiple Regression
Since we have only two variables–one outcome and one predictor—in simple regression, it is easy to visualize the result in a 2D figure. However, since we have more than two variables–including the outcome–in multiple regression, it is little–only little– harder to visuzlize the results.
There are a few different approaches.
(1) Residual plots
One reason we use multiple regression is to epxalin the part of the variation that is not explained by an explnatory variable.
Suppose we try to explain the outcome \(v\) by the explanatory variable \(x\), but we thinkg that theare are other imoportant causes of \(v\), \(m\). By adding the second epxlanatory variable \(m\) to the model, we can estimate
- how much our prediction would improve if we take into account \(m\), given we already know \(x\) and
- how much our prediction would improve if we take into account \(x\), given we already know \(m\).
The effect of an explanatory variable is the effect that remains after we remove the effect of the other explanatory variables. So let’s remove the effect of an explanatory variable and display the remining effect.
First, we focus on the effect of experience, so we will remove the effect of campaign spending. To do so, we regress experience on spending and removre the effect of spending on experience. Here is a model.
\[x_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = a_1 + b_1 m\] \[a_1 \sim \mbox{Normal}(0, 100)\] \[b_1 \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]
Estimate \(a_1\), \(b_1\), and \(\sigma\).
res_model1_1 <- alist(
experience ~ dnorm(mu, sigma),
mu <- a1 + b1*exp_m_c,
a1 ~ dnorm(0, 100),
b1 ~ dnorm(0, 100),
sigma ~ dunif(0, 50)
)
res_fit1_1 <- map(res_model1_1, data = HR09_imp)Using this result, calculate the predicted vote share (\(\mu\)) of each candidate.
mu <- coef(res_fit1_1)[1] + coef(res_fit1_1)[2]*HR09_imp$exp_m_cUsing these predicted values, calculate the residuals. The residual is the difference between the observd value and the predicted value.
HR09_imp$experience_res <- HR09_imp$experience - muNext, regress the outcome \(v\) on the residual. Here is a model. \[v_i \sim \mbox{Normal}(x_i^{\mathrm{res}}, \sigma)\] \[\mu_i = \rho + \beta x_i^{\mathrm{res}}\] \[\rho \sim \mbox{Normal}(0, 100)\] \[\beta \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]
Estimate parameters.
res_model1_2 <- alist(
voteshare ~ dnorm(mu, sigma),
mu <- rho + beta*experience_res,
rho ~ dnorm(0, 100),
beta ~ dnorm(0, 100),
sigma ~ dunif(0, 50)
)
res_fit1_2 <- map(res_model1_2, data = HR09_imp)
precis(res_fit1_2)## Mean StdDev 5.5% 94.5%
## rho 26.34 0.63 25.33 27.34
## beta 18.00 1.76 15.19 20.81
## sigma 21.23 0.44 20.52 21.95The estimate of beta here matches that of the multiple regression above (up to rounding error)! This is not a coincidence. As stated above, multiple regression is just a simple way to remove the effect of other variables by a sigle shot. We can reach the same results with a series of simple regerssions.
SInce the second step of this procedure is simple regression, we can visualize the results as we did for simple regression.
post <- mvrnorm(10000, mu = coef(res_fit1_2), Sigma = vcov(res_fit1_2)) %>%
as.data.frame()
x <- seq(min(HR09_imp$experience_res), max(HR09_imp$experience_res),
length.out = 100)
mu <- sapply(x, function(x) post$rho + post$beta*x)
mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI) ## 89% HPDI
df <- data.frame(experience_res = x, mu_mean,
lb = mu_HPDI[1,], ub = mu_HPDI[2,])
plt_res1 <- ggplot(HR09_imp, aes(x = experience_res)) +
geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
geom_point(aes(y = voteshare), color = 'royalblue') +
geom_line(data = df, aes(y = mu_mean)) +
xlab('Experinece residuals') +
geom_vline(xintercept = 0, linetype = 'dashed')
print(plt_res1)
In this figure, the vertical dashed line shows where the residual is zero. On the left side are the candidates whose experience is lower than that expected by the spending; on the right are the candidates whose experience is higher than that expected by the spending. The line in this figure shows that the experince increases the vote share.
You should be able to create a similar figure whose horizontal axis is campaign spending.
(2) Marginal Effects
Next, let’s create a figure displaying the marginal effect of an explanatory variable, where we fix the value of other explanatory varialbes at some specific values.
We will show the marginal effect of experience. The variable experience takes either 0 or 1, but we compute the predicted values for the values between 0 and 1 to make a figure.
df <- data.frame(experience = seq(0, 1, length.out = 100))Next, we need to decide the value at which we fix the other predictor. Here, let’s fix the spending at its mean.
df$money_mean <- with(HR09_imp, mean(exp_m_c)) ## this must be 0Using these values, calculate the predicted value \(\mu\).
post <- mvrnorm(10000, mu = coef(fit5), Sigma = vcov(fit5)) %>%
as.data.frame()
mu <- sapply(df$experience,
function(x) post$alpha + post$beta*x + post$gamma*df$money_mean)
# Because money_mean = 0, you can run the following instead.
# mu <- sapply(df$x, function(x) post$alpha + post$beta*x)
df$mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI) # 89% HPDI
df$lb <- mu_HPDI[1,]
df$ub <- mu_HPDI[2,]
plt_x <- ggplot(HR09_imp, aes(x = experience)) +
geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
geom_jitter(aes(y = voteshare), color = 'royalblue',
position = position_jitter(width = 0.05)) +
geom_line(data = df, aes(y = mu_mean)) +
scale_x_continuous(breaks = c(0, 1))
print(plt_x)
We can fix the spending at a difference value. For instance, we can use the minimum and the maximum in addition to the mean.
df$money_min <- with(HR09_imp, min(exp_m_c))
df$money_max <- with(HR09_imp, max(exp_m_c))
mu_min <- sapply(df$experience,
function(x) post$alpha + post$beta*x + post$gamma*df$money_min)
mu_max <- sapply(df$experience,
function(x) post$alpha + post$beta*x + post$gamma*df$money_max)
df$mu_min_mean <- apply(mu_min, 2, mean)
df$mu_max_mean <- apply(mu_max, 2, mean)
mu_min_HPDI <- apply(mu_min, 2, HPDI) # 89% HPDI
mu_max_HPDI <- apply(mu_max, 2, HPDI) # 89% HPDI
df$lb_min <- mu_min_HPDI[1,]
df$ub_min <- mu_min_HPDI[2,]
df$lb_max <- mu_max_HPDI[1,]
df$ub_max <- mu_max_HPDI[2,]
plt_x <- ggplot(HR09_imp, aes(x = experience)) +
geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
geom_ribbon(data = df, aes(ymin = lb_min, ymax = ub_min), fill = 'gray') +
geom_ribbon(data = df, aes(ymin = lb_max, ymax = ub_max), fill = 'gray') +
geom_jitter(aes(y = voteshare), position = position_jitter(width = 0.05)) +
geom_line(data = df, aes(y = mu_mean)) +
geom_line(data = df, aes(y = mu_min_mean), color = 'blue') +
geom_line(data = df, aes(y = mu_max_mean), color = 'red') +
scale_x_continuous(breaks = c(0, 1))
print(plt_x)
In this figure, the black, red, and blue lines presents the predicted vote share when the spending is at the mean, the maximum, and the minimum, respectively.
You should be able to make a similar figure to show the marginal effect of campaign spending.
(3) Posterior Predictive Plots
Finally, let’s examine if our model fits the ovserved data. To do so, we use a scatter plot of posterior predictions versus the observed outvome.
First, we need to calculate (simulate) posterior predictive values using the obseved values of predictors.
predictor <- list()
for(i in 1:nrow(HR09_imp)) {
row_i <- HR09_imp[i,]
predictor[[i]] <- c(row_i$experience, row_i$exp_m_c)
}
mu <- sapply(predictor,
function(x) post$alpha + post$beta*x[1] + post$gamma*x[2])
HR09_imp$mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI) # 89% HPDI
HR09_imp$mu_lb <- mu_HPDI[1,]
HR09_imp$mu_ub <- mu_HPDI[2,]Then, visualize this.
pred_obs <- ggplot(HR09_imp, aes(voteshare, mu_mean)) +
geom_point(size = .5) +
geom_pointrange(aes(ymin = mu_lb, ymax = mu_ub)) +
xlim(0, 100) + ylim(0, 100) +
geom_abline(intercept = 0, slope = 1, linetype = 'dashed') +
labs(x = 'observed vote share', y = 'predicted vote share')
print(pred_obs)
If our model predicted the observations perfectly–which we do not want for reasons explained late, the points will be on the dashed line.