PDF of Example | SAS/STAT Focus Area Examples

---

Overview

This SAS Web Example demonstrates how to fit one-, two-, three-, and four-parameter (1PL, 2PL, 3PL, and 4PL) unidimensional binary item response theory (IRT) models by using the MCMC procedure. The Analysis section briefly presents mathematical representations of the models. The Examples section presents one example for each of the four models. The discussion focuses on Bayesian model specification and PROC MCMC syntax. There is also sample SAS code for producing plots of item characteristic curves (ICC), item information curves (IIC), test characteristic curves (TCC), and test information curves (TIC).

Analysis

In unidimensional binary item response theory models, an instrument (test) consists of a number of items (questions) that require a binary response (for example, yes/no, true/false, or agree/disagree). The purpose of the instrument is to measure a single latent trait (or factor) of the subjects. The latent trait is assumed to be measurable and to have a range that encompasses the real line. An individual’s location within this range, , is assumed to be a continuous random variable. IRT models enable you to estimate the probability of a correct response for each item, to estimate the levels of the latent traits of the subjects, and to evaluate how well the items, individually and collectively, measure the latent trait.

In the models that follow, the probability that a subject will answer an item correctly (or affirmatively) is assumed to be a function of the subject’s location () on the range of possible values of the latent trait:

It is also assumed that a subject’s response to an item is independent of the subject’s response to any other item, conditional on . The function is known as an item response function (IRF) and is assumed to be a sigmoid, resembling a normal or logistic cumulative distribution. Choosing an IRT model amounts to fully specifying the functional form for . The models that are presented in this web example are all variations of a logistic function.

Note: The parameter estimates from a logistic model and a normal model that are fitted to the same data set differ by an approximately linear transformation; the parameter estimates from a normal model are approximately 1.7 times the estimates from a logistic model. You need to modify only one statement in a PROC MCMC model specification to switch between a normal model and a logistic IRT model.

One-Parameter Logistic (1PL) Model

In a 1PL model, the probability that a subject will respond correctly or affirmatively to an item is a function of and an item parameter, , which measures an item’s location on a latent difficulty scale. An item’s difficulty is defined as the level of that a random test subject must possess, such that the probability of a correct response is . The basic premise of the 1PL model is that the greater the distance between and , the greater the certainty in how a subject is expected to respond to an item. As the distance approaches zero, the more likely that there is a 50:50 chance that a subject will correctly respond to an item (De Ayala 2009). Mathematically, the 1PL model is written as

where i indexes subjects and j indexes items. The parameter is known as the item discrimination parameter. It is related to the IRF’s slope and reflects how well an item discriminates among individuals located at different points along the continuum (De Ayala 2009). In the 1PL model, is a singular parameter (it has no subscript). The fact that is common to all the items implies that all the IRFs are parallel. If you remove from the model, implying that it has a value of 1, you get what is known as the Rasch model.

The item information function for the 1PL model is

Two-Parameter Logistic (2PL) Model

The 2PL model has the following form:

It is identical to the 1PL model, except that now has a subscript, indicating that each item now has its own discrimination parameter. This means that the IRFs of the items are no longer constrained to be parallel. The cost of this flexibility is that if you have items, you now have to estimate additional parameters compared to the 1PL model.

The item information function for the 2PL model is

Three-Parameter Logistic (3PL) Model

In both the 1PL and 2PL models, the IRFs have a lower asymptote of 0, implying that as approaches , the probability of a correct or affirmative response approaches 0. However, if a subject randomly selects a response (guesses), the probability of a correct or affirmative response would be greater than 0 regardless of the subject’s or the difficulty of the item. The 3PL model includes an additional parameter, , for each item, enabling you to estimate the lower asymptote. The 3PL model has the following form:

For some items, it might be that guessing behavior is unlikely. In such cases, you can easily constrain selected parameters to be 0.

The item information function for the 3PL model is

Four-Parameter Logistic (4PL) Model

Although the 3PL model allows for the possibility that the lower asymptote is nonzero, the upper asymptote is still 1. This means that as approaches , the probability of a correct or affirmative response approaches 1. However, subjects with a very large can still make clerical errors. The 4PL model includes a ceiling parameter, , for each item, enabling you to estimate the upper asymptote. The 4PL model has the following form:

The item information function for the 4PL model is

ICC, TIC, IIC, and TIC Plots

After you choose a particular model and sample the posterior distributions of the model parameters, you can produce plots of versus for each item. Such plots are known as item characteristic curves. If you compute the sum of over all the items, you get what is known as the test characteristic function. A test characteristic function provides the expected number of correct or affirmative answers as a function of . A plot of the test characteristic function versus is known as a test characteristic curve. You can also assess the informational value of the items, individually and collectively. You do this by first computing the individual item information functions, which are defined as the negative of the expected value of the second derivative of the log-likelihood function. A plot of an item information function versus is known as an item information curve. If you then compute the sum of all the item information functions, you get what is known as the test information function. A plot of the test information function versus is known as a test information curve.

Examples

