Add partial credit (#14)

* Adding partial credit model.

* Added multidimensional pcm model to the main class.

* Unittested partial credit model.

README.md

@@ -40,11 +40,13 @@ pip install . -t $PYTHONPATH --upgrade
 * 2PL Model
 * 3PL Model
 * Graded Response Model
+* Partial Credit Model
 
 **Multi-dimensional**
 
 * 2PL Model
 * Graded Response Model
+* Partial Credit Model
 
 # Usage
 

girth_mcmc/dichotomous/onepl_model.py

@@ -1,6 +1,6 @@
 import pymc3 as pm
 
-from girth_mcmc.utils import Rayleigh
+from girth_mcmc.distributions import Rayleigh
 
 
 __all__ = ['onepl_model', 'onepl_parameters']

girth_mcmc/dichotomous/threepl_model.py

@@ -1,6 +1,6 @@
 import pymc3 as pm
 
-from girth_mcmc.utils import Rayleigh
+from girth_mcmc.distributions import Rayleigh
 
 
 __all__ = ["threepl_model", "threepl_parameters"]

girth_mcmc/dichotomous/twopl_model.py

@@ -1,6 +1,6 @@
 import pymc3 as pm
 
-from girth_mcmc.utils import Rayleigh
+from girth_mcmc.distributions import Rayleigh
 
 
 __all__ = ["twopl_model", "twopl_parameters"]

girth_mcmc/distributions/__init__.py

@@ -0,0 +1,2 @@
+from .rayleigh import *
+from .partial_credit import *

girth_mcmc/distributions/partial_credit.py

@@ -0,0 +1,31 @@
+import numpy as np
+import theano.tensor as tt
+
+from theano.tensor.nnet import softmax
+from theano.tensor.extra_ops import cumsum
+from pymc3.theanof import floatX
+
+from pymc3.distributions.discrete import Categorical
+
+
+__all__ = ['PartialCredit']
+
+
+class PartialCredit(Categorical):
+    """Computed the probability for the partial credit model given a set of
+    cutpoints and observations.
+    """
+
+    def __init__(self, eta, cutpoints, *args, **kwargs):
+        eta = tt.as_tensor_variable(floatX(eta))
+        cutpoints = tt.concatenate(
+            [
+                tt.as_tensor_variable([0.0]),
+                tt.as_tensor_variable(cutpoints)
+            ])
+        cutpoints = tt.shape_padaxis(cutpoints, 0)
+        eta = tt.shape_padaxis(eta, 1)
+
+        p = softmax(cumsum(eta - cutpoints, axis=1))
+
+        super().__init__(p=p, *args, **kwargs)

girth_mcmc/utils/rayleigh.py → girth_mcmc/distributions/rayleigh.py

@@ -2,6 +2,7 @@
 
 import pymc3 as pm
 import theano.tensor as tt
+
 from pymc3.distributions.distribution import draw_values, generate_samples
 from pymc3.distributions.dist_math import bound
 

girth_mcmc/girth_class.py

@@ -3,23 +3,33 @@
 import pymc3 as pm
 
 from girth_mcmc.utils import validate_mcmc_options
-from girth_mcmc.dichotomous import (rasch_model, rasch_parameters,
-                                    onepl_model, onepl_parameters,
-                                    twopl_model, twopl_parameters,
-                                    multidimensional_twopl_model, multidimensional_twopl_parameters,
-                                    multidimensional_twopl_initial_guess, threepl_model,
-                                    threepl_parameters)
-from girth_mcmc.polytomous import (graded_response_model, graded_response_parameters,
-                                   multidimensional_graded_model,
-                                   multidimensional_graded_parameters)
+from girth_mcmc.dichotomous import (
+    rasch_model, rasch_parameters,
+    onepl_model, onepl_parameters,
+    twopl_model, twopl_parameters,
+    multidimensional_twopl_model, multidimensional_twopl_parameters,
+    multidimensional_twopl_initial_guess, threepl_model,
+    threepl_parameters)
+
+from girth_mcmc.polytomous import (
+    graded_response_model, 
+    graded_response_parameters,
+    multidimensional_graded_model,
+    multidimensional_graded_parameters,
+    partial_credit_model,
+    multidimensional_credit_model
+    )
 
 
 class GirthMCMC(object):
     """GIRTH MCMC class to run estimation models using PyMC3.
 
     Parameters:
-        model: (string) ['Rasch', '1PL', '2PL', '3PL', 'GRM', '2PL_md', 'GRM_md'] 
-            which model to run
+        model: (string) which model to run
+                        ['Rasch', '1PL', '2PL', '3PL', 
+                         'GRM', 'PCM', 
+                         '2PL_md', 'GRM_md', 'PCM_md'] 
+            
         model_args: (tuple) tuple of arguments to pass to model
         options: (dict) mcmc options dictionary
     
@@ -33,8 +43,8 @@ class GirthMCMC(object):
     Notes:
         'GRM' requires setting the number of levels
         '2PL_md' requires setting the number of factors
+        'GRM_md' and 'PCM_md' require setting the number of categories and factors
     """
-
     def __init__(self, model, model_args=None, options=None):
         """Constructor method to run markov models."""
         self.options = validate_mcmc_options(options)
@@ -49,12 +59,16 @@ def __init__(self, model, model_args=None, options=None):
             '2pl': (twopl_model, twopl_parameters, None),
             '3pl': (threepl_model, threepl_parameters, None),
             'grm': (graded_response_model, graded_response_parameters, None),
+            'pcm': (partial_credit_model, graded_response_parameters, None),
 
             # Multidimensional Models
             '2pl_md': (multidimensional_twopl_model, 
                        multidimensional_twopl_parameters,
                        multidimensional_twopl_initial_guess),
             'grm_md': (multidimensional_graded_model, 
+                       multidimensional_graded_parameters,
+                       lambda x, *y: multidimensional_twopl_initial_guess(x, y[1])),
+            'pcm_md': (multidimensional_credit_model, 
                        multidimensional_graded_parameters,
                        lambda x, *y: multidimensional_twopl_initial_guess(x, y[1]))
         }[model.lower()]

