---
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.
  
  

