01_AlphaAndOmega.R

# Clear all existing objects in the workspace
rm(list=ls())

# Load required library for statistical analysis
library(psych) # used for PCA, Cronbach's Alpha, Guttman's Lambda6, McDonald's Omega

# Set seed for reproducibility
set.seed(1828)

# Configure options to suppress significance stars in outputs
options(show.signif.stars=FALSE)

# Load the dataset
load('pcDat.rda')

###################
### Proficiency ###
###################

# Calculate Cronbach's Alpha for reliability analysis of selected variables
# Automatically reverses the items if needed
psych::alpha(pcDat[,3:8], check.keys=TRUE)
# Reliability analysis
# Call: psych::alpha(x = pcDat[, 3:8], check.keys = TRUE)
#
#   raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
#       0.78      0.78    0.83      0.37 3.5 0.056  2.5 0.63     0.28
#
#     95% confidence boundaries
#          lower alpha upper
# Feldt     0.65  0.78  0.87
# Duhachek  0.67  0.78  0.89
#
#  Reliability if an item is dropped:
#                       raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
# proficiency_listening      0.72      0.72    0.78      0.34 2.6    0.072 0.076  0.24
# proficiency_writing        0.71      0.71    0.75      0.33 2.5    0.076 0.059  0.24
# proficiency_speaking       0.73      0.73    0.79      0.35 2.7    0.071 0.076  0.24
# proficiency_reading        0.68      0.69    0.74      0.30 2.2    0.083 0.071  0.12
# proficiency_grammar-       0.80      0.81    0.83      0.46 4.2    0.052 0.073  0.60
# proficiency_vocab.z-       0.80      0.80    0.82      0.44 3.9    0.047 0.085  0.60

# Compute McDonald's Omega to measure internal consistency reliability
# Specify 2 factors based on theoretical assumptions or EFA results
omega(pcDat[,3:8], nfactors=2)
# Omega 
# Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
#     digits = digits, title = title, sl = sl, labels = labels, 
#     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
#     covar = covar)
# Alpha:                 0.78 
# G.6:                   0.83 
# Omega Hierarchical:    0.29 
# Omega H asymptotic:    0.34 
# Omega Total            0.87 
# 
# Schmid Leiman Factor loadings greater than  0.2 
#                          g   F1*   F2*   h2   u2   p2
# proficiency_listening 0.36  0.68       0.59 0.41 0.22
# proficiency_writing   0.39  0.79       0.77 0.23 0.19
# proficiency_speaking  0.33  0.67       0.56 0.44 0.20
# proficiency_reading   0.46  0.74       0.77 0.23 0.28
# proficiency_grammar-  0.31       -0.58 0.43 0.57 0.22
# proficiency_vocab.z-  0.38       -0.70 0.64 0.36 0.23
# 
# With Sums of squares  of:
#    g  F1*  F2* 
# 0.85 2.07 0.86 
# 
# general/max  0.41   max/min =   2.41
# mean percent general =  0.22    with sd =  0.03 and cv of  0.14 
# Explained Common Variance of the general factor =  0.22 
# 
# The degrees of freedom are 4  and the fit is  0.17 
# The number of observations was  40  with Chi Square =  5.99  with prob <  0.2
# The root mean square of the residuals is  0.03 
# The df corrected root mean square of the residuals is  0.07
# RMSEA index =  0.108  and the 10 % confidence intervals are  0 0.286
# BIC =  -8.77
# 
# Compare this with the adequacy of just a general factor and no group factors
# The degrees of freedom for just the general factor are 9  and the fit is  1.92 
# The number of observations was  40  with Chi Square =  68.33  with prob <  3.2e-11
# The root mean square of the residuals is  0.34 
# The df corrected root mean square of the residuals is  0.45 
# 
# RMSEA index =  0.405  and the 10 % confidence intervals are  0.323 0.505
# BIC =  35.13 
# 
# Measures of factor score adequacy             
#                                                   g  F1*  F2*
# Correlation of scores with factors             0.56 0.85 0.76
# Multiple R square of scores with factors       0.31 0.72 0.58
# Minimum correlation of factor score estimates -0.38 0.44 0.17
# 
#  Total, General and Subset omega for each subset
#                                                  g  F1*  F2*
# Omega total for total scores and subscales    0.87 0.89 0.70
# Omega general for total scores and subscales  0.29 0.20 0.16
# Omega group for total scores and subscales    0.58 0.69 0.54

