### A simple example of using RNetica
### Russell Almond, 2015-03-04
## The goal of this exercise is to run though the basic commands to
## build, compile and draw inferences from a Netica network.
## To do this, will we will answer Problem 6.17 from Almond, Mislevy,
## Steinberg, Yan and Williamson (2015), Bayesian Models in
## Educational Assessment, Springer.
## Homework problem 6.17
###########################################################################
## Create Network and variables
### Create the network and set metadata
irt3 <- CreateNetwork("IRT3")
NetworkTitle(irt3) <- "3-item IRT network, four levels for theta"
NetworkComment(irt3) <- "Almond et al (2015) problem 6.17"
### Create nodes
## Create the latent ability node
theta <- NewDiscreteNode(irt3,"theta",c("VeryLow","Low","High","VeryHigh"))
NodeStateTitles(theta) <- paste(c(-1.5,-0.5,0.5,1.5))
NodeLevels(theta) <- c(-1.5,-0.5,0.5,1.5) # Associate numeric value with level
## This makes the node numeric, so we can later calculate an expected
## value.
## We can create multiple nodes at once, handy when they have a common
## pattern. Here create Item_1, Item_2, and Item_3
items <- NewDiscreteNode(irt3,paste("Item",1:3,sep="_"),c("Right","Wrong"))
### Add links
## There are three ways we can create links
AddLink(theta,items[[1]])
NodeParents(items[[2]]) <- list(theta) #directly set the parents:
#Note: this can also be used
#to change the order of the parents
## Note, we can often loop over item sets
for(i in 3:3) {
AddLink(theta,items[[i]])
}
### Add Conditional Probability Tables (CPTs)
## Most of the heavy lifting is done by the "[" operator, which when
## used with a NeticaNode, accesses its CPT. This is actually really
## powerful and has a lot of different options and modes. For more
## information, see:
help("Extract.NeticaNode")
## Simple node without any parents, just build a triangle distribution.
theta[] <- c(1,3,3,1)/8
## The book gives recommended values for the CPTs for the items. For
## Item_1 the book gives probabilities .378, .622, .818, .924 for Very
## Low, Low, High & Very High respectively
## Access parent state by name
items[[1]]["VeryLow"]<-.378 #Note that when number of
#elements is one less than
#number of states, the last
#value is calculated to
#complete the normalization
## Access parent state combinations by order
items[[1]][2:4]<- c(0.622,0.818,0.924)
## Secret: Values in the table are from Rasch model IRT. Item 1 has
## difficulty -1, Item 2 has difficulty 0 and Item 3 has difficulty
## 1. We can use the Rasch equations (basically inverse logit
## (theta-difficulty) to calculate the CPT.
invlogit <- function (x) exp(x)/(1+exp(x))
items[[2]][1:4] <- invlogit(NodeLevels(theta)-0)
items[[3]][1:4] <- invlogit(NodeLevels(theta)-1)
## The CPTtools package has more sophisticated tools for building CPTs
## which are set to be inserted into Netica.
## Save the network
## Probably want to save the network in its uncompiled state here.
#WriteNetworks(irt3,"T:/My Documents/BayesNets/irt3.dne")
############################################
## Compilation and inference
## Before doing inference, we need to compile network.
CompileNetwork(irt3)
## Now we can access prior beliefs
NodeBeliefs(theta)
NodeExpectedValue(theta) #Numeric nodes only
## NodeExpectedValue puts the standard deviation of expected value as
## an attribute, to access it use:
attr(NodeExpectedValue(theta),"std_dev")
### Enter Findings
## We will illustrate the inference functions by verifying that the
## net is monotonic, that is as the student gets more items right, the
## probabilities in theta (and hence the expected value) shift to
## higher values.
## Start by setting all of the observables to wrong
for( i in 1:3) {
NodeFinding(items[[i]])<-"Wrong"
}
## Now query the key node
NodeBeliefs(theta)
NodeExpectedValue(theta)
## Next, one at a time change the nodes findings to Right and recheck
## the expected value.
for (i in 1:3) {
NodeFinding(items[[i]])<-"Right"
print(NodeExpectedValue(theta))
}
### Retract Findings
## If we want to reset everything, we can retract findings
RetractNetFindings(irt3) #RetractNodeFindings affects
#just one node
NodeBeliefs(theta)
### Save the network
## Saving the network at this point, will also save all findings entered
#WriteNetworks(irt3,"T:/My Documents/BayesNets/irt3.dne")
## Don't need pathname as it is the same as the one above.