C:/programs/etirm/src/Start3PL.h

Go to the documentation of this file.

/*! \file Start3PL.h
 
  \brief
  Compute item parameter estimate starting values for one, two, or
  three-parameter logistic model using nonlinear regression for each item
  based on an estimated ability for each examinee. Ability estimates
  are obtained using the PROX method.

  Estimation Toolkit for Item Response Models (ETIRM)
  http://www.smallwaters.com/software/cpp/etirm.html

  Author(s): 
  Werner Wothke, maintenance (http://www.smallwaters.com)
  Brad Hanson (http://www.b-a-h.com/)
  See the file LICENSE for information on usage and redistribution.

  Copyright (C) 2008, Werner Wothke
  Copyright (c) 2000-2001, Bradley A. Hanson
 */

#ifndef ETIRM_START3PL_H_
#define ETIRM_START3PL_H_

#ifdef ETIRM_NO_DIR_PREFIX
#include "etirmtypes.h"
#else
#include "etirm/etirmtypes.h"
#endif

#include <set>
#include <cmath> // for log, sqrt, fabs
// for compilers which do not put C library functions in std namespace
#ifdef BOOST_NO_STDC_NAMESPACE
namespace std
{ using ::sqrt; using ::log; using ::fabs;}
#endif

namespace etirm
{

  const Real invalidLogit = 10000.0; 
  //!< value to indicate invalid logit (e.g., logit(1))

  /*!
    \brief
    Compute mean and standard deviation of logit item difficulty over items 
    responded to by an examinee.
   
    \section template_args Template Parameters
   
    \param RI Type of iterator over examinee responses (iterator returned 
        by responses_begin() for an Examinee object)
    \param II Type of iterator over items (iterator over pointers to Item objects)

    \section function_args Function Parameters
   
    \param[in] resp Iterator pointing to first response of examinee.
    \param[in] iitem  Iterator pointing to first item.
    \param[in] itemdiff Item difficulties.
    \param[in] notPresentedResponse Response which indicates an item was not 
        presented to an examinee. Assumed to be the same for all items.
    \param[out] mean  Mean of abilities for persons responding to the item.
    \param[out] variance  Variance of abilities for persons responding to the item.
   */
  template <class RI, class II> void PersonMoments(RI iresp, II iitem, RealVector &itemdiff,
      Response notPresentedResponse, Real *mean, Real *variance)
  {

    /* Compute logit of examinee proportion correct scores */
    RealVector::iterator idiff = itemdiff.begin();
    int total = 0;
    *mean = 0.0;
    *variance = 0.0;
    for (int i = itemdiff.size(); i--; ++iitem, ++idiff)
    {
      Response resp = iresp[(*iitem)->Index()];
      if (resp != notPresentedResponse && *idiff != invalidLogit)
      {
        *mean += *idiff;
        *variance += *idiff * *idiff;
        ++total;
      }
    }

    *mean /= total;
    *variance /= total;
    *variance -= *mean * *mean;
  }

  /*!
    \brief
    Compute mean and standard deviation of person ability over examinees
    who reponded to an item.
   
    \section template_args Template Parameters
   
    \param EI  Type of iterator over examinees (iterator over pointers to Examinee objects)
    \param RI Type of iterator over examinee responses (iterator returned by responses_begin() for an Examinee object)
   
    \section function_args Function Parameters
   
    \param[in] num_examinees  Number of examinees.
    \param[in] begin_examinee Iterator to pointer to first examinee.
    \param[in] item Index of item in response string to compute moments for (0 offset).
    \param[in] &abilities  Address of ability vector for all persons.
    \param[in] notPresentedResponse Response which indicates an item was not presented
       to an examinee. Assumed to be the same for all items.
    \param[out] mean  Mean of abilities for persons resonding to the item.
    \param[out] variance  Variance of abilities for persons resonding to the item.
   */
  template <class IE, class RI> void ItemMoments(int num_examinees, IE begin_examinee, int item,
      RealVector &abilities, Response notPresentedResponse, double *mean, double *variance)
  {
    /* Compute logit of examinee proportion correct scores */
    RealVector::iterator iabil = abilities.begin();
    IE ie = begin_examinee;
    Real total = 0;
    *mean = 0.0;
    *variance = 0.0;
    for (int i = num_examinees; i--; ++ie, ++iabil)
    {
      RI resp = (*ie)->responses_begin() + item;
      Real count = (*ie)->Count();

      if (*resp != notPresentedResponse && *iabil != invalidLogit)
      {
        *mean += *iabil * count;
        *variance += *iabil * *iabil * count;
        total += count;
      }
    }

    *mean /= total;
    *variance /= total;
    *variance -= *mean * *mean;
  }

