Repository: eribean/girth

Issue:
Polytomous Missing Values Handling
Allow missing values in polytomous datasets for parameter estimation.

Discussion:

Patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
index 6179386..783e86b 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -6,7 +6,7 @@ orbs:
 jobs:
   build-and-test:
     docker:
-    - image: circleci/python:3.7
+    - image: circleci/python:3.8
     steps:
       - checkout
       - restore_cache:
diff --git a/README.md b/README.md
index 39c12f6..ce49f63 100644
--- a/README.md
+++ b/README.md
@@ -16,9 +16,13 @@ Girth is a python package for estimating item response theory (IRT) parameters.
    * Joint Maximum Likelihood
    * Conditional Likelihood
    * Marginal Maximum Likelihood
-2. One / Two Parameter Logistic Models
+2. One Parameter Logistic Models
    * Joint Maximum Likelihood
    * Marginal Maximum Likelihood
+3. Two Parameter Logistic Models
+   * Joint Maximum Likelihood
+   * Marginal Maximum Likelihood
+   * Mixed Expected A Prior / Marginal Maximum Likelihood
 3. Three Parameter Logistic Models
    * Marginal Maximum Likelihood (No Optimization and Minimal Support)
 
@@ -26,6 +30,7 @@ Girth is a python package for estimating item response theory (IRT) parameters.
 1. Graded Response Model
    * Joint Maximum Likelihood
    * Marginal Maximum Likelihood
+   * Mixed Expected A Prior / Marginal Maximum Likelihood
 2. Partial Credit Model
    * Joint Maximum Likelihood
    * Marginal Maximum Likelihood
diff --git a/girth/__init__.py b/girth/__init__.py
index e3b39b4..6ef4f2f 100644
--- a/girth/__init__.py
+++ b/girth/__init__.py
@@ -2,7 +2,8 @@ from .numba_functions import numba_expit, _compute_partial_integral, _array_LUT
 from .utils import (validate_estimation_options, irt_evaluation,
                     trim_response_set_and_counts,
                     convert_responses_to_kernel_sign,
-                    get_true_false_counts, mml_approx, create_beta_LUT)
+                    get_true_false_counts, mml_approx, create_beta_LUT,
+                    tag_missing_data, INVALID_RESPONSE)
 from .polytomous_utils import condition_polytomous_response
 from .synthetic import (create_correlated_abilities,
                         create_synthetic_irt_dichotomous,
diff --git a/girth/jml_methods.py b/girth/jml_methods.py
index f09bc98..aa1b8aa 100644
--- a/girth/jml_methods.py
+++ b/girth/jml_methods.py
@@ -220,7 +220,9 @@ def grm_jml(dataset, options=None):
     """
     options = validate_estimation_options(options)
 
-    responses, item_counts = condition_polytomous_response(dataset)
+    cpr_result = condition_polytomous_response(dataset)
+    responses, item_counts, valid_response_mask = cpr_result
+    invalid_response_mask = ~valid_response_mask    
     n_items, n_takers = responses.shape
 
     # Set initial parameter estimates to default
@@ -241,6 +243,11 @@ def grm_jml(dataset, options=None):
     betas_roll = np.roll(betas, -1)
     betas_roll[cumulative_item_counts-1] = 10000
 
+    # Set invalid index to zero, this allows minimal
+    # changes for invalid data and it is corrected
+    # during integration
+    responses[invalid_response_mask] = 0
+    
     for iteration in range(options['max_iteration']):
         previous_betas = betas.copy()
 
@@ -256,7 +263,7 @@ def grm_jml(dataset, options=None):
                                irt_evaluation(betas_roll, discrimination, theta))
 
                 values = graded_prob[responses[:, ndx]]
-                return -np.log(values).sum()
+                return -np.log(values[valid_response_mask[:, ndx]] + 1e-313).sum()
 
             thetas[ndx] = fminbound(_theta_min, -6, 6)
 
@@ -284,9 +291,9 @@ def grm_jml(dataset, options=None):
                                irt_evaluation(betas_roll, discrimination, thetas))
 
                 values = np.take_along_axis(
-                    graded_prob, responses[None, ndx], axis=0)
+                    graded_prob, responses[None, ndx], axis=0).squeeze()
                 np.clip(values, 1e-23, np.inf, out=values)
-                return -np.log(values).sum()
+                return -np.log(values[valid_response_mask[ndx]]).sum()
 
             # Solves jointly for parameters using numerical derivatives
             initial_guess = np.concatenate(([discrimination[start_ndx]],
@@ -333,8 +340,9 @@ def pcm_jml(dataset, options=None):
     """
     options = validate_estimation_options(options)
 
-    responses, item_counts = condition_polytomous_response(
-        dataset, _reference=0.0)
+    cpr_result = condition_polytomous_response(dataset, _reference=0.0)
+    responses, item_counts, valid_response_mask = cpr_result
+    invalid_response_mask = ~valid_response_mask    
     n_items, n_takers = responses.shape
 
     # Set initial parameter estimates to default
@@ -349,6 +357,11 @@ def pcm_jml(dataset, options=None):
         item_length = item_counts[ndx] - 1
         betas[ndx, :item_length] = np.linspace(-1, 1, item_length)
 
+    # Set invalid index to zero, this allows minimal
+    # changes for invalid data and it is corrected
+    # during integration
+    responses[invalid_response_mask] = 0        
+
     for iteration in range(options['max_iteration']):
         previous_discrimination = discrimination.copy()
 
@@ -374,8 +387,8 @@ def pcm_jml(dataset, options=None):
 
                 # Probability associated with response
                 values = np.take_along_axis(
-                    scratch, response_set[:, None], axis=1)
-                return -np.log(values + 1e-23).sum()
+                    scratch, response_set[:, None], axis=1).squeeze()
+                return -np.log(values[valid_response_mask[:, ndx]] + 1e-313).sum()
 
             thetas[ndx] = fminbound(_theta_min, -6, 6, args=(scratch,))
 