There are several different ways to specify IRT models in PROC MCMC. The method that the following examples use requires the data to be in long form (also known as hierarchical form, panel data form, or longitudinal form). When the data are in long form, the rows of the data set are indexed by both person and item, and there is a single response variable. Having the data in this form enables you to specify most of the model parameters as random effects by using the RANDOM statement in PROC MCMC. Treating the parameters as random effects is possible because of the conditional independence assumption stated in the Analysis section. This approach simplifies model specification in PROC MCMC, and experimentation shows that using this approach can reduce execution times significantly. In the following examples, the reduction in execution times ranges from 18% for the 4PL model to 37% for the 1PL model compared to alternative model specifications.

Preparing the Data

The following examples use the IrtBinary data set, which is used in the Getting Started example in the IRT procedure chapter in the SAS/STAT User's Guide. There are 50 subjects, and each subject responds to 10 items. These 10 items have binary responses: 1 indicates correct and 0 indicates incorrect. The following DATA step reads the data and creates the variable Person, which indexes the rows. The data set is in wide form.

data IrtBinary;
   input item1-item10 @@;
   person=_N_;
   datalines;
1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1
0 0 0 0 1 0 1 0 0 1 1 1 1 1 1 0 0 0 0 0
1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 1 0 1
1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 1 1

   ... more lines ...

0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
0 1 0 1 0 0 0 0 0 0 1 1 1 0 0 1 0 1 1 1
1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1
;

The following SAS statements use PROC TRANSPOSE to reshape IrtBinary into long form. PROC TRANSPOSE saves the output data set IrtBinaryRE. IrtBinaryRE contains three variables: Person, which indexes the subjects; Item, which indexes the items; and Response, which records whether the subjects’ responses to each item are correct. IrtBinaryRE contains 10 rows of data for each person; IrtBinary has only 1 row per person.

proc transpose data=IrtBinary out=IrtBinaryRE(rename=(col1=response _NAME_=item));
  by person;
run;

proc print data=IrtBinaryRE(obs=10);
run;

Output 1 shows the first 10 observations of the output data set IrtBinaryRE.

 

Output 1: First 10 Observations of IrtBinaryRE

 

Example 1: 1PL Model

Recall that the probability of a correct or affirmative response in the 1PL model is

To perform Bayesian estimation of this model, you must specify prior distributions for , , and . The prior for is usually specified as a standard normal distribution; is also often assigned a normal prior. However, unless you actually have prior information that you want to incorporate into the model, it is common practice to assign a large variance to the prior for so that it is an uninformative or diffuse prior. The theoretical range of encompasses the real line, so you might consider assigning a diffuse normal distribution as the prior for . However, because of the indeterminacy of metric that is inherent in IRT models (De Ayala 2009, p. 41), you sometimes get a negative estimate for , which would be interpreted to mean that subjects with lower s have a higher probability of a correct response than subjects with higher s. As a result, it is common practice for modelers to assign prior distributions such as a lognormal or a truncated normal distribution that restrict to be nonnegative.

You can also treat the parameters of the prior distributions for , , or as hyperparameters and assign hyperprior distributions. This example assigns a standard normal distribution to and a normal distribution with a large variance to . Both and are treated as random effects, and their prior distributions are specified by using RANDOM statements in PROC MCMC; is assigned a truncated normal prior with a mean of 1 (the value implied by the Rasch model), a large variance, and a lower truncation boundary of 0.

The following statements specify the 1PL model by using PROC MCMC:

ods graphics on;
proc mcmc data=IrtBinaryRE nmc=20000 outpost=out1p  seed=1000 nthreads=-1 ;
  parms alpha;
  prior alpha ~normal(1, var=10, lower=0);
  random theta~normal(0, var=1) subject=person namesuffix=position nooutpost;
  random delta~normal(0, var=10) subject=item monitor=(delta) namesuffix=position;
  p=logistic(alpha*(theta-delta));
  model response ~ binary(p);
run;

The OUTPOST= option in the PROC MCMC statement saves the MCMC samples in a data set named Out1p. The NTHREADS=–1 option sets the number of available threads to the number of hyperthreaded cores available on your system.

The PARMS statement declares the parameter , and the PRIOR statement specifies the prior distribution for as a truncated normal distribution with a mean of 1, a variance of 10, and a lower truncation boundary of 0.

The first RANDOM statement specifies the prior distribution for as a standard normal distribution. The SUBJECT= option specifies that the variable Person identifies the subjects. The NOOUTPOST option suppresses the output of the posterior samples of the random-effects parameters to the OUTPOST= data set; this reduces the execution time. However, if you want to perform analysis on the posterior samples of , you can omit this option. The NAMESUFFIX= option specifies how the names of the random-effects parameters are internally created from the SUBJECT= variable; when you specify NAMESUFFIX=POSITION, the MCMC procedure constructs the parameter names by appending the numbers 1, 2, 3, and so on, where the number indicates the order in which the SUBJECT= variable appears in the data set.

