examples.Rmd

---
title: "Examples"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Examples}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

I provide three examples on how the MMMM can be employed to study spatial, network, and aggregation problems

- [Example 1: The effect of air quality on home values](#e1) (spatial regression)
- [Example 2: All friends or just your best friend?](#e2) (network regression)
- [Example 3: The effect of political parties' financial dependency on the survival coalition governments](#e3) (aggregation regression)

### The effect of air quality on home values {#e1}

The R file of this example can be found [here](https://benrosche.github.io/rmm/articles/examples_files/mmmm_spatial.Rmd).

This **spatial analysis** example is based on [Harrison & Rubenfield 1978](https://www.sciencedirect.com/science/article/abs/pii/0095069678900062),
who study the effect of air quality on home values, and a re-analysis by [Bivand 2017](https://openjournals.wu.ac.at/region/paper_107/107.html).

The study employs census tract data from the Boston Standard Metropolitan Statistical Area in 1970. With tracts containing no housing units or comprised entirely of institutions excluded, the sample contains 506 census tracts. Air quality is measured by the concentration of nitric oxides in the air, which is obtained from a meteorological model (Transportation and Air Shed Simulation Model). I refer to their paper for more information on data and operationalization. 

Let us load in the data and plot the home values across town:

```{r, message = FALSE, warning = FALSE}

library(spData)     # spatial datasets
library(sf)         # read spatial datasets
library(spdep)      # create spatial weights
library(spatialreg) # spatial regression models
library(dplyr)
library(ggplot2)
library(rmm)

# Load spatial data from Boston
boston <- 
  read_sf(system.file("shapes/boston_tracts.shp", package = "spData")) %>%
  select(CMEDV, NOX, CRIM, RM, DIS, AGE, LSTAT, geometry) %>% 
  st_transform(crs = 5070) %>% # use Albers equal-area conic projection 
  mutate(tid = row_number(), lnCMEDV=log(CMEDV), across(c(NOX, CRIM, RM, DIS, AGE, LSTAT), scale)) %>% 
  relocate(tid, CMEDV, lnCMEDV) 

# Dependent variable:
# lnCMEDV = ln(median home value in $1000)

# Plot median home values across Boston
ggplot(boston, aes(fill = CMEDV)) +
  geom_sf(color = NA) +
  labs(fill = "Median home value") +
  scale_fill_viridis_b() +
  theme(legend.position = "bottom")
# Explanatory variables (all have been standardized):
# NOX     = nitric oxides concentration 
# CRIM    = per capita crime
# RM      = avg. number of rooms per dwelling
# DIS     = weighted distance to five Boston employment centers
# AGE     = proportion of units built prior 1940 
# LSTAT   = percentage working-class population
```

Let us first consider three canonical spatial regression models:

1. **Residual spatial effect**: $Y = XB + \lambda Wu + \epsilon$. This model is labeled *spatial error model* because the residuals of other spatial units enter the regression equation.
2. **Exogenous spatial effects**: $Y = XB + WXB + \epsilon$. This model is labeled the *spatial lag-x model* because the covariates of other spatial units enter the regression equation.
3. **Endogenous spatial effect**: $Y = \rho WY + XB + \epsilon$. This model is labeled the *spatial autoregressive model* because the outcomes of other spatial units enter the regression equation.

More information on those models and combinations of them can be found in [Gibbons, Overman & Patacchini 2015](https://www.sciencedirect.com/science/article/pii/B9780444595171000039). 

To specify any of those models, we have to create a *spatial weight matrix* $W$, which imposes a structure in terms of what are the neighbors for each location. $W$ is a NxN table, where N equals the number of neighborhoods. Each weight $w_{ij}$ represents the relationship between location $i$ and $j$ and, by convention, $w_{ij}=0$ for diagonal elements. The relationship can either be based on direct contiguity and therefore binary (1=neighbor, 0=otherwise), or based on distance and therefore continuous. Continuous weights are often row-standardized so that the weights of all neighbors j of location i sum to 1.

The R package `spatialreg` takes the weight matrix as list:

```{r}
# Create row-standardized weight matrix
boston_nb <- poly2nb(as_Spatial(boston), row.names = boston$tid) # from polygon list to neighbor list

boston_wmat  <- nb2mat(boston_nb, zero.policy = TRUE) %>% as.matrix() # weight matrix
boston_wlist <- nb2listw(boston_nb, style = "W") # weight matrix as list
```

Now let us estimate the three models:

```{r}

# Residual spatial effect (spatial error model):
mod1 <- errorsarlm(lnCMEDV ~ NOX + CRIM + RM + DIS + AGE,
                   data = boston,
                   listw = boston_wlist)
summary(mod1)

# Exogenous spatial effect (spatial lag-x model):  
mod2 <- lmSLX(lnCMEDV ~ NOX + CRIM + RM + DIS + AGE,
              data = boston,
              listw = boston_wlist,
              Durbin =  ~ NOX + CRIM + RM + DIS + AGE) 
summary(mod2)

# Endogenous spatial effect (spatial autoregressive model):
mod3 <- lagsarlm(lnCMEDV ~ NOX + CRIM + RM + DIS + AGE,
                 data = boston,
                 listw = boston_wlist)
summary(mod3)

```

**Results**:

1. The spatial error model estimates that after conditioning on those five covariates, the residual is still spatially correlated, as $\lambda=`r round(mod1$lambda,2)`$.

2. The spatial lag-x model estimates an exogenous spatial effect for each of the considered covariates. The effect of air quality of neighboring locations, for instance, is estimated to be $\beta=`r round(mod2$coefficients['lag.NOX'],2)`$. That is, the home values of a given neighborhood increases if the air quality of surrounding neighborhoods decreases.

3. The spatial autoregessive (SAR) model estimates a endogenous spatial effect of $\rho=`r round(mod3$rho,2)`$. That is, the SAR summarizes the spatial dependency in one coefficient.

#### The MMMM for spatial analysis

Let $y_{i}$ be the outcome of location $i$. Using the MMMM, we can model this outcome in terms of (i) effects of the location's own features
$x_{i}^{\intercal}\beta+\epsilon_{i}$ and (ii) effects its neighbor's features $\sum_{j \in n(i)}w_{j}(z_{j}^{\intercal}\gamma+u_{j})$, where $j$ indexes the neighbors, $n(i)$ is the set neighbors of location $i$, $z_{j}$ represents the observed features of neighbor $j$, and $u_{j}$ represents the combined influence of unobserved features.

This model is almost identical to the combination of a spatial lag and a spatial error model. 
The only difference is that the error of each location is split into a random effect for its role as focal location and a random effect for its role as neighbor. 

To me, the combination of exogenous spatial effects and spatial error makes the most sense as the entire right-hand side of a neighbor affects a location. 
The endogenous spatial effect model is more difficult to interpret and runs into identification problems 
when both endogenous and exogenous effects are included (spatial Durbin model).

To estimate a spatial MMMM, the neighbors of each location must be included in the dataframe as individual rows:

```{r}
# Neighbor list to data.frame
nb2df <- function(nb) {
  return(
    data.frame(tid = unlist( mapply(rep, 1:length(nb), sapply(nb, length), SIMPLIFY = FALSE) ), tid_nb = unlist(nb) )
  )
}

boston_df <- 
  nb2df(boston_nb) %>% 
  group_by(tid) %>% 
    mutate(n=n()) %>% 
  ungroup() %>% 
  inner_join(boston, by=c("tid")) %>% # own features
  inner_join(                         # neighbor features
    as.data.frame(boston) %>% 
      select(-CMEDV,-lnCMEDV, -geometry) %>% 
      rename_with(~paste0(.,"_nb")), 
    by=c("tid_nb"))

head(boston_df %>% select(tid, tid_nb, NOX, CRIM, NOX_nb, CRIM_nb))

```
We can see that, for each tract `tid`, we have one row for each of its neighbors `tid_nb`. 
These rows include the covariates of `tid`, which don't change across its neighbors, and covariates of the neighbors themselves `NOX_nb, CRIM_nb, ...`.

Now we are ready to estimate the MMMM:

```{r}

# Spatial random effect:
mod.rmm1 <-
  rmm(lnCMEDV ~
        NOX + CRIM + RM + DIS + AGE +
        mm(
          id(tid, tid_nb),
          mmc(),
          mmw(w ~ 1/offset(n), constraint=1)
        ),
      iter = 1000, burnin = 100, seed=1, monitor = T,
      data = boston_df)

names(mod.rmm1)

mod.rmm1$reg.table

# Spatial fixed effects + spatial random effect:
mod.rmm2 <-
  rmm(lnCMEDV ~
        NOX + CRIM + RM + DIS + AGE +
        mm(
          id(tid, tid_nb),
          mmc(NOX_nb + CRIM_nb + RM_nb + DIS_nb + AGE_nb),
          mmw(w ~ 1/offset(n), constraint=1)
        ),
      iter = 1000, burnin = 100, seed=1, monitor = T,
      data = boston_df)
mod.rmm2$reg.table

# Calculate spatial correlation in the residual
getLambda <- function(x) {
  s.l1 <- x$reg.table["sigma.l1", "coefficients"]
  s.l2 <- x$reg.table["sigma.l2", "coefficients"]
  return(s.l1^2/(s.l1^2+s.l2^2))
}

(lambda1 <- getLambda(mod.rmm1))
(lambda2 <- getLambda(mod.rmm2))

```

Lets look at the `rmm()` function in more detail by typing `?rmm`:

![](C:/Users/benja/OneDrive - Cornell University/GitHub/rmm/docs/articles/examples_files/rmm-help.jpg){height=70%; width=70%}

**The formula object**

The most important part is the `formula` object, which is our case looks like this:
`lnCMEDV ~ NOX + CRIM + RM + DIS + AGE + mm(id(tid, tid_nb), mmc(NOX_nb + CRIM_nb + RM_nb + DIS_nb + AGE_nb), mmw(w ~ 1/offset(n), constraint=1))`

The only change compared to a `lm()` formula is the `mm()` container. Within this container, we specify 3 containers, `ids()` for the ids, `mmc()` for the considered covariates, and `mmw()` to endogenize the weight function. Here `w ~ 1/offset(n)` is specified, which implements the row-standardized weight.

The combination of spatial lag and spatial error model can also be estimated in the spatial regression framework.
Let us do that and compare the estimates:

```{r}

# Exogenous + residual spatial effect (combination of spatial lag model and spatial error model):  
mod2 <- errorsarlm(lnCMEDV ~ NOX + CRIM + RM + DIS + AGE,
              data = boston,
              listw = boston_wlist,
              Durbin =  ~ NOX + CRIM + RM + DIS + AGE) 
summary(mod2)

# Plot coefficients next to each other
coefs <- 
  mod.rmm2$reg.table[,c(1,2,4,5)] %>% 
  mutate(model="MMMM") %>% 
  add_row(
    data.frame(model="SLM", coefficients=coef(mod2)[-1], confint(mod2, level=0.95)[-1,]) %>% 
      rename(lb=X2.5.., ub=X97.5..) %>% 
      rownames_to_column("variable") %>% 
      mutate(
        variable=
          case_when(
            startsWith(variable, "lag.") ~ paste0(sub("lag.", "", variable), "_nb"), 
            variable=="(Intercept)" ~ "X0",
            TRUE ~ variable))
  ) %>% 
  filter(!variable %in% c("X0", "sigma.l1", "sigma.l2", "DIC")) 

ggplot(coefs, aes(x=variable, y=coefficients, color=model))+
  geom_point(position=position_dodge(width=0.3))+
  geom_pointrange(aes(ymin = lb, ymax = ub), position=position_dodge(width=0.3))+
  labs(title = "Coefficients", x = "Variables", y="") + coord_flip()

```

The estimates are similar but not identical. This is due to the different estimation algorithm (maximum likelihood vs Bayesian MCMC) and because the models are not 100% identical.

The advantage of `errorsarlm` is that estimation is faster and the errors are endogenous. So if this model is what you want, 
it's probably best to estimate it with the `spatialreg` package.

**Weight function regression**

The MMMM, however, allows to make the weights a function of covariates. That is, instead of assuming a specific weighting regime,
we can estimate whether the influence of a neighbor on a location depends on some covariates. 

For this example, I hypothesize that the weight of a neighbor in the total neighborhood effect depends on the similarity between neighbor and focal location,
which is an idea I borrow from social network theory.

I calculate the *similarity* as the average absolute difference between a neighbor and the focal location on the six considered covariates (i.e., the opposite of similarity). I call the variable $DIFF$: 

```{r}
# Homogeneity of focal location to its neighbors 
boston_df2 <- boston_df %>% mutate(DIFF=1/6*(abs(NOX-NOX_nb)+abs(CRIM-CRIM_nb)+abs(RM-RM_nb)+abs(DIS-DIS_nb)+abs(AGE-AGE_nb)+abs(LSTAT-LSTAT_nb)))
```

The weight function is specified in the `mmw()` container. Any function that produces weights that are bounded can be specified here. 
Functions, such as *min, max, sum* etc., can also be included in the weight function. Here I consider the following functional form:
$w=\frac{1}{n^{exp(X\beta)}}$, where $n$ is the number of neighbors of location $i$, and $X\beta$ are the covariates on which to base the weights on.
The benefit of this function is that it is bounded between 0 and 1 and if $X\beta=\boldsymbol{0}$, the weights reduce to $w=\frac{1}{n}$, which are the default row-standardized weights.

**The issue of scaling**:

The sum of all weights with row-standardized weights is $\sum_{i=1}^{506}\sum_{j}w_{ij}=\sum_{i=1}^{506}1=506$. 
We need to make sure that this overall sum does not change.
Remember that aggregation of neighbor effects is $\sum_{j \in n(i)}w_{j}(z_{j}^{\intercal}\gamma+u_{j})$.
Changes of the sum of weights will rescale the regression coefficients $\gamma$. To avoid that, two
constraints can be specified:

- `constraint=1`: constrains the weights of neighbors to sum 1 for each focal location (but allows them to differ within each location)
- `constraint=2`: constrains the all weights to sum up to the number of neighborhoods (allows them to differ within and across location)

Both constrains identify the model but have different substantive interpretations. 

Here, I allow the weights to differ within locations:

```{r}
# Spatial fixed effects + spatial random effect with "endogenized" weights:
mod.rmm3 <-
  rmm(lnCMEDV ~
        NOX + CRIM + RM + DIS + AGE +
        mm(
          id(tid, tid_nb),
          mmc(NOX_nb + CRIM_nb + RM_nb + DIS_nb + AGE_nb),
          mmw(w ~ 1/offset(n)^exp(-DIFF), constraint=1)
        ),
      iter = 10000, burnin = 2500, seed=1, monitor = T,
      data = boston_df2)

mod.rmm3$reg.table

monetPlot(mod.rmm3, "b.w[2]", "Difference")

```

We find that the dissimilarity of location and neighbors significantly impacts the weights. However,
rather than more similar neighbors having a stronger effect, we find that less similar neighbors are more impactful. 
That is, the covariates of those neighbors that are very different to the focal location in terms of air quality, level of crime, etc. 
have a stronger influence on the home values in the focal neighborhood. 

Such a insight is impossible to generate with the conventional spatial regression package!

Finally, the `rmm()` function allows us to have a look at the weight estimates:

```{r}

head(data.frame(mod.rmm3$w, sum=rowSums(mod.rmm3$w, na.rm = T)))

sum(rowSums(mod.rmm3$w, na.rm = T))

```
The weights differ across the neighbors of each location but they sum up to 1 (and thus the total of all weights equals 506).

To do:

- Introduce predict() function
- Compare predictive performance

### All your friends or just your best friend? {#e2}

A popular peer effect model in economics is the linear-in-means model, which assumes that all peers in a peer group have the same effect on an individual. In sociology, by contrast, the best-friend model is considered more often, which assumes that mainly a person's best friend will influence that individual. In this example, I use the MMMM to empirically test which of those models fits the data better.

Tbd.

### The effect of political parties' financial dependency on the survival coalition governments {#e3}

Here, I go through the example considered in the paper and analyze whether the survival of coalition governments is influences by parties’ financial dependencies, and whether this influence is homogenous across coalition parties.

Tbd.