SachaEpskamp/global.R

## global.R
library("bootnet")
library("qgraph")
library("dplyr")
library("ggplot2")
library("knitr")

nNodes <- 12
neighborhood <- 2

## mathNotes2.Rmd
---
output:
  word_document: default
  html_document:
    mathjax: local
    self_contained: yes
---


Suppose that the score of person $p$ on a set of variables $\pmb{y}$ at measurement $i$ is a linear combination of trait level variation stable over all measurements, $\pmb{\mathrm{t}}_{p}$, and state level variation in measurement $i$, $\pmb{\mathrm{s}}_{pi}$:
\[
\pmb{y}_{pi} = \pmb{\mathrm{t}}_{p} + \pmb{\mathrm{s}}_{pi}.
\]
For simplicity, assume that the variance-covariance matrix of state-level variation is homogeneous across people and measurements and has a mean of $\pmb{0}$ (using capitalized subscripts to indicate over which subscript the variable is random):
\[
\mathrm{Var}(\pmb{\mathrm{s}}_{pI}) = \pmb{\Theta}.
\]
Because of the assumption of homogeneity, the variance-covariance structure of state-level variation is identical across people in a fixed measurement occasion $i$ over people:
\[
\mathrm{Var}(\pmb{\mathrm{s}}_{Pi}) = \pmb{\Theta}.
\]
Assume furthermore that the traits have a different variance-covariance matrix across people:
\[
\mathrm{Var}(\pmb{\mathrm{t}}_P) = \pmb{\Omega}.
\]
Because of the assumption of homogeneity and all state-level means being $\pmb{0}$ , it follows that:
\[
\mathrm{Cov}(\pmb{\mathrm{s}}_{Pi}, \pmb{\mathrm{t}}_P) = \pmb{O},
\]
as otherwise a high trait level one some variable for person $p$ would indicate a non-zero state-level mean for that person. From this it follows that the variance-covariance structure of an \emph{cross-sectional} study at time point $i$ becomes:
\[
\mathrm{Var}\left(\pmb{y}_{Pi}\right) = \pmb{\Omega} + \pmb{\Theta},
\]
and that for any two measurements $i$ and $j$:
\[
\begin{aligned}
\mathrm{Cov}\left(\pmb{y}_{Pi}, \pmb{y}_{Pj} \right) &= \mathrm{Cov}\left(\pmb{\mathrm{t}}_{P} + \pmb{\mathrm{s}}_{Pi},\pmb{\mathrm{t}}_{P} + \pmb{\mathrm{s}}_{Pj} \right) \\\\
&= \pmb{\Omega} +
\mathrm{Cov}\left( \pmb{\mathrm{s}}_{Pi}, \pmb{\mathrm{s}}_{Pj} \right),
\end{aligned}
\]
in which $\mathrm{Cov}\left( \pmb{\mathrm{s}}_{Pi}, \pmb{\mathrm{s}}_{Pj} \right)$ denotes the within-person auto-correlations and cross-lagged effects. As null-model we will here assume that there is no within-person auto-correlation or cross-lagged effects and thus  $\mathrm{Cov}\left( \pmb{\mathrm{s}}_{Pi}, \pmb{\mathrm{s}}_{Pj} \right) = \pmb{O}$. It then follows that the covariances between two measurement occasions exactly equal the trait-level covariance structure, $\mathrm{Cov}\left(\pmb{y}_{Pi}, \pmb{y}_{Pj} \right) = \pmb{\Omega}$.

