Skip to content

Commit 94e9cb2

Browse files
jjackson-ecojjackson-eco
committedOct 12, 2020
Additions to markdown 13-14
1 parent c89def6 commit 94e9cb2

File tree

1 file changed

+212
-0
lines changed

1 file changed

+212
-0
lines changed
 

‎Statistical_Rethinking_13_14.Rmd

+212
Original file line numberDiff line numberDiff line change
@@ -1012,6 +1012,218 @@ $$
10121012
$$
10131013
We have two linear models here, one for E and one for W, with which we estimate our effect of Education on Wages, and the effect of Quarter on Education. Then, in the covariance matrix **S** we account for our impact of U on both i.e. their correlation with each other.
10141014

1015+
## Social relations model
1016+
1017+
Context here is that we are interested in dyadic interactions between units, like pairs of individuals in a social network. How can we separate behaviour from specific dyadic relationships. For this we need to specify these dyadic relationships between units, but these are given by data in the model we want (not some additional covariance matrix like phylogenetic). This requires a specific covariance structure. Often called social relation models.
1018+
1019+
Example from Nicaraguan Households and the gifts that are given between households. Usually gifts of meat. 25 households, 300 dyads. Each row is a dyad and the gifts that flow between that Dyad. Gifts flow in either direction. Just raw correlation is 0.25 between A->B and B->A. But this is not the way to measure reciprocity. This includes other dyadic specific effects as well as just reciprocity such as kinship, distance etc. Number of gifts given follows a poisson distribution, we can think of the model for gifts from A to B as.
1020+
1021+
$$
1022+
\begin{aligned}
1023+
y_{A\rightarrow B} &\sim \text{Poisson}(\lambda_{AB}) \\
1024+
\log \lambda_{AB} &= \alpha + g_{A} + r_{A} + d_{AB}
1025+
1026+
\end{aligned}
1027+
$$
1028+
1029+
Where the log rate of gift giving from A to B is given by a linear model, which has an average gift giving $\alpha$, a giving offset for household A $g$, a receiving offset for household A $r$, and dyadic interactions between A and B $d$. The giving offset is individual houses givingness, some households are really giving. The $g$ and $r$ are varying effects for each household A. Then our **reciprocity** is the dyadic interactions between A and B. We have the same thing repeated for our effect the other way i.e. B to A. But these two things are related to each other. So there are **Two linear models**.
1030+
1031+
We have two covariance matrices in this model, which link together the two linear models. The first is for the giving and receiving for the household A (denoted with $i$) (those give more probably receive more generally). No means because we have the average gift giving term already ($\alpha$). The second is the dyadic covariance between household A $i$ and household B $j$. The covariance matrix here is special because it only has one variance term in it i.e. $\sigma_d$, which is because the matrix is symmetrical. The labels A and B are arbitrary, because it is just the dyadic correlation, same in both directions.
1032+
1033+
1034+
$$
1035+
\begin{aligned}
1036+
\begin{bmatrix} g_i \\ r_i \end{bmatrix} &\sim \text{MVNormal}\left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{pmatrix}
1037+
\sigma^{2}_g & \sigma_g\sigma_r\rho_{gr} \\
1038+
\sigma_g\sigma_r\rho_{gr} & \sigma^{2}_r
1039+
\end{pmatrix}\right) \\
1040+
\begin{bmatrix} d_{ij} \\ d_{ji} \end{bmatrix} &\sim \text{MVNormal}\left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{pmatrix}
1041+
\sigma^{2}_d & \sigma^{2}_d\rho_{d} \\
1042+
\sigma^{2}_d\rho_{d} & \sigma^{2}_d
1043+
\end{pmatrix}\right) \\
1044+
\end{aligned}
1045+
$$
1046+
Remember that we have 300 dyads for each direction (A-B and B-A), and because we have $d_{AB}$ and $d_{BA}$, we have estimates of all of these dyad interactions. So the model output just for the dyad terms has 600 parameters, as well as all of our other terms for giving and receiving for each A. But Hamiltonian Monte Carlo manages this no problem at all. Richards inspection of the correlation matrices shows us that for `Rho_gr` (giving and receiving), there is a negative correlation. This is because there are a few important/wealthy households that provide subsistence for others and don't take much back. For the `Rho_d` there is a really strong positive correlation i.e. in general if A gives a lot to B, B gives a lot to A. These are usually to maintain kinship bonds etc. Super interesting stuff with this model.
1047+
1048+
# Lecture 19 - Gaussian Processess
1049+
1050+
### Continuous covariance - Spatial autocorrelation (phylo gets its own markdown)
1051+
1052+
Starting with a story on American politics. Birth year with republican vote share for different elections. As you get later in time i.e. more recent elections, more younger ages (later birth years) can vote. Interesting phenomenon is that people born from similar birth years tend to vote the same way. People born in the mid 60s vote for republican more, regardless of election year. Wouldn't see this with age. There are cohort effects for political preference, and this seems to last for life.
1053+
1054+
Doesn't mean something at birth lol but something happens to cohorts at some point in their life, and seems to be happening in peoples early voting (18) that sets their voting off for life. Who is in power, and a are they popular at that age for people seems to predict voting preference. Nice paper looking at this fit a different weighting effect for how current president popularity predicts on each age. Complex non-linear pattern with age an weighting of the president popularity.
1055+
1056+
But age isn't really a category, its a continuous variable. 18 more similar to 19 and 17 than to 25. We need to deal with continuous categories like age.
1057+
1058+
Thus far we have been treating our clusters in varying effects models as categories. But can't just model them in a linear way because the relationships are complex. But not really categories. So they like ordered categories but just a lot of them. Location, phylogenetic distance are proxies for other things in nature (evolutionary distance and similar environment) that act in this way.
1059+
1060+
The set of tools to deal with these are **Gaussian Process Regression**. Pooling for infinite numbers of Clusters. Not very informative title, but normal distribution in there.
1061+
1062+
## Spatial Autocorrelation
1063+
1064+
Return to our tool use in pacific islands example. Want to predict the number of tools on each island as a function of each islands connectivity to other islands. Had a category for this before. But we can use spatial distance as a proxy for this and this is one of these continuous variable, and a proxy for for similar geology/ecology/connectivity between islands.
1065+
1066+
Lets remind ourselves of the data and plot the spatial distances, which Richard has also put in the package
1067+
1068+
```{r Kline tool plots, fig.width= 13, fig.height= 6}
1069+
data(Kline2)
1070+
data(islandsDistMatrix)
1071+
Kline <- mutate(Kline2, pop_std = as.numeric(scale(log(population))),
1072+
society = 1:10)
1073+
1074+
## Data cleaning for later
1075+
kline_list <- list(tools = Kline$total_tools,
1076+
population = Kline$population,
1077+
society = Kline$society,
1078+
Dmat = islandsDistMatrix)
1079+
1080+
a <- ggplot(Kline,aes(x = pop_std, y = total_tools)) +
1081+
geom_point(size = 4, colour = "green" ) +
1082+
geom_text(aes(label = culture, colour = NULL),
1083+
nudge_y = 2.3, nudge_x = 0.08, show.legend = F) +
1084+
labs(x = "Standardised population size", y = "Total number of tools") +
1085+
theme_bw(base_size = 14) + theme(panel.grid = element_blank())
1086+
1087+
b <- tibble(island1 = rownames(islandsDistMatrix), as_tibble(islandsDistMatrix)) %>%
1088+
pivot_longer(-island1) %>%
1089+
ggplot(aes(x = island1, y = name)) +
1090+
geom_tile(aes(fill = value)) +
1091+
geom_text(aes(label = round(value,2))) +
1092+
scale_fill_viridis_c(option = "D",begin = 0.1,
1093+
end = 0.8, name = "Distance\n(1000 km)") +
1094+
scale_x_discrete(expand = c(0,0)) +
1095+
scale_y_discrete(expand = c(0,0)) +
1096+
labs(x = NULL, y = NULL) +
1097+
theme_bw(base_size = 14)
1098+
1099+
grid.arrange(a,b, ncol = 2, widths = c(4,5))
1100+
1101+
```
1102+
1103+
We can think of this like a map of confound threat - We think that tool use is going to be more similar in islands that are close together.
1104+
1105+
We're going to use the scientific tools model, where we are modelling our expected number of tools $\lambda$ as a Poisson, and with our innovation $\alpha$, our diminishing returns $\beta$ and our loss $\gamma$
1106+
1107+
$$
1108+
\begin{aligned}
1109+
\text{TOOL}_i &\sim \text{Poisson}(\lambda_i) \\
1110+
\lambda_i &= \alpha P^{\beta}_{i}/\gamma
1111+
\end{aligned}
1112+
$$
1113+
We get the Gaussian process of spatial autocorrelation in here by adding in a factor for our society. We add in a term $k$ which is a varying effect of the society at each observation $i$. We exponentiate it because it has to be positive, following $\lambda$ which is the rate. It as a factor, which means that it adjusts our expectation for the whole equation. So, when or value of $k = 0$, indicating no effect of society, then $\exp0 = 1$ so the equation doesn't change. If $k$ is negative, then $\lamda$ is worse proportionally than expected, and if greater than 0 then it is better than expected.
1114+
1115+
$$
1116+
\begin{aligned}
1117+
\text{TOOL}_i &\sim \text{Poisson}(\lambda_i) \\
1118+
\lambda_i &= \exp({k_{\text{SOCIETY}[i]}})\alpha P^{\beta}_{i}/\gamma
1119+
\end{aligned}
1120+
$$
1121+
1122+
We are going to estimate $k$ using the matrix of spatial distances for all the islands. Not distinct category memberships like before.
1123+
1124+
## The Gaussian process
1125+
1126+
At the heart of the Gaussian process regression is the multivariate gaussian prior for the intercepts with a covariance matrix **K**. We also define our covariance matrix, which stores values between society $i$ and society $j$
1127+
1128+
$$
1129+
\begin{aligned}
1130+
\begin{bmatrix} k_1 \\ k_2 \\ k_3 \\ \cdots \\ k_{10} \end{bmatrix} &\sim \text{MVNormal}\left(\begin{bmatrix} 0 \\ 0 \\ 0 \\ \cdots \\ 0 \end{bmatrix}, \textbf{K}\right) \\ \\
1131+
K_{ij} &= \eta^2 + \exp(-\rho^2 D^2_{ij}) + \delta_{ij}\sigma^2
1132+
\end{aligned}
1133+
$$
1134+
So this time we're building the covariance matrix not just with info on the variance and covariance of our two clusters. We parameterise it from the distance matrix. This is sometimes called the *L2norm* formulation, a common way of incorporating spatial autocorrelation. It has just three parameters, $\eta$, $\rho$ and $\sigma$.
1135+
1136+
So **K** is our covariance between islands, $\eta^2$ is the maximum covariance, $\rho^2$ is a rate of decline with distance, $D^2_{ij}$ is our distance squared. *This bit with $e$ to the power of minus something squared is out explicit Gaussian Part of the Gaussian process.* It gives a bell curve in shape. Then our final part $\delta_{ij}\sigma^2$ is sometimes called the jitter, where the $\delta$ is just an identity function telling you in binary whether i and j are equal ($\delta_{ij} = 1$) or different ($\delta_{ij} = 0$). This turns the variance, $\sigma^2$ of an island with itself on and off. So gives the variance of an island with itself when you have multiple observations in your cluster (important for you). This isn't needed here really, so we just set our sigma value at something arbritrary above 0. The other two parameters need priors though.
1137+
1138+
The part in the middle is the key part, creating this Gaussian decline in our covariance with distance. If it were linear, it would be much more exponential in decline, but because we square it it creates a nice smooth decline with distance.
1139+
1140+
Our full joint model is
1141+
1142+
$$
1143+
\begin{aligned}
1144+
\text{TOOL}_i &\sim \text{Poisson}(\lambda_i) \\
1145+
\lambda_i &= \exp({k_{\text{SOCIETY}[i]}})\alpha P^{\beta}_{i}/\gamma \\
1146+
\textbf{k} &\sim \text{MVNormal}((0,0,\ldots,0), \textbf{K}) \\
1147+
\textbf{K}_{ij} &= \eta^2 + \exp(-\rho^2 D^2_{ij}) + \delta_{ij}\sigma^2 \\
1148+
\alpha &\sim \text{Exponential}(1) \\
1149+
\beta &\sim \text{Exponential}(1) \\
1150+
\eta^2 &\sim \text{Exponential}(2) \\
1151+
\rho^2 &\sim \text{Exponential}(0.5)
1152+
\end{aligned}
1153+
$$
1154+
And that is it. This time our **K** matrix is specified with a special prior which does the calculations for us called `cov_GPL2` in `rethinking``. Richard did some prior predictive simulation for the rates of rho and eta in the covariance matrix. These both imply a pretty large drop-off in spatial covariance as you move away. This will get shifted a bit with our posterior.
1155+
1156+
## Fitting a Gaussian Process L2 distance model
1157+
1158+
Now lets fit it.
1159+
1160+
```{r tools spatial ulam, error = F, results = F, mesage = F}
1161+
m14.8 <- ulam(alist(
1162+
# model
1163+
tools ~ dpois(lambda),
1164+
lambda <- (a*population^b/g)*exp(k[society]),
1165+
1166+
# Gaussian process - our multivariate varying effect and distance cov matrix
1167+
vector[10]:k ~ multi_normal(0, SIGMA),
1168+
matrix[10,10]:SIGMA <- cov_GPL2(Dmat, etasq, rhosq, 0.01),
1169+
1170+
# fixed priors
1171+
c(a,b,g) ~ dexp(1),
1172+
etasq ~ dexp(2),
1173+
rhosq ~ dexp(0.5)),
1174+
data = kline_list, chains = 4, cores = 4, iter = 2000)
1175+
```
1176+
1177+
And now we want to take a look at this model. We're going to pull out some posterior predictions for the median relationship between population size and tools, excluding our distance measure. Then, we're going to compute the spatial covariance at the median from our matrix **K** by pushing our rho and eta through our equation again, and see how the spatial covariance actually looks between populations. Richard also shows the posterior covariance with distance and the matrix.
1178+
1179+
```{r kline spatial post, fig.width = 7, fig.height= 7}
1180+
precis(m14.8, depth = 2)
1181+
1182+
kpost <- extract.samples(m14.8)
1183+
1184+
## Compute the covariance matrix from the posterior
1185+
K <- matrix(0,nrow = 10, ncol = 10)
1186+
for(i in 1:10){
1187+
for(j in 1:10){
1188+
K[i,j] <- median(kpost$etasq)*
1189+
exp(-median(kpost$rhosq)*islandsDistMatrix[i,j]^2)}}
1190+
1191+
diag(K) <- median(kpost$etasq) + 0.01
1192+
1193+
# convert to correlation matrix
1194+
Rho <- round(cov2cor(K), 2)
1195+
1196+
# convert to data
1197+
Rho_dat <- tibble(as_tibble(Rho), society1 = 1:10) %>%
1198+
pivot_longer(-society1, names_to = "society2", values_to = "Corr",
1199+
names_prefix = "V") %>%
1200+
mutate(society2 = as.integer(society2)) %>%
1201+
left_join(x = ., y = Kline[,c(11,4,9)], by = c("society1" = "society")) %>%
1202+
left_join(x = ., y = Kline[,c(11,4,9)], by = c("society2" = "society"),
1203+
suffix = c("1", "2"))
1204+
1205+
# Posterior predictions excluding d
1206+
logpop.seq <- seq(6,14,length.out = 30)
1207+
lambda <- sapply(logpop.seq, function(lp) exp(kpost$a + kpost$b*lp - kpost$g))
1208+
preddat <- tibble(logpop = logpop.seq,
1209+
md = apply(lambda, 2, median),
1210+
lwr = apply(lambda, 2, PI, prob = 0.8)[1,],
1211+
upr = apply(lambda, 2, PI, prob = 0.8)[2,])
1212+
1213+
# The plot
1214+
ggplot(Kline,aes(x = logpop, y = total_tools)) +
1215+
geom_segment(data = Rho_dat,
1216+
aes(x = logpop1, xend = logpop2,
1217+
y = total_tools1, yend = total_tools2, alpha = Corr)) +
1218+
geom_point(aes(size = logpop), show.legend = F) +
1219+
geom_line(data = preddat, aes(x = logpop, y = md),
1220+
colour = "cornflowerblue") +
1221+
geom_text(aes(label = culture, colour = NULL),
1222+
nudge_y = 3, nudge_x = 0.08) +
1223+
labs(x = "log Population size", y = "Number of tools") +
1224+
theme_bw(base_size = 14) + theme(panel.grid = element_blank())
1225+
```
1226+
10151227

10161228

10171229

0 commit comments

Comments
 (0)
Please sign in to comment.