Diagnostics code for review of the NetworkChange package

diagnostics.Rmd

---
title: "Diagnostics"
output: html_notebook
---


```{r}
library(NetworkChange)
require(sna)
```

```{r}
set.seed(11173)
n <- 10 ## number of nodes in each cluster
Y <- MakeBlockNetworkChange(n=n, break.point = .5,
                            base.prob=.05, block.prob=.7,
                            T=20, type ="split")
```


```{r}
set.seed(52253)
G <- 100
Yout <- list()
for(i in 1:8) {
  Yout[[i]] <-  NetworkChange(Y, R=2, m=1, mcmc=G, burnin=G, verbose=0)
}

```

```{r}
get_U_samples <- function(list_of_fits, period, index) {
  n_iterations <- dim(attr(list_of_fits[[1]], "Umat")[[period]])[1]
  s <- array(NA_real_, c(n_iterations, length(list_of_fits)))
  for(i in 1:length(list_of_fits)) {
    s[,i] <- attr(list_of_fits[[i]], "Umat")[[period]][, index]
  }
  s
}

get_V_samples <- function(list_of_fits, index) {
  n_iterations <- dim(attr(list_of_fits[[1]], "Vmat"))[1]
  s <- array(NA_real_, c(n_iterations, length(list_of_fits)))
  for(i in 1:length(list_of_fits)) {
    s[,i] <- attr(list_of_fits[[i]], "Vmat")[, index]
  }
  s
}

get_s2_samples <- function(list_of_fits, period, index = 1) {
  n_iterations <- dim(attr(list_of_fits[[1]], "s2mat")[[period]])[1]
  s <- array(NA_real_, c(n_iterations, length(list_of_fits)))
  for(i in 1:length(list_of_fits)) {
    s[,i] <- attr(list_of_fits[[i]], "s2mat")[[period]][, index]
  }
  s
}


to_diagnose <- list(
  U1_1 = get_U_samples(Yout, 1, 1),
  U1_15 = get_U_samples(Yout, 1, 15),
  U2_40 = get_U_samples(Yout, 2, 40),
  U2_60 = get_U_samples(Yout, 2, 40),
  V_1 = get_V_samples(Yout, 1), 
  V_15 = get_V_samples(Yout, 15), 
  V_33 = get_V_samples(Yout, 33), 
  s2_1 = get_s2_samples(Yout, 1),
  s2_2 = get_s2_samples(Yout, 2)
)

purrr::imap_dfr(to_diagnose, function(s, name) { data.frame(par = name, Rhat = posterior::rhat(s), ess_bulk = posterior::ess_bulk(s)) })
```
```{r}
matplot(to_diagnose$U1_15, type = "l")
matplot(to_diagnose$U2_40, type = "l")
matplot(to_diagnose$s2_1, type = "l")
matplot(to_diagnose$s2_2, type = "l")
```

```{r}
for(i in 1:8) {
  plotV(Yout[[i]], cex=2)
}

```

```{r}
for(i in 1:8) {
  Ydraw <- drawPostAnalysis(Yout[[i]], Y, n.cluster=c(2,3))
  multiplot(plotlist=Ydraw, cols=2)
}
```


```{r}
data(MajorAlly)
Yally <- MajorAlly
time <- dim(Y)[3]
drop.state <- c(which(colnames(Yally) == "USA"), which(colnames(Yally) == "CHN"))
newY <- Yally[-drop.state, -drop.state, 1:62]


G <- 100
set.seed(1990)
test.run <- NetworkStatic(newY, R=2, mcmc=G, burnin=G, verbose=0,
                          v0=10, v1=time*2)
V <- attr(test.run, "V")
sigma.mu = abs(mean(apply(V, 2, mean)))
sigma.var = 10*mean(apply(V, 2, var))
v0 <- 4 + 2 * (sigma.mu^2/sigma.var)
v1 <- 2 * sigma.mu * (v0/2 - 1)

G <- 100
K <- dim(newY)
m <- 2
initial.s <- sort(rep(1:(m+1), length=K[[3]]))
set.seed(11223);
resAlly <- list()
for(i in 1:8) {
  resAlly[[i]] <- NetworkChange(newY, R=2, m=m, mcmc=G, initial.s = initial.s,
                     burnin=G, verbose=0, v0=v0, v1=v1)

}
```

```{r}
to_diagnose_ally <- list(
  U1_1 = get_U_samples(resAlly, 1, 1),
  U1_5 = get_U_samples(resAlly, 1, 5),
  U2_3 = get_U_samples(resAlly, 2, 3),
  U2_10 = get_U_samples(resAlly, 2, 10),
  U3_7 = get_U_samples(resAlly, 3, 7),
  U3_11 = get_U_samples(resAlly, 3, 11),
  V_1 = get_V_samples(resAlly, 1), 
  V_33 = get_V_samples(resAlly, 33), 
  s2_1 = get_s2_samples(resAlly, 1),
  s2_2 = get_s2_samples(resAlly, 2),
  s2_3 = get_s2_samples(resAlly, 3)
)

purrr::imap_dfr(to_diagnose_ally, function(s, name) { data.frame(par = name, Rhat = posterior::rhat(s), ess_bulk = posterior::ess_bulk(s)) })
```

```{r}
matplot(to_diagnose_ally$s2_1, type = "l")
matplot(to_diagnose_ally$s2_2, type = "l")
matplot(to_diagnose_ally$s2_3, type = "l")

```

```{r}
for(i in 1:8) {
  fit <- resAlly[[i]]
  attr(fit, "y") <- 1:K[[3]]
  plotState(fit, start=1)
}
```

```{r}
for(i in 1:8) {
  fit <- resAlly[[i]]
  p.list <- drawPostAnalysis(fit, newY, n.cluster=c(4, 4, 3))
  multiplot(plotlist = p.list, cols=3)
}
```