Let $\pmb{c}_p = \pmb{y}_{pi} - \pmb{y}_{pj}$ denote a vector of change scores between $\pmb{y}_{pi}$ and $\pmb{y}_{pj}$, two cross-sectional observation of the same person $p$. We then obtain:
\[
\begin{aligned}
\mathrm{Var}( \pmb{c}_P) &=\mathrm{Cov} (\pmb{c}_P, \pmb{c}_P) \\\\
&=\mathrm{Cov}(\pmb{y}_{Pi} - \pmb{y}_{Pj}, \pmb{y}_{Pi} - \pmb{y}_{Pj}) \\\\
&=\mathrm{Var}(\pmb{y}_{Pi}) + \mathrm{Var}(\pmb{y}_{Pj})  - \mathrm{Cov}(\pmb{y}_{Pi} ,  \pmb{y}_{Pj}) - \mathrm{Cov}(\pmb{y}_{Pj}, \pmb{y}_{Pi} ) \\\\
&= 2(\pmb{\Omega} + \pmb{\Theta}) - 2 \pmb{\Omega}\\\\
&= 2 \pmb{\Theta}.
\end{aligned}
\]
Now let $\pmb{w}_1$ and $\pmb{w}_2$ be a parameter vectors of the same length as $\pmb{c}_p$, such that $\pmb{w}_1^\top \pmb{c}_p$ denotes a weighted sumscore of change-scores and $\pmb{w}_2^\top \pmb{c}_p$ denotes a different weighted sumscore. We can, for example, specify:
\[
\pmb{w}_1 = \begin{bmatrix} 1 \\\\ 0 \\\\0 \\\\ \vdots \\\\ 0 \end{bmatrix},
\pmb{w}_2 = \begin{bmatrix} 0 \\\\ 1 \\\\ 1 \\\\ \vdots \\\\ 1 \end{bmatrix},
\]
such that $\pmb{w}_1^\top \pmb{c}_p$ denotes person $p$'s change in the first variable and $\pmb{w}_2^\top \pmb{c}_p$ the change in all other variables (corresponding to the main analysis of the manuscript). We can then see that:
\[
\begin{aligned}
\mathrm{Cov}(\pmb{w}_1^\top \pmb{c}_P, \pmb{w}_2^\top \pmb{c}_P) &= 2 \pmb{w}_1^\top \pmb{\Theta}\pmb{w}_2,
\end{aligned}
\]
which equals twice the sum of state-level covariances between variable 1 and all other variables.

The cross-sectional partial correlation network at measurement $i$ can be obtained from the inverse variance--covariance matrix:
\[
\pmb{K} = \left( \pmb{\Omega} + \pmb{\Theta} \right)^{-1},
\]
which we can further standardize as:
\[
\mathrm{pcor} = \frac{\kappa_{kl}}{\sqrt{\kappa_{kl}}\sqrt{\kappa_{kl}}}
\]
We would expect variables that correlate high with other variables to show high centrality in the partial correlation network. As a result, we would expect that centrality predicts $\mathrm{Cov}(\pmb{w}_1^\top \pmb{c}_P, \pmb{w}_2^\top \pmb{c}_P)$ very well.

<!-- \[ -->
<!-- \begin{aligned} -->
<!-- \mathrm{Cov}(\pmb{y}_1, \pmb{c}) &= \mathrm{Cov}(\pmb{y}_1, \pmb{y}_2 - \pmb{y}_1) \\\\ -->
<!-- &= \mathrm{Cov}(\pmb{y}_1, \pmb{y}_2)  - \mathrm{Cov}(\pmb{y}_1, \pmb{y}_1)  \\\\ -->
<!-- &= - \pmb{\Sigma}_1 -->
<!-- \end{aligned} -->
<!-- \] -->

## server.R
#
# This is the server logic of a Shiny web application. You can run the
# application by clicking 'Run App' above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#

library(shiny)