The second RANDOM statement specifies the prior distribution for as a normal distribution with a mean of 0 and a variance of 10. The SUBJECT= option specifies that the variable Item identifies the subjects. The MONITOR= option requests that the random-effects parameters for be monitored and reported in the output.

The following programming statement generates the variable P, which contains the probability of a correct response. If you prefer a normal model instead of a logistic model, just replace the LOGISTIC function with the PROBNORM function. That is, specify the following:

p=probnorm(alpha*(theta-delta));

The MODEL statement specifies the conditional distribution of the data, given the parameters (the likelihood function). In this example, it specifies that Response has a binary (Bernoulli) distribution with a probability of success equal to P.

Output 2 shows the posterior summaries and intervals table for the 1PL model. Difficulty parameters measure the difficulties of the items. As the value of the difficulty parameter increases, the item becomes more difficult. You can observe that all the difficulty parameters are less than 0, which suggests that all the items in this example are relatively easy.

 

Output 2: Posterior Summaries and Intervals

 

 

Producing ICC, IIC, TCC, and TIC Plots

To produce ICC, IIC, TCC, and TIC plots, you use the means of the posterior samples that are saved in the OUTPOST= data set Out1p to compute the following information: < /p>

• the probability of a correct or affirmative response for each item over a range of values of (ICC)

• the sum of the probabilities across all items (TCC)

• the item information functions for each item over a range of values of (IIC)

• the sum of the item information functions across all items (TIC)

The following statements create the data set Plots1PL, from which you can produce ICC, TCC, IIC, and TIC plots by using PROC SGPLOT:

data plots;
  set out1p;
  array d[10] delta_1-delta_10;
  array p[10] p1-p10;
  array i[10] i1-i10;
  do theta=-6 to 6 by .25;
    do j=1 to 10;
      p[j]=logistic(alpha*(theta-d[j]));
	  i[j]=alpha**2*p[j]*(1-p[j]);
	end;
	tcf=sum(of p1-p10);
	tif=sum(of i1-i10);
	output;
  end;
run;

proc sort data=plots;
  by theta;
run;

proc means data=plots noprint;
  var p1-p10 i1-i10 tcf tif;
  by theta;
  output out=plots1pl
    mean=
    p5=lb_p1-lb_p10 lb_i1-lb_i10 lb_tcf lb_tif
    p95=ub_p1-ub_p10 ub_i1-ub_i10 ub_tcf ub_tif;
run;

The following statements use the data set Plots1PL and PROC SGPLOT to produce a plot of the item characteristic curve and a 95% credible interval for the first item:

proc sgplot data=plots1pl;
  title "Item Characteristic Curve";
  title2 "item=1";
  band x=theta upper=ub_p1 lower=lb_p1 / modelname="icc" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=p1 / name="icc";
  refline .5 / axis=y;
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Probability";
run;
title;

Figure 1 displays a plot of the ICC and a 95% credible interval from the 1PL model for the first item. The curve shows the how the probability of a correct answer varies as the value of the latent trait varies. The difficult parameter for the first item is –1.1553; the ICC curve intersects the horizontal line at this point. This indicates that a person with a below average trait value of –1.1553 has a 50% chance of answering the question correctly.

 

Figure 1: Item Characteristic Curve

 

The following statements produce a plot of the test characteristic curve that shows you the expected number of correct responses as a function of :

proc sgplot data=plots1pl;
  title "Test Characteristic Curve";
  band x=theta upper=ub_tcf lower=lb_tcf / modelname="tcc" transparency=.5
    legendlabel="95% Credible Interval";
  series x=theta y=tcf / name="tcc";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Expected Number of Correct Answers";
run;
title;

Figure 2 displays a plot of the TCC and a 95% credible interval from the 1PL model. The TCC curve shows how the expected number of correct answers varies as the value of the latent trait varies. The graph indicates that a person with an average trait value of 0 is expected to answer approximately eight of the ten questions correctly.

 

Figure 2: Test Characteristic Curve

 

The following statements produce a plot of the item information curve for the first item from the 1PL model:

proc sgplot data=plots1pl;
  title "Item Information Curve";
  title2 "item=1";
  band x=theta upper=ub_i1 lower=lb_i1 / modelname="iic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=i1 / name="iic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;

Figure 3 displays a plot of the IIC and a 95% credible interval of the first item from the 1PL model. The curve shows that an item provides its maximum information at the point where the latent trait value equals the item’s difficulty.

 

Figure 3: Item Information Curve

 

The following statements produce a plot of the test information curve from the 1PL model:

proc sgplot data=plots1pl;
  title "Test Information Curve";
  band x=theta upper=ub_tif lower=lb_tif / modelname="tic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=tif / name="tic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;

Figure 4 displays a plot of the TIC and a 95% credible interval from the 1PL model. The curve shows that the instrument provides its maximum information for estimating the latent trait at a value of approximately –1.

 

Figure 4: Test Information Curve

 

Example 2: 2PL Model

The 2PL model extends the 1PL model by relaxing the assumption that the IRFs are parallel. This is accomplished by introducing the discrimination parameter for each item. The probability of a correct or affirmative response is

