Search code examples
rcontinuous-integrationalphapsych

Access general confidence boundaries of cronbach's alpha in R psych


I have a dataframe on which I use psych::alpha. In the output there are general confidence boundaries around a general cronbach's alpha value. I want to access those but they don't appear in the results when I save the output as a variable. In the documentation they're called itemboot.ci but that doesn't exist in the alpha object.enter code here

library(psych)
set.seed(1)
foo <- data.frame(X = rnorm(25), Y = rnorm(25), Z = rnorm(25))
bar <- alpha(foo)

alpha(foo) and bar both give this output:

Reliability analysis   
Call: alpha(x = foo)

  raw_alpha std.alpha G6(smc) average_r      S/N  ase mean   sd median_r
     0.023  -0.00091  0.0049    -3e-04 -0.00091 0.32 0.12 0.54    0.022

 lower alpha upper     95% confidence boundaries
-0.61 0.02 0.66 

 Reliability if an item is dropped:
  raw_alpha std.alpha G6(smc) average_r    S/N alpha se var.r  med.r
X     0.039     0.043   0.022     0.022  0.045     0.35    NA  0.022
Y     0.093     0.094   0.049     0.049  0.103     0.36    NA  0.049
Z    -0.148    -0.155  -0.072    -0.072 -0.134     0.44    NA -0.072

 Item statistics 
   n raw.r std.r r.cor  r.drop  mean   sd
X 25  0.58  0.56 -0.12  0.0027 0.169 0.95
Y 25  0.41  0.55 -0.36 -0.0296 0.032 0.71
Z 25  0.72  0.62  0.61  0.0543 0.165 1.11

Where

lower alpha upper     95% confidence boundaries
-0.61 0.02 0.66

are the values I want to access and save as a variable. bar$itemboot.ci gives NULL even though it's mentioned in the documentation and even looping through a doesn't show the confidence boundaries:

for(i in 1:length(bar)){
    print(bar[i])
}

gives:

$total
  raw_alpha    std.alpha     G6(smc)    average_r           S/N      ase      mean       sd   median_r
 0.02297052 -0.000914273 0.004926574 -0.000304572 -0.0009134379 0.323334 0.1221019 0.544037 0.02180445

$alpha.drop
    raw_alpha   std.alpha     G6(smc)   average_r         S/N  alpha se var.r       med.r
X  0.03880703  0.04267832  0.02180445  0.02180445  0.04458096 0.3488970    NA  0.02180445
Y  0.09264777  0.09365915  0.04913033  0.04913033  0.10333767 0.3588215    NA  0.04913033
Z -0.14776806 -0.15482062 -0.07184849 -0.07184849 -0.13406464 0.4395377    NA -0.07184849

$item.stats
   n     raw.r     std.r      r.cor       r.drop       mean        sd
X 25 0.5843116 0.5644059 -0.1232943  0.002682803 0.16866521 0.9501080
Y 25 0.4056801 0.5486244 -0.3639981 -0.029593160 0.03223135 0.7062806
Z 25 0.7151930 0.6184929  0.6088273  0.054340331 0.16540907 1.1051952

$response.freq
NULL

$keys
X Y Z 
1 1 1 

$scores
 [1] -0.094825557 -0.194726192 -0.655087102 -0.004077450  0.726824345  0.839537022  0.005806616  0.027287230  0.363898646
[10] -0.605834183  1.166134789 -0.014562240  0.003061795 -0.362224119  0.381611153 -0.006888302 -0.691503524  1.035451483
[19]  0.510379244  0.692585766  0.228997248  0.145590611  0.484608087 -1.011931847  0.082433358

$nvar
[1] 3

$boot.ci
NULL

$boot
NULL

$Unidim
$Unidim$Unidim
[1] -0.0006057996


$var.r
[1] 0.004025575

$Fit
$Fit$Fit.off
[1] 0.9973251


$call
alpha(x = foo)

$title
NULL

Solution

  • When you print an object, either by using print or by sending it to the R console, some extra processing may happen. Objects often have their own print and in this case you can see that the print.psych method (called behind the scenes instead of print on any psych package object) is doing the following with your object of (sub)class alpha:

      }, alpha = {
        cat("\nReliability analysis ", x$title, " \n")
        cat("Call: ")
        print(x$call)
        cat("\n ")
        print(x$total, digits = digits)
        if (!is.null(x$total$ase)) {
          cat("\n lower alpha upper     95% confidence boundaries\n")
          cat(round(c(x$total$raw_alpha - 1.96 * x$total$ase, 
            x$total$raw_alpha, x$total$raw_alpha + 1.96 * 
              x$total$ase), digits = digits), "\n")
        }
        if (!is.null(x$boot.ci)) {
          cat("\n lower median upper bootstrapped confidence intervals\n", 
            round(x$boot.ci, digits = digits))
        }
    

    If you do the same thing with your bar output variable; you get those values:

    library(psych)
    set.seed(1)
    
    foo <- data.frame(bar = rnorm(25), Y = rnorm(25), Z = rnorm(25))
    
    bar <- alpha(foo)
    
    
    cat(round(c(bar$total$raw_alpha - 1.96 * bar$total$ase, 
                bar$total$raw_alpha, bar$total$raw_alpha + 1.96 * 
                  bar$total$ase), digits = 2), "\n")
    #> -0.61 0.02 0.66
    

    So, in summary, those values are obtained by "manually" calculating confidence intervals with the coefficient estimate raw_alpha plus or minus the z-score for 95% CI (1.96) times the standard errors (ase).