girth_mcmc/polytomous/__init__.py

@@ -1,2 +1,4 @@
 from .graded_response_model import *
-from .multidimensional_grm import *
+from .multidimensional_grm import *
+from .partial_credit_model import *
+from .multidimensional_pcm import *

girth_mcmc/polytomous/graded_response_model.py

@@ -1,7 +1,7 @@
 import pymc3 as pm
 from numpy import linspace, zeros, unique
 
-from girth_mcmc.utils import Rayleigh
+from girth_mcmc.distributions import Rayleigh
 
 
 __all__ = ["graded_response_model", "graded_response_parameters"]

girth_mcmc/polytomous/multidimensional_pcm.py

@@ -0,0 +1,72 @@
+import pymc3 as pm
+from numpy import linspace, zeros, unique
+
+import theano
+from theano import tensor as tt
+
+from girth.multidimensional import initial_guess_md
+from girth_mcmc.utils import get_discrimination_indices
+from girth_mcmc.distributions import PartialCredit
+
+
+__all__= ["multidimensional_credit_model"]
+
+
+def multidimensional_credit_model(dataset, n_categories, n_factors):
+    """Defines the mcmc model for the multidimensional partial credit model.
+    
+    Args:
+        dataset: [n_items, n_participants] 2d array of measured responses
+        n_categories: (int) number of polytomous values (i.e. Number of Likert Levels)
+        n_factors: (int) number of factors to extract
+
+    Returns:
+        model: PyMC3 model to run
+    """
+    if n_factors < 2:
+        raise AssertionError(f"Multidimensional GRM model requires "
+                             f"two or more factors specified!")
+
+    n_items, n_people = dataset.shape
+    n_levels = n_categories - 1
+
+    # Need small deviation in offset to
+    # fit into pymc framework
+    mu_value = linspace(-0.05, 0.05, n_levels)
+
+    # Run through 0, K - 1
+    observed = dataset - dataset.min()
+
+    diagonal_indices, lower_indices = get_discrimination_indices(n_items, n_factors)
+    lower_length = lower_indices[0].shape[0]
+
+    graded_mcmc_model = pm.Model()
+
+    with graded_mcmc_model:
+        # Ability Parameters
+        ability = pm.Normal("Ability", mu=0, sigma=1, shape=(n_factors, n_people))
+
+        # Multidimensional Discrimination
+        discrimination = tt.zeros((n_items, n_factors), dtype=theano.config.floatX)
+        diagonal_discrimination = pm.Lognormal('Diagonal Discrimination', mu=0, 
+                                               sigma=0.25, shape=n_factors)
+        lower_discrimination = pm.Normal('Lower Discrimination', sigma=1, 
+                                          shape=lower_length)
+        discrimination = tt.set_subtensor(discrimination[diagonal_indices], 
+                                          diagonal_discrimination)
+
+        discrimination = tt.set_subtensor(discrimination[lower_indices], 
+                                          lower_discrimination)
+
+        # Threshold multilevel prior
+        sigma_difficulty = pm.HalfNormal('Difficulty_SD', sigma=1, shape=1)
+        for ndx in range(n_items):
+            thresholds = pm.Normal(f"Thresholds{ndx}", mu=mu_value, 
+                                   sigma=sigma_difficulty, shape=n_levels)
+
+            # Compute the log likelihood
+            kernel = pm.math.dot(discrimination[ndx], ability)
+            probabilities = PartialCredit(f'Log_Likelihood{ndx}', cutpoints=thresholds, 
+                                          eta=kernel, observed=observed[ndx])
+
+    return graded_mcmc_model