The following statements specify the 2PL model by using PROC MCMC:

proc mcmc data=IrtBinaryRE nmc=20000 outpost=out2pl  seed=1000 nthreads=-1;
  random theta~normal(0, var=1) subject=person namesuffix=position nooutpost;
  random delta~normal(0, var=10) subject=item monitor=(delta) namesuffix=position;
  random alpha~normal(1, var=10) subject=item monitor=(alpha) namesuffix=position;
  p=logistic(alpha*(theta-delta));
  model response ~ binary(p);
run;

The only change from the 1PL model is that is now treated as a random effect, so you no longer have a PARMS statement and a PRIOR statement for ; instead, you now have a RANDOM statement. In the new RANDOM statement, the prior distribution for is specified as normal with a mean of 1 and a variance of 10. The SUBJECT= option specifies that the variable Item identifies the subjects, and the MONITOR= option requests that the random-effects parameters for be monitored and reported in the output. The statement that generates the variable P is changed to reflect the 2PL model’s probability of a correct response.

Output 3 shows the posterior summaries and intervals table for the 2PL model.

 

Output 3: Posterior Summaries and Intervals

 

 

The estimates for are quite variable, indicating that the parallel IRF assumption of the 1PL model is likely to be too restrictive.

The SAS program that produces the ICC, IIC, TCC, and TIC plots is essentially the same as for the 1PL model except for slight modifications to the code that creates the Plots2PL data set. The modifications are needed to accommodate the new parameters, and the probability and information function formulas are updated for the 2PL model. Only the DATA step that creates the data set Plots2PL is shown here; everything else is unchanged and is not shown to save space. The complete SAS program is included (for all the models) in the downloadable sample SAS program file that accompanies this document.

data plots;
  set out2pl;
  array alpha[10] alpha_1-alpha_10;
  array d[10] delta_1-delta_10;
  array p[10] p1-p10;
  array i[10] i1-i10;
  do theta=-6 to 6 by .25;
    do j=1 to 10;
      p[j]=logistic(alpha[j]*(theta-d[j]));
      i[j]=alpha[j]**2*p[j]*(1-p[j]);
    end;
    tcf=sum(of p1-p10);
    tif=sum(of i1-i10);
    output;
  end;
run;

Example 3: 3PL Model

The 3PL model extends the 2PL model by adding the parameters , which are used to estimate the lower asymptote of the IRF. This addresses the effect that guessing has on estimating the probability of a correct or affirmative response. The probability of a correct response for the 3PL model is

The following statements specify the 3PL model by using PROC MCMC:

proc mcmc data=IrtBinaryRE nmc=20000 outpost=out3pl  seed=1000 nthreads=-1;
  random theta~normal(0, var=1) subject=person namesuffix=position nooutpost;
  random delta~normal(0, var=10) subject=item monitor=(delta) namesuffix=position;
  random alpha~normal(1, var=10) subject=item monitor=(alpha) namesuffix=position;
  random g~beta(2,2)  subject=item monitor=(g) namesuffix=position;
  p=g + (1 - g)*logistic(alpha*(theta-delta));
  model response ~ binary(p);
run;

There are only two changes from the 2PL model. The first is a new RANDOM statement for the parameters. The prior distribution for is specified as a beta distribution with both shape parameters equal to 2. The SUBJECT= option specifies that the variable Item identifies the subjects, and the MONITOR= option requests that the random-effects parameters for be monitored and reported in the output. The second change is in the statement that generates the variable P. It is changed to reflect the 3PL model’s probability of a correct response.

Output 4 shows the posterior summaries and intervals table for the 3PL model.

 

Output 4: Posterior Summaries and Intervals

 

 

The strictly positive 95% HPD intervals for the parameters indicate that the lower asymptote is not likely to be 0, as implied by the 1PL and 2PL models. The inclusion of in the model has a significant effect on the estimates for the and parameters.

To create the ICC, IIC, TCC, and TIC plots, you must modify the DATA step that creates the Plots3PL data set to accommodate the new parameters and to update the probability and information function formulas for the 3PL model. The modified statements are as follows:

data plots;
  set out3pl;
  array alpha[10] alpha_1-alpha_10;
  array d[10] delta_1-delta_10;
  array p[10] p1-p10;
  array i[10] i1-i10;
  array g[10] g_1-g_10;
  do theta=-6 to 6 by .25;
    do j=1 to 10;
      p[j]=g[j] + (1 - g[j])*logistic(alpha[j]*(theta-d[j]));
      i[j]=alpha[j]**2*((p[j]-g[j])**2/(1-g[j])**2)*((1-p[j])/p[j]);
    end;
    tcf=sum(of p1-p10);
    tif=sum(of i1-i10);
    output;
  end;
run;

Example 4: Four-Parameter Logistic IRT Model

The 4PL model extends the 3PL model by adding the parameters , which are used to estimate the upper asymptote of the IRF. This addresses the effect that clerical errors have on estimating the probability of a correct or affirmative response. The probability of a correct response for the 4PL model is

