--- title: "Regression Distributions" author: "Russell Almond" date: "2/14/2022" output: html_document runtime: shiny --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) library(shiny) library(shinyjs) ``` This demonstration shows how to build and set the parameters for a regression distribution using the DiBello framework. The DiBello method has three steps: 1. Mapping the input variables to continuous values (effective thetas), using the `PnodeStateValues` property of the `PNode` object. 2. Combining the parent thetas into a single \dQuote{effective theta} for the child variable using the rule specified in the `PnodeRule` property of the child node. 3. Mapping the effective theta back into a conditional probability table using the `PnodeLink` function. ## Preliminaries These gadgets use code from the [Peanut](https://pluto.coe.fsu.edu/RNetica/Peanut.hmtl) package, but that requires a Bayes net implementation to support the code. The implementation we will use is `PNetica`. Loading `PNetica` will load `Peanut`, `RNetica` and `CPTtools` (which does most of the computations). If you want to try tweaking this, the Rmd source can be found at https://pluto.coe.fsu.edu/svn/common/rgroup-shiny/Peanut/CompensatoryGadget.Rmd ```{r library} library(PNetica) ``` Next we need to start a Netica session to store our networks. ```{r session} sess <- RNetica::NeticaSession() RNetica::startSession(sess) ``` We are only building a really small net, so the demo version is fine. ## Build a small network. Create a network with two parent variables (`theta1`, `theta2`) and one child (`partialCredit`). ```{r structure} tNet <- RNetica::CreateNetwork("TestNet",sess) theta1 <- RNetica::NewDiscreteNode(tNet,"theta1", c("High","Mid","Low")) theta2 <- RNetica::NewDiscreteNode(tNet,"theta2", c("High","Mid","Low")) PnodeParents(theta2) <- list(theta1) partial3 <- RNetica::NewDiscreteNode(tNet,"partial3", c("FullCredit","PartialCredit","NoCredit")) PnodeParents(partial3) <- list(theta1,theta2) ``` ## Effective Thetas The key idea of the DiBello model is that associated with each observable outcome variable (in other words, a child of several skill variables) is an effective dimension, $\theta$, in which students who are higher in this dimension will perform better than students who are lower. DiBello suggested mapping each configuration of the parent variables (corresponding to the rows in the conditional probability table) into an \dQuote{effective theta}. This is done in two steps: first, each parent state is mapped to a real number; and second, a combination rule is used to combine the parent values. In Psychometrics, it is customary to think of $\theta$ and a unit normal variable, corresponding to the population distribution of interest. As the scale and location are arbitrary, the mean and standard deviation is set to one. Therefore 0, corresponds to the median of the population, 1 to the `r round(100*qnorm(1),1)` percentile and -1 to the `r round(100*qnorm(-1),1)` percentile. For a number of reasons (see Almond, et al., 2015), equally spaced quantiles of a normal distribution are a good starting point for setting these values. The function `CPTtools::effectiveThetas(n)` calculates these quantiles for a variable with $n$ categories. Thus, the code sample below is the recommended practice. ```{r effectiveTheta} effectiveThetas(3) PnodeStateValues(theta1) <- effectiveThetas(PnodeNumStates(theta1)) PnodeStateValues(theta2) <- effectiveThetas(PnodeNumStates(theta2)) ``` The setter form of `PnodeStateValues(node)` should take a numeric vector of the same length as the number of states of `node`. Analysts are free to experiment with other mappings, but not that the mapping is defined at the level of the parent node, so that any child that has `theta1` as a parent will use the same state values. ## Compensatory Combitnation Rule The compensatory combination rule is simply a weighted average of the effective thetas. Let $r$ be an index for the row of the conditional probability table, and let $\tilde {\theta_k(r)}$ be the effective theta value for the $k$th parent in that row. Then the effective theta values for that row is: $$ \tilde{\theta(r)} = \frac{1}{\sqrt{K}} \sum_k \alpha_{k},\tilde{\theta_{k}(r) - \beta_v$$ Here $k$ is an index over the parent variables, and $v$ is an index over the states of the child variable. By convention, the lowest value child variable is given an effective theta of $-\infty$, so it does not need to be specified. The factor $1/\sqrt{K}$ is a variance stabilization term. (Without it the varaiance of $\tilde{\theta(r)} would increase as more and more parents are added). We have two types of parameters: * _discriminations_ (also called slopes), the ${\alpha_k}$ values. These are inspected/set with the `PnodeAlphas()` function. The value should be a vector corresponding to the parent variables (so of length $K$). In most psychometric applications, the alphas are constrained to be strictly positive, so for the purposes of model fitting, $\log({\alpha_k})$ is taken as the parameter. The function `PnodeLnAlphas()` works with this transformed parameter space. * _demands_ (also called difficulties or negative intercepts), the $\beta_v$ values are related to how prevelent each state is for a person of median ability. The exact meaning depends on the link function. This value should be a __list__ with one fewer elements than the number of states (by convention, the lowest state is omitted). The convention used in both `CPTtools` and `Peanut` is that if the parameter is indexed by parent variables, it is a `vector`, and if it is indexed by states of the child variable, it is a `list`. This allows some more complex expressions where parameters vary by both parents and child states, see [DPCGadget](DPCGadget.Rmd). The function `CPTtools::eThetaFrame()` is used to show the effective theta calculations. It takes four arguments: `skillLevels` which is a list of the `ParentStateValues()` for each parent node, `lnAlphas` which is the vector of log alphas, `beta` which is the demand, should be a scalar, and the `rule` which is a special function. In this case, `CPTtools::Compensatory`. ```{r compensatory eTheta, echo=TRUE} eThetaFrame(ParentStates(partial3),log(c(Skill1=1.2,Skill2=.8)),0, Compensatory) ``` The first two numeric columns in this data frame show the effective thetas for the parent variables. The final column, shows the result of applying the `Compensatory()` function with the specified parameters. Note that the combination rule is itself a parameter (set with the `PnodeRule()` function). Any function with the same arguments would work. The `CPTtools` package supplies two, `Conjunctive()` (based on minimum) and `Disjunctive()` (based on max). However, two variants of these `OffsetConjuctive()` and `OffsetDisjunctive()`, which use a slightly different parameterization are more useful. See [OffsetGadget](OffsetGadget.Rmd) to look at those. ## Link Function Note that the first step goes from the discrete scale to the continuous scale by using quantiles of the normal distribution. Its inverse, should also use the normal distribution. Use the same mapping as before, to divide the normal distribution into $K$ groups, where $K$ is the number of states of child variable. Next, create a new normal curve to represent the information from the parent variables. Applying the combination function to the parent effective thetas produces the mean of that normal distribution. The standard deviation is an additional parameter, called the _link scale_. (It is set with the `PnodeLinkScale()` function.) The conditional probabilities are now the area under curve between cut points. ```{r normal offset, echo=FALSE} states <- c("Low","Medium","High") M <- length(states) c <- qnorm((1L:(M-1L))/M,0,1) thetas <- qnorm((1L:M -.5)/M,0,1) curve(dnorm(x),xlim=c(-3,3),xaxt="n",yaxt="n",ylim=c(0,.5),col="red", ylab="",xlab="Effective Theta") segments(c,-.1,c,.6) text(thetas,.45,states) axis(1,c, do.call(expression, lapply(1:(M-1), function (m) substitute(c[m],list(m=m))))) theta <- .5 sig <- .8 curve(dnorm(x,theta,sig),col="sienna",add=TRUE) cc <- seq(c[1],c[2],.025) polygon(c(c[1],cc,c[2]),c(0,dnorm(cc,theta,sig),0),angle=45,col="sienna2") cplus <- c(-Inf,c,Inf) pvals <- diff(pnorm(cplus,theta,sig)) names(pvals) <- states pvals ``` ## Nodes and Pnodes A `Pnode` is an ordinary node wrapped in some extra functions so it can store information about the parameterization and parameters. The package `PNetica` (the first and currently only implementation of the `Peanut` protocol) adds implementation for the function in the `Peanut` protocol. The functions are: * `is.Pnode`, `as.Pnode`, `Pnode` --- class definition objects. * `PnodeNet` --- fetches the network. * `PnodeName`, `PnodeTitle`, `PnodeDescription`, `PnodeLabels`, `PnodeStates`, `PnodeStateTitles`, `PnodeStateDescriptions`, `PnodeStateBounds`, `isPnodeContinuous`, `PnodeProbs` --- These are just generic names for common opperations supported by most Bayesian network packages. * `PnodeParents`, `PnodeParentNames`, `PnodeNumParents` --- These functions are related to the structure of the network. The functions listed above are all generic. Any implementation of the `Peanut` protocol will supply methods for these. Therefore, if a programmer uses these functions, a different Bayesian network engine can be easily substituted for `PNetica`. These next functions are important to the `Peanut` protocol. * `PnodeStateValues` --- This associates each state with a real value (see effective thetas above). `PnodeParentTvals` fetches the state values for all of the parent nodes. * `PnodeRules` --- The combination rule, e.g., `CPTtools::Compensatory` * `PnodeLink` --- The link function, e.g., `CPTtools::partialCredit` or `CPTtools::gradedResponse`. The link function `CPTtools::normalLink` needs a scale parameter, this is accessed with `PnodeLinkScale`. * `PnodeAlphas`, `PnodeLnAlphas` --- These access the the slope or discrimination parameters. The two functions access the same set of parameters, it is just that one has the parameters on the log scale and one on the natural scale. For the `Compensatory` rule, the length of the alpha vector should match the number of parents. The function `defaultAlphas` will return a vector of log alphas the appropriate length. * `PnodeBetas` -- These access the demand or difficulty (negative intercept) parameters. For the `partialCredit` and `gradedResponse` vectors, these should be a __list__ whose length is one less than the number of states. The function `defaultAlphas` will return a vector of betas the appropriate length. * `PnodeQ` -- With the `partialCredit` link function, only a subset of the parent variables need to be included in each transtion. See the [DPCGadget](DPCGadget.Rmd) for details. * `PnodePriorWeight`, `PnodePostWeight`, `GetPriorWeight` --- These functions are related to doing learning, they are weights (equivalent to numbers of observations) which should be given to the prior distribution when fitting the models. * `BuildTable` --- This is the workhorse function which takes the parameters desribed above. The \dQuote{constructor} `Pnode()` is not a real constructor function, rather it takes a node object, and sets the default properties so that it functions as a `Pnode`. This will set the rule and link properties (what is shown below is the default). ```{r defaultParameters} # Usual way to set rules is in constructor partial3 <- Pnode(partial3,rules="Compensatory", link="normalLink") ``` We can set either the alphas, or the log of the alphas. The default is to set the log alphas to 0, which makes the alphas 1. These are set as a vector whose names should match the parents ```{r setAlpha} PnodeAlphas(partial3) PnodeLnAlphas(partial3) PnodeParentNames(partial3) PnodeAlphas(partial3) <- c(theta1=1.2,theta2=.8) PnodeLnAlphas(partial3) ``` The betas are set as a list. This should correspond to the states of the node, with the last (smallest value) omitted. ```{r setBeta} PnodeBetas(partial3) PnodeStates(partial3) PnodeBetas(partial3) <- list(FullCredit=1.0,PartialCredit=0.0) PnodeBetas(partial3) ``` The function `BuildTable()` will produce the conditional probability table. Although we do not use the prior weights unless we are calibrating the model, the `BuildTable` function checks that they are set, so we need to do that before building the table. The `PnodeProbs()` function returns an array; the function `CPTtools::as.CPF()` converts it into the more compact data frame representations. ## Link Scale Parameter The link scale parameter is the residual standard deviation. In general, as the standard deviation of the reference distribution is 1, values around one will make the distribution flat (although not uniform unless the expected value is zero). Values above one will favor the extremes, and values less than one will concentrate the distribution in the middle. The `PnodeLinkScale` function will set the link scale parameter. ```{r linkScale} PnodeLinkScale(partial3) <- .5 ``` This is generally a difficult parameter for experts to work with. ## Building. ```{r buildit} PnodePriorWeight(partial3) <- 10 BuildTable(partial3) as.CPF(PnodeProbs(partial3)) ``` The function `barchart.CPF` (which extends the lattice function `barchart`) will build a visualization of the CPF. The `baseCol` argument can be any R color specification, it is then used as the base color for the graph. ```{r barchart} barchart.CPF(as.CPF(PnodeProbs(partial3)),baseCol="royalblue") ``` ## The Gadget This is a lot to remember. The `Peanut` package provides a shiny gadget which allows an analyst to adjust the parameters of the model, and then save them back to the node object if the analyst finds a set of values that the analyst likes. The `RegressionGadget` makes the following assumptions about the node: * The rule used has the same arguments as `Compensatory`, that is it expects one alpha per parent, and a scaler beta. That is, it is one of `Compensatory`, `Conjunctive` or `Disjunctive`. * The link function is `normalLink`, and a link scale parameter is needed. * The number of relevant alphas is the same for each state of the child variable, and the alpha values are the same. Adjusting the parameter causes the display of the conditional probability table to adjust itself to the new parameters. The tabs allow the view to switch to the table or a table of the effective thetas. Pressing the "OK" button will cause the values of the node being edited to be adjusted to the new values. Pressing "Cancel" will raise a cancel error, and the node parameters will not change. This is invoked using the following command: ```{r gadget, eval=FALSE} partial3 <- RegressionGadget(partial3) ``` [As Pnode objects are reference (R6) objects not functional objects, `partial3` would be modified without the assignment. However, putting the assignment is good style as it reminds R programers that `partial3` is modified.] The function `CompensatoryGadget` opens the gadget in a web browser or as a special web page. The function `MakeCompensatoryGadget` returns the guts of the gadget and allows it to be embedded in a Rmarkdown document. So you can play with the compensatory gadget below. ```{css, echo = FALSE} .shiny-frame{height: 1000px;} ``` ```{r shinyapp, fig.height=8, out.height=8, echo=FALSE} gadget <- MakeRegressionGadget(partial3) shinyApp(gadget$ui,gadget$server,options(height="2000px")) ``` ## Other Gadgets * [CompensatoryGadget](CompensatoryGadget.Rmd) --- Simple applications of the `Compensatory` rule where there is a single discrimination (alpha) parameter for each parent variable. * [OffsetGadget](OffsetGadget.Rmd) --- Simple applications of the `OffsetConjunctive` and `OffsetDisjunctive` rules, where there is a different demand (beta) parameter for each parent variable. * [RegressionGadget](RegressionGadget.Rmd) --- For situation in which the `normalLink` is used. * [DPCGadget](DPCGadget.Rmd) --- Fully featured gadget which allows for different alphas and betas for all parent--child state combinations.