Skip to content

Instantly share code, notes, and snippets.

@SachaEpskamp
Created November 17, 2020 16:27
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save SachaEpskamp/16af6e6dd4a967361a1dec9e71873da2 to your computer and use it in GitHub Desktop.
Save SachaEpskamp/16af6e6dd4a967361a1dec9e71873da2 to your computer and use it in GitHub Desktop.
Adaptive Ising app
library("shiny")
library("psychonetrics")
library("IsingFit")
library("IsingSampler")
library("qgraph")
# Load networks:
trueNetwork <- read.csv('http://sachaepskamp.com/files/weiadj.csv')[,-1]
trueNetwork <- as.matrix(trueNetwork)
Symptoms <- rownames(trueNetwork) <- colnames(trueNetwork)
Thresholds <- read.csv('http://sachaepskamp.com/files/thr.csv')[,-1]
# Transform to -1 and 1:
trans <- LinTransform(trueNetwork, Thresholds, from = c(0,1),to=c(-1,1))
# Estimate the distribution:
startdistribution <- CURDIST <- IsingLikelihood(trans$graph, trans$thresholds, beta = 1, responses = c(-1,1))
# Compute layout:
Layout <- qgraph(trueNetwork, layout = "spring", DoNotPlot=TRUE)$layout
# Current state (-1 is off, 0 is unknown, 1 is on)
CUR_STATE <- rep(0, ncol(trueNetwork))
# Compute last means:
P <- startdistribution$Probability
nodes <- as.matrix(startdistribution[,-1])
nodes[nodes==-1] <- 0
LAST_MEANS <- as.vector(t(nodes) %*% P)
shinyServer(function(input, output, session) {
dist <- reactive({
lastX <- input$click$x
lastY <- input$click$y
if (!is.null(lastX) && !is.null(lastY)){
# Detect which node:
distances <- apply(Layout,1,function(x)
sqrt((x[1]-lastX)^2 + (x[2]-lastY)^2))
node <- which.min(distances)
# If distance is more than 0.1, do nothing
if (distances[node] < 0.1){
# What is the current probability of this node being on?
curMean <- sum((as.matrix(CURDIST[,-1])[,node]==1) * CURDIST$Probability)
# Make a copy:
CUR_STATE_COPY <- CUR_STATE
# If the node is unknown, switch on:
if (CUR_STATE[node] == 0){
curMean <- sum((as.matrix(CURDIST[,-1])[,node]==1) * CURDIST$Probability)
# Overwrite last means:
LAST_MEANS_COPY <- LAST_MEANS
LAST_MEANS_COPY[node] <- curMean
LAST_MEANS <<- LAST_MEANS_COPY
CUR_STATE_COPY[node] <- ifelse(curMean < 0.5, -1, 1)
CUR_STATE <<- CUR_STATE_COPY
} else {
# If cur state is in line with current prediction, switch to other state, else switch off:
if (CUR_STATE[node] == ifelse(LAST_MEANS[node] < 0.5, -1, 1)){
CUR_STATE_COPY[node] <- ifelse(LAST_MEANS[node] < 0.5, 1, -1)
CUR_STATE <<- CUR_STATE_COPY
} else {
CUR_STATE_COPY[node] <- 0
CUR_STATE <<- CUR_STATE_COPY
}
}
}
}
# Retrieve the start distribution:
dist <- startdistribution
# Matrix version:
matdist <- as.matrix(dist[,-1])
# Which rows to cut?
cut <- which(colMeans((CUR_STATE != 0) & t(matdist) != CUR_STATE) > 0)
# Cut rows:
if (length(cut) > 0){
dist <- dist[-cut,]
dist$Probability <- dist$Probability / sum(dist$Probability)
}
# Write out cur dist (to avoid some recursive stuff):
CURDIST <<- dist
return(dist)
})
output$plot= renderPlot({
# Color:
distTable <- dist()
P <- distTable$Probability
nodes <- as.matrix(distTable[,-1])
nodes[nodes==-1] <- 0
means <- t(nodes) %*% P
colvals <- round(2*means-1,5)
col <- ifelse(colvals > 0,
qgraph:::Fade("darkblue",abs(colvals)),
qgraph:::Fade("red",abs(colvals))
)
# Border:
known <- CUR_STATE != 0
bwidth <- ifelse(known, 4,1)
qgraph(trueNetwork, cut = 0, theme = "colorblind", layout = Layout,
labels = Symptoms, color = col, rescale = FALSE,
border.width = bwidth, vsize = 10)
})
})
library(shiny)
fluidPage(
h4("Click on nodes to make a node known (wide border). Click again to change the state (first click) or to turn the node unknown (second click)."),
# # sliderInput("mywidth", "width of the pencil", min=1, max=30, step=1, value=10),
plotOutput("plot", width = "800px", height = "500px",
click="click")
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment