Classical Test Theory Reliability Functions

The data used to demonstrate these functions is read in first.

data_demo <- read.table("C:\\Users\\grant_morgan\\Desktop\\Teaching\\EDP 6337 - Psychometric Theory\\Spring 2014\\Data\\data_demo.DAT", 
    sep = ",", header = FALSE)

The code below contains 3 functions that can aid in classical theory analysis.

1) The alphafunction returns Cronbach’s alpha and Guttman’s lambda2.

2) The sem function returns the standard error of measurement based on Cronbach’s alpha and Guttman’s lambda2.

3) The ctt function 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 output 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:
a) x = observed score,
b) treg = the t-hat estimated using the regression method for that X,
c) count = the number of examinees with that observed score,
d) losem and hisem = the confidence interval for T using the form X+zSEM,
e) loreg and hireg = the confidence interval for T using the form T-hat+zSEE.

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)
}

The alpha function and its output are shown below.

alpha(data_demo)

## $alpha
## [1] 0.6788
## 
## $lambda2
## [1] 0.684

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)
}

The sem function and its output are shown below.

sem(data_demo)

## $sema
## [1] 1.008
## 
## $seml2
## [1] 1

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 ctt function and its output are shown below.

ctt(data_demo)

## $avgx
## [1] 8.48
## 
## $varx
## [1] 3.165
## 
## $sdx
## [1] 1.779
## 
## $summary
##        x  treg count    losem hisem loreg  hireg
##  [1,]  0 2.679     1 -1.95998  1.96 1.058  4.300
##  [2,]  1 3.363     3 -0.95998  2.96 1.742  4.984
##  [3,]  2 4.047     4  0.04002  3.96 2.426  5.668
##  [4,]  3 4.731     5  1.04002  4.96 3.110  6.352
##  [5,]  4 5.415     6  2.04002  5.96 3.794  7.037
##  [6,]  5 6.100    15  3.04002  6.96 4.478  7.721
##  [7,]  6 6.784    26  4.04002  7.96 5.163  8.405
##  [8,]  7 7.468    44  5.04002  8.96 5.847  9.089
##  [9,]  8 8.152    83  6.04002  9.96 6.531  9.773
## [10,]  9 8.836   137  7.04002 10.96 7.215 10.457
## [11,] 10 9.520   172  8.04002 11.96 7.899 11.141