  /*!
    \brief
    Compute Rasch item difficulties and examinee abilities using the Normal
    Approximation Estimation Algorithm (PROX). 
    
    The scale of the latent variable is set such that the sum of the item 
    difficulties is zero.
    
    For information on PROX see
   
    Linacre, J. M. (1994). PROX with missing data. Rasch Measurement
        Transactions, 8(3), 378.
   
    Cohen, L. (1979). Approximate expressions for parameter estimates in the Rasch model.
        British Journal of Mathematical and Statistical Psychology, 32(1), 113-120.
   
    \section template_args Template Parameters
       
    \param  EI  Type of iterator over examinees (iterator over pointers to Examinee 
        objects).
    \param  RI  Type of iterator over examinee responses (iterator returned by 
        responses_begin() for an Examinee object).
    \param  II  Type of iterator over items (iterator over pointers to Item objects).

    \section function_args Function Parameters
   
    \param[in]  numItems  Number of items.
    \param[in]  useAll  If true calculate PROX estimates for all examinees and items, 
        even examinees who get all items right or all items wrong, and even for items 
        answered all correctly or all incorrectly.
    \param[in]  begin_examinee  Iterator to pointer to first examinees.
    \param[in]  end_examinee  Iterator to pointer to one past last examinees.
    \param[in]  begin_item  Iterator to pointers to first item.
    \param[in]  notPresentedResponse  Response which indicates an item was not 
        presented to an examinee. Assumed to be the same for all items.
    \param[out] itemdiff  Item difficulties.
    \param[out] abilities Examinee abilities.
   */
  template <class EI, class RI, class II> void PROX(int numItems, bool useAll, EI begin_examinee,
      EI end_examinee, II begin_item, Response notPresentedResponse, RealVector &itemdiff,
      RealVector &abilities)
  {
    const int maxIter = 10; // maximum number of iterations for PROX algorithm
    const double convergenceCrit = 0.01; // if maximum difference in examinee abilities is less then this
    // over consecutive iterations then stop
    int i, iitem;
    int num_examinees = end_examinee - begin_examinee;
    RealVector logitDiff(numItems);
    RealVector logitAbility(num_examinees);

    /* Compute logit of item difficulties */
    RealVector::iterator ilogit = logitDiff.begin();
    for (i=0; i<numItems; ++i, ++ilogit)
    {
      Real total = 0.0;
      Real sum = 0.0;
      EI ie = begin_examinee;
      int nvalid = 0;
      iitem = begin_item[i]->Index();
      Response correct = begin_item[i]->CorrectResponse();
      while (ie != end_examinee)
      {
        RI ir = (*ie)->responses_begin() + iitem;
        if (*ir != notPresentedResponse)
        {
          if (*ir == correct)
            sum += (*ie)->Count();
          total += (*ie)->Count();
          ++nvalid;
        }
        ++ie;
      }

      if (total == 0.0)
        throw RuntimeError("No valid responses for an item", "PROX");

      if (sum == 0.0 || sum == total)
      {
        if (useAll)
        {
          Real halfCount = total / (2*nvalid); // 1/2 average count per examinee (usually 0.5)
          Real maxLogit = std::log(halfCount / (total-halfCount));

          /* Use sum = halfcount if sum == 0.0 and use sum = total - halfCount if sum == total */
          *ilogit = (sum == 0.0) ? maxLogit : (1.0 / maxLogit);
        }
        else
        {
          *ilogit = invalidLogit;
        }
      }
      else
      {
        *ilogit = std::log(sum / (total-sum));
      }
    }

    /* Compute logit of examinee proportion correct scores */
    ilogit = logitAbility.begin();
    for (EI ie = begin_examinee; ie != end_examinee; ++ie, ++ilogit)
    {
      Real total = 0.0;
      Real sum = 0.0;
      II item = begin_item;
      RI bresp = (*ie)->responses_begin();
      for (i=numItems; i--; ++item)
      {
        Response resp = bresp[(*item)->Index()];
        if (resp != notPresentedResponse)
        {
          if (resp == (*item)->CorrectResponse())
            ++sum;
          ++total;
        }
      }

      if (total == 0.0)
      {
        *ilogit = invalidLogit;
      }
      else if (sum == 0.0 || sum == total)
      {
        if (useAll)
        {
          /* Use sum = 0.5 if sum == 0.0 and use sum = total - 0.5 if sum == total */
          *ilogit = (sum == 0.0) ? (std::log(0.5 / (total-0.5))) : std::log((total - 0.5) / 0.5);
        }
        else
        {
          *ilogit = invalidLogit;
        }
      }
      else
      {
        *ilogit = std::log(sum / (total-sum));
      }
    }

    int iter;
    double mean = 0.0;
    double variance = 1.0;
    RealVector prevAbility(abilities); // ability estimates at previous iteration - used to check for convergence
    for (iter=0; iter<maxIter; ++iter)
    {
      /* Compute difficulties for each item */
      RealVector::iterator idiff = itemdiff.begin();
      ilogit = logitDiff.begin();
      Real sum = 0.0;
      int nvalid = 0;
      for (i=0; i<numItems; ++i, ++ilogit, ++idiff)
      {
        if (*ilogit != invalidLogit)
        {
          /* Compute mean and variance of abilities for persons responding to item i */
          if (iter > 0) // mean and variance of item difficulties set to 0 and 1 for first iteration
          {
            iitem = begin_item[i]->Index();
            ItemMoments<EI, RI>(num_examinees, begin_examinee, iitem, abilities,
                notPresentedResponse, &mean, &variance);
          }

          /* Compute new item difficulty */
          *idiff = mean - std::sqrt(1 + variance/2.9) * *ilogit;
          sum += *idiff;

          ++nvalid;
        }
      }

      /* adjust item difficulties to sum to zero */
      sum /= nvalid;
      idiff = itemdiff.begin();
      for (i = numItems; i--; ++idiff)
      {
        if (*idiff != invalidLogit)
          *idiff -= sum;
      }

      /* Compute abilities for each examinee */
      RealVector::iterator iabil = abilities.begin();
      ilogit = logitAbility.begin();
      for (EI ie = begin_examinee; ie != end_examinee; ++ie, ++iabil, ++ilogit)
      {
        if (*ilogit != invalidLogit)
        {
          /* Compute mean and variance of item difficulties for items person responded to */
          PersonMoments<RI, II>((*ie)->responses_begin(), begin_item, itemdiff, notPresentedResponse, &mean, &variance);

          /* Compute new person ability */
          *iabil = mean + std::sqrt(1 + variance/2.9) * *ilogit;
        }
      }

      /* Check if abilities have converged */
      if (iter > 0)
      {
        Real maxdiff = 0.0;
        for (i=0; i<num_examinees; ++i)
        {
          if (abilities[i] != invalidLogit)
          {
            Real diff = std::fabs(abilities[i] - prevAbility[i]);
            if (diff > maxdiff)
              maxdiff = diff;
          }
        }
        if (maxdiff < convergenceCrit)
          break;
      }
      prevAbility = abilities;
    }
  }

