---
title: "DIY modeling of a binary network outcome"
author: "Peter Hoff"
date: "`r substring(Sys.time(),1,10)`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{DIY modeling of a binary network outcome}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

The `amen` package provides wrapper functions `ame` and `ame_rep`
for fitting a wide variety of AME models. However, 
wrapper functions can be opaque, and even wrapper functions 
with many options are limited by the energy and imagination of their 
authors. Therefore, some users may prefer to use some of the 
low-level `amen` functions to build their own custom 
Gibbs samplers for AME model fitting. 

In this vignette, 
I show how to build such a sampler for a binary AME model. 
The example data we will use is the international relations 
data `IR90s` that comes with the `amen` package. We will 
model the binary indicator of conflict initiation between 
pairs of countries. Note that this dyadic variable is 
a count of the number of conflicts. For the purpose of this example, 
we dichotomize it:
```{r,message=FALSE}
library(amen)

Y<-1*(IR90s$dyadvars[,,"conflicts"] >0 )

netplot(Y,plot.iso=FALSE,plotnames=TRUE)
```

Which pairs of countries have conflicts? 
Well, poor countries that are far away from 
each other don't have conflicts. With this in mind, 
let's fit a model of international conflict 
that has as predictors 
the GDPs of the initiator and target countries, 
as well as the 
geographic distance between the pair (both on a log scale):


```{r}
Xn<-log(IR90s$nodevars[,2])
Xd<-log(IR90s$dyadvars[,,3]+1) 
```

The binary  probit AME model we will fit to `Y`, `Xn` and `Xd` is 
\[
\begin{align*}
z_{i,j} & = \beta_0 + \beta_r \text{gdp}_i + \beta_c \text{gdp}_j + \beta_d \text{dist}_{i,j}  + u_i^\top v_j + a_i +b_j + \epsilon_{i,j}  \\
y_{i,j} & = 1( z_{i,j}> 0 ) 
\end{align*} 
\]
where $\text{Var}{ (\begin{smallmatrix} a_i \\ b_i \end{smallmatrix} ) } = \Sigma_{ab}$,
$\text{Var}{  (\begin{smallmatrix} u_i \\ v_i \end{smallmatrix} ) } = \Sigma_{uv}$, and
\[ \text{Var}{ ( \begin{smallmatrix} \epsilon_{i,j} \\ \epsilon_{j,i}\end{smallmatrix}) } = \Sigma_\epsilon =   \sigma^2  \begin{pmatrix}  1 & \rho \\ \rho & 1 \end{pmatrix}.
\]


We will construct a Gibbs sampler to approximate 
the joint posterior distribution of all unknown 
quantities, which includes 

* the global parameters $\beta, \Sigma_{ab}, \Sigma_{uv}, \rho$; 

* the latent nodal effects $\{ (a_i, b_i, u_i, v_i): i=1,\ldots, n\}$; 

* the unobserved latent dyadic variables $z_{i,j}$. 


Here we go. 
First note that the "linear regression" part of the model 
may be written as $\beta^\top x_{i,j}$ where 
$x_{i,j} = ( 1, \text{gdp}_i ,  \text{gdp}_j ,  \text{dist}_{i,j} )$. 
and $\beta = ( \beta_0 , \beta_r, \beta_c \beta_d)$. 
We construct the  "design array", the $n\times n \times p = 130 \times 130 \times 4$:

```{r}
X<-design_array(Xrow=Xn,Xcol=Xn,Xdyad=Xd,n=nrow(Y)) 
```

If you are not sure what this is doing, examine `X` yourself or 
study the following output:
```{r}
dim(X) 

dimnames(X)[[3]] 

X[1:3,1:3,1]

X[1:3,1:3,2]

X[1:3,1:3,3]

X[1:3,1:3,4]
```


Next we need starting values for all (but one) of our 
unknown quantities. 
Because this is just a tutorial  I will set these lazily:

```{r}
Sab<-diag(2) 

rho<-0 

R<-2 

Suv<-diag(2*R)

U<-V<-matrix(0,nrow(Y),R) 

# just need Z[i,j]<0 if Y[i,j]=0, Z[i,j]>0 if Y[i,j]=1  
Z<-matrix(zscores(Y,ties.method="random"),nrow(Y),nrow(Y))
Z<-Z-max(Z[!is.na(Y) & Y==0 ] ) 
diag(Z)<-0
```