# Perform Principal Component Analysis (PCA) on proficiency variables
# Using Varimax rotation for clearer factor structure
pca1 = principal(pcDat[,3:8], nfactors=2, rotate='varimax')
# Principal Components Analysis
# Call: principal(r = pcDat[, 3:8], nfactors = 2, rotate = "varimax")
# Standardized loadings (pattern matrix) based upon correlation matrix
#                         RC1   RC2   h2   u2 com
# proficiency_listening  0.84 -0.06 0.71 0.29 1.0
# proficiency_writing    0.90 -0.02 0.81 0.19 1.0
# proficiency_speaking   0.83 -0.01 0.69 0.31 1.0
# proficiency_reading    0.87 -0.22 0.80 0.20 1.1
# proficiency_grammar   -0.05  0.87 0.76 0.24 1.0
# proficiency_vocab.z   -0.11  0.87 0.76 0.24 1.0
# 
#                        RC1  RC2
# SS loadings           2.97 1.56
# Proportion Var        0.50 0.26
# Cumulative Var        0.50 0.76
# Proportion Explained  0.66 0.34
# Cumulative Proportion 0.66 1.00
# 
# Mean item complexity =  1
# Test of the hypothesis that 2 components are sufficient.
# 
# The root mean square of the residuals (RMSR) is  0.09 
#  with the empirical chi square  9.48  with prob <  0.05 
# 
# Fit based upon off diagonal values = 0.96

# Add the resulting principal components as new variables in the dataset
pcDat$Proficiency.PC1 = pca1$scores[,1] # First principal component
pcDat$Proficiency.PC2 = pca1$scores[,2] # Second principal component

# Analyze the relationship between years of immersion and the principal components
cor.test(pcDat$immersion_years, pcDat$Proficiency.PC1)
# 	Pearson's product-moment correlation
#
# data:  pcDat$immersion_years and pcDat$Proficiency.PC1
# t = 1.5157, df = 38, p-value = 0.1379
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.07858242  0.51218506
# sample estimates:
#       cor
# 0.2387719

cor.test(pcDat$immersion_years, pcDat$Proficiency.PC2)
# 	Pearson's product-moment correlation
#
# data:  pcDat$immersion_years and pcDat$Proficiency.PC2
# t = 0.016475, df = 38, p-value = 0.9869
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.3090940  0.3139207
# sample estimates:
#         cor
# 0.002672663

# Check mutual correlation between the two principal components
cor.test(pcDat$Proficiency.PC1, pcDat$Proficiency.PC2)
# 	Pearson's product-moment correlation
#
# data:  pcDat$Proficiency.PC1 and pcDat$Proficiency.PC2
# t = 1.6841e-15, df = 38, p-value = 1
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.3115093  0.3115093
# sample estimates:
#          cor
# 2.731901e-16

################
### Accuracy ###
################

# Create a correlation matrix for accuracy-related variables
# Set upper triangular matrix (including diagonal) to NA for cleaner output
# Round the correlation coefficients to two decimal places
cormtx = round(cor(pcDat[,9:12]), 2)
cormtx[upper.tri(cormtx, diag=TRUE)] = NA
print(cormtx, na.print='')
#                    accuracy_correct accuracy_incorrect accuracy_article accuracy_tense
# accuracy_correct                                                                      
# accuracy_incorrect             0.66                                                   
# accuracy_article               0.87               0.85                                
# accuracy_tense                 0.89               0.84             0.81               

# omega(pcDat[,9:12], nfactors=2)
# In factor.scores, the correlation matrix is singular, the pseudo inverse is  used

# Perform Principal Component Analysis (PCA) on accuracy variables
# Specify a single component based on theoretical or data-driven rationale
pca2 = principal(pcDat[,9:12], nfactors=1)
pca2
# Principal Components Analysis
# Call: principal(r = pcDat[, 9:12], nfactors = 1)
# Standardized loadings (pattern matrix) based upon correlation matrix
#                     PC1   h2    u2 com
# accuracy_correct   0.92 0.85 0.153   1
# accuracy_incorrect 0.90 0.81 0.191   1
# accuracy_article   0.95 0.90 0.099   1
# accuracy_tense     0.95 0.90 0.096   1
# 
#                 PC1
# SS loadings    3.46
# Proportion Var 0.87
# 
# Mean item complexity =  1
# Test of the hypothesis that 1 component is sufficient.
# 
# The root mean square of the residuals (RMSR) is  0.08 
#  with the empirical chi square  3.13  with prob <  0.21 
# 
# Fit based upon off diagonal values = 0.99

# Add the resulting principal component as a new variable in the dataset
pcDat$Accuracy.PC1 = pca2$scores[,1]

# Analyze the relationship between a combined accuracy score and the principal component
cor.test(pcDat$accuracy_all, pcDat$Accuracy.PC1)
# 	Pearson's product-moment correlation
#
# data:  pcDat$accuracy_all and pcDat$Accuracy.PC1
# t = 329.24, df = 38, p-value < 2.2e-16
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  0.9996662 0.9999080
# sample estimates:
#       cor
# 0.9998248