@@ -403,8 +416,8 @@ def pcm_jml(dataset, options=None):
                 kernel /= np.nansum(kernel, axis=1)[:, None]
                 # Probability associated with response
                 values = np.take_along_axis(
-                    kernel, response_set[:, None], axis=1)
-                return -np.log(values).sum()
+                    kernel, response_set[:, None], axis=1).squeeze()
+                return -np.log(values[valid_response_mask[ndx]]).sum()
 
             # Solves jointly for parameters using numerical derivatives
             initial_guess = np.concatenate(([discrimination[ndx]],
diff --git a/girth/mml_eap_methods.py b/girth/mml_eap_methods.py
index 8d55e4a..f72c4ba 100644
--- a/girth/mml_eap_methods.py
+++ b/girth/mml_eap_methods.py
@@ -77,7 +77,10 @@ def grm_mml_eap(dataset, options=None):
     """
     options = validate_estimation_options(options)
 
-    responses, item_counts = condition_polytomous_response(dataset, trim_ends=False)
+    cpr_result = condition_polytomous_response(dataset, trim_ends=False)
+    responses, item_counts, valid_response_mask = cpr_result
+    invalid_response_mask = ~valid_response_mask
+
     n_items = responses.shape[0]
     
     # Only use LUT
@@ -91,7 +94,7 @@ def grm_mml_eap(dataset, options=None):
     # Compute the values needed for integral equations
     integral_counts = list()
     for ndx in range(n_items):
-        temp_output = _solve_for_constants(responses[ndx])
+        temp_output = _solve_for_constants(responses[ndx, valid_response_mask[ndx]])
         integral_counts.append(temp_output)
     
     # Initialize difficulty parameters for estimation
@@ -109,6 +112,11 @@ def grm_mml_eap(dataset, options=None):
     betas_roll = np.roll(betas, -1)
     betas_roll[cumulative_item_counts-1] = 10000
     
+    # Set invalid index to zero, this allows minimal
+    # changes for invalid data and it is corrected
+    # during integration
+    responses[invalid_response_mask] = 0    
+
     # Prior Parameters
     ray_scale = 1.0
     eap_options = {'distribution': stats.rayleigh(loc=.25, scale=ray_scale).pdf,
@@ -137,7 +145,8 @@ def grm_mml_eap(dataset, options=None):
         for item_ndx in range(n_items):
             partial_int *= _graded_partial_integral(theta, betas, betas_roll,
                                                     discrimination,
-                                                    responses[item_ndx])
+                                                    responses[item_ndx],
+                                                    invalid_response_mask[item_ndx])
         
         # Estimate the distribution if requested
         distribution_x_weight = latent_pdf(partial_int, iteration)
@@ -168,7 +177,8 @@ def grm_mml_eap(dataset, options=None):
             old_values = _graded_partial_integral(theta, previous_betas,
                                                   previous_betas_roll,
                                                   previous_discrimination,
-                                                  responses[item_ndx])
+                                                  responses[item_ndx],
+                                                  invalid_response_mask[item_ndx])
             partial_int /= old_values
 
             def _local_min_func(estimate):
@@ -183,7 +193,8 @@ def grm_mml_eap(dataset, options=None):
 
                 new_values = _graded_partial_integral(theta, betas, betas_roll,
                                                       discrimination,
-                                                      responses[item_ndx])
+                                                      responses[item_ndx],
+                                                      invalid_response_mask[item_ndx])
 
                 new_values *= partial_int
                 otpt = np.sum(new_values, axis=1)
@@ -206,7 +217,8 @@ def grm_mml_eap(dataset, options=None):
      
             new_values = _graded_partial_integral(theta, betas, betas_roll,
                                                   discrimination,
-                                                  responses[item_ndx])
+                                                  responses[item_ndx],
+                                                  invalid_response_mask[item_ndx])
 
             partial_int *= new_values
             
@@ -230,7 +242,8 @@ def grm_mml_eap(dataset, options=None):
     for item_ndx in range(n_items):
         partial_int *= _graded_partial_integral(theta, betas, betas_roll,
                                                 discrimination,
-                                                responses[item_ndx])
+                                                responses[item_ndx],
+                                                invalid_response_mask[item_ndx])
 
     # Trim difficulties to conform to standard output
     # TODO:  look where missing values are and place NAN there instead
diff --git a/girth/mml_full_methods.py b/girth/mml_full_methods.py
index ab06168..3b153d7 100644
--- a/girth/mml_full_methods.py
+++ b/girth/mml_full_methods.py
@@ -243,8 +243,10 @@ def pcm_mml(dataset, options=None):
     """
     options = validate_estimation_options(options)
 
-    responses, item_counts = condition_polytomous_response(dataset, trim_ends=False,
-                                                           _reference=0.0)
+    cpr_result =  condition_polytomous_response(dataset, trim_ends=False, _reference=0.0)
+    responses, item_counts, valid_response_mask = cpr_result
+    invalid_response_mask = ~valid_response_mask
+
     n_items = responses.shape[0]
 
     # Quadrature Locations
@@ -262,6 +264,11 @@ def pcm_mml(dataset, options=None):
     for ndx in range(n_items):
         betas[ndx, 1:item_counts[ndx]] = np.linspace(-1, 1, item_counts[ndx]-1)
 
+    # Set invalid index to zero, this allows minimal
+    # changes for invalid data and it is corrected
+    # during integration
+    responses[invalid_response_mask] = 0
+
     #############
     # 1. Start the iteration loop
     # 2. Estimate Dicriminatin/Difficulty Jointly
@@ -279,7 +286,8 @@ def pcm_mml(dataset, options=None):
         for item_ndx in range(n_items):
             partial_int *= _credit_partial_integral(theta, betas[item_ndx],
                                                     discrimination[item_ndx],
-                                                    responses[item_ndx])
+                                                    responses[item_ndx],
+                                                    invalid_response_mask[item_ndx])
         # Estimate the distribution if requested
         distribution_x_weight = latent_pdf(partial_int, iteration)
         partial_int *= distribution_x_weight        
@@ -293,14 +301,16 @@ def pcm_mml(dataset, options=None):
             # Remove the previous output
             old_values = _credit_partial_integral(theta, previous_betas[item_ndx],
                                                   previous_discrimination[item_ndx],
-                                                  responses[item_ndx])
+                                                  responses[item_ndx],
+                                                  invalid_response_mask[item_ndx])
             partial_int /= old_values
 
             def _local_min_func(estimate):
                 new_betas[1:] = estimate[1:]
                 new_values = _credit_partial_integral(theta, new_betas,
                                                       estimate[0],
-                                                      responses[item_ndx])
+                                                      responses[item_ndx],
+                                                      invalid_response_mask[item_ndx])
                 new_values *= partial_int
                 otpt = np.sum(new_values, axis=1)
                 return -np.log(otpt).sum()