Finally, we specify some aspects of the prior distribution. 
The functions in the `amen` 
package use 
a $g$-prior for $\beta$ by default (but this can be changed)
and an arc-sine prior or uniform prior for $\rho$. The priors for 
$\Sigma_{ab}$ and $\Sigma_{uv}$ are

* $\Sigma_{ab}^{-1} \sim  W( \Sigma_{ab0}^{-1}, \eta_0 )$; 
* $\Sigma_{uv}^{-1} \sim W( \Sigma_{uv0}^{-1} , \kappa_0)$. 

Sparse network datasets may not provide much information 
with which to estimate $\Sigma_{ab}$ and $\Sigma_{uv}$. 
For such data, estimation stability can be improved by 
concentrating the scale of these covariance matrices 
near some empirical Bayes estimate, where 
the concentration is decreasing with the amount of information 
in the data. 
The following calculations provide one way to adaptively 
choose the priors for $\Sigma_{ab}$ and $\Sigma_{uv}$ 
based on the information in the data:
```{r}
ybar<-mean(Y,na.rm=TRUE) ; mu<-qnorm(ybar)
E<- (Y - ybar)/dnorm(qnorm(ybar)) ; diag(E)<-0
a<-rowMeans(E,na.rm=TRUE)  ; a[is.na(a)]<-0
b<-colMeans(E,na.rm=TRUE)  ; b[is.na(b)]<-0

vscale<-mean(diag(cov(cbind(a,b))))
PHAT<-pnorm(mu+outer(a,b,"+"))
vdfmlt<-.25/mean(PHAT*(1-PHAT))

Sab0<-diag(2)*vscale 
eta0<-round(4*vdfmlt)       

Suv0<-diag(2*R)*vscale 
kappa0<-round((2*R+2)*vdfmlt)  
``` 


Now all we need is a for-loop that iteratively updates 
the unknown quantities from their full conditional distributions:


```{r,results="hide",cache=TRUE,fig.keep='last'}  
BETA<-NULL ; RHO<-NULL ; UV<-matrix(0,nrow(Y),nrow(Y))

for(s in 1:10000)
{ 
  # update beta, a and b  
  tmp<-rbeta_ab_fc(Z,Sab,rho,X,offset=U%*%t(V))
  beta<-tmp$beta ; a<-tmp$a ; b<-tmp$b

  # update UV 
  tmp<-rUV_fc(Z,U,V,Suv,rho,offset=Xbeta(X,beta) + outer(a,b,"+")) 
  U<-tmp$U ; V<-tmp$V

  # update Suv
  Suv<-rSuv_fc(U,V,Suv0=Suv0,kappa0=kappa0) 

  # update Sab
  Sab<-rSab_fc(a,b,Sab0=Sab0,eta0=eta0)

  # update rho 
  rho<-rrho_mh(Z,rho,offset=Xbeta(X,beta)+outer(a,b,"+")+U%*%t(V)) 

  # update Z 
  Z<-rZ_bin_fc(Z, Xbeta(X,beta)+outer(a,b,"+")+U%*%t(V), rho, Y)  

  # periodically save some output   
  if(s%%25==0 & s>500) 
  { 
    cat(round(100*s/10000),"% complete\n") 
    BETA<-rbind(BETA,beta) 
    RHO<-c(RHO,rho) 
    UV<-UV+U%*%t(V)  
    matplot(BETA,type="l") 
  }
}
```	

The mixing is a bit slow, so you would probably want to 
run the Markov chain longer to obtain some final results. 


Here is the "wrapper" version, fitting the same model:
```{r,results='hide',fig.keep='last',cache=TRUE}
fit<-ame(Y,Xrow=Xn,Xcol=Xn,Xdyad=Xd,R=2,family="bin")
```

Now compare posterior density approximations:
```{r,echo=c(-1)}
par(mfrow=c(2,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
for(j in 1:4)
{
  plot(density(BETA[,j] ),main="",xlab="")
  lines(density(fit$BETA[,j] ),col="green") 
}
```

The longer you run the Markov chains, the closer the 
two approximations should get to each other (and to the
actual posterior distribution). 





