R, like the commercial S-Plus, is based on the programming language S. It can be downloaded for from www.r-project.org [At home, after reading the instructions!, choose: CRAN, then Windows (95 and later), then base, then R-2.4.0-win32.exe]. This page also has links to FAQs and other information about R.
One of the nice things about R (besides the fact that it's powerful and free) is that you can save everything you've done after you've finished each time. If you are using R on a lab PC in LeConte, you will need to manually load and save the .Rdata file in your Z drive each time you want to start and end the program, respectively. (This is because you don't have access to the C drive.) This wouldn't be a problem if you are using it on your own PC.
At the heart of R are the various objects that you enter. At any time (when you are in R) you can type in the function objects() to see which ones you have entered previously. Presumably there aren't any there right now. R also has many built in objects, most of which are functions. Because these objects are always there, it will not list them when you type objects().
A list of the various functions for R can be found in the manuals that are installed with the program. The various manuals can be found under the Manuals submenu of the Help menu. In particular, the R Reference Manual lists not only all the functions in the base package of R, but also some of those in the supplementary packages as well. Because this manual is nearly 1000 pages long, you probably shouldn't print out a copy!!!!
Once you know the name of a function, you can get help about it simply by typing help(functionname). In R, parentheses indicate either mathematical grouping, like they usually do, or that you are applying a function. Braces [ ] indicate that you are finding an element in a string or array.
This brief example shows some of R's flexibility.
x<-c(5,4,3,2)
The <- is the command that assigns what ever is on the left to the name that is on the right. The c indicates that an array or vector is going to be entered. Simply typing x now will return those four values. Typing x[3] will return the third value, mean(x) will return the mean, and objects() will show you that you now have an object called x.
R also has a variety of graphical functions and statistical functions built in. hist(x) will construct a histogram, and t.test(x,alternative="greater",mu=5) will test the null hypothesis mu=5 vs. the alternate that mu > 5. (If you happen to mistype a command and get an error, just hit the up arrow to bring it up again for editing).
Besides just getting the result of a function, we can store the results somewhere. For example x.t<-t.test(x,alternative="greater",mu=5) will store the result of the t-test in x.t. attributes(x.t) will then show you the various properties of x.t. Thus, x.t$p.value will return the p-value.
R can also be used to write your own functions. For example, if you wanted a function that returned sin(x)+x2 we could write:
weird<-function(x){
y<-sin(x)+x^2
return(y)}
Try weird(5) and weird(x). What is it doing?
Loading in Data: Two common formats for testing data are simply big matrices, either with or without spaces. Two examples can be found at: http://www.stat.sc.edu/~habing/courses/data/mdatab.txt and http://www.stat.sc.edu/~habing/courses/data/mdatabspace.txt. To load these into R as an object called math we could do:
math<-read.fwf("http://www.stat.sc.edu/~habing/courses/data/mdatab.txt",widths=rep(1,32))
or
math<-read.fwf("http://www.stat.sc.edu/~habing/courses/data/mdatabspace.txt",widths=rep(2,32))
The widths=rep(x,32) tells r that there will be 32 columns, each of width x. Notice that just typing math won't be very interesting. However, if you type
dim(math)
you will see that math is a matrix with 2642 rows and 32 columns. Using
apply(math,2,mean)
will give you the average of each column (the 2 is for column, 1 would be for row). Finally, we could get a histogram of the sums of the rows
hist(apply(math,1,sum),breaks=c((0:33)-0.5))
Classical Test Theory Statistics
The three functions below give some of the basic Classical Test Theory statistics that we talked about in class. It is probably best to copy all three of the functions over at one time since sem uses alpha and ctt uses sem. Copy them over just as they are, you don't need to make any changes.
Once you have copied the three functions over, you can
run them by simply using the form function(arguments).
Thus, to get the measures of reliability for the math data set we could
use alpha(math). To get the 60% confidence intervals for true score
using alpha (instead of lambda2) you would type
ctt(math,rtype="alpha",conf.level=0.6)
Brief Descriptions of the Three Functions
alpha - Calculates both the estimated Cronbach's alpha and Guttman's lambda2 statistics.
sem - Calculates the estimated sem using both alpha and lambda2 to estimate the reliability. sem1 uses alpha, and and sem2 uses lambda2.
ctt - Gives the frequency distribution of observed scores, the regression estimate of true score, the sem based confidence interval for T, and the regression based confidence interval for T. It gives you the option of using either alpha or lambda2 (lambda 2 is the default) and your choice of confidence level (0.95 is the default). At the top of the ouput it gives you the sample mean, variance and standard deviation of the observed scores (X). Below that it gives a table. The columns of the table are:
The Functions Themselves
alpha<-function(testdata){
n<-ncol(testdata)
nexmn<-nrow(testdata)
x<-apply(testdata,1,sum)
s2y<-diag(var(testdata))*(nexmn-1)/nexmn
s2x<-var(x)*(nexmn-1)/nexmn
alpha<-(n/(n-1))*(1-sum(s2y)/s2x)
s2yy<-(((var(testdata)-diag(diag(var(testdata))))*(nexmn-1)/nexmn))^2
lambda2<-1-sum(s2y)/s2x+sqrt((n/(n-1))*sum(s2yy))/s2x
list(alpha=alpha,lambda2=lambda2)}
sem<-function(testdata){
nexmn<-nrow(testdata)
x<-apply(testdata,1,sum)
rhat<-alpha(testdata)
sx<-sqrt(var(x)*(nexmn-1)/nexmn)
sema<-sx*sqrt(1-rhat$alpha)
seml2<-sx*sqrt(1-rhat$lambda2)
list(sema=sema,seml2=seml2)}
ctt<-function(testdata,rtype="lambda2",conf.level=0.95){
nexmn<-nrow(testdata)
x<-apply(testdata,1,sum)
avgx<-mean(x)
varx<-var(x)*(nexmn-1)/nexmn
sdx<-sqrt(varx)
maxx<-max(x)
xrange<-c(0:maxx)
count<-rep(0,maxx)
for (i in 0:maxx){
count[i+1]<-sum(rep(1,nexmn)[x==i])}
x<-c(0:maxx)
z<--1*qnorm((1-conf.level)/2)
if (rtype=="alpha") {
rhat<-alpha(testdata)$alpha
sem.est<-sem(testdata)$sema
}
else {
rhat<-alpha(testdata)$lambda2
sem.est<-sem(testdata)$seml2
}
losem<-x-z*sem.est
hisem<-x+z*sem.est
treg<-rhat*x+(1-rhat)*avgx
see.est<-sem.est*sqrt(rhat)
loreg<-treg-z*see.est
hireg<-treg+z*see.est
list(avgx=avgx,varx=varx,sdx=sdx,summary=cbind(x,treg,count,losem,hisem,loreg,hireg))}
The function items will calculate the p-value, the variance, the D based on a user specified percent of the examinees (default=27%), the point biserial, and the biserial. If correct=T (the default) it will remove the item in question from the sum score, otherwise it will include that item. dpct tells it what percent of the examinees to use at each end for calculating D. digits tells it how many digits to report after the decimal (the default is four).
Thus, items(math) will analyze
math at the default settings, while
items(math,correct=F,dpct=0.5,digits=3) will analyze it
without removing the item in question, using 50% instead of 27%, and
using only 3 digits after the decimal place.
items<-function(testdata,correct=T,dpct=0.27,digits=4){
n<-ncol(testdata)
nexmn<-nrow(testdata)
x<-apply(testdata,1,sum)
medx<-median(x)
lowx<-sort(x)[floor((nexmn+1)*dpct)]
highx<-sort(x)[ceiling((nexmn+1)*(1-dpct))]
pvalue<-as.numeric(apply(testdata,2,mean))
var<-as.numeric(apply(testdata,2,var)*(nexmn-1)/nexmn)
D<-rep(0,n)
PBis<-rep(0,n)
Bis<-rep(0,n)
for (i in 1:n){
if (correct==T){
xx<-x-testdata[,i]
lowx<-sort(xx)[floor((nexmn+1)*dpct)]
highx<-sort(xx)[ceiling((nexmn+1)*(1-dpct))]
}
else {xx<-x}
pu<-mean(testdata[xx >= highx,i])
pl<-mean(testdata[xx <= lowx,i])
D[i]<-pu-pl
PBis[i]<-cor(xx,testdata[,i])
Bis[i]<-sqrt(pvalue[i]*(1-pvalue[i]))*PBis[i]/dnorm(qnorm(pvalue[i]))
}
round((10^digits)*cbind(pvalue,var,D,PBis,Bis))/(10^digits)}
The three functions below are:
irf - given theta, a, b, c, and the type of model ("3pl" or "norm"), this function returns the probability the examinee will get the item correct. By default the guessing parameter cc is set to 0. The following would give the values to plot for item 1 in problem 1 on page 28 of the text (answer on page 30).
irf(theta=c(-3,-2,-1,0,1,2,3),a=1.8,b=1.0,cc=0.0,type="3pl") round(irf(theta=c(-3,-2,-1,0,1,2,3),a=1.8,b=1.0,cc=0.0,type="3pl"),3)Notice that you are able to enter multiple theta values as a vector and it will return the value of the irf at each separate theta. If you enter multiple items at the same time it will give you the sum of the item response functions at those theta values. For example, the following will give column 6 of table 5.2 on page 87.
irf(theta=c(-3,-2,-1,0,1,2,3),a=c(0.8,1,1.2,1.5,2),b=c(-2,-1,0,1,2),cc=c(0,0,0.1,0.15,0.2))
irfplot - will give a plot of the specified item response function for abilities ranting from -3.5 to +3.5. It allows you to label the graph, and make it dashed if you would like. The following will construct a graph of the same parameters for both the normal ogive (dashed) and logistic (not-dashed) models. 3PL and not-dashed are both the default settings.
irfplot(a=1.8,b=1.0,cc=0.0,label="item 1 comparison") par(new=T) irfplot(a=1.8,b=1.0,cc=0.0,type="norm",dashed=T)The command par(new=T) allows you to overlay the first and second plot (otherwise it would have erased the first one). If vectors are entered for a, b, and c the expected total score will be plotted for the entire set of items.
irtgen - this function will simulate data from an exam that follows either the normal-ogive or logistic models. It requires the number of examinees (default=10); the vector of a, b, and c parameters; the type of item ("3PL" or "norm"); and the mean and standard deviation of the examinee abilities (which are simulated to follow the normal distribution). The following sample code simulates an exam consisting of 10 items and 100 examinees, where each item had a=1, b=0, and c=1, where it followed the normal-ogive model, and where the distribution of examinee abilities was standard normal. It calls the simulated data set sampdat and runs items on it.
a<-rep(1,10) b<-rep(0,10) cc<-rep(0,10) sampdat<-irtgen(nexmn=100,avec=a,bvec=b,cvec=cc,type="norm",mthet=0,sdthet=1) sampdat items(sampdat)The rep(1,10) is telling R to make a vector that repeats 1 ten times. Type a to see thats what it did. Also, note that we had to use cc instead of c for the name of that vector. This is because R reserves c to be the command for making a vector (e.g. c(1,2,3). Finally, make sure to enter all three functions because both irfplot and irtgen use irf.
The Functions
irf<-function(theta=0,a=1,b=0,cc=0,type="3pl"){
n<-length(a)
prob<-0
if (type=="3pl"){
for (i in 1:n){
prob<-prob+cc[i]+(1-cc[i])/(1+exp(-1.7*a[i]*(theta-b[i])))
}
}
else{
for (i in 1:n){
prob<-prob+cc[i]+(1-cc)[i]*pnorm(a[i]*(theta-b[i]))
}
}
return(prob)
}
irfplot<-function(a=1,b=0,cc=0,type="3pl",label="",dashed=F){
prob<-rep(0,71)
theta<-rep(0,71)
lim<-length(a)
if (lim==1){lab="Probability"}
else {lab="Expected Score"}
for (i in 1:71){
theta[i]<-(i-36)/10
prob[i]<-irf(theta[i],a,b,cc,type)
}
if (dashed==F){
plot(theta,prob,xlab="Ability",ylab=lab,type="l",
xlim=c(-3.5,3.5),ylim=c(0,lim),main=label)
}
if (dashed==T){
plot(theta,prob,xlab="Ability",ylab=lab,type="l",
xlim=c(-3.5,3.5),ylim=c(0,lim),main=label,lty=2)
}
}
irtgen<-function(nexmn=10,avec=c(1),bvec=c(0),cvec=c(0),type="3pl",
mthet=0,sdthet=1){
n<-length(avec)
data<-matrix(0,nrow=nexmn,ncol=n)
ability<-rnorm(nexmn,mean=mthet,sd=sdthet)
for (j in 1:nexmn){
for (i in 1:n){
try<-runif(1,0,1)
if (try < irf(ability[j],avec[i],bvec[i],cvec[i],type)){
data[j,i]<-1
}
}
}
return(data)
}
The following two functions concern the information functions of items and tests. As the second function relies on the first, it is probably best to enter them as a pair. In each a single item can be analyzed by entering its particular a, b, and c parameters. If vectors of parameters are entered then the analysis will be done for a test consisting of those items. For example a=c(1,1.5,1), b=c(0.5,1,1), cc=c(0,0,0.2) would be used to get the information for a three item test with parameters: a1=1, b1=0.5, c1=0, a2=1.5, b2=1, c2=0, a3=1, b3=1, and c3=0.2.
irfinfo - returns the information for the specified item or test at the theta values requested. It will function for either a logistic or ogive form using the options type="3pl" or type="norm".
The following commands will generate part of the answers to question 3 on page 97 of the text.
irfinfo(theta=c(-2,-1,0,1,2),a=c(1.25,1.5,1.25),b=c(-0.5,0,1), cc=c(0,0,0),type="3pl")
irfinfoplot - returns the plotted information function, with the same options as irfplot above. The only difference in the options is that the user can specify the maximum value of the y-axis for the plots using maxy=. The default is 1.1 times the maximum information for that item. This is important for overlaying the plots of several information functions
The following commands will generate the information curves for the two items (and the two item test) where the parameters are a1=1.5, b1=-0.5, c1=0.2, a2=1, b2=1, c2=0.25.
par(mfrow=c(2,1)) irfinfoplot(a=1.5,b=-0.5,cc=0.2,type="3pl",label="two items information", dashed="T",maxy=1.25) par(new=T) irfinfoplot(a=1,b=1,c=0.25,type="3pl",label="",dashed="F",maxy=1.25) irfinfoplot(a=1.5,b=-0.5,cc=0.2,label="test info for two items", dashed="T",maxy=1.25) par(new=T) irfinfoplot(a=1,b=1,cc=0.25,dashed="T",maxy=1.25) par(new=T) irfinfoplot(a=c(1.5,1),b=c(-0.5,1),cc=c(0.2,0.25),maxy=1.25)
Recall that the par(mfrow=c(2,1)) tells R to plot the graphics so that there are two rows and one column in the plot window, and that par(new=T) means to plot over the current figure. Also notice that the default settings are: type="3pl", dashed="F", and label="", so we don't need to include those.
The code for the two functions is:
#-----------------------------------------------------------
irfinfo<-function(theta=0,a=1,b=0,cc=0,type="3pl"){
n<-length(a)
info<-0
if (type=="3pl"){
for (i in 1:n){
info<-info+2.89*a[i]^2*(1-cc[i])/(
(cc[i]+exp(1.7*a[i]*(theta-b[i])))*
((1+exp(-1.7*a[i]*(theta-b[i])))^2)
)
}
}
else{
for (i in 1:n){
info<-info+
a[i]^2*(1-cc[i])^2*dnorm(a[i]*(theta-b[i]))^2/
(irf(theta=theta,a=a[i],b=b[i],cc=cc[i],type="norm")*
(1-irf(theta=theta,a=a[i],b=b[i],cc=cc[i],type="norm")))
}
}
return(info)
}
irfinfoplot<-function(a=1,b=0,cc=0,type="3pl",label="",dashed="F",maxy=0){
info<-rep(0,71)
theta<-rep(0,71)
for (i in 1:71){
theta[i]<-(i-36)/10
info[i]<-irfinfo(theta=theta[i],a=a,b=b,cc=cc,type=type)
}
if (maxy==0){maxy<-1.1*max(info)}
if (dashed=="F"){
plot(theta,info,xlab="theta",ylab="information",type="l",
xlim=c(-3.5,3.5),ylim=c(0,maxy),main=label)
}
if (dashed=="T"){
plot(theta,info,xlab="theta",ylab="information",type="l",
xlim=c(-3.5,3.5),ylim=c(0,maxy),main=label,lty=2)
}
}
#-----------------------------------------------------------
irtsave - will save an IRT data matrix with several different options. The default is with an id number, but no spaces between the response. To remove the id number use id=F and to add spaces between the responses use sep=T. The default is for the id number to be at least five digits long (it will automatically adjust be long enough to have a separate number for each examinee). The following examples show what form the file will take for the various options. (It is not limited to dichotomous items.)
--------------------------------------------------- irtsave(math,"d:/newmath.txt") 00001 11011111110001100000001110000010 00002 11111100110000001010001000100000 00003 00001011110000000000000010001000 etc... --------------------------------------------------- irtsave(math,"d:/newmath.txt",sep=T) 00001 1 1 0 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 0 00002 1 1 1 1 1 1 0 0 1 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 00003 0 0 0 0 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 etc... --------------------------------------------------- irtsave(math,"d:/newmath.txt",id=F) 11011111110001100000001110000010 11111100110000001010001000100000 00001011110000000000000010001000 etc... --------------------------------------------------- irtsave(math,"d:/newmath.txt",ndigits=8) 00000001 11011111110001100000001110000010 00000002 11111100110000001010001000100000 00000003 00001011110000000000000010001000 etc... ---------------------------------------------------The following code creates the data set discussed in homework 3 question 7 where each group of 4 items on the math test are combined into a single polytomous item. It then saves it as the file polymath.dat.
polymath<-cbind(apply(math[,1:4],1,sum), apply(math[,5:8],1,sum),apply(math[,9:12],1,sum), apply(math[,13:16],1,sum),apply(math[,17:20],1,sum), apply(math[,21:24],1,sum),apply(math[,25:28],1,sum), apply(math[,29:32],1,sum)) irtsave(polymath,"z:/polymath.dat")
The function itself:
#---------------------------------------------------------
irtsave<-function(datamatrix,flnm,id=T,sep=F,ndigits=5){
nitem<-length(datamatrix[1,])
nexmn<-length(datamatrix[,1])
ndigits<-max(ndigits,ceiling(log10(nexmn+1)))
datamatrix<-matrix(as.character(as.matrix(datamatrix)),ncol=nitem,nrow=nexmn,byrow=F)
if (sep==F) {
s<-""
}
else {
s<-" "
}
datamatrix<-apply(datamatrix,1,paste,collapse=s)
datamatrix<-as.matrix(datamatrix)
if (id==T){
idvec<-NULL
for (j in 1:nexmn){
idvec<-rbind(idvec,paste(
paste(rep("0",(ndigits-ceiling(log10(j+1)))),sep="",collapse=""),
as.character(j),sep="",collapse=""))
}
datamatrix<-apply(cbind(idvec,datamatrix),1,paste,collapse=" ")
}
write(as.table(as.matrix(datamatrix)),file=flnm)
}
#---------------------------------------------------------
Reading In BILOG and PARSCALE Output
BILOG: The following code will read in the estimated item parameters and ability estimates produced by BILOG for the data set mdataid.dat (as produced by the basics of BILOG example).
mbscores<-read.table("z:/MATH.SCO",head=F,fill=T,skip=2)[(1:2642)*2,6:7]
colnames(mbscores)<-list("estimate","se")
rownames(mbscores)<-1:2642
The 1:2642 refers to the number of examinees.
mbpars<-read.table("z:/MATH.PAR",head=F,fill=T,skip=4)[,c(4,5,6,7,10,11)]
colnames(mbpars)<-list("a","se.a","b","se.b","c","se.c")
PARSCALE: The following code will read in the estimated item parameter and ability
estimates produced by PARSCALE for the data set polymath.dat (as produced by
the basics
of PARSCALE example).
pmpscores<-read.table("z:/POLYMATH.SCO",head=F,fill=T)[(1:2642)*2,7:8]
colnames(pmpscores)<-list("estimate","se")
rownames(pmpscores)<-1:2642
The 1:2,642 refers to the number of examinees
pmpbase<-read.table("z:/POLYMATH.PAR",head=F,fill=T,skip=5)
pmppars<-cbind(pmpbase[(1:8)*3-2,c(3,5)],pmpbase[(1:8)*3-1,2:5],
pmpbase[(1:8)*3-2,c(4,6)],pmpbase[(1:8)*3,2:5])
colnames(pmppars)=list("a","b","d1","d2","d3","d4",
"se.a","se.b","se.d1","se.d2","se.d3","se.d4")
rownames(pmppars)=1:8
The 1:8 refers to the number of items. The 2:5 is because there were five response categories.
A set of functions to produce residual plots for dichotomous IRT models as suggested in chapter 4 of the text can be found at 778residS06.txt and loaded into R using:
source("http://www.stat.sc.edu/~habing/courses/778residS06.txt").
A brief description of the functions are given below, along with examples of their use with the math data set, mbscores, and mbpars estimates loaded in previously.
Several inputs appear in the various functions:
The functions available are:
irfresidplot - produces graphs like Figure 4.5 to 4.7 in the text. For example:
irfresidplot(2,math,mbscores,mbpars,label="Observed Proportions for Item 2")
itemresidplot - produces graphs like Figures 4.8 to 4.10 in the text. For example:
itemresidplot(2,math,mbscores,mbpars,
label="Standardized Residuals for Item 2",ylim=c(-5,5))
bisresidplot - produces a graph like Figures 4.14 to 4.16 in the text. For example:
bisresidplot(math,mbscores,mbpars,label="|StResid| vs. Biserial Correlation")
pvalresidplot - produces a graph of the absolute standardized residuals versus the item p-values (see Figures 9-8 and 9-9 in Hambleton and Swaminathan, 1985). For example:
pvalresidplot(math,mbscores,mbpars,label="|StResid| vs. % Correct")
residsummary - produces a histogram of the standardzied residuals (half of Figures 4.11 to 4.13 in the text). It also outputs the summary of the absolute standardized residuals needed to construct a table like 9-4 in Hambleton and Swaminathan (1985). For example:
residsummary(math,mbscores,mbpars)
Markov Chain Monte Carlo estimation can be performed in R using the packages BRugs, lattice, and coda. saved BRugs zip file. This requires having at least R 2.1.1 installed on your computer. The lattice packages should automatically be available with any of the semi-recent versions of R and can be loaded by simply typing: library(lattice). To download the others, go to Packages and select Install Packages. Choose one of the USA sites listed and then choose both the BRugs and coda packages for downloading. Once they have finished downloading, use the commands library(BRugs) and library(coda) to make them active.
BRugs requires three saved files in order to perform its estimation: the model file, the data file, and the inital values files. All three of these must be in a very specific format. A 3PL model file can be downloaded from model3PL.txt and saved in the directory of your choice. The other two files can be created using the irtbugdata and irtbuginit functions that can be loaded from 778mcmcS06.txt using source("http://www.stat.sc.edu/~habing/courses/778mcmcS06.txt") . irtbugdata takes the data matrix and the file you want it saved in for its arguments. irtbuginit takes those also, but also can take the number of categories for the initial estimate of the c parameter (default=5).
Assuming the three libraries have been loaded in, the model file saved, and the two functions sourced (not to mention that you've loaded in your data set!), the following will perform example of the estimation method using a portion of the math data set:
#Take only the first 100 examinees and 5 items to make it a fast example.
# You wouldn't do this for really analzying the data, you would use the
# whole thing and take a lot longer!
minimath<-math[1:100,1:5]
#Save the file, in this case on the Z-drive
irtbugdata(minimath,"z:/minimath.txt")
#Check the model, which was previously saved on the Z-drive
modelCheck("z:/model3PL.txt")
#Load in and check the data we previously saved above
modelData("z:/minimath.txt")
#Compile the model. This step can take 5 minutes or more for large
# data sets. Note that it is generally unadvisable to use only a
# single chain as it doesn't let you check for convergence. For purposes
# of the exam that is fine though.
modelCompile(numChains=1)
#Write the file of initial values. Starting with values reasonably close to
# the true values should help convergence. The function sets the inital
# a and b parameters using the formulas relating the classical item
# statistics to the normal ogive model, c to be 1/#categories, and theta
# to the normal quantile based on the rank of the examinees total score
# divided by the total number of examinees +1.
irtbuginit(minimath,"z:/minimathinits.txt",cats=5)
#Read in the inital values.
modelInits("z:/minimathinits.txt")
#Perform a 1,000 long burn-in so that the chain will hopefully
# converge. This takes about 20 seconds on my older
# PC for this example. For a bigger example with almost
# 2,000 examinees, 25 items, and only 500 updates my newer PC
# takes about 10 minutes for this step (about 1 sec. per update).
# It may even look like the program has crashed, but it hasn't.
# The length of the burn in should not be arbitrarily set
# and you should examine how the chains behave to see if
# it appears to have converged. In some applications this
# may need to be set to 10,000 or more.
modelUpdate(1000)
#Start recording the values of the parameters
samplesSet(c("a","b","c","theta"))
#Get 1,000 more observations to use for estimation. This is the
# same process as the modelUpdate two steps ago. In some
# applications you may want to make this 100,000 or more and
# will also want to investigate if thinning is needed.
modelUpdate(1000)
#View numerical summaries of the a's, b's, and c's and the
# theta estimates of the first three examinees. The summary
# statistics include the posterior mean and sd, the MC_error is
# the estimated standard error for the mean, the next three values
# are the 2.5%, 50%, and 97.5% values from the posterior, and finally
# the first value that was included in these calculations and
# how many values were included total.
samplesStats("a")
samplesStats("b")
samplesStats("c")
samplesStats("theta[1:3]")
#A large number of additional options are discussed at:
# http://cran.r-project.org/doc/packages/BRugs.pdf
# Of particular note, are functions such as samplesDensity,
# and samplesHistory. samplesSample is used similarly
# and gives you the actual values. The following code
# would give you the histogram of the posterior for the
# first a parameter and the first examinee's ability.
hist(samplesSample("a[1]"))
hist(samplesSample("theta[1]"))
CCOV Analysis and Mokken Scaling
CCOV Analysis: The code for running the conditional covariance analysis has functions in both R as well as in the form of a fortran based .dll file. To load the functions first use
source("http://www.stat.sc.edu/~habing/courses/778/mirtr.txt")
to load in the R functions. Next, right click on the link: ccprox.dll and choose save as. Save this file on a drive you have access to. Then run code similar to:
dyn.unload("z:/ccprox.dll")
dyn.load("z:/ccprox.dll")
where the directory is where you actually have saved the file (not necessarily z:/). Note that the first line may give an error and you should not worry about this.
The two main functions are:
hcaccprox - performs the HCA/CCPROX hierarchical cluster analysis (see Roussos, Stout, and Marden - 1998 JEM) on either dichotomous or polytomous data. It could be performed on the math data set using hcaccprox(math). The function has several options that may be specified:
mdsccprox - performs a two-dimensional isocmetric multidimensional scaling of the items using the CCPROX distances. Use help(isoMDS) for references. The options are the same as for the hcaccprox function except that method= is not applicable. It would be run on the math data set using mdsccprox(math).
Mokken Scaling: Mokken scaling is based around scalability coefficients,
which are ratios of observed covariances to the maximum possible covariances
given the item p-values. There are three types of scalability coefficients:
Hil is the scalability for item pair (i,l); Hi is the
scalability that item i has with the other items in the cluster; and H is
a measure of the total scalability for the cluster. A set of items is said
to form a Mokken scale if all of the Hil for that scale are positive
and if all of the Hi for items in that scale with the rest of
the scale are above some constant c. If all of the Hi (and therefore
H) are between 0.3 and 0.4 it is said to be a weak scale, between 0.4 and 0.5
a moderate scale, and a strong scale if at least 0.5.
One method of using Mokken scaling is to have a search algorithm that will refine a large set of items into smaller sets that satisfy the properties of a Mokken scale for a given cut-value c. Code to perform Mokken Scaling cluster construction can be loaded using: source("http://www.stat.sc.edu/~habing/courses/778/mokken.txt") and can be run on the math data using the default settings by: mokken(math).
The output from the Mokken scaling function is a list of the clusters that
formed using the chosen options. Each item is listed with its corresponding
Hi (how strongly it is related to the other items) as well
as the total H for the entire cluster, the smallest Hi for those
items chosen, and the estimated reliability using alpha.
Running Yen's Q3 requires having the original data, the estimated abilities, and the parameter estimates, just like getting the residual plots would. Q3 is simply the correlation between the residuals for two of the items (with the suggested bias correction of adding 1/(#items-1)). The test statistic and p-value are calculated using Fisher's R-to-Z transformation x=invtanh(r) is approximately normal with variance 1/(#examinees-3) when the true correlation is zero.
Assuming you have read in the data, the estimated abilities, and the estimated item parameters, the following code will run the procedure on items 1 and 2 of the math test we have been studying.
source("http://www.stat.sc.edu/~habing/courses/778residS06.txt")
q3(1,2,math,mbscores,mbpars
The relationship is positive (Q3=0.0512) indicating they are on the
"same side", and the p-value is 0.0084 for the two sided hypothesis test.
The code for the function itself is:
#----Q3 function--------------------------
q3<-function(item1,item2,dat,scores,pars){
n<-ncol(dat)
N<-nrow(dat)
r<-resids(dat,scores,pars,ngroup=0)$resid
q3<-cor(r[,item1],r[,item2])+1/(n-1)
zstat<-atanh(q3)*sqrt(N-3)
pvalgreater<-1-pnorm(zstat)
pvalless<-pnorm(zstat)
pvalnoteq<-2*(min(pvalgreater,pvalless))
out<-c(q3,zstat,pvalgreater,pvalnoteq,pvalless)
names(out)<-c("Q3","Z","HA:Q3>0","HA:Q3!=0","HA:Q3<0")
round(out,4)}
#-------------------------------------------
R features a built in procedure for calculating the Mantel-Haenszel test statistic for the null hypothesis that the common odds ratio is 1, and for constructing a 95% confidence interval for the common odds ratio. The code given below organizes the reference and focal groups appropriately and calculates the odds ratio estimate and confidence interval on the ETS delta scale.
Note that it is not set to give the p-value from the hypothesis test, but instead the test results can be gotten by using the confidence interval.
The following code examines the first 10 items on the data sets m1 and m2 used on homework assignment 2, treating m1 as the reference group.
m1<-read.fwf("http://www.stat.sc.edu/~habing/courses/data/m1.txt",widths=rep(1,32))
m2<-read.fwf("http://www.stat.sc.edu/~habing/courses/data/m2.txt",widths=rep(1,32))
MHDIF(m1,m2,1:10)
Note that the confidence intervals for 1-6 and 8-10 all contain 0, and so all of those items are A items (they are not significantly different from zero). Item 7 is significantly different from 0 (the CI does not contain it), but the |delta| is less than 1, and so it is still classified as an A item. An item with 1.6, (1.1, 2.4) would be a C item favoring m1, -1.4, (-2.2,-1.1) would be an B item favoring m2.
#--- function for performing MH DIF analysis--------------------
MHDIF<-function(refdat,focdat,items=(1:ncol(refdat))){
fulldat<-rbind(refdat,focdat)
group<-c(rep(1,nrow(refdat)),rep(2,nrow(focdat)))
s<-apply(fulldat,1,sum)
bad<-rep(0,(ncol(refdat)+1))
for (i in 1:(ncol(refdat)+1)){
if (sum(s[group==1]==(i-1))==0){bad[i]=1}
if (sum(s[group==2]==(i-1))==0){bad[i]=1}
}
bad<-(1:(ncol(refdat)+1))[bad==1]-1
for (i in bad){
fulldat<-fulldat[s!=i,]
group<-group[s!=i]
s<-s[s!=i]
}
out<-matrix(0,ncol=4,nrow=length(items))
for (i in items){
temp<-mantelhaen.test(as.factor(group),as.factor(fulldat[,i]),as.factor(s))
out[i,2]<-temp$estimate
out[i,4:3]<-temp$conf.int
}
out<--2.35*log(out)
out[,1]<-items
colnames(out)<-c("item","delta","lower-ci","upper-ci")
round(out,3)
}
#-----------------------------------------------------------------
The function eqyx below carries out equalpercentile equating as defined on pages 43-47 of Kolen and Brennan (2004). The required input is:
x<-c(0,0,1,1,1,2,2,3,3,4)
y<-c(0,1,1,2,2,3,3,3,4,4)
The following code would give the form y equalpercentile equivalent to an x score of 2, as well as the entire conversion table.
eqyx(2,x,y)
eqyx(0:4,x,y)
One way of constructing a graphical display would be to construct a scatter plot showing how form x raw scores convert to form y scores and vice-versa. Based on the piece-wise nature of the functions this is best done by observing how the scores from -0.5, 0.5, to max-score+0.5 convert, and this is the procedure implemented using the function: eqplot(x,y).
The required functions are:
#-------------------------------------
Fx<-function(x,pop){
max(0,min(1,mean(pop<=x)))}
Px<-function(x,pop){
xstar<-floor(x+.5)
max(0,min(100, 100*(Fx(xstar-1,pop)+(x-(xstar-.5))*(Fx(xstar,pop)-Fx(xstar-1,pop)))))}
xstaru<-function(pstar,pop){
val<-pop[max(1,floor((pstar/100)*length(pop)))]
while (Fx(val,pop)<=(pstar/100)){
val<-val+1}
val}
xstarl<-function(pstar,pop){
val<-pop[min(length(pop),ceiling((pstar/100)*length(pop)))]
while (Fx(val,pop)>=(pstar/100)){
val<-val-1}
val}
Pinv<-function(pstar,pop,Kx=max(pop)){
if (pstar==100){xupstar<-Kx+.5}
else{
xsu<-xstaru(pstar,pop)
xupstar<-(pstar/100-Fx(xsu-1,pop))/(Fx(xsu,pop)-Fx(xsu-1,pop))+(xsu-.5)
}
if (pstar==0){xlpstar<--.5}
else{
xsl<-xstarl(pstar,pop)
xlpstar<-(pstar/100-Fx(xsl,pop))/(Fx(xsl+1,pop)-Fx(xsl,pop))+(xsl+.5)
}
mean(c(xlpstar,xupstar))}
eqyx<-function(xval,xpop,ypop,Ky=max(ypop)){
xpop<-sort(xpop)
ypop<-sort(ypop)
if (length(xval)==1){
Pinv(Px(xval,xpop),ypop,Ky)}
else{
temp<-sapply(xval,Px,pop=xpop)
out<-rbind(xval,sapply(temp,Pinv,pop=ypop))
rownames(out)<-list("score","equated to")
t(out)}
}
eqplot<-function(xpop,ypop,Kx=max(xpop),Ky=max(ypop)){
x<-(-1:Kx)+.5
y<-(-1:Ky)+.5
temp1<-eqyx(x,xpop,ypop,Ky)
temp2<-eqyx(y,ypop,xpop,Kx)
form1<-sort(c(temp1[,1],temp2[,2]))
form2<-sort(c(temp2[,1],temp1[,2]))
plot(form1,form2)
lines(form1,form2)
}
#------------------------------------------------