girth_mcmc/polytomous/partial_credit_model.py

@@ -0,0 +1,54 @@
+import pymc3 as pm
+from numpy import linspace, zeros, unique
+
+from girth_mcmc.distributions import PartialCredit, Rayleigh
+
+
+__all__ = ["partial_credit_model"]
+
+
+def partial_credit_model(dataset, n_categories):
+    """Defines the mcmc model for the partial credit model.
+    
+    Args:
+        dataset: [n_items, n_participants] 2d array of measured responses
+        n_categories: number of polytomous values (i.e. Number of Likert Levels)
+
+    Returns:
+        model: PyMC3 model to run
+    """
+    n_items, n_people = dataset.shape
+    n_levels = n_categories - 1
+
+    # Need small dither in offset to
+    # fit into pymc framework
+    mu_value = linspace(-0.05, 0.05, n_levels)
+
+    # Run through 0, K - 1
+    observed = dataset - dataset.min()
+
+    partial_mcmc_model = pm.Model()
+
+    with partial_mcmc_model:
+        # Ability Parameters
+        ability = pm.Normal("Ability", mu=0, sigma=1, shape=n_people)
+
+        # Discrimination multilevel prior
+        rayleigh_scale = pm.Lognormal("Rayleigh_Scale", mu=0, sigma=1/4, shape=1)
+        discrimination = pm.Bound(Rayleigh, lower=0.25)(name='Discrimination', 
+                                  beta=rayleigh_scale, offset=0.25, shape=n_items)
+
+        # Threshold multilevel prior
+        sigma_difficulty = pm.HalfNormal('Difficulty_SD', sigma=1, shape=1)
+
+        # Possible Unorderd Categories
+        for ndx in range(n_items):
+            thresholds = pm.Normal(f"Thresholds{ndx}", mu=mu_value, 
+                                   sigma=sigma_difficulty, shape=n_levels)
+
+            # Compute the log likelihood
+            kernel = discrimination[ndx] * ability
+            probabilities = PartialCredit(f'Log_Likelihood{ndx}', cutpoints=thresholds, 
+                                          eta=kernel, observed=observed[ndx])
+
+    return partial_mcmc_model

girth_mcmc/utils/__init__.py

@@ -1,4 +1,3 @@
 from .options import *
 from .multidimensional_utils import *
 from .missing_data import *
-from .rayleigh import *

requirements.txt

@@ -1,4 +1,4 @@
 numpy>=1.18.2
 scipy>=1.4.1
 pymc3>=3.10.0
-girth>=0.5.0
+girth>=0.7.0

setup.py

@@ -10,9 +10,9 @@
     setup(
         name="girth_mcmc", 
         packages=['girth_mcmc', 'girth_mcmc.dichotomous', 'girth_mcmc.polytomous', 
-                  'girth_mcmc.utils'],
+                  'girth_mcmc.utils', 'girth_mcmc.distributions'],
         package_dir={'girth_mcmc': 'girth_mcmc'},
-        version="0.5.0",
+        version="0.5.1",
         license="MIT",
         description="Bayesian Item Response Theory Estimation.",
         long_description=long_description.replace('<ins>','').replace('</ins>',''),