@@ -318,7 +328,8 @@ def pcm_mml(dataset, options=None):
 
             new_values = _credit_partial_integral(theta, betas[item_ndx],
                                                   discrimination[item_ndx],
-                                                  responses[item_ndx])
+                                                  responses[item_ndx],
+                                                  invalid_response_mask[item_ndx])
 
             partial_int *= new_values
 
@@ -330,7 +341,8 @@ def pcm_mml(dataset, options=None):
     for item_ndx in range(n_items):
         partial_int *= _credit_partial_integral(theta, betas[item_ndx],
                                                 discrimination[item_ndx],
-                                                responses[item_ndx])
+                                                responses[item_ndx],
+                                                invalid_response_mask[item_ndx])
 
     # TODO:  look where missing values are and place NAN there instead
     # of appending them to the end
@@ -382,8 +394,10 @@ def gum_mml(dataset, delta_sign=(0, 1), options=None):
     """
     options = validate_estimation_options(options)
 
-    responses, item_counts = condition_polytomous_response(dataset, trim_ends=False,
-                                                           _reference=0.0)
+    cpr_result = condition_polytomous_response(dataset, trim_ends=False, _reference=0.0)
+    responses, item_counts, valid_response_mask = cpr_result
+    invalid_response_mask = ~valid_response_mask
+
     n_items = responses.shape[0]
 
     # Interpolation Locations
@@ -410,6 +424,11 @@ def gum_mml(dataset, delta_sign=(0, 1), options=None):
     delta_ndx = delta_sign[0]
     delta_multiplier = np.sign(delta_sign[1])
 
+    # Set invalid index to zero, this allows minimal
+    # changes for invalid data and it is corrected
+    # during integration
+    responses[invalid_response_mask] = 0    
+
     #############
     # 1. Start the iteration loop
     # 2. Estimate Dicriminatin/Difficulty Jointly
@@ -430,7 +449,9 @@ def gum_mml(dataset, delta_sign=(0, 1), options=None):
                                                     betas[item_ndx],
                                                     discrimination[item_ndx],
                                                     fold_span[item_ndx],
-                                                    responses[item_ndx])
+                                                    responses[item_ndx],
+                                                    invalid_response_mask[item_ndx])
+
         # Estimate the distribution if requested
         distribution_x_weight = latent_pdf(partial_int, iteration)
         partial_int *= distribution_x_weight
@@ -445,7 +466,8 @@ def gum_mml(dataset, delta_sign=(0, 1), options=None):
                                                   previous_betas[item_ndx],
                                                   previous_discrimination[item_ndx],
                                                   fold_span[item_ndx],
-                                                  responses[item_ndx])
+                                                  responses[item_ndx],
+                                                  invalid_response_mask[item_ndx])
             partial_int /= old_values
             
             new_betas = np.full((betas.shape[1],), np.nan)
@@ -454,7 +476,8 @@ def gum_mml(dataset, delta_sign=(0, 1), options=None):
                 new_values = _unfold_partial_integral(theta, estimate[1],
                                                       new_betas,
                                                       estimate[0], fold_span[item_ndx],
-                                                      responses[item_ndx])
+                                                      responses[item_ndx], 
+                                                      invalid_response_mask[item_ndx])
 
                 new_values *= partial_int
                 otpt = np.sum(new_values, axis=1)
@@ -477,7 +500,8 @@ def gum_mml(dataset, delta_sign=(0, 1), options=None):
                                                   betas[item_ndx],
                                                   discrimination[item_ndx],
                                                   fold_span[item_ndx],
-                                                  responses[item_ndx])
+                                                  responses[item_ndx],
+                                                  invalid_response_mask[item_ndx])
 
             partial_int *= new_values
 
@@ -494,7 +518,8 @@ def gum_mml(dataset, delta_sign=(0, 1), options=None):
                                                 betas[item_ndx],
                                                 discrimination[item_ndx],
                                                 fold_span[item_ndx],
-                                                responses[item_ndx])    
+                                                responses[item_ndx],
+                                                invalid_response_mask[item_ndx])    
     # Compute statistics for final iteration
     null_metrics = latent_pdf.compute_metrics(partial_int, latent_pdf.null_distribution * 
                                               latent_pdf.weights, 0)
diff --git a/girth/mml_methods.py b/girth/mml_methods.py
index a14295a..6b25bcd 100644
--- a/girth/mml_methods.py
+++ b/girth/mml_methods.py
@@ -236,7 +236,10 @@ def grm_mml(dataset, options=None):
     """
     options = validate_estimation_options(options)
 