# Define server logic required to draw a histogram
shinyServer(function(input, output) {

  # Trait network:
  traitNet <- reactive({
    genGGM(nNodes,propPositive = 1,nei = neighborhood)
  })
  traitCor <- reactive({
    cov2cor(solve(diag(nNodes) - traitNet()))
  })

  # State network:
  stateNet <- reactive({
    genGGM(nNodes,input$stateRewire,propPositive = 1,nei = neighborhood)
  })
  stateCor <-  reactive({
    cov2cor(solve(diag(nNodes) - stateNet()))
  })

  # Cross sectional cor matrix:
  crossSectionCor <- reactive({
    input$traitProp*traitCor() + (1-input$traitProp)*stateCor()
  })

  # Cross sectional network:
  crossSectionPcor <- reactive({
    g <- qgraph(crossSectionCor(), graph = "pcor", DoNotPlot=TRUE)
    getWmat(g)
  })


  output$networks <- renderPlot({
    trait <- traitNet()
    state <- stateNet()
    cross <- crossSectionPcor()

    max <- max(c(trait[upper.tri(trait)], state[upper.tri(state)], cross[upper.tri(cross)]))

    layout(t(1:3))
    qgraph(trait, layout = "circle", theme = "colorblind", labels = FALSE,
           title = "Trait-level network", maximum = max, cut = 0)
    qgraph(state, layout = "circle", theme = "colorblind", labels = FALSE,
           title = "State-level network", maximum = max, cut = 0)
    qgraph(cross, layout = "circle", theme = "colorblind", labels = FALSE,
           title = "Cross-sectional network", maximum = max, cut = 0)
  })

  output$centrality <- renderPlot({
    trait <- traitNet()
    state <- stateNet()
    cross <- crossSectionPcor()

    centralityPlot(list(Trait=trait,State=state,CrossSectional=cross),include =
                     c("Strength","Closeness","Betweenness","ExpectedInfluence"))

  })


  output$changecor <- renderPlot({
    cent <- centrality(crossSectionPcor())

    centNames <- c("Strength","Closeness","Betweenness","ExpectedInfluence")

   # for (coef in seq_along(centNames)){
     # if (centNames[coef] == "Strength"){
     #   centIndex <- "OutDegree"
     # }
     #  if (centNames[coef] == "ExpectedInfluence"){
     #    centIndex <- "OutExpectedInfluence"
     #  }


    # Expected:
    expCov <- expCor <- numeric(nNodes)
    for (i in 1:nNodes){
      w1 <- w2 <- rep(0, nNodes)
      w1[i] <- 1
      w2[-i] <- 1
      expCov[i] <- 2 * t(w1) %*% stateCor() %*% w2
      expCor[i] <- expCov[i] / (sqrt(2 * t(w1) %*% stateCor() %*% w1) * sqrt(2 * t(w2) %*% stateCor() %*% w2))
    }

    # Expected change correlation:
    expChangecor <- sapply(seq_along(centNames),
                           function(coef){
                             centIndex <- centNames[coef]
                             if (centNames[coef] == "Strength"){
                               centIndex <- "OutDegree"
                             }
                             if (centNames[coef] == "ExpectedInfluence"){
                               centIndex <- "OutExpectedInfluence"
                             }

                             cor(cent[[centIndex]], expCor)
                           })
    expChangecorDF <- data.frame(index = centNames, value = expChangecor)

    # 100 replications (n = sampleSize:
    obsList <- replicate(100,{
      Trait <- ggmGenerator()(input$sampleSize,traitNet())
      State1 <- ggmGenerator()(input$sampleSize,stateNet())
      State2 <- ggmGenerator()(input$sampleSize,stateNet())

      Change <- (Trait + State1) - (Trait + State2)

      obsCor <- numeric(nNodes)
      for (i in 1:nNodes){
        obsCor[i] <- cor(Change[,i],rowSums(Change[,-i]))
      }

      obsChangecor <- sapply(seq_along(centNames),
                             function(coef){
                               centIndex <- centNames[coef]
                               if (centNames[coef] == "Strength"){
                                 centIndex <- "OutDegree"
                               }
                               if (centNames[coef] == "ExpectedInfluence"){
                                 centIndex <- "OutExpectedInfluence"
                               }

                               cor(cent[[centIndex]], obsCor)
                             })
      data.frame(index = centNames, value = obsChangecor)
    }, simplify = FALSE)
    obsDF <- bind_rows(obsList)


    # Plot:
    ggplot(obsDF, aes(x = index, y = value)) +
      geom_boxplot() +
      ylim(0,1) +
      theme_bw() +
      geom_point(data = expChangecorDF, color = "red",
                 cex = 3) +
      xlab("") + ylab("")

  })


  output$mathMD <- renderUI({
    HTML(markdown::markdownToHTML(knitr::knit('mathNotes2.Rmd', quiet = TRUE)))
  })


})

## ui.R
#
# This is the user-interface definition of a Shiny web application. You can
# run the application by clicking 'Run App' above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#
#
# # Input:
# nNodes <- 12
# stateRewire <- 0.2
# neighborhood <- 2
# # proportion of variance due to trait:
# traitProp <- 0.5

library(shiny)