The following statements specify the 4PL model by using PROC MCMC:

proc mcmc data=IrtBinaryRE nmc=20000 outpost=out4pl  seed=1000 nthreads=-1 ;
  random theta~normal(0, var=1) subject=person namesuffix=position nooutpost;
  random delta~normal(0, var=10) subject=item monitor=(delta) namesuffix=position;
  random alpha~normal(1, var=10) subject=item monitor=(alpha) namesuffix=position;
  random g~beta(2,2)  subject=item monitor=(g) namesuffix=position;
  llc=lpdfbeta(c,2,2,g);
  random c~general(llc)  subject=item monitor=(c) initial=.999 namesuffix=position;
  p=g + (c - g)*logistic(alpha*(theta-delta));
  model response ~ binary(p);
run;

There are three changes from the 3PL model. Like the parameters of the 3PL model, the parameters are assigned a beta prior distribution. However, you need to restrict the parameters, which are the upper asymptotes of the IRF, to be greater than the corresponding parameters, which are the lower asymptotes of the IRF. You can accomplish this by using a truncated beta distribution and specifying as the lower boundary for the distribution of . However, the RANDOM statement does not directly support the truncated beta distribution. Therefore, you must use an indirect method that is documented in the section "Using Density Functions in the Programming Statements" of the chapter "The MCMC Procedure" in the SAS/STAT User's Guide. The LPDFBETA function computes the density of a beta distribution, and it permits both upper and lower truncation boundaries. So you use a programming statement to create a variable Llc and set it equal to the LPDFBETA function with appropriate parameters as follows:

llc=lpdfbeta(c,2,2,g);

You then specify a RANDOM statement for , specify the prior distribution as a general distribution, and specify the variable Llc as the general distribution’s argument as follows:

random c~general(llc)  subject=item monitor=(c) initial=.999 namesuffix=position;

When you use the general distribution in a RANDOM statement, PROC MCMC requires you to specify the INITIAL= option to provide initial values for the random-effects parameters.

Finally, you need to update the programming statement that specifies the formula for P to reflect the appropriate probability for the 4PL model.

Output 5 shows the posterior summaries and intervals table for the 4PL model.

 

Output 5: Posterior Summaries and Intervals

 

 

The estimated upper asymptotes are consistently less than 1. The inclusion of the parameters also significantly affects the estimates of the and parameters.

To create the ICC, IIC, TCC, and TIC plots, you must modify the DATA step that creates the Plots4PL data set to accommodate the new parameters and to update the probability and information function formulas for the 4PL model. The modified DATA step is as follows:

data plots;
  set out4pl;
  array alpha[10] alpha_1-alpha_10;
  array d[10] delta_1-delta_10;
  array p[10] p1-p10;
  array i[10] i1-i10;
  array g[10] g_1-g_10;
  array c[10] c_1-c_10;
  do theta=-6 to 6 by .25;
    do j=1 to 10;
      p[j]=g[j] + (c[j] - g[j])*logistic(alpha[j]*(theta-d[j]));
      i[j]=(alpha[j]**2*(p[j]-c[j])**2*(g[j]-p[j])**2)/((g[j]-c[j])**2*p[j]*(1-p[j]));
    end;
    tcf=sum(of p1-p10);
    tif=sum(of i1-i10);
    output;
  end;
run;

References

• De Ayala, R.J. (2009). The Theory and Practice of Item Response Theory. New York: Guilford Press.

---

These sample files and code examples are provided by SAS Institute Inc. "as is" without warranty of any kind, either express or implied, including but not limited to the implied warranties of merchantability and fitness for a particular purpose. Recipients acknowledge and agree that SAS Institute shall not be liable for any damages whatsoever arising out of their use of this material. In addition, SAS Institute will provide no support for the materials contained herein.

---

data IrtBinary;
   input item1-item10 @@;
   person=_N_;
   datalines;