-    responses, item_counts = condition_polytomous_response(dataset, trim_ends=False)
+    cpr_result = condition_polytomous_response(dataset, trim_ends=False)
+    responses, item_counts, valid_response_mask = cpr_result
+    invalid_response_mask = ~valid_response_mask
+    
     n_items = responses.shape[0]
     
     # Should we use the LUT
@@ -253,7 +256,7 @@ def grm_mml(dataset, options=None):
     # Compute the values needed for integral equations
     integral_counts = list()
     for ndx in range(n_items):
-        temp_output = _solve_for_constants(responses[ndx])
+        temp_output = _solve_for_constants(responses[ndx, valid_response_mask[ndx]])
         integral_counts.append(temp_output)
     
     # Initialize difficulty parameters for estimation
@@ -270,6 +273,11 @@ def grm_mml(dataset, options=None):
                                                item_counts[ndx] - 1)
     betas_roll = np.roll(betas, -1)
     betas_roll[cumulative_item_counts-1] = 10000
+
+    # Set invalid index to zero, this allows minimal
+    # changes for invalid data and it is corrected
+    # during integration
+    responses[invalid_response_mask] = 0
     
     #############
     # 1. Start the iteration loop
@@ -288,7 +296,8 @@ def grm_mml(dataset, options=None):
         for item_ndx in range(n_items):
             partial_int *= _graded_partial_integral(theta, betas, betas_roll,
                                                     discrimination,
-                                                    responses[item_ndx])
+                                                    responses[item_ndx],
+                                                    invalid_response_mask[item_ndx])
         
         # Estimate the distribution if requested
         distribution_x_weight = latent_pdf(partial_int, iteration)
@@ -314,7 +323,8 @@ def grm_mml(dataset, options=None):
             old_values = _graded_partial_integral(theta, previous_betas,
                                                   previous_betas_roll,
                                                   previous_discrimination,
-                                                  responses[item_ndx])
+                                                  responses[item_ndx],
+                                                  invalid_response_mask[item_ndx])
             partial_int /= old_values
 
             def _local_min_func(estimate):
@@ -329,7 +339,8 @@ def grm_mml(dataset, options=None):
 
                 new_values = _graded_partial_integral(theta, betas, betas_roll,
                                                       discrimination,
-                                                      responses[item_ndx])
+                                                      responses[item_ndx],
+                                                      invalid_response_mask[item_ndx])
 
                 new_values *= partial_int
                 otpt = np.sum(new_values, axis=1)
@@ -341,7 +352,8 @@ def grm_mml(dataset, options=None):
 
             new_values = _graded_partial_integral(theta, betas, betas_roll,
                                                   discrimination,
-                                                  responses[item_ndx])
+                                                  responses[item_ndx],
+                                                  invalid_response_mask[item_ndx])
 
             partial_int *= new_values
 
@@ -353,7 +365,8 @@ def grm_mml(dataset, options=None):
     for item_ndx in range(n_items):
         partial_int *= _graded_partial_integral(theta, betas, betas_roll,
                                                 discrimination,
-                                                responses[item_ndx])
+                                                responses[item_ndx],
+                                                invalid_response_mask[item_ndx])
 
     # Trim difficulties to conform to standard output
     # TODO:  look where missing values are and place NAN there instead
diff --git a/girth/polytomous_utils.py b/girth/polytomous_utils.py
index 8c5565a..75e4620 100644
--- a/girth/polytomous_utils.py
+++ b/girth/polytomous_utils.py
@@ -2,6 +2,7 @@ import numpy as np
 from scipy import integrate
 from scipy.optimize import fminbound
 
+from girth import INVALID_RESPONSE
 from girth.numba_functions import numba_expit
 
 
@@ -19,16 +20,24 @@ def condition_polytomous_response(dataset, trim_ends=True, _reference=1.0):
     Returns:
         dataset: (2d array) ordinal values converted to linear indices
         beta_length: (1d array) number of unique values per item