# Define UI for application that draws a histogram
shinyUI(fluidPage(

  # Application title
  titlePanel("Cross-sectional Network Change Score App"),

  # Sidebar with a slider input for number of bins
  sidebarLayout(
    sidebarPanel(
      sliderInput("traitProp",
                  "Proportion trait variance",
                  min = 0,
                  max = 1,
                  value = 0.5,step = 0.01),
      sliderInput("stateRewire",
                  "State network rewire probability",
                  min = 0,
                  max = 1,
                  value = 0.5,step = 0.01),
      sliderInput("sampleSize",
                 "Sample size replications",
                 min = 100,
                 max = 10000,
                 value = 1000,step = 100)
    ),

    # Show a plot of the generated distribution
    mainPanel(
      tabsetPanel(
        tabPanel("Networks",
                 plotOutput("networks", width = "900px", height = "300px"),
                 p("Explanation: Two networks are generated using bootnet::genGGM(). The first is a chain graph for the between-subject (trait) level and the second is a rewired chain graph for the within-subject (state) level. The cross-sectional network shown is the implied cross-sectional network given both between- and within-subject network structures. As can be seen in the mathematical notes section, the higher the proportion of within-subject variance, the more we expect centrality of the cross-sectional network to predict correlation between change in a variable compared to the sum of change on all other variables, given a null model of no within-person autocorrelation (all auto-correlation due to between-subject effects). The panel 'Change-score correlation' shows the theoretically derived value as well as observed correlations between centrality (based on true network) compared to correlations between changescores.")


        ),
        tabPanel("Centrality",
                 plotOutput("centrality", width = "900px", height = "500px")
        ),

        tabPanel("Change-score correlation",
                 plotOutput("changecor", width = "900px", height = "500px"),
                 p("Correlation between (a) correlation between change score in one node and sum of all other nodes, and (b) centrality index of node. Red dot indicates mathematically expected value given null model.")

        ),
        tabPanel("Math notes",
                 uiOutput("mathMD")
        )

      )
    )
  )
)
)
	library("bootnet")
	library("qgraph")
	library("dplyr")
	library("ggplot2")
	library("knitr")

	nNodes <- 12
	neighborhood <- 2
	---
	output:
	word_document: default
	html_document:
	mathjax: local
	self_contained: yes
	---


	Suppose that the score of person $p$ on a set of variables $\pmb{y}$ at measurement $i$ is a linear combination of trait level variation stable over all measurements, $\pmb{\mathrm{t}}_{p}$, and state level variation in measurement $i$, $\pmb{\mathrm{s}}_{pi}$:
	\[
	\pmb{y}_{pi} = \pmb{\mathrm{t}}_{p} + \pmb{\mathrm{s}}_{pi}.
	\]
	For simplicity, assume that the variance-covariance matrix of state-level variation is homogeneous across people and measurements and has a mean of $\pmb{0}$ (using capitalized subscripts to indicate over which subscript the variable is random):
	\[
	\mathrm{Var}(\pmb{\mathrm{s}}_{pI}) = \pmb{\Theta}.
	\]
	Because of the assumption of homogeneity, the variance-covariance structure of state-level variation is identical across people in a fixed measurement occasion $i$ over people:
	\[
	\mathrm{Var}(\pmb{\mathrm{s}}_{Pi}) = \pmb{\Theta}.
	\]
	Assume furthermore that the traits have a different variance-covariance matrix across people:
	\[
	\mathrm{Var}(\pmb{\mathrm{t}}_P) = \pmb{\Omega}.
	\]
	Because of the assumption of homogeneity and all state-level means being $\pmb{0}$ , it follows that:
	\[
	\mathrm{Cov}(\pmb{\mathrm{s}}_{Pi}, \pmb{\mathrm{t}}_P) = \pmb{O},
	\]
	as otherwise a high trait level one some variable for person $p$ would indicate a non-zero state-level mean for that person. From this it follows that the variance-covariance structure of an \emph{cross-sectional} study at time point $i$ becomes:
	\[
	\mathrm{Var}\left(\pmb{y}_{Pi}\right) = \pmb{\Omega} + \pmb{\Theta},
	\]
	and that for any two measurements $i$ and $j$:
	\[
	\begin{aligned}
	\mathrm{Cov}\left(\pmb{y}_{Pi}, \pmb{y}_{Pj} \right) &= \mathrm{Cov}\left(\pmb{\mathrm{t}}_{P} + \pmb{\mathrm{s}}_{Pi},\pmb{\mathrm{t}}_{P} + \pmb{\mathrm{s}}_{Pj} \right) \\\\
	&= \pmb{\Omega} +
	\mathrm{Cov}\left( \pmb{\mathrm{s}}_{Pi}, \pmb{\mathrm{s}}_{Pj} \right),
	\end{aligned}
	\]
	in which $\mathrm{Cov}\left( \pmb{\mathrm{s}}_{Pi}, \pmb{\mathrm{s}}_{Pj} \right)$ denotes the within-person auto-correlations and cross-lagged effects. As null-model we will here assume that there is no within-person auto-correlation or cross-lagged effects and thus $\mathrm{Cov}\left( \pmb{\mathrm{s}}_{Pi}, \pmb{\mathrm{s}}_{Pj} \right) = \pmb{O}$. It then follows that the covariances between two measurement occasions exactly equal the trait-level covariance structure, $\mathrm{Cov}\left(\pmb{y}_{Pi}, \pmb{y}_{Pj} \right) = \pmb{\Omega}$.

	Let $\pmb{c}_p = \pmb{y}_{pi} - \pmb{y}_{pj}$ denote a vector of change scores between $\pmb{y}_{pi}$ and $\pmb{y}_{pj}$, two cross-sectional observation of the same person $p$. We then obtain:
	\[
	\begin{aligned}
	\mathrm{Var}( \pmb{c}_P) &=\mathrm{Cov} (\pmb{c}_P, \pmb{c}_P) \\\\
	&=\mathrm{Cov}(\pmb{y}_{Pi} - \pmb{y}_{Pj}, \pmb{y}_{Pi} - \pmb{y}_{Pj}) \\\\
	&=\mathrm{Var}(\pmb{y}_{Pi}) + \mathrm{Var}(\pmb{y}_{Pj}) - \mathrm{Cov}(\pmb{y}_{Pi} , \pmb{y}_{Pj}) - \mathrm{Cov}(\pmb{y}_{Pj}, \pmb{y}_{Pi} ) \\\\
	&= 2(\pmb{\Omega} + \pmb{\Theta}) - 2 \pmb{\Omega}\\\\
	&= 2 \pmb{\Theta}.
	\end{aligned}
	\]
	Now let $\pmb{w}_1$ and $\pmb{w}_2$ be a parameter vectors of the same length as $\pmb{c}_p$, such that $\pmb{w}_1^\top \pmb{c}_p$ denotes a weighted sumscore of change-scores and $\pmb{w}_2^\top \pmb{c}_p$ denotes a different weighted sumscore. We can, for example, specify:
	\[
	\pmb{w}_1 = \begin{bmatrix} 1 \\\\ 0 \\\\0 \\\\ \vdots \\\\ 0 \end{bmatrix},
	\pmb{w}_2 = \begin{bmatrix} 0 \\\\ 1 \\\\ 1 \\\\ \vdots \\\\ 1 \end{bmatrix},
	\]
	such that $\pmb{w}_1^\top \pmb{c}_p$ denotes person $p$'s change in the first variable and $\pmb{w}_2^\top \pmb{c}_p$ the change in all other variables (corresponding to the main analysis of the manuscript). We can then see that:
	\[
	\begin{aligned}
	\mathrm{Cov}(\pmb{w}_1^\top \pmb{c}_P, \pmb{w}_2^\top \pmb{c}_P) &= 2 \pmb{w}_1^\top \pmb{\Theta}\pmb{w}_2,
	\end{aligned}
	\]
	which equals twice the sum of state-level covariances between variable 1 and all other variables.

	The cross-sectional partial correlation network at measurement $i$ can be obtained from the inverse variance--covariance matrix:
	\[
	\pmb{K} = \left( \pmb{\Omega} + \pmb{\Theta} \right)^{-1},
	\]
	which we can further standardize as:
	\[
	\mathrm{pcor} = \frac{\kappa_{kl}}{\sqrt{\kappa_{kl}}\sqrt{\kappa_{kl}}}
	\]
	We would expect variables that correlate high with other variables to show high centrality in the partial correlation network. As a result, we would expect that centrality predicts $\mathrm{Cov}(\pmb{w}_1^\top \pmb{c}_P, \pmb{w}_2^\top \pmb{c}_P)$ very well.

	<!-- \[ -->
	<!-- \begin{aligned} -->
	<!-- \mathrm{Cov}(\pmb{y}_1, \pmb{c}) &= \mathrm{Cov}(\pmb{y}_1, \pmb{y}_2 - \pmb{y}_1) \\\\ -->
	<!-- &= \mathrm{Cov}(\pmb{y}_1, \pmb{y}_2) - \mathrm{Cov}(\pmb{y}_1, \pmb{y}_1) \\\\ -->
	<!-- &= - \pmb{\Sigma}_1 -->
	<!-- \end{aligned} -->
	<!-- \] -->
	#
	# This is the server logic of a Shiny web application. You can run the
	# application by clicking 'Run App' above.
	#
	# Find out more about building applications with Shiny here:
	#
	# http://shiny.rstudio.com/
	#

	library(shiny)

	# Define server logic required to draw a histogram
	shinyServer(function(input, output) {

	# Trait network:
	traitNet <- reactive({
	genGGM(nNodes,propPositive = 1,nei = neighborhood)
	})
	traitCor <- reactive({
	cov2cor(solve(diag(nNodes) - traitNet()))
	})

	# State network:
	stateNet <- reactive({
	genGGM(nNodes,input$stateRewire,propPositive = 1,nei = neighborhood)
	})
	stateCor <- reactive({
	cov2cor(solve(diag(nNodes) - stateNet()))
	})

	# Cross sectional cor matrix:
	crossSectionCor <- reactive({
	input$traitProptraitCor() + (1-input$traitProp)stateCor()
	})

	# Cross sectional network:
	crossSectionPcor <- reactive({
	g <- qgraph(crossSectionCor(), graph = "pcor", DoNotPlot=TRUE)
	getWmat(g)
	})


	output$networks <- renderPlot({
	trait <- traitNet()
	state <- stateNet()
	cross <- crossSectionPcor()

	max <- max(c(trait[upper.tri(trait)], state[upper.tri(state)], cross[upper.tri(cross)]))

	layout(t(1:3))
	qgraph(trait, layout = "circle", theme = "colorblind", labels = FALSE,
	title = "Trait-level network", maximum = max, cut = 0)
	qgraph(state, layout = "circle", theme = "colorblind", labels = FALSE,
	title = "State-level network", maximum = max, cut = 0)
	qgraph(cross, layout = "circle", theme = "colorblind", labels = FALSE,
	title = "Cross-sectional network", maximum = max, cut = 0)
	})

	output$centrality <- renderPlot({
	trait <- traitNet()
	state <- stateNet()
	cross <- crossSectionPcor()

	centralityPlot(list(Trait=trait,State=state,CrossSectional=cross),include =
	c("Strength","Closeness","Betweenness","ExpectedInfluence"))

	})


	output$changecor <- renderPlot({
	cent <- centrality(crossSectionPcor())

	centNames <- c("Strength","Closeness","Betweenness","ExpectedInfluence")

	# for (coef in seq_along(centNames)){
	# if (centNames[coef] == "Strength"){
	# centIndex <- "OutDegree"
	# }
	# if (centNames[coef] == "ExpectedInfluence"){
	# centIndex <- "OutExpectedInfluence"
	# }


	# Expected:
	expCov <- expCor <- numeric(nNodes)
	for (i in 1:nNodes){
	w1 <- w2 <- rep(0, nNodes)
	w1[i] <- 1
	w2[-i] <- 1
	expCov[i] <- 2 * t(w1) %% stateCor() %% w2
	expCor[i] <- expCov[i] / (sqrt(2 * t(w1) %% stateCor() %% w1) * sqrt(2 * t(w2) %% stateCor() %% w2))
	}

	# Expected change correlation:
	expChangecor <- sapply(seq_along(centNames),
	function(coef){
	centIndex <- centNames[coef]
	if (centNames[coef] == "Strength"){
	centIndex <- "OutDegree"
	}
	if (centNames[coef] == "ExpectedInfluence"){
	centIndex <- "OutExpectedInfluence"
	}

	cor(cent[[centIndex]], expCor)
	})
	expChangecorDF <- data.frame(index = centNames, value = expChangecor)

	# 100 replications (n = sampleSize:
	obsList <- replicate(100,{
	Trait <- ggmGenerator()(input$sampleSize,traitNet())
	State1 <- ggmGenerator()(input$sampleSize,stateNet())
	State2 <- ggmGenerator()(input$sampleSize,stateNet())

	Change <- (Trait + State1) - (Trait + State2)

	obsCor <- numeric(nNodes)
	for (i in 1:nNodes){
	obsCor[i] <- cor(Change[,i],rowSums(Change[,-i]))
	}

	obsChangecor <- sapply(seq_along(centNames),
	function(coef){
	centIndex <- centNames[coef]
	if (centNames[coef] == "Strength"){
	centIndex <- "OutDegree"
	}
	if (centNames[coef] == "ExpectedInfluence"){
	centIndex <- "OutExpectedInfluence"
	}

	cor(cent[[centIndex]], obsCor)
	})
	data.frame(index = centNames, value = obsChangecor)
	}, simplify = FALSE)
	obsDF <- bind_rows(obsList)


	# Plot:
	ggplot(obsDF, aes(x = index, y = value)) +
	geom_boxplot() +
	ylim(0,1) +
	theme_bw() +
	geom_point(data = expChangecorDF, color = "red",
	cex = 3) +
	xlab("") + ylab("")

	})


	output$mathMD <- renderUI({
	HTML(markdown::markdownToHTML(knitr::knit('mathNotes2.Rmd', quiet = TRUE)))
	})




	})
	#
	# This is the user-interface definition of a Shiny web application. You can
	# run the application by clicking 'Run App' above.
	#
	# Find out more about building applications with Shiny here:
	#
	# http://shiny.rstudio.com/
	#
	#
	# # Input:
	# nNodes <- 12
	# stateRewire <- 0.2
	# neighborhood <- 2
	# # proportion of variance due to trait:
	# traitProp <- 0.5

	library(shiny)

	# Define UI for application that draws a histogram
	shinyUI(fluidPage(

	# Application title
	titlePanel("Cross-sectional Network Change Score App"),

	# Sidebar with a slider input for number of bins
	sidebarLayout(
	sidebarPanel(
	sliderInput("traitProp",
	"Proportion trait variance",
	min = 0,
	max = 1,
	value = 0.5,step = 0.01),
	sliderInput("stateRewire",
	"State network rewire probability",
	min = 0,
	max = 1,
	value = 0.5,step = 0.01),
	sliderInput("sampleSize",
	"Sample size replications",
	min = 100,
	max = 10000,
	value = 1000,step = 100)
	),

	# Show a plot of the generated distribution
	mainPanel(
	tabsetPanel(
	tabPanel("Networks",
	plotOutput("networks", width = "900px", height = "300px"),
	p("Explanation: Two networks are generated using bootnet::genGGM(). The first is a chain graph for the between-subject (trait) level and the second is a rewired chain graph for the within-subject (state) level. The cross-sectional network shown is the implied cross-sectional network given both between- and within-subject network structures. As can be seen in the mathematical notes section, the higher the proportion of within-subject variance, the more we expect centrality of the cross-sectional network to predict correlation between change in a variable compared to the sum of change on all other variables, given a null model of no within-person autocorrelation (all auto-correlation due to between-subject effects). The panel 'Change-score correlation' shows the theoretically derived value as well as observed correlations between centrality (based on true network) compared to correlations between changescores.")


	),
	tabPanel("Centrality",
	plotOutput("centrality", width = "900px", height = "500px")
	),

	tabPanel("Change-score correlation",
	plotOutput("changecor", width = "900px", height = "500px"),
	p("Correlation between (a) correlation between change score in one node and sum of all other nodes, and (b) centrality index of node. Red dot indicates mathematically expected value given null model.")

	),
	tabPanel("Math notes",
	uiOutput("mathMD")
	)

	)
	)
	)
	)
	)