  /*! 
    \brief
    Assigns individual examinee thetas into discrete latent variable
    categories and updates expected counts n and r for item M-step.
   
    Given category values of the discrete latent variable as 
    c_1, c_2, ..., c_ncat, classifies individual thetas into categories 
    using ncat-1 cutpoints (c_1+c_2)/2, (c_2+c_3)/2,..., where the 
    examinee is classified into the first category if the theta value 
    is less than (c_1+c_2)/2, into the second category if their theta 
    value is greater than or equal to (c_1+c_2)/2 and less than 
    (c_2+c_3)/2, etc.
    
    \section template_args Template Parameter
   
    \param I  Item type.
   
    \section function_args Function Parameters
   
    \param[in,out]  n Number of examinees who respond to the item on input, 0 on output.
    \param[in]  itemTheta Vector of thetas for examinees who reponded to an item.
    \param[in,out]  nVec  Iterator to vector of expected item response counts for each theta value?
    \param[in,out]  rVec  Iterator to vector of expected number correct counts for each theta value?.
    \param[in,out]  *item Pointer to item object.
   */
  template <class I> void AssignItemNR(int n, RealVector::iterator itemTheta,
      RealVector::iterator nVec, RealVector::iterator rVec, I *item)
  {
    int i;

    int ncat = item->NumLatentVarCat();

    if (ncat == 0)
      throw RuntimeError("Values of discrete latent variable categories are missing for an item",
          "AssignItemNR");

    /* Compute cutpoints for categorizing individual latent variable values */
    RealVector cutpoints(ncat-1);
    RealVector::iterator ic = cutpoints.begin();
    typename I::point_iterator ip = item->GetLatentVarPoints(1); // use points for first group
    for (i=ncat-1; i--; ++ic, ++ip)
    {
      *ic = (ip[0] + ip[1]) / 2.0;
    }

    item->InitializeNR();
    typename I::n_iterator in = item->NVector();
    typename I::r_iterator ir = item->RVector(item->CorrectResponse());
    for (; n--; ++itemTheta, ++nVec, ++rVec)
    {
      /* find category for theta using inefficient linear search */
      ic = cutpoints.begin();
      for (i=0; i<ncat-1; ++i, ++ic)
      {
        if (*itemTheta < *ic)
          break;
      }

      in[i] += *nVec;
      ir[i] += *rVec;
    }
  }