+        valid_response_mask: (2d_array) boolean mask of valid responses
+                             that is the same size as dataset
     """
     # Remove all no / yes endorsements
-    min_value, max_value = dataset.min(), dataset.max()
+    valid_response_mask = dataset != INVALID_RESPONSE
+    min_value = np.min(dataset[valid_response_mask])
+    max_value = np.max(dataset[valid_response_mask])
     n_items = dataset.shape[0]
 
     if trim_ends:
-        raw_score_sums = dataset.sum(0)
-        mask = ((raw_score_sums == (n_items * min_value)) | 
-                (raw_score_sums == (n_items * max_value)))
+        raw_score_sums = np.sum(dataset, where=valid_response_mask, axis=0)
+        valid_score_sums = np.count_nonzero(valid_response_mask, axis=0)
+        ratio_score_sums = raw_score_sums / valid_score_sums
+
+        mask = ((ratio_score_sums == min_value) | 
+                (ratio_score_sums == max_value))
         dataset = dataset[:, ~mask]
+        valid_response_mask = dataset != INVALID_RESPONSE
     
     betas_length = np.zeros((n_items,), dtype='int')
     the_output = dataset.copy()
@@ -38,7 +47,8 @@ def condition_polytomous_response(dataset, trim_ends=True, _reference=1.0):
     # responses, and replace with linear indexing
     cnt = 0
     for ndx, item in enumerate(the_output):
-        values, indices = np.unique(item, return_inverse=True)
+        values, indices = np.unique(item[valid_response_mask[ndx]], 
+                                    return_inverse=True)
         betas_length[ndx] = max(values.size, 2)
 
         # Recode from zero to N-1
@@ -46,12 +56,15 @@ def condition_polytomous_response(dataset, trim_ends=True, _reference=1.0):
         if values.size == 1:
             recode = np.array([values[0] != 0], dtype='int')
         values = recode + cnt * _reference
-        the_output[ndx] = values[indices]
+        the_output[ndx, valid_response_mask[ndx]] = values[indices]
 
         # Update linear index
         cnt += betas_length[ndx]
+    
+    # Reset the invalid responses
+    the_output[~valid_response_mask] = INVALID_RESPONSE
 
-    return the_output, betas_length
+    return the_output, betas_length, valid_response_mask
 
 
 def _solve_for_constants(item_responses):
@@ -71,14 +84,19 @@ def _solve_for_constants(item_responses):
     
 
 def _graded_partial_integral(theta, betas, betas_roll,
-                             discrimination, responses):
+                             discrimination, response_set,
+                             invalid_response_mask):
     """Computes the partial integral for the graded response."""
     temp1 = (betas[:, None] - theta) * discrimination[:, None]
     temp2 = (betas_roll[:, None] - theta) * discrimination[:, None]
     graded_prob = numba_expit(temp1) 
     graded_prob -= numba_expit(temp2)
 
-    return graded_prob[responses, :]
+    # Set all the responses and fix afterward
+    temp_output = graded_prob[response_set, :]
+    temp_output[invalid_response_mask] = 1.0
+
+    return temp_output
 
 
 def _solve_integral_equations_LUT(discrimination, ratio, _, __, interpolate_function):
@@ -104,7 +122,7 @@ def _solve_integral_equations(discrimination, ratio, distribution, theta, _):
 
 
 def _credit_partial_integral(theta, betas, discrimination, 
-                             response_set):
+                             response_set, invalid_response_mask):
     """Computes the partial integral for the partial credit model."""
     # Creates a 2d array [beta x thetas]
     kernel = theta[None, :] - betas[:, None]
@@ -121,12 +139,16 @@ def _credit_partial_integral(theta, betas, discrimination,
     # Normalize probability to equal one
     kernel /= np.nansum(kernel, axis=0)[None, :]
 
-    return kernel[response_set, :]
+    # Set all the responses and fix afterward
+    temp_output = kernel[response_set, :]
+    temp_output[invalid_response_mask] = 1.0
+
+    return temp_output
 
 
 def _unfold_partial_integral(theta, delta, betas, 
                              discrimination, fold_span,
-                             response_set):
+                             response_set, invalid_response_mask):
     """Computes the partial integral for the _GGUM model."""
     # Unfolding_Model
     thresholds = np.exp(-discrimination * np.cumsum(betas))
@@ -137,4 +159,8 @@ def _unfold_partial_integral(theta, delta, betas,
     kernel[1:, :] *= thresholds[:, None]
     kernel /= np.nansum(kernel, axis=0)
 
-    return kernel[response_set, :]
\ No newline at end of file
+    # Set all the responses and fix afterward
+    temp_output = kernel[response_set, :]
+    temp_output[invalid_response_mask] = 1.0
+
+    return temp_output
\ No newline at end of file
diff --git a/girth/synthetic.py b/girth/synthetic.py
index 31e8ce2..cf3fa83 100644
--- a/girth/synthetic.py
+++ b/girth/synthetic.py
@@ -58,7 +58,7 @@ def create_synthetic_irt_dichotomous(difficulty, discrimination, thetas,
     # convert to binary based on probability
     random_compare = np.random.rand(*continuous_output.shape)
 
-    return random_compare <= continuous_output
+    return (random_compare <= continuous_output).astype('int')
 
 
 def create_synthetic_mirt_dichotomous(difficulty, discrimination, thetas,
@@ -103,7 +103,7 @@ def create_synthetic_mirt_dichotomous(difficulty, discrimination, thetas,
     # convert to binary based on probability
     random_compare = np.random.rand(*continuous_output.shape)
 
-    return random_compare <= continuous_output
+    return (random_compare <= continuous_output).astype('int')
 
 
 # Private functions for polytomous outputs
diff --git a/girth/utils.py b/girth/utils.py
index 6846eed..7cb94c6 100644
--- a/girth/utils.py
+++ b/girth/utils.py
@@ -6,6 +6,10 @@ from scipy.special import roots_legendre
 from girth import _array_LUT
 
 
+# There is no int nan, so use this as a placeholder
+INVALID_RESPONSE = -99999
+
+
 def default_options():
     """ Dictionary of options used in Girth.
 
@@ -85,6 +89,24 @@ def validate_estimation_options(options_dict=None):
     return full_options
 
 