1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1
0 0 0 0 1 0 1 0 0 1 1 1 1 1 1 0 0 0 0 0
1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 1 0 1
1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 1 1
1 1 1 1 1 0 1 1 0 1 0 0 0 0 0 0 0 1 0 0
1 1 1 1 0 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0
1 1 0 1 1 1 1 1 1 1 1 0 0 0 1 0 1 0 0 1
0 1 1 1 0 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1
0 1 1 0 0 0 1 0 0 1 0 1 1 1 0 0 1 0 0 0
1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 0
0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 0 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1
1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 0 1
0 1 0 0 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1
1 0 0 1 0 1 0 0 0 0 1 1 1 1 1 1 1 0 1 1
1 1 1 1 1 0 1 1 0 1 1 1 1 1 0 1 1 1 0 1
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
1 1 0 1 1 1 1 0 1 0 0 0 1 1 0 1 0 0 0 0
1 0 1 1 1 0 1 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 0 0
1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 0 1
1 1 1 1 1 1 0 0 1 0 1 1 1 1 1 0 0 1 1 1
1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 0 1 1 1 1
1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 0 0 1 1 1 1 1 1 0 1 0 0 0 0 1
1 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 1 1 0 0
0 1 1 1 1 0 0 1 0 1 1 1 1 0 1 1 1 1 0 0
1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0
0 0 0 0 1 0 0 0 0 1 1 1 1 1 0 1 1 1 1 1
1 1 0 0 0 1 0 0 1 1 1 1 1 1 0 0 0 1 1 0
0 0 1 0 1 0 0 0 1 0 1 1 0 0 1 0 0 0 1 1
1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 1 1 0
1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 0 0 1 1
0 0 0 1 1 0 0 0 1 0 1 1 1 1 0 1 1 1 0 1
0 0 1 1 1 0 0 1 1 0 0 1 1 0 0 1 1 1 1 0
1 1 1 0 0 0 1 0 1 0 0 1 0 1 1 0 1 0 0 0
1 1 1 0 1 1 1 0 1 1 0 0 0 0 0 1 1 1 0 0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 0
1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 1 1 0 1
1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1
1 1 1 0 1 1 0 1 0 1 1 1 1 1 1 1 1 0 1 1
0 1 0 0 1 1 0 0 0 0 1 1 1 1 1 1 1 0 0 0
0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
0 1 0 1 0 0 0 0 0 0 1 1 1 0 0 1 0 1 1 1
1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1
;


proc transpose data=IrtBinary out=IrtBinaryRE(rename=(col1=response _NAME_=item));
  by person;
run;


proc print data=IrtBinaryRE(obs=10);
run;


ods graphics on;
proc mcmc data=IrtBinaryRE nmc=20000 outpost=out1p  seed=1000 nthreads=-1 ;
  parms alpha;
  prior alpha ~normal(1, var=10, lower=0);
  random theta~normal(0, var=1) subject=person namesuffix=position nooutpost;
  random delta~normal(0, var=10) subject=item monitor=(delta) namesuffix=position;
  p=logistic(alpha*(theta-delta));
  model response ~ binary(p);
run;


data plots;
  set out1p;
  array d[10] delta_1-delta_10;
  array p[10] p1-p10;
  array i[10] i1-i10;
  do theta=-6 to 6 by .25;
    do j=1 to 10;
      p[j]=logistic(alpha*(theta-d[j]));
	  i[j]=alpha**2*p[j]*(1-p[j]);
	end;
	tcf=sum(of p1-p10);
	tif=sum(of i1-i10);
	output;
  end;
run;


proc sort data=plots;
  by theta;
run;

proc means data=plots noprint;
  var p1-p10 i1-i10 tcf tif;
  by theta;
  output out=plots1pl
    mean=
    p5=lb_p1-lb_p10 lb_i1-lb_i10 lb_tcf lb_tif
    p95=ub_p1-ub_p10 ub_i1-ub_i10 ub_tcf ub_tif;
run;


proc sgplot data=plots1pl;
  title "Item Characteristic Curve";
  title2 "item=1";
  band x=theta upper=ub_p1 lower=lb_p1 / modelname="icc" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=p1 / name="icc";
  refline .5 / axis=y;
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Probability";
run;
title;


proc sgplot data=plots1pl;
  title "Test Characteristic Curve";
  band x=theta upper=ub_tcf lower=lb_tcf / modelname="tcc" transparency=.5
    legendlabel="95% Credible Interval";
  series x=theta y=tcf / name="tcc";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Expected Number of Correct Answers";
run;
title;


proc sgplot data=plots1pl;
  title "Item Information Curve";
  title2 "item=1";
  band x=theta upper=ub_i1 lower=lb_i1 / modelname="iic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=i1 / name="iic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;


proc sgplot data=plots1pl;
  title "Test Information Curve";
  band x=theta upper=ub_tif lower=lb_tif / modelname="tic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=tif / name="tic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;


proc mcmc data=IrtBinaryRE nmc=20000 outpost=out2pl  seed=1000 nthreads=-1;
  random theta~normal(0, var=1) subject=person namesuffix=position nooutpost;
  random delta~normal(0, var=10) subject=item monitor=(delta) namesuffix=position;
  random alpha~normal(1, var=10) subject=item monitor=(alpha) namesuffix=position;
  p=logistic(alpha*(theta-delta));
  model response ~ binary(p);
run;


data plots;
  set out2pl;
  array alpha[10] alpha_1-alpha_10;
  array d[10] delta_1-delta_10;
  array p[10] p1-p10;
  array i[10] i1-i10;
  do theta=-6 to 6 by .25;
    do j=1 to 10;
      p[j]=logistic(alpha[j]*(theta-d[j]));
      i[j]=alpha[j]**2*p[j]*(1-p[j]);
    end;
    tcf=sum(of p1-p10);
    tif=sum(of i1-i10);
    output;
  end;
run;


proc sort data=plots;
  by theta;
run;

proc means data=plots noprint;
  var p1-p10 i1-i10 tcf tif;
  by theta;
  output out=plots2pl
    mean=
    p5=lb_p1-lb_p10 lb_i1-lb_i10 lb_tcf lb_tif
    p95=ub_p1-ub_p10 ub_i1-ub_i10 ub_tcf ub_tif;