  /*!
    \brief
    Computes starting values for a sequence of items modeled with the
    three-parameter logistic model. 
    
    Nonlinear regression is employed for each item, treating the theta estimate for 
    each person as the independent variable and the item responses as the dependent
    variable. Theta estimates and initial b parameters for the nonlinear regression 
    are calculated using the PROX procedure. Theta estimates for each person are 
    categorized using latent variable categories defined for each item (using the
    interface of class ItemNR).
   
    Returns the number of times minimization procedure failed while using nonlinear 
    regression. Nonconvergence after the maximum number of iterations has been 
    exceeded is not considered a failure.
   
    \section template_args Template Parameters

    \param RI Type of iterator over examinee responses (iterator returned by 
        responses_begin() for an Examinee object).
    \param MI Type of iterator over minimization routines used for each item.
    \param EI Type of iterator over examinees (iterator over pointers to 
        Examinee objects).
    \param I  Item type.
    \param II Type of iterator over items (iterator over pointers to Item objects).
   
    \section function_args Function Parameters
   
    \param[in]  begin_min Iterator to pointers to minimization routines used for items.
         These should correspond to the items given by begin_item and end_item that 
         are in the set calc_items.
    \param[in]  begin_examinee  Iterator to pointer to the first examinee.
    \param[in]  end_examinee    Iterator to pointer to one past the last examinee. 
    \param[in,out]  begin_item  Iterator to pointer to the first item.
    \param[in,out]  end_item  Iterator to pointer to one past the last item.   
        All these items are used in computing PROX estimates of examinee 
        abilities. Starting values are computed for any of these items that are 
        also in the set calc_items. On output the item parameter estimates for each 
        item also present in calc_items are set to the starting values calculated 
        in this function. 
    \param[in]  notPresentedResponse  Response which indicates an item was not 
        presented to an examinee. Assumed to be the same for all items.
    \param[in]  useAll  If true calculate PROX estimates for all examinees and 
        items, even examinees who get all items right or all items wrong, and 
        even for items answered all correctly or all incorrectly.
    \param[in]  baseGroup Number of group to use as base group where mean and 
        s.d. of latent variable distribution are zero and one.
    \param[in]  calc_items  Set containing items for which starting values are 
        computed. All items given by begin_item and end_item are used to compute 
        PROX estimates of examinee ability, but starting values are only 
        calculated for items in calc_items. If calc_items is null then starting 
        values are computed for all items given by begin_item and end_item.
  */
  template <class RI, class MI, class EI, class I, class II> int StartingValues3PL(MI begin_min,
      EI begin_examinee, EI end_examinee, II begin_item, II end_item,
      Response notPresentedResponse, bool useAll = false, int baseGroup = 1,
      std::set<I *> *calc_items = 0)
  {

    int numExaminees = end_examinee - begin_examinee;
    int numItems = end_item - begin_item;

    RealVector thetas(numExaminees, invalidLogit);
    RealVector itemDiff(numItems, invalidLogit);

    typedef std::vector<Response> cvector;

    /* Compute PROX estimates of examinee ability and item difficulty */
    PROX<EI, RI, II>(numItems, useAll, begin_examinee, end_examinee, begin_item,
        notPresentedResponse, itemDiff, thetas);

    /* Compute mean of PROX estimates of examinee ability in base group and assign case weights */
    int i;
    RealVector nVec(numExaminees); // store case weights to use in regression
    RealVector::iterator in = nVec.begin();
    RealVector::iterator itheta = thetas.begin();
    Real mean = 0.0;
    Real nwt = 0.0;
    EI ie = begin_examinee;
    for (i = numExaminees; i--; ++in, ++itheta, ++ie)
    {
      *in = (*ie)->Count();
      if ((*ie)->Group() == baseGroup && *itheta != invalidLogit)
      {
        mean += *itheta * *in;
        nwt += *in;
      }
    }
    if (nwt == 0)
      throw RuntimeError("No examinees with valid item responses in base group",
          "StartingValues3PL");

    mean /= nwt;

    /* Compute standard deviation of PROX estimates of examinee ability in base group */
    Real sd = 0.0;
    itheta = thetas.begin();
    in = nVec.begin();
    ie = begin_examinee;
    for (i = numExaminees; i--; ++in, ++itheta, ++ie)
    {
      if ((*ie)->Group() == baseGroup && *itheta != invalidLogit)
      {
        Real diff = *itheta - mean;
        sd += diff * diff * *in;
      }
    }
    sd /= nwt;
    sd = std::sqrt(sd);

    /* Put thetas on latent variable scale so mean and s.d. of thetas for examinees in base group are 0 and 1 */
    itheta = thetas.begin();
    for (i=numExaminees; i--; ++itheta)
    {
      if (*itheta != invalidLogit)
      {
        *itheta -= mean;
        *itheta /= sd;
      }
    }

    /* Put b parameters on latent variable scale so mean and s.d. of thetas for examinees in base group are 0 and 1 */
    RealVector::iterator idiff = itemDiff.begin();
    for (i = numItems; i--; ++idiff)
    {
      if (*idiff != invalidLogit)
      {
        *idiff -= mean;
        *idiff /= sd;
      }
    }

    RealVector itemTheta(numExaminees); // vector of thetas for examinees who reponded to an item
    RealVector rVec(numExaminees); // count of number of correct answers for each theta value 
    // (r vector in M-step computation)
    // nVec holds count for each theta value (n vector in M-step computation)

    /* Compute starting values for each item */
    int minError = 0; // Count of number of time minimization failed in computing regression
    for (int item = 0; item < numItems; ++item, ++begin_item)
    {
      if (calc_items)
      {
        if (calc_items->find(*begin_item) == calc_items->end())
          continue;
      }

      /* Create vectors of item responses and thetas for item */
      EI ie = begin_examinee;
      RealVector::iterator it = itemTheta.begin();
      in = nVec.begin();
      RealVector::iterator ir = rVec.begin();
      RealVector::iterator alltheta = thetas.begin();
      int n = 0; // number of examinees who respond to the item
      nwt = 0.0; // weighted number of examinees who respond to the item
      Real pvalue = 0.0;
      Response correct = (*begin_item)->CorrectResponse();
      while (ie != end_examinee)
      {
        RI iresp = (*ie)->responses_begin() + (*begin_item)->Index();
        if (*iresp != notPresentedResponse && *alltheta != invalidLogit)
        {
          *it = *alltheta;
          *in = (*ie)->Count();
          nwt += *in;
          if (*iresp == correct)
          {
            pvalue += *in;
            *ir = *in;
          }
          else
          {
            *ir = 0.0;
          }
          ++n;
          ++it;
          ++in;
          ++ir;
        }
        ++alltheta;

        ++ie;
      }
      pvalue /= nwt;

      /* Compute starting values for item */
      int numParameters = (*begin_item)->NumParameters();
      RealVector start(numParameters);
      Real startb = itemDiff[item]; // use b from PROX as starting value

      if (numParameters > 1)
      {
        // PROX estimate assumes D=1.0, so if D for item is different from
        // 1.0 adjust the a-parameter.
        start(1) = 1.0 / (*begin_item)->NormalizingConstant();
        start(2) = startb;
      }
      else
      {
        start(1) = startb;
      }

      /* Make starting c parameter slightly larger than zero since prior density of c
       is zero when c is equal to zero */
      if (numParameters > 2)
        start(3) = 0.05;

      if (pvalue == 1.0 || pvalue == 0.0)
      {
        if (!useAll)
        {
          /* If all responses to an item are correct or incorrect and useAll is false then PROX
           estimate of b will not exist, so use logit of number correct score 1/2 different from zero
           or perfect as b starting value. */
          Real maxlogit = std::log((numExaminees - 0.5) / 0.5);
          startb = (pvalue == 1.0) ? maxlogit : 1.0/maxlogit;
          if (numParameters > 1)
            start(2) = startb;
          else
            start(1) = startb;
        }
        (*begin_item)->NonZeroPriors(start);
        (*begin_item)->SetParameters(start);
      }
      else
      {
        /* Use nonlinear regression to find starting values */

        /* Assign expected counts n and r for item */
        AssignItemNR<I>(n, itemTheta.begin(), nVec.begin(), rVec.begin(), *begin_item);

        (*begin_min)->SetFunction(*begin_item);

        RealVector xpls(numParameters), gpls(numParameters);
        double fpls;

        /* Modify starting values which for which prior density is zero */
        (*begin_item)->NonZeroPriors(start);

        int result = (*begin_min)->Minimize(start, xpls, fpls, gpls);

        // assign values computed if minimization succeeded 
        // or failed because of the maximum number of iterations
        // being exceeded without convergence
        if (!result || (*begin_min)->GetMessage() == 4)
        {
          (*begin_item)->SetParameters(xpls);
        }
        else // if minimization failed use initial starting value for b
        {
          start(1) = 1.0;
          start(3) = 0.1;
          (*begin_item)->SetParameters(start);
          minError++;
        }
      }
      ++begin_min; // only increment for items for which starting values are computed
    }
    return minError;
  }

#endif // ETIRM_START3PL_H_
} // namespace etirm

---