+def tag_missing_data(dataset, valid_responses):
+    """Checks the data for valid responses.
+    
+    Args:
+        dataset: (array) array to validate
+        valid_responses: (array-like) list of valid responses
+        
+    Returns:
+        updated_dataset: (array) data that holds only valid_responses and
+                         invalid_fill
+    """
+    mask = np.isin(dataset, valid_responses)
+    output = dataset.copy()
+    output[~mask] = INVALID_RESPONSE
+    
+    return output
+
+
 def get_true_false_counts(responses):
     """ Returns the number of true and false for each item.
 


Test patch:
diff --git a/testing/test_missing_data.py b/testing/test_missing_data.py
new file mode 100644
index 0000000..9141ec8
--- /dev/null
+++ b/testing/test_missing_data.py
@@ -0,0 +1,252 @@
+import unittest  # pylint: disable=cyclic-import
+
+import numpy as np
+
+from girth import tag_missing_data
+
+from girth import create_synthetic_irt_dichotomous, create_synthetic_irt_polytomous
+from girth import rasch_jml, onepl_jml, twopl_jml
+from girth import rasch_mml, onepl_mml, twopl_mml
+from girth import rasch_full, onepl_full, twopl_full
+
+from girth import grm_jml, pcm_jml
+from girth import grm_mml_eap, grm_mml, pcm_mml, gum_mml
+
+
+def _create_missing_data(dataset, seed, threshold):
+    """Creates missing data for unittests."""
+    np.random.seed(seed)
+    mask = np.random.rand(*dataset.shape) < threshold
+    dataset2 = dataset.copy()
+    dataset2[mask] = 2222
+    return dataset2
+
+
+def _rmse(dataset1, dataset2):
+    return np.sqrt(np.square(dataset1 - dataset2).mean())
+
+
+class TestMissingDichotomous(unittest.TestCase):
+    """Testing Dichotomous methods for missing data."""
+
+    ## SMOKE / REGRESSION TESTS with human-in-loop spot check
+
+    def test_rasch_jml_mml(self):
+        """Testing rasch mml/jml for missing data."""
+        np.random.seed(7321)
+        n_items = 10
+        n_people = 200
+
+        difficulty = np.random.randn(n_items)
+        theta = np.random.randn(n_people)
+        syn_data = create_synthetic_irt_dichotomous(difficulty, 1.0, theta)
+        syn_data_tagged = _create_missing_data(syn_data, 73435, 0.1)
+        syn_data_missing = tag_missing_data(syn_data_tagged, [0, 1])
+
+        # JML
+        result_all_good = rasch_jml(syn_data)
+        result_missing = rasch_jml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.100458, 4)
+
+        # MML
+        result_all_good = rasch_mml(syn_data)
+        result_missing = rasch_mml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.054726, 4)
+
+        # MML FULL
+        result_all_good = rasch_full(syn_data)
+        result_missing = rasch_full(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.0590626, 4)        
+
+    def test_onepl_jml_mml(self):
+        """Testing oneple mml/jml for missing data."""
+        np.random.seed(789595)
+        n_items = 10
+        n_people = 200
+
+        discrimination = 1.53
+        difficulty = np.random.randn(n_items)
+        theta = np.random.randn(n_people)
+        syn_data = create_synthetic_irt_dichotomous(difficulty, discrimination, theta)
+        syn_data_tagged = _create_missing_data(syn_data, 66877, 0.1)
+        syn_data_missing = tag_missing_data(syn_data_tagged, [0, 1])
+
+        # JML (This is weirdly bad)
+        result_all_good = onepl_jml(syn_data)
+        result_missing = onepl_jml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.862138914, 4)        
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 2.99675, 4)        
+
+        # MML
+        result_all_good = onepl_mml(syn_data)
+        result_missing = onepl_mml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.031971, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.045726, 4)
+
+        # MML FULL
+        result_all_good = onepl_full(syn_data)
+        result_missing = onepl_full(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.038550, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.052542, 4)
+
+    def test_twopl_jml_mml(self):
+        """Testing rasch mml/jml for missing data."""
+        np.random.seed(15335)
+        n_items = 10
+        n_people = 200
+
+        discrimination = np.sqrt(-2 * np.log(np.random.rand(n_items)))
+        difficulty = np.random.randn(n_items)
+        theta = np.random.randn(n_people)
+
+        syn_data = create_synthetic_irt_dichotomous(difficulty, discrimination, theta)
+        syn_data_tagged = _create_missing_data(syn_data, 84433, 0.1)
+        syn_data_missing = tag_missing_data(syn_data_tagged, [0, 1])
+
+        # JML (This is weirdly bad)
+        result_all_good = twopl_jml(syn_data)
+        result_missing = twopl_jml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.636390, 4)        
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 2.027410, 4)        
+
+        # MML
+        result_all_good = twopl_mml(syn_data)
+        result_missing = twopl_mml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.064496, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.155868, 4)
+
+        # MML FULL
+        result_all_good = twopl_full(syn_data)
+        result_missing = twopl_full(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'], result_missing['Difficulty'])
+        self.assertAlmostEqual(difference_rmse, 0.066202, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.166226, 4)
+
+
+class TestMissingPolytomous(unittest.TestCase):
+    """Testing Polytomous methods for missing data."""    
+
+    def test_grm_jml_mml(self):
+        """Testing graded response mml/jml for missing data."""
+        np.random.seed(69871)
+        n_items = 10
+        n_people = 300
+
+        discrimination = np.sqrt(-2 * np.log(np.random.rand(n_items)))
+        difficulty = np.random.randn(n_items, 3)
+        difficulty = np.sort(difficulty, 1)
+        theta = np.random.randn(n_people)
+
+        syn_data = create_synthetic_irt_polytomous(difficulty, discrimination, 
+                                                   theta, model='grm')
+        syn_data_tagged = _create_missing_data(syn_data, 879858, 0.1)
+        syn_data_missing = tag_missing_data(syn_data_tagged, [1, 2, 3, 4])
+
+        # JML
+        result_all_good = grm_jml(syn_data)
+        result_missing = grm_jml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'].ravel(), 
+                                result_missing['Difficulty'].ravel())
+        self.assertAlmostEqual(difference_rmse, 0.097251, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.112861, 4)
+
+        # MML
+        result_all_good = grm_mml(syn_data)
+        result_missing = grm_mml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'].ravel(), 
+                                result_missing['Difficulty'].ravel())
+        self.assertAlmostEqual(difference_rmse, 0.088626, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.065834, 4)
+
+        # EAP/MML
+        result_all_good = grm_mml_eap(syn_data)
+        result_missing = grm_mml_eap(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'].ravel(), 
+                                result_missing['Difficulty'].ravel())
+        self.assertAlmostEqual(difference_rmse, 0.073556, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.062288, 4)
+
+
+    def test_pcm_jml_mml(self):
+        """Testing partial response mml/jml for missing data."""
+        np.random.seed(499867)
+        n_items = 10
+        n_people = 300
+
+        discrimination = np.sqrt(-2 * np.log(np.random.rand(n_items)))
+        difficulty = np.random.randn(n_items, 3)
+        difficulty = np.sort(difficulty, 1)
+        theta = np.random.randn(n_people)
+
+        syn_data = create_synthetic_irt_polytomous(difficulty, discrimination, 
+                                                   theta, model='pcm')
+        syn_data_tagged = _create_missing_data(syn_data, 1986778, 0.1)
+        syn_data_missing = tag_missing_data(syn_data_tagged, [1, 2, 3, 4])
+
+        # JML
+        result_all_good = pcm_jml(syn_data)
+        result_missing = pcm_jml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'].ravel(), 
+                                result_missing['Difficulty'].ravel())
+        self.assertAlmostEqual(difference_rmse, 0.078015, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.147213, 4)
+        
+        # MML
+        result_all_good = pcm_mml(syn_data)
+        result_missing = pcm_mml(syn_data_missing)
+        difference_rmse = _rmse(result_all_good['Difficulty'].ravel(), 
+                                result_missing['Difficulty'].ravel())
+        self.assertAlmostEqual(difference_rmse, 0.085013, 4)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.079640, 4)
+
+
+    def test_gum_mml(self):
+        """Testing unfolding response mml/jml for missing data."""
+        np.random.seed(7382)
+        n_items = 15
+        n_people = 600
+
+        discrimination = np.sqrt(-2 * np.log(np.random.rand(n_items)))
+        delta = np.random.randn(n_items, 1)
+        difficulty = np.random.randn(n_items, 3)
+        difficulty = np.sort(difficulty, 1)
+        difficulty = np.c_[-difficulty, np.zeros((n_items, 1)), difficulty[:, ::-1]] + delta
+        theta = np.random.randn(n_people)
+
+        syn_data = create_synthetic_irt_polytomous(difficulty, discrimination, 
+                                                   theta, model='gum')
+        syn_data_tagged = _create_missing_data(syn_data, 112, 0.1)
+        syn_data_missing = tag_missing_data(syn_data_tagged, [1, 2, 3, 4])
+
+        # MML
+        result_all_good = gum_mml(syn_data)
+        result_missing = gum_mml(syn_data_missing)
+        print(result_all_good['Discrimination'], 
+              result_missing['Discrimination'])
+
+        difference_rmse = _rmse(result_all_good['Delta'], result_missing['Delta'])
+        self.assertAlmostEqual(difference_rmse, 0.0664, 3)
+        difference_rmse = _rmse(result_all_good['Discrimination'], result_missing['Discrimination'])
+        self.assertAlmostEqual(difference_rmse, 0.1382, 2)
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file
diff --git a/testing/test_utils.py b/testing/test_utils.py
index a89c2ba..27eb1f0 100644
--- a/testing/test_utils.py
+++ b/testing/test_utils.py
@@ -10,7 +10,8 @@ from girth import (create_synthetic_irt_dichotomous,
                    convert_responses_to_kernel_sign,
                    validate_estimation_options, condition_polytomous_response,
                    get_true_false_counts, mml_approx)
-from girth.utils import _get_quadrature_points, default_options, create_beta_LUT
+from girth.utils import (_get_quadrature_points, default_options, 
+                         create_beta_LUT, tag_missing_data, INVALID_RESPONSE)
 from girth.numba_functions import _compute_partial_integral
 from girth.polytomous_utils import (_graded_partial_integral, _solve_for_constants,
                                     _solve_integral_equations, _credit_partial_integral,
@@ -163,6 +164,19 @@ class TestUtilitiesMethods(unittest.TestCase):
         np.testing.assert_array_almost_equal(expected, result)
         np.testing.assert_array_almost_equal(expected, result2)
 
+    def test_tag_missing_values(self):
+        """Testing the tagging missing values."""
+        np.random.seed(148549)
+        synthetic_data = np.random.randint(1, 6, (20, 650))
+        mask = np.random.rand(*synthetic_data.shape) < .15
+        synthetic_data[mask] = 5555
+        
+        updated_data = tag_missing_data(synthetic_data, list(range(1, 6)))
+
+        new_mask = updated_data == INVALID_RESPONSE
+
+        np.testing.assert_equal(new_mask, mask)
+
 
 class TestPolytomousUtilities(unittest.TestCase):
     """Tests the polytomous utilities"""
@@ -194,7 +208,7 @@ class TestPolytomousUtilities(unittest.TestCase):
         self.assertTrue(output[0].std(axis=0)[0] != 0)
 
         # Trim First/Last Column but not last
-        dataset[:, -1] = 1
+        dataset[:, -1] = 5
         output = condition_polytomous_response(dataset, trim_ends=True)
         self.assertTupleEqual(
             output[0].shape, (dataset.shape[0], dataset.shape[1]-2))
@@ -248,7 +262,6 @@ class TestPolytomousUtilities(unittest.TestCase):
         output = _solve_for_constants(dataset[1])        
         self.assertTrue(output == 1)
 
-
     def test_integral_equations(self):
         """Tests solving for integral given a ratio."""
         np.random.seed(786)
@@ -259,8 +272,8 @@ class TestPolytomousUtilities(unittest.TestCase):
         # Compare against dichotomous data
         syn_data = create_synthetic_irt_dichotomous(
             difficulty, discrimination, theta)
-        n0 = np.count_nonzero(~syn_data, axis=1)
-        n1 = np.count_nonzero(syn_data, axis=1)
+        n0 = np.count_nonzero(syn_data == 0, axis=1)
+        n1 = np.count_nonzero(syn_data == 1, axis=1)
         ratio = n1 / (n1 + n0)
 
         theta, weights = _get_quadrature_points(61, -5, 5)
@@ -276,11 +289,13 @@ class TestPolytomousUtilities(unittest.TestCase):
         betas = np.array([-10000, -.3, 0.1, 1.2])
         betas_roll = np.roll(betas, -1)
         betas_roll[-1] = 10000
+        invalid_response_mask = np.zeros_like(responses, dtype='bool')
 
         output = np.ones((responses.shape[1], theta.size))
         for ndx in range(responses.shape[0]):
             output *= _graded_partial_integral(theta, betas, betas_roll,
-                                               np.array([1,]), responses[ndx])
+                                               np.array([1,]), responses[ndx],
+                                               invalid_response_mask[ndx])
 
         # Compare to hand calculations
         hand_calc = list()
@@ -295,12 +310,27 @@ class TestPolytomousUtilities(unittest.TestCase):
 
         np.testing.assert_array_equal(hand_calc, output)
 
+        # Test invalid response
+        invalid_response_mask[0, 1] = True
+        invalid_response_mask[0, 7] = True
+        output = _graded_partial_integral(theta, betas, betas_roll,
+                                          np.array([1,]), responses[0],
+                                          invalid_response_mask[0])
+        
+        np.testing.assert_equal(output[1], np.ones(61,))
+        np.testing.assert_equal(output[7], np.ones(61,))
+
+        with np.testing.assert_raises(AssertionError):
+            for ndx in [0, 2, 3, 4, 5, 6, 8, 9]:
+                np.testing.assert_equal(output[ndx], np.ones(61,))
+
     def test_credit_partial_integration(self):
         """Testing the partial integral in the graded model."""
         theta, _ = _get_quadrature_points(61, -5, 5)
         response_set = np.array([0, 1, 2, 2, 1, 0, 3, 1, 3, 2, 2, 2])
         betas = np.array([0, -0.4, 0.94, -.37])
         discrimination = 1.42
+        invalid_response_mask = np.zeros_like(response_set, dtype='bool')
 
         # Hand calculations
         offsets = np.cumsum(betas)[1:]
@@ -317,16 +347,30 @@ class TestPolytomousUtilities(unittest.TestCase):
             expected[ndx] = probability_values[response]
 
         result = _credit_partial_integral(theta, betas, discrimination,
-                                          response_set)
+                                          response_set, invalid_response_mask)
 
         np.testing.assert_array_almost_equal(result, expected)
 
+        invalid_response_mask[1] = True
+        invalid_response_mask[7] = True
+        result = _credit_partial_integral(theta, betas, discrimination,
+                                          response_set, invalid_response_mask)
+
+        np.testing.assert_equal(result[1], np.ones(61,))
+        np.testing.assert_equal(result[7], np.ones(61,))
+
+        with np.testing.assert_raises(AssertionError):
+            for ndx in [0, 2, 3, 4, 5, 6, 8, 9]:
+                np.testing.assert_equal(result[ndx], np.ones(61,))
+
     def test_unfold_partial_integration(self):
         """Testing the unfolding integral."""
         theta, _ = _get_quadrature_points(61, -5, 5)
         response_set = np.array([0, 1, 2, 2, 1, 0, 3, 1, 3, 2, 2, 2])
         betas = np.array([-1.3, -.4, 0.2])
         delta = -0.76
+        invalid_response_mask = np.zeros_like(response_set, dtype='bool')
+
         # (2N -1) / 2 - n
         folding = 3.5 - np.arange(4)
         discrimination = 1.42
@@ -343,9 +387,22 @@ class TestPolytomousUtilities(unittest.TestCase):
 
         result = _unfold_partial_integral(theta, delta, betas,
                                           discrimination, folding,
-                                          response_set)
+                                          response_set, invalid_response_mask)
         np.testing.assert_array_almost_equal(result, expected)
 
+        invalid_response_mask[1] = True
+        invalid_response_mask[7] = True
+        result = _unfold_partial_integral(theta, delta, betas,
+                                          discrimination, folding,
+                                          response_set, invalid_response_mask)
+
+        np.testing.assert_equal(result[1], np.ones(61,))
+        np.testing.assert_equal(result[7], np.ones(61,))
+
+        with np.testing.assert_raises(AssertionError):
+            for ndx in [0, 2, 3, 4, 5, 6, 8, 9]:
+                np.testing.assert_equal(result[ndx], np.ones(61,))
+
     def test_lut_creation(self):
         """Test the lookup table creation function."""
         lut_func = create_beta_LUT((0.5, 2, 500), (-3, 3, 500))