run;

proc sgplot data=plots2pl;
  title "Item Characteristic Curve";
  title2 "item=1";
  band x=theta upper=ub_p1 lower=lb_p1 / modelname="icc" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=p1 / name="icc";
  refline .5 / axis=y;
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Probability";
run;
title;

proc sgplot data=plots2pl;
  title "Test Characteristic Curve";
  band x=theta upper=ub_tcf lower=lb_tcf / modelname="tcc" transparency=.5
    legendlabel="95% Credible Interval";
  series x=theta y=tcf / name="tcc";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Expected Number of Correct Answers";
run;
title;

proc sgplot data=plots2pl;
  title "Item Information Curve";
  title2 "item=1";
  band x=theta upper=ub_i1 lower=lb_i1 / modelname="iic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=i1 / name="iic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;

proc sgplot data=plots2pl;
  title "Test Information Curve";
  band x=theta upper=ub_tif lower=lb_tif / modelname="tic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=tif / name="tic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;


proc mcmc data=IrtBinaryRE nmc=20000 outpost=out3pl  seed=1000 nthreads=-1;
  random theta~normal(0, var=1) subject=person namesuffix=position nooutpost;
  random delta~normal(0, var=10) subject=item monitor=(delta) namesuffix=position;
  random alpha~normal(1, var=10) subject=item monitor=(alpha) namesuffix=position;
  random g~beta(2,2)  subject=item monitor=(g) namesuffix=position;
  p=g + (1 - g)*logistic(alpha*(theta-delta));
  model response ~ binary(p);
run;


data plots;
  set out3pl;
  array alpha[10] alpha_1-alpha_10;
  array d[10] delta_1-delta_10;
  array p[10] p1-p10;
  array i[10] i1-i10;
  array g[10] g_1-g_10;
  do theta=-6 to 6 by .25;
    do j=1 to 10;
      p[j]=g[j] + (1 - g[j])*logistic(alpha[j]*(theta-d[j]));
      i[j]=alpha[j]**2*((p[j]-g[j])**2/(1-g[j])**2)*((1-p[j])/p[j]);
    end;
    tcf=sum(of p1-p10);
    tif=sum(of i1-i10);
    output;
  end;
run;


proc sort data=plots;
  by theta;
run;

proc means data=plots noprint;
  var p1-p10 i1-i10 tcf tif;
  by theta;
  output out=plots3pl
    mean=
    p5=lb_p1-lb_p10 lb_i1-lb_i10 lb_tcf lb_tif
    p95=ub_p1-ub_p10 ub_i1-ub_i10 ub_tcf ub_tif;
run;


proc sgplot data=plots3pl;
  title "Item Characteristic Curve";
  title2 "item=1";
  band x=theta upper=ub_p1 lower=lb_p1 / modelname="icc" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=p1 / name="icc";
  refline .5 / axis=y;
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Probability";
run;
title;

proc sgplot data=plots3pl;
  title "Test Characteristic Curve";
  band x=theta upper=ub_tcf lower=lb_tcf / modelname="tcc" transparency=.5
    legendlabel="95% Credible Interval";
  series x=theta y=tcf / name="tcc";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Expected Number of Correct Answers";
run;
title;

proc sgplot data=plots3pl;
  title "Item Information Curve";
  title2 "item=1";
  band x=theta upper=ub_i1 lower=lb_i1 / modelname="iic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=i1 / name="iic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;

proc sgplot data=plots3pl;
  title "Test Information Curve";
  band x=theta upper=ub_tif lower=lb_tif / modelname="tic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=tif / name="tic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;


proc mcmc data=IrtBinaryRE nmc=20000 outpost=out4pl  seed=1000 nthreads=-1 ;
  random theta~normal(0, var=1) subject=person namesuffix=position nooutpost;
  random delta~normal(0, var=10) subject=item monitor=(delta) namesuffix=position;
  random alpha~normal(1, var=10) subject=item monitor=(alpha) namesuffix=position;
  random g~beta(2,2)  subject=item monitor=(g) namesuffix=position;
  llc=lpdfbeta(c,2,2,g);
  random c~general(llc)  subject=item monitor=(c) initial=.999 namesuffix=position;
  p=g + (c - g)*logistic(alpha*(theta-delta));
  model response ~ binary(p);
run;


data plots;
  set out4pl;
  array alpha[10] alpha_1-alpha_10;
  array d[10] delta_1-delta_10;
  array p[10] p1-p10;
  array i[10] i1-i10;
  array g[10] g_1-g_10;
  array c[10] c_1-c_10;
  do theta=-6 to 6 by .25;
    do j=1 to 10;
      p[j]=g[j] + (c[j] - g[j])*logistic(alpha[j]*(theta-d[j]));
      i[j]=(alpha[j]**2*(p[j]-c[j])**2*(g[j]-p[j])**2)/((g[j]-c[j])**2*p[j]*(1-p[j]));
    end;
    tcf=sum(of p1-p10);
    tif=sum(of i1-i10);
    output;
  end;
