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;