run;


proc sort data=plots;
  by theta;
run;

proc means data=plots noprint;
  var p1-p10 i1-i10 tcf tif;
  by theta;
  output out=plots4pl
    mean=
    p5=lb_p1-lb_p10 lb_i1-lb_i10 lb_tcf lb_tif
    p95=ub_p1-ub_p10 ub_i1-ub_i10 ub_tcf ub_tif;
run;


proc sgplot data=plots4pl;
  title "Item Characteristic Curve";
  title2 "item=1";
  band x=theta upper=ub_p1 lower=lb_p1 / modelname="icc" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=p1 / name="icc";
  refline .5 / axis=y;
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Probability";
run;
title;

proc sgplot data=plots4pl;
  title "Test Characteristic Curve";
  band x=theta upper=ub_tcf lower=lb_tcf / modelname="tcc" transparency=.5
    legendlabel="95% Credible Interval";
  series x=theta y=tcf / name="tcc";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Expected Number of Correct Answers";
run;
title;

proc sgplot data=plots4pl;
  title "Item Information Curve";
  title2 "item=1";
  band x=theta upper=ub_i1 lower=lb_i1 / modelname="iic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=i1 / name="iic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;

proc sgplot data=plots4pl;
  title "Test Information Curve";
  band x=theta upper=ub_tif lower=lb_tif / modelname="tic" transparency=.5
       legendlabel="95% Credible Interval";
  series x=theta y=tif / name="tic";
  xaxis label="Trait ((*ESC*){unicode theta})";
  yaxis label="Information";
run;
title;

---

These sample files and code examples are provided by SAS Institute Inc. "as is" without warranty of any kind, either express or implied, including but not limited to the implied warranties of merchantability and fitness for a particular purpose. Recipients acknowledge and agree that SAS Institute shall not be liable for any damages whatsoever arising out of their use of this material. In addition, SAS Institute will provide no support for the materials contained herein.

This sample demonstrates how to fit one-, two-, three-, and four-parameter (1PL, 2PL, 3PL, and 4PL) unidimensional binary item response theory (IRT) models by using the MCMC procedure.

Type:	Sample
Topic:	Analytics ==> Bayesian Analysis SAS Reference ==> Procedures ==> MCMC

Date Modified:	2017-01-12 16:08:34
Date Created:	2016-10-25 15:18:00

Operating System and Release Information

Product Family	Product	Host	SAS Release
			Starting	Ending
SAS System	SAS/STAT	z/OS
		z/OS 64-bit
		OpenVMS VAX
		Microsoft® Windows® for 64-Bit Itanium-based Systems
		Microsoft Windows Server 2003 Datacenter 64-bit Edition
		Microsoft Windows Server 2003 Enterprise 64-bit Edition
		Microsoft Windows XP 64-bit Edition
		Microsoft® Windows® for x64
		OS/2
		Microsoft Windows 8 Enterprise 32-bit
		Microsoft Windows 8 Enterprise x64
		Microsoft Windows 8 Pro 32-bit
		Microsoft Windows 8 Pro x64
		Microsoft Windows 8.1 Enterprise 32-bit
		Microsoft Windows 8.1 Enterprise x64
		Microsoft Windows 8.1 Pro 32-bit
		Microsoft Windows 8.1 Pro x64
		Microsoft Windows 10
		Microsoft Windows 95/98
		Microsoft Windows 2000 Advanced Server
		Microsoft Windows 2000 Datacenter Server
		Microsoft Windows 2000 Server
		Microsoft Windows 2000 Professional
		Microsoft Windows NT Workstation
		Microsoft Windows Server 2003 Datacenter Edition
		Microsoft Windows Server 2003 Enterprise Edition
		Microsoft Windows Server 2003 Standard Edition
		Microsoft Windows Server 2003 for x64
		Microsoft Windows Server 2008
		Microsoft Windows Server 2008 R2
		Microsoft Windows Server 2008 for x64
		Microsoft Windows Server 2012 Datacenter
		Microsoft Windows Server 2012 R2 Datacenter
		Microsoft Windows Server 2012 R2 Std
		Microsoft Windows Server 2012 Std
		Microsoft Windows XP Professional
		Windows 7 Enterprise 32 bit
		Windows 7 Enterprise x64
		Windows 7 Home Premium 32 bit
		Windows 7 Home Premium x64
		Windows 7 Professional 32 bit
		Windows 7 Professional x64
		Windows 7 Ultimate 32 bit
		Windows 7 Ultimate x64
		Windows Millennium Edition (Me)
		Windows Vista
		Windows Vista for x64
		64-bit Enabled AIX
		64-bit Enabled HP-UX
		64-bit Enabled Solaris
		ABI+ for Intel Architecture
		AIX
		HP-UX
		HP-UX IPF
		IRIX
		Linux
		Linux for x64
		Linux on Itanium
		OpenVMS Alpha
		OpenVMS on HP Integrity
		Solaris
		Solaris for x64
		Tru64 UNIX