Tetrachoric Corr

1 Pre-Process Data

Midterms <- cohorts |> filter(exam == "Midterm")

M_lst <- Midterms |> group_by(term) |> group_split()

names(M_lst) <- c("Fall", "Spring", "Winter")

Fall_binary <- M_lst$Fall |> select(contains("Q"))
Winter_binary <- M_lst$Winter |> select(contains("Q"))
Winter_binary <- Winter_binary |> select(Q1:Q33) # only has 33 Q's
Spring_binary <- M_lst$Spring |> select(contains("Q"))

2 Calculation

library(psych) # where tetrachoric() is from

Fall_tetra <- tetrachoric(Fall_binary, smooth = TRUE)

## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done

F_corr <- Fall_tetra$rho

Winter_tetra <- tetrachoric(Winter_binary, smooth = TRUE)

## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done

W_corr <- Winter_tetra$rho

Spring_tetra <- tetrachoric(Spring_binary, smooth = TRUE)

## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done

S_corr <- Spring_tetra$rho


?tetrachoric

• Notice each is not postive definite (x’Ax ≥ 0) so according to literature, its suggested we apply some smoothing.

2.1 Fall tetrachoric Vis.

2.2 Based on question order

library(corrplot)

## corrplot 0.95 loaded

corrplot(
  F_corr,
  method = "color",
  type = "lower",
  tl.col = "black",
  tl.cex = 0.75,
  diag = FALSE
)

• we notice nearly all are positive

2.2.1 Based on clustering

corrplot(
  F_corr,
  method = "color",
  order = "hclust",
  hclust.method = "complete",
  type = "lower",
  tl.col = "black",
  tl.cex = 0.75,
  diag = FALSE
)

?corrplot

• looks

F_hclust <- hclust(
  as.dist(1 - F_corr),
  method = "complete"
)

plot(F_hclust)

3 Winter tetrachoric Vis.

corrplot(
  W_corr,
  method = "color",
  type = "lower",
  tl.col = "black",
  tl.cex = 0.75,
  diag = FALSE
)

• nearly all positive

3.0.1 Based on clustering

corrplot(
  W_corr,
  method = "color",
  order = "hclust",
  hclust.method = "complete",
  type = "lower",
  tl.col = "black",
  tl.cex = 0.75,
  diag = FALSE
)

W_hclust <- hclust(
  as.dist(1 - W_corr),
  method = "complete"
)

plot(W_hclust)

4 Spring tetrachoric Vis.

corrplot(
  S_corr,
  method = "color",
  type = "lower",
  tl.col = "black",
  tl.cex = 0.75,
  diag = FALSE
)

• appears to have weaker correlations

• more negative instances

4.0.1 Based on clustering

corrplot(
  S_corr,
  method = "color",
  order = "hclust",
  hclust.method = "complete",
  type = "lower",
  tl.col = "black",
  tl.cex = 0.75,
  diag = FALSE
)

S_hclust <- hclust(
  as.dist(1 - S_corr),
  method = "complete"
)

plot(S_hclust)

5 Together Global

5.1 Correlation plts

par(mfrow = c(1, 3))

corrplot(F_corr, method = "color", type = "lower",
         title = "Fall", mar = c(0, 0, 2, 0))

corrplot(W_corr, method = "color", type = "lower",
         title = "Winter", mar = c(0, 0, 2, 0))

corrplot(S_corr, method = "color", type = "lower",
         title = "Spring", mar = c(0, 0, 2, 0))

par(mfrow = c(1, 1))

• we must remember that the exams arent constant, they change as discussed in the following data below

# F23 example
Midterms_lst$F23_Midterm |> 
  select(question, Question_Order, difference)

5.2 Hierarchy plots

plot(F_hclust, main = "Fall", xlab = "", sub = "",
     cex = 0.65)

plot(W_hclust, main = "Winter", xlab = "", sub = "",
     cex = 0.65)

plot(S_hclust, main = "Spring", xlab = "", sub = "",
     cex = 0.65)

# reset plotting
par(mfrow = c(1, 1))

• notice that globally we see a similar trend across cohorts.

• it looks like for each theres 2 main clusters with 2 sub-cluster associated with each main clusters

• Spring and winter are more similar, which makes sense because Fall is the most dis-similar exam

6 Re-Plot

library(dendextend) # suggested by bio textbook

# forgot the names lol
grep("clust", ls(), value = TRUE)

## [1] "F_hclust" "S_hclust" "W_hclust"

# make dendrogram objects
F_dendrogram <- as.dendrogram(F_hclust)
W_dendrogram <- as.dendrogram(W_hclust)
S_dendrogram <- as.dendrogram(S_hclust)

print("F_dendrogram")

## [1] "F_dendrogram"

F_dendrogram

## 'dendrogram' with 2 branches and 34 members total, at height 1.19648

print("W_dendrogram")

## [1] "W_dendrogram"

W_dendrogram

## 'dendrogram' with 2 branches and 33 members total, at height 1.189564

print("S_dendrogram")

## [1] "S_dendrogram"

S_dendrogram

## 'dendrogram' with 2 branches and 34 members total, at height 1.283387

• notice they all have 2 branches!

?plot.dendrogram

plot(color_branches(F_dendrogram, k=2),leaflab="none")

plot(color_branches(W_dendrogram, k=2),leaflab="none")

plot(color_branches(S_dendrogram, k=2),leaflab="none")

• this plot makes those 2 main branches more clear

library(dendextend)

color_tree <- function(dend) {
  dend |>
    color_branches(
      k = 4,
      col = c("red", "blue", "green", "purple")
    ) |>
    set(
      "labels_col",
      value = c("red", "blue", "green", "purple"),
      k = 4
    ) |>
    set("labels_cex", 0.7)
}

plot(color_tree(F_dendrogram), main = "Fall")

plot(color_tree(W_dendrogram), main = "Winter")

plot(color_tree(S_dendrogram), main = "Spring")

• 2 main branches: Preserved at the highest split of each dendrogram.

• 4 subgroups: Defined by cutting the hierarchy into four clusters.

• these are not the same as the number of latent variables

7 Scree Plot : PCA

• Technically speaking we shouldn’t be using PCA

  • Our data is binary, not continous

# PCA from tetrachoric correlation matrices
F_pca <- princomp(covmat = F_corr)
W_pca <- princomp(covmat = W_corr)
S_pca <- princomp(covmat = S_corr)

7.1 Fall

screeplot(F_pca, type = "lines", main = "Fall")
abline(h = 1, col = "red", lty = 2) 

• suggests 1 strong factor; maybe 2-3 (based on elbow)

7.2 Winter

screeplot(W_pca, type = "lines", main = "Winter")
abline(h = 1, col = "red", lty = 2) 

• suggests 1 strong factor; maybe 2-3 (based on elbow)

7.3 Spring

screeplot(S_pca, type = "lines", main = "Spring")
abline(h = 1, col = "red", lty = 2) 

• suggests 1 strong factor; maybe 2-3 (based on elbow)

7.4 Global

# Three scree plots
par(mfrow = c(1, 3))

screeplot(F_pca, type = "lines", main = "Fall")
screeplot(W_pca, type = "lines", main = "Winter")
screeplot(S_pca, type = "lines", main = "Spring")

par(mfrow = c(1, 1))

• notice globally a similar shape

• It looks like for all, 1 dominant factor with 2-3 weaker-factors

  • similar to what dendrograms suggest

7.4.1 Kaiser criterion

print("Fall")

## [1] "Fall"

sum(F_pca$sdev^2 > 1)

## [1] 9

print("Winter")

## [1] "Winter"

sum(W_pca$sdev^2 > 1)

## [1] 11

print("Spring")

## [1] "Spring"

sum(S_pca$sdev^2 > 1)

## [1] 12

• here we compare it to just 1 standardized unit; which is a bit arbitrary

8 Parallel Analysis

• notice unlike before where we used PCA, here we use traditional factor analysis

• Comparing randomly generated data results to our actual data

8.1 Fall

?fa.parallel

set.seed(123)
par(cex = 0.7)

F_parallel <- fa.parallel(
  Fall_binary,
  fa = "both",
  cor = "tet", # tetrachoric corr
  n.iter = 300,
  quant = 0.95,
  main = "Fall"
)

## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done

## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0

## Parallel analysis suggests that the number of factors =  12  and the number of components =  8

abline(h = 1, col = "black") 

par(cex = 1)

• Suggests 12 factors

• 8 components

• 1-3 based on elbow test

8.2 Winter

set.seed(123)
par(cex = 0.7) # adjust legend size
W_parallel <- fa.parallel(
  Winter_binary,
  fa = "both",
  cor = "tet",
  n.iter = 300,
  quant = 0.95, # .95th quantile? 
  main = "Winter"
)

## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done

## Parallel analysis suggests that the number of factors =  14  and the number of components =  7

abline(h = 1, col = "black") 

par(cex = 1) # re-adjust legend size

• suggests 14 factor

• 7 components

• 1-3 based on elbow test

8.3 Spring

set.seed(123)
par(cex = 0.7)
S_parallel <- fa.parallel(
  Spring_binary,
  fa = "both",
  cor = "tet",
  n.iter = 300,
  quant = 0.95,
  main = "Spring"
)

## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done

## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0

## Warning in matrix(unlist(templist), ncol = n.iter): data length [40702] is not
## a sub-multiple or multiple of the number of rows [136]

## Parallel analysis suggests that the number of factors =  2  and the number of components =  4

abline(h = 1, col = "black") 

par(cex = 1)

• suggest 1 factor–very different

• number of components 4

• Based on elbow test 3

8.3.1 Examining PCAs more closely

pca_lst <- list(
  Fall = F_pca,
  Winter = W_pca,
  Spring = S_pca
)

lst <-
  lapply(pca_lst, function(pca) {
  eigenvalues <- pca$sdev^2

  c(
    PC1_variance = eigenvalues[1] / sum(eigenvalues),
    PC2_variance = eigenvalues[2] / sum(eigenvalues),
    PC1_PC2_ratio = eigenvalues[1] / eigenvalues[2]
  )
})

lst

## $Fall
##  PC1_variance.Comp.1  PC2_variance.Comp.2 PC1_PC2_ratio.Comp.1 
##           0.36948820           0.06996483           5.28105639 
## 
## $Winter
##  PC1_variance.Comp.1  PC2_variance.Comp.2 PC1_PC2_ratio.Comp.1 
##           0.33419720           0.05745207           5.81697430 
## 
## $Spring
##  PC1_variance.Comp.1  PC2_variance.Comp.2 PC1_PC2_ratio.Comp.1 
##           0.24018436           0.08635198           2.78145751

• Notice that this suggests :

  • Fall

    • ~37% of var explained by comp.1, ~7% comp.2

  • Winter

    • ~33% of var explained by comp.1, ~6% comp.2

  • Spring

    • ~24% of var explained by comp.1, ~9% comp.2

So overall, Fall & winter suggest 1 strong component & spring shows more spread out variation. It could be that the strength of measuring the common factor may vary across exam, not that its mostly just 1 factor.

9 Extract (and rotate) factors

• now we will investigate the factor structure. We will assume at most 2 factors

• here we use oblimin because we assume factors are correlated

  • more likely true

• here we use equamax because we assume factors are uncorrelated

  • probably not true given its a stats test

9.1 Fall

9.1.1 dependennt (correlated) factors

# Dependent factors (Oblimin); 1-3 factors
F_fa_dependent_lst <- lapply(1:3, function(k) {
  fa(
    F_corr,
    nfactors = k,
    fm = "pa",
    rotate = if (k == 1) "none" else "oblimin",
    n.obs = nrow(Fall_binary)
  )
})

## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0

## Loading required namespace: GPArotation

## In smc, smcs < 0 were set to .0

names(F_fa_dependent_lst) <- paste0("Dim", 1:3)

9.1.1.1 Factor Loadings

print(F_fa_dependent_lst$Dim1$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1  
## Q1  0.326
## Q2  0.683
## Q3  0.643
## Q4  0.710
## Q5  0.645
## Q6  0.723
## Q7  0.471
## Q8  0.684
## Q9  0.496
## Q10      
## Q11 0.810
## Q12      
## Q13 0.535
## Q14 0.517
## Q15 0.555
## Q16 0.584
## Q17 0.560
## Q18 0.629
## Q19 0.589
## Q20 0.621
## Q21 0.568
## Q22 0.660
## Q23      
## Q24 0.454
## Q25 0.708
## Q26 0.627
## Q27 0.760
## Q28 0.698
## Q29 0.706
## Q30 0.504
## Q31 0.670
## Q32 0.621
## Q33 0.688
## Q34 0.475
## 
##                   PA1
## SS loadings    11.979
## Proportion Var  0.352

print(F_fa_dependent_lst$Dim2$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA2   
## Q1               
## Q2   0.565       
## Q3          0.534
## Q4          0.765
## Q5   0.491       
## Q6   0.593       
## Q7          0.924
## Q8          0.653
## Q9   0.398       
## Q10  0.423 -0.309
## Q11  0.618  0.315
## Q12              
## Q13         0.363
## Q14  0.438       
## Q15  0.512       
## Q16  0.413       
## Q17  0.525       
## Q18  0.854       
## Q19  0.603       
## Q20  0.448       
## Q21  0.752       
## Q22  0.373  0.438
## Q23              
## Q24  0.592       
## Q25  0.812       
## Q26  0.504       
## Q27  0.496  0.411
## Q28  0.696       
## Q29  0.700       
## Q30  0.494       
## Q31  0.524       
## Q32  0.441       
## Q33  0.596       
## Q34  0.570       
## 
##                  PA1   PA2
## SS loadings    8.804 3.623
## Proportion Var 0.259 0.107
## Cumulative Var 0.259 0.366

print(F_fa_dependent_lst$Dim3$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA3    PA2   
## Q1                      
## Q2   0.680              
## Q3          0.367  0.452
## Q4          0.405  0.663
## Q5   0.332  0.303       
## Q6   0.643              
## Q7                 0.906
## Q8                 0.594
## Q9          0.331       
## Q10  0.532              
## Q11  0.505              
## Q12                     
## Q13  0.436         0.375
## Q14  0.462              
## Q15         0.401       
## Q16  0.655         0.308
## Q17  0.310  0.348       
## Q18  0.536  0.434       
## Q19  0.529              
## Q20  0.438              
## Q21  0.627              
## Q22         0.314  0.364
## Q23         0.628       
## Q24  0.494              
## Q25  0.889              
## Q26  0.431              
## Q27         0.806       
## Q28  0.482  0.349       
## Q29         0.832       
## Q30  0.415              
## Q31  0.442              
## Q32         0.412       
## Q33  0.470              
## Q34  0.336  0.334       
## 
##                  PA1   PA3   PA2
## SS loadings    6.171 3.648 2.998
## Proportion Var 0.181 0.107 0.088
## Cumulative Var 0.181 0.289 0.377

9.1.1.2 Factor correlations :

F_fa_dependent_lst$Dim2$Phi

##           PA1       PA2
## PA1 1.0000000 0.4435671
## PA2 0.4435671 1.0000000

F_fa_dependent_lst$Dim3$Phi

##           PA1       PA3       PA2
## PA1 1.0000000 0.4999271 0.3280693
## PA3 0.4999271 1.0000000 0.2791827
## PA2 0.3280693 0.2791827 1.0000000

• moderate positive correlations between the latent factors; supporting dependent structure

• Looks like 2 factor model is more convincing

9.1.2 independent factors

# Independent factors (Equamax); 1-3 factors
F_fa_independent_lst <- lapply(1:3, function(k) {
  fa(
    F_corr,
    nfactors = k,
    fm = "pa",
    rotate = if (k == 1) "none" else "equamax",
    n.obs = nrow(Fall_binary)
  )
})

## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0

names(F_fa_independent_lst) <- paste0("Dim", 1:3)

9.1.2.1 Factor Loadings

print(F_fa_independent_lst$Dim1$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1  
## Q1  0.326
## Q2  0.683
## Q3  0.643
## Q4  0.710
## Q5  0.645
## Q6  0.723
## Q7  0.471
## Q8  0.684
## Q9  0.496
## Q10      
## Q11 0.810
## Q12      
## Q13 0.535
## Q14 0.517
## Q15 0.555
## Q16 0.584
## Q17 0.560
## Q18 0.629
## Q19 0.589
## Q20 0.621
## Q21 0.568
## Q22 0.660
## Q23      
## Q24 0.454
## Q25 0.708
## Q26 0.627
## Q27 0.760
## Q28 0.698
## Q29 0.706
## Q30 0.504
## Q31 0.670
## Q32 0.621
## Q33 0.688
## Q34 0.475
## 
##                   PA1
## SS loadings    11.979
## Proportion Var  0.352

print(F_fa_independent_lst$Dim2$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA2   
## Q1               
## Q2   0.531  0.428
## Q3          0.646
## Q4          0.842
## Q5   0.466  0.446
## Q6   0.557  0.460
## Q7          0.858
## Q8          0.749
## Q9   0.375  0.323
## Q10  0.372       
## Q11  0.585  0.557
## Q12              
## Q13         0.479
## Q14  0.410  0.315
## Q15  0.476       
## Q16  0.394  0.433
## Q17  0.488       
## Q18  0.771       
## Q19  0.556       
## Q20  0.427  0.453
## Q21  0.680       
## Q22  0.367  0.584
## Q23              
## Q24  0.536       
## Q25  0.741       
## Q26  0.475  0.408
## Q27  0.479  0.605
## Q28  0.642  0.327
## Q29  0.647  0.334
## Q30  0.457       
## Q31  0.495  0.449
## Q32  0.421  0.460
## Q33  0.557  0.405
## Q34  0.519       
## 
##                  PA1   PA2
## SS loadings    7.625 6.339
## Proportion Var 0.224 0.186
## Cumulative Var 0.224 0.411

print(F_fa_independent_lst$Dim3$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA3    PA2   
## Q1                      
## Q2   0.607         0.402
## Q3          0.408  0.572
## Q4          0.428  0.767
## Q5   0.343  0.406  0.365
## Q6   0.586         0.419
## Q7                 0.895
## Q8                 0.706
## Q9          0.386       
## Q10  0.451              
## Q11  0.495  0.426  0.476
## Q12                     
## Q13  0.387         0.483
## Q14  0.422              
## Q15         0.469       
## Q16  0.564         0.450
## Q17  0.327  0.435       
## Q18  0.536  0.576       
## Q19  0.495  0.309       
## Q20  0.413         0.404
## Q21  0.584  0.361       
## Q22         0.392  0.507
## Q23         0.529       
## Q24  0.461              
## Q25  0.793              
## Q26  0.415  0.316  0.349
## Q27         0.783  0.471
## Q28  0.480  0.489       
## Q29         0.846       
## Q30  0.395              
## Q31  0.428  0.343  0.383
## Q32         0.470  0.366
## Q33  0.459  0.396  0.325
## Q34  0.345  0.423       
## 
##                  PA1   PA3   PA2
## SS loadings    5.485 5.084 4.986
## Proportion Var 0.161 0.150 0.147
## Cumulative Var 0.161 0.311 0.458

9.1.2.1.1 summary :

• One-factor model: Most questions load moderately to strongly onto a common factor, explaining approximately 35.2% of total item variance. Q10, Q12, and Q23 have weak loadings.

• Two-factor model: Reveals a potential secondary construct, particularly associated with Q7, Q4, Q8, and Q3.

• Three-factor model: Identifies more differentiated question groupings, with Q23, Q27, and Q29 strongly associated with an additional factor.

• Rotation comparison: Oblimin produces a clearer loading structure than equamax, suggesting that allowing correlated factors may be more appropriate

  • We will keep depednt version

9.1.2.2 Visualize Models

par(mfrow = c(1, 2))

fa.diagram(
  F_fa_dependent_lst$Dim1,
  main = "1-Factor Model",
  cut = 0.3
)

fa.diagram(
  F_fa_dependent_lst$Dim2,
  main = "2-Factor Model",
  cut = 0.3
)

par(mfrow = c(1, 1))

9.2 Winter

9.2.1 dependent (correlated) factors

# Dependent factors (Oblimin); 1-3 factors
W_fa_dependent_lst <- lapply(1:3, function(k) {
  fa(
    W_corr,
    nfactors = k,
    fm = "pa",
    rotate = if (k == 1) "none" else "oblimin",
    n.obs = nrow(Winter_binary)
  )
})

names(W_fa_dependent_lst) <- paste0("Dim", 1:3)

9.2.1.1 Factor Loadings

print(W_fa_dependent_lst$Dim1$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1  
## Q1       
## Q2  0.559
## Q3  0.535
## Q4  0.674
## Q5  0.555
## Q6       
## Q7  0.702
## Q8  0.660
## Q9  0.709
## Q10 0.595
## Q11 0.388
## Q12 0.553
## Q13 0.735
## Q14 0.305
## Q15      
## Q16 0.518
## Q17 0.330
## Q18 0.664
## Q19 0.672
## Q20 0.544
## Q21 0.629
## Q22 0.397
## Q23 0.624
## Q24      
## Q25 0.568
## Q26 0.493
## Q27 0.493
## Q28 0.685
## Q29 0.684
## Q30 0.662
## Q31 0.612
## Q32 0.622
## Q33 0.640
## 
##                   PA1
## SS loadings    10.407
## Proportion Var  0.315

print(W_fa_dependent_lst$Dim2$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA2   
## Q1               
## Q2   0.454       
## Q3          0.512
## Q4   0.532       
## Q5          0.514
## Q6               
## Q7   0.784       
## Q8   0.570       
## Q9   0.913       
## Q10  0.385       
## Q11  0.482       
## Q12  0.444       
## Q13  0.627       
## Q14         0.406
## Q15              
## Q16         0.582
## Q17              
## Q18         0.615
## Q19  0.503       
## Q20  0.418       
## Q21  0.454       
## Q22         0.379
## Q23  0.580       
## Q24         0.598
## Q25  0.327  0.334
## Q26  0.508       
## Q27  0.705       
## Q28  0.559       
## Q29  0.628       
## Q30  0.596       
## Q31  0.402       
## Q32  0.539       
## Q33  0.581       
## 
##                  PA1   PA2
## SS loadings    7.139 2.932
## Proportion Var 0.216 0.089
## Cumulative Var 0.216 0.305

print(W_fa_dependent_lst$Dim3$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA3    PA2   
## Q1                      
## Q2   0.661              
## Q3   0.474              
## Q4   0.700              
## Q5   0.606              
## Q6                      
## Q7          0.904       
## Q8          0.575       
## Q9          0.736       
## Q10  0.642              
## Q11         0.488       
## Q12  0.592              
## Q13  0.614              
## Q14                0.336
## Q15  0.337              
## Q16                0.746
## Q17                0.339
## Q18  0.470         0.433
## Q19  0.578              
## Q20  0.390              
## Q21  0.767              
## Q22  0.462              
## Q23  0.485              
## Q24                0.526
## Q25  0.478              
## Q26  0.336              
## Q27  0.378  0.346       
## Q28  0.497              
## Q29  0.639              
## Q30  0.523              
## Q31  0.506              
## Q32  0.462              
## Q33  0.358  0.334       
## 
##                  PA1   PA3   PA2
## SS loadings    6.685 2.763 1.728
## Proportion Var 0.203 0.084 0.052
## Cumulative Var 0.203 0.286 0.339

9.2.1.2 Factor correlations

W_fa_dependent_lst$Dim2$Phi

##           PA1       PA2
## PA1 1.0000000 0.5379647
## PA2 0.5379647 1.0000000

W_fa_dependent_lst$Dim3$Phi

##           PA1       PA3       PA2
## PA1 1.0000000 0.5652099 0.3430508
## PA3 0.5652099 1.0000000 0.2527580
## PA2 0.3430508 0.2527580 1.0000000

9.2.2 independent factors

# Independent factors (Equamax); 1-3 factors
W_fa_independent_lst <- lapply(1:3, function(k) {
  fa(
    W_corr,
    nfactors = k,
    fm = "pa",
    rotate = if (k == 1) "none" else "equamax",
    n.obs = nrow(Winter_binary)
  )
})

names(W_fa_independent_lst) <- paste0("Dim", 1:3)

9.2.2.1 Factor Loadings

print(W_fa_independent_lst$Dim1$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1  
## Q1       
## Q2  0.559
## Q3  0.535
## Q4  0.674
## Q5  0.555
## Q6       
## Q7  0.702
## Q8  0.660
## Q9  0.709
## Q10 0.595
## Q11 0.388
## Q12 0.553
## Q13 0.735
## Q14 0.305
## Q15      
## Q16 0.518
## Q17 0.330
## Q18 0.664
## Q19 0.672
## Q20 0.544
## Q21 0.629
## Q22 0.397
## Q23 0.624
## Q24      
## Q25 0.568
## Q26 0.493
## Q27 0.493
## Q28 0.685
## Q29 0.684
## Q30 0.662
## Q31 0.612
## Q32 0.622
## Q33 0.640
## 
##                   PA1
## SS loadings    10.407
## Proportion Var  0.315

print(W_fa_independent_lst$Dim2$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA2   
## Q1               
## Q2   0.488       
## Q3   0.305  0.527
## Q4   0.580  0.341
## Q5   0.324  0.534
## Q6               
## Q7   0.741       
## Q8   0.597       
## Q9   0.822       
## Q10  0.462  0.381
## Q11  0.437       
## Q12  0.480       
## Q13  0.660  0.329
## Q14         0.386
## Q15              
## Q16         0.576
## Q17              
## Q18  0.390  0.639
## Q19  0.563  0.366
## Q20  0.461       
## Q21  0.517  0.358
## Q22         0.390
## Q23  0.589       
## Q24         0.526
## Q25  0.418  0.402
## Q26  0.494       
## Q27  0.612       
## Q28  0.600  0.331
## Q29  0.641       
## Q30  0.613       
## Q31  0.479  0.386
## Q32  0.563       
## Q33  0.596       
## 
##                  PA1   PA2
## SS loadings    7.841 3.867
## Proportion Var 0.238 0.117
## Cumulative Var 0.238 0.355

print(W_fa_independent_lst$Dim3$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA3    PA2   
## Q1                      
## Q2   0.561              
## Q3   0.394         0.422
## Q4   0.605  0.322       
## Q5   0.491         0.424
## Q6                      
## Q7          0.888       
## Q8          0.635  0.380
## Q9          0.804       
## Q10  0.543              
## Q11         0.484       
## Q12  0.510              
## Q13  0.552  0.449       
## Q14                0.381
## Q15                     
## Q16         0.302  0.776
## Q17                0.371
## Q18  0.400         0.578
## Q19  0.510  0.360       
## Q20  0.355  0.337       
## Q21  0.645              
## Q22  0.375              
## Q23  0.445  0.420       
## Q24                0.544
## Q25  0.414         0.309
## Q26  0.321  0.381       
## Q27  0.371  0.469       
## Q28  0.455  0.440       
## Q29  0.570  0.418       
## Q30  0.477  0.433       
## Q31  0.444  0.308       
## Q32  0.425  0.415       
## Q33  0.349  0.484       
## 
##                  PA1   PA3   PA2
## SS loadings    5.181 4.802 3.051
## Proportion Var 0.157 0.146 0.092
## Cumulative Var 0.157 0.303 0.395

9.2.2.1.1 summary

• One-factor model: Most questions load moderately to strongly onto a common factor, explaining approximately 31.5% of total item variance. Q1, Q6, Q15, and Q24 exhibit weak loadings (below 0.30).

• Two-factor model: Identifies a possible secondary construct, particularly associated with Q3, Q5, Q16, Q18, and Q24. The factors are moderately correlated (r=0.538).

• Three-factor model: Reveals more differentiated question groupings, notably Q7/Q8/Q9 and Q16/Q18/Q24. Factor correlations range from 0.253 to 0.565.

• Rotation comparison: Oblimin permits correlated factors and produces relatively distinct loading patterns, while equamax imposes independence and exhibits more cross-loadings

Similar result to last model

9.2.2.2 Visualize Models

par(mfrow = c(1, 2))

fa.diagram(
  W_fa_dependent_lst$Dim1,
  main = "1-Factor Model",
  cut = 0.3
)


fa.diagram(
  W_fa_dependent_lst$Dim2,
  main = "2-Factor Model",
  cut = 0.3
)

par(mfrow = c(1, 1))

9.3 Spring

9.3.1 dependent (correlated) factors

# Dependent factors (Oblimin); 1-3 factors
S_fa_dependent_lst <- lapply(1:3, function(k) {
  fa(
    S_corr,
    nfactors = k,
    fm = "pa",
    rotate = if (k == 1) "none" else "oblimin",
    n.obs = nrow(Spring_binary)
  )
})

## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0

names(S_fa_dependent_lst) <- paste0("Dim", 1:3)

9.3.1.1 Factor Loadings

print(S_fa_dependent_lst$Dim1$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1   
## Q1        
## Q2   0.424
## Q3   0.509
## Q4   0.334
## Q5   0.463
## Q6   0.465
## Q7   0.357
## Q8   0.523
## Q9   0.535
## Q10  0.676
## Q11  0.408
## Q12  0.614
## Q13  0.473
## Q14       
## Q15       
## Q16  0.454
## Q17       
## Q18  0.614
## Q19  0.593
## Q20  0.736
## Q21  0.425
## Q22       
## Q23  0.503
## Q24       
## Q25  0.533
## Q26  0.348
## Q27  0.701
## Q28  0.362
## Q29  0.461
## Q30  0.418
## Q31  0.614
## Q32  0.356
## Q33  0.517
## Q34  0.464
## 
##                  PA1
## SS loadings    7.462
## Proportion Var 0.219

print(S_fa_dependent_lst$Dim2$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA2   
## Q1          0.554
## Q2          0.798
## Q3   0.321  0.343
## Q4               
## Q5   0.462       
## Q6   0.523       
## Q7               
## Q8   0.387       
## Q9   0.345  0.348
## Q10  0.372  0.547
## Q11  0.391       
## Q12  0.491       
## Q13  0.679       
## Q14              
## Q15              
## Q16  0.447       
## Q17         0.558
## Q18  0.722       
## Q19  0.691       
## Q20  0.823       
## Q21  0.431       
## Q22              
## Q23  0.513       
## Q24              
## Q25  0.628       
## Q26  0.439       
## Q27  0.606       
## Q28         0.331
## Q29  0.349       
## Q30         0.808
## Q31  0.475       
## Q32              
## Q33  0.437       
## Q34  0.345       
## 
##                  PA1   PA2
## SS loadings    6.053 3.338
## Proportion Var 0.178 0.098
## Cumulative Var 0.178 0.276

print(S_fa_dependent_lst$Dim3$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA2    PA3   
## Q1          0.540       
## Q2          0.782       
## Q3          0.331       
## Q4   0.373              
## Q5   0.341              
## Q6   0.543              
## Q7                 0.344
## Q8                 0.318
## Q9          0.332  0.558
## Q10  0.330  0.532       
## Q11                0.772
## Q12  0.440              
## Q13  0.678              
## Q14                     
## Q15                     
## Q16                0.410
## Q17         0.567       
## Q18  0.658              
## Q19  0.638              
## Q20  0.596         0.512
## Q21  0.422              
## Q22                     
## Q23  0.631              
## Q24                0.354
## Q25  0.609              
## Q26                0.356
## Q27  0.598              
## Q28         0.334       
## Q29  0.411              
## Q30         0.825       
## Q31  0.446              
## Q32  0.395        -0.320
## Q33  0.616              
## Q34  0.408              
## 
##                  PA1   PA2   PA3
## SS loadings    5.370 3.291 2.470
## Proportion Var 0.158 0.097 0.073
## Cumulative Var 0.158 0.255 0.327

9.3.1.2 Factor correlations

S_fa_dependent_lst$Dim2$Phi

##           PA1       PA2
## PA1 1.0000000 0.2822044
## PA2 0.2822044 1.0000000

S_fa_dependent_lst$Dim3$Phi

##           PA1       PA2       PA3
## PA1 1.0000000 0.2496154 0.2747172
## PA2 0.2496154 1.0000000 0.1133091
## PA3 0.2747172 0.1133091 1.0000000

9.3.2 independent factors

# Independent factors (Equamax); 1-3 factors
S_fa_independent_lst <- lapply(1:3, function(k) {
  fa(
    S_corr,
    nfactors = k,
    fm = "pa",
    rotate = if (k == 1) "none" else "equamax",
    n.obs = nrow(Spring_binary)
  )
})

## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0
## In smc, smcs < 0 were set to .0

names(S_fa_independent_lst) <- paste0("Dim", 1:3)

9.3.2.1 Factor Loadings

print(S_fa_independent_lst$Dim1$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1   
## Q1        
## Q2   0.424
## Q3   0.509
## Q4   0.334
## Q5   0.463
## Q6   0.465
## Q7   0.357
## Q8   0.523
## Q9   0.535
## Q10  0.676
## Q11  0.408
## Q12  0.614
## Q13  0.473
## Q14       
## Q15       
## Q16  0.454
## Q17       
## Q18  0.614
## Q19  0.593
## Q20  0.736
## Q21  0.425
## Q22       
## Q23  0.503
## Q24       
## Q25  0.533
## Q26  0.348
## Q27  0.701
## Q28  0.362
## Q29  0.461
## Q30  0.418
## Q31  0.614
## Q32  0.356
## Q33  0.517
## Q34  0.464
## 
##                  PA1
## SS loadings    7.462
## Proportion Var 0.219

print(S_fa_independent_lst$Dim2$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA2   
## Q1          0.541
## Q2          0.785
## Q3   0.354  0.397
## Q4               
## Q5   0.461       
## Q6   0.511       
## Q7          0.304
## Q8   0.410  0.329
## Q9   0.378  0.406
## Q10  0.427  0.609
## Q11  0.393       
## Q12  0.511  0.334
## Q13  0.642       
## Q14              
## Q15              
## Q16  0.446       
## Q17         0.542
## Q18  0.701       
## Q19  0.672       
## Q20  0.807       
## Q21  0.429       
## Q22              
## Q23  0.510       
## Q24              
## Q25  0.610       
## Q26  0.421       
## Q27  0.621  0.323
## Q28         0.359
## Q29  0.368       
## Q30         0.792
## Q31  0.498  0.356
## Q32              
## Q33  0.449       
## Q34  0.365       
## 
##                  PA1   PA2
## SS loadings    6.106 3.796
## Proportion Var 0.180 0.112
## Cumulative Var 0.180 0.291

print(S_fa_independent_lst$Dim3$loadings, cutoff = 0.3)

## 
## Loadings:
##     PA1    PA2    PA3   
## Q1          0.520       
## Q2          0.759       
## Q3          0.370  0.301
## Q4   0.366              
## Q5   0.309         0.378
## Q6   0.507              
## Q7                 0.383
## Q8                 0.420
## Q9          0.353  0.616
## Q10  0.338  0.582  0.305
## Q11                0.784
## Q12  0.421  0.308  0.314
## Q13  0.617              
## Q14                     
## Q15                     
## Q16                0.494
## Q17         0.566       
## Q18  0.603         0.353
## Q19  0.586         0.329
## Q20  0.529         0.693
## Q21  0.395              
## Q22                     
## Q23  0.602              
## Q24                     
## Q25  0.561              
## Q26                0.429
## Q27  0.571  0.312       
## Q28         0.365       
## Q29  0.403              
## Q30         0.823       
## Q31  0.431  0.338       
## Q32  0.404  0.343       
## Q33  0.603              
## Q34  0.401              
## 
##                  PA1   PA2   PA3
## SS loadings    4.771 3.637 3.510
## Proportion Var 0.140 0.107 0.103
## Cumulative Var 0.140 0.247 0.351

9.3.2.1.1 summary :

• One-factor model: Explains approximately 21.9% of item variance, with Q20, Q27, and Q10 showing the strongest loadings. Six questions have loadings below 0.30.

• Two-factor model: Identifies a potential secondary dimension associated with Q30, Q2, Q17, Q1, and Q10. The factors are weakly correlated (r=0.282).

• Three-factor model: Reveals additional question groupings, particularly Q13/Q18/Q19 and Q11/Q9/Q20. Factor correlations range from 0.113 to 0.275.

• Rotation comparison: Oblimin and equamax produce broadly similar structures, with weaker factor correlations making orthogonal rotation more defensible than in Fall or Winter.

9.3.2.2 Visualize Models

par(mfrow = c(1, 3))

fa.diagram(S_fa_dependent_lst$Dim1, main = "1 Factor", cut = 0.3)
fa.diagram(S_fa_dependent_lst$Dim2, main = "2 Factors", cut = 0.3)
fa.diagram(S_fa_dependent_lst$Dim3, main = "3 Factors", cut = 0.3)

par(mfrow = c(1, 1))

10 Exploratory IRT :

• Here we continue from yesterdays work based upon the findings.

• The purpose is to :

  • Fit Model(s)

  • Understand fit-indices

  • Understand Information Criteria

  • Compare Models to choose “best”

  • and, inspect candidate model closely.

Aside :

• We will fit 1PL, 2PL : 1 Dim IRT and,

• 2PL : 1, 2 m IRT

  • Correlated Structure

  • Uncorrelated Structure

grep("binary", ls(), value=T)

## [1] "Fall_binary"   "Spring_binary" "Winter_binary"

library(mirt)

## Warning: package 'mirt' was built under R version 4.4.3

## Loading required package: stats4

## Loading required package: lattice

## Registered S3 method overwritten by 'vegan':
##   method     from      
##   rev.hclust dendextend

?mirt

10.1 Fall

library(mirt)

library(mirt)

# Rasch: 1 Dimension
F_Rasch <- mirt(Fall_binary, 1, itemtype = "Rasch")
extract.mirt(F_Rasch, "converged") # good!

## [1] TRUE

# 2PL: 1 Dimension
F_Dim1 <- mirt(Fall_binary, 1, itemtype = "2PL")
extract.mirt(F_Dim1, "converged") # good!

## [1] TRUE

# 2PL: 2 Dimensions : Uncorrelated (orthogonal)
F_Dim2 <- 
  mirt(
  Fall_binary, 2,
  itemtype = "2PL"
)

## Warning: EM cycles terminated after 500 iterations.

extract.mirt(F_Dim2, "converged") # BAD!

## [1] FALSE

• notice our model didnt converge, maybe lets allow it more iterations : Expectation-Maximization (EM-Alg)

F_Dim2 <- 
  mirt(
  Fall_binary, 2,
  itemtype = "2PL",
  technical = list(NCYCLES = 10000)
)

extract.mirt(F_Dim2, "converged") # nice!

## [1] TRUE

• here we see it did converge, but after a quite substantial amount of iterations 6757

• Further notice our Log-Lik for our simpler model (F_Dim1) is -6653.186 ; and F_Dim2 is -6587.884 so, our model is only doing slightly better despite taking way more iterations and being a more complex model.

  • Suspect

-6587.884 - -6653.186 

## [1] 65.302

• very small improvement

• this isnt proof of dimensionality but rather something to think about

# Store fitted models
Fall_models <- list(
  Rasch = F_Rasch,
  Dim1 = F_Dim1,
  Dim2 = F_Dim2
)

10.1.1 Examine Rotations

• ie, is it a correlated or uncorrelated structure?

• Notice we use the same rotations as in EFA

# Orthogonal (uncorrelated)
summary(Fall_models$Dim2, rotate = "varimax")

## 
## Rotation:  varimax 
## 
## Rotated factor loadings: 
## 
##        F1     F2    h2
## Q1  0.377  0.105 0.153
## Q2  0.629  0.409 0.563
## Q3  0.155  0.810 0.680
## Q4  0.174  0.969 0.969
## Q5  0.486  0.465 0.453
## Q6  0.635  0.584 0.744
## Q7  0.026  0.881 0.777
## Q8  0.464  0.635 0.618
## Q9  0.289  0.358 0.212
## Q10 0.597 -0.349 0.478
## Q11 0.579  0.636 0.740
## Q12 0.041  0.105 0.013
## Q13 0.223  0.560 0.364
## Q14 0.444  0.299 0.287
## Q15 0.322  0.512 0.366
## Q16 0.527  0.467 0.496
## Q17 0.526  0.379 0.420
## Q18 0.694  0.270 0.555
## Q19 0.508  0.400 0.418
## Q20 0.480  0.390 0.382
## Q21 0.699  0.170 0.518
## Q22 0.464  0.517 0.482
## Q23 0.146  0.206 0.064
## Q24 0.493  0.226 0.294
## Q25 0.859  0.216 0.784
## Q26 0.432  0.427 0.368
## Q27 0.399  0.762 0.741
## Q28 0.579  0.326 0.441
## Q29 0.573  0.588 0.674
## Q30 0.402  0.372 0.300
## Q31 0.538  0.504 0.544
## Q32 0.434  0.459 0.400
## Q33 0.586  0.484 0.578
## Q34 0.375  0.289 0.224
## 
## Rotated SS loadings:  7.938 8.16 
## 
## Factor correlations: 
## 
##    F1 F2
## F1  1   
## F2  0  1

• here H represents the Communality; or pct of response variation explained

• Notice we see the correlation structure is the \(I_2\) matrix as we asked for

# Oblique (correlated)
summary(Fall_models$Dim2, rotate = "oblimin")

## 
## Rotation:  oblimin 
## 
## Rotated factor loadings: 
## 
##        F1     F2    h2
## Q1  0.423 -0.095 0.153
## Q2  0.718  0.068 0.563
## Q3  0.216  0.708 0.680
## Q4  0.246  0.853 0.969
## Q5  0.564  0.198 0.453
## Q6  0.734  0.237 0.744
## Q7  0.078  0.845 0.777
## Q8  0.549  0.375 0.618
## Q9  0.340  0.197 0.212
## Q10 0.640 -0.653 0.478
## Q11 0.676  0.316 0.740
## Q12 0.051  0.081 0.013
## Q13 0.278  0.429 0.364
## Q14 0.508  0.058 0.287
## Q15 0.385  0.330 0.366
## Q16 0.609  0.178 0.496
## Q17 0.602  0.094 0.420
## Q18 0.783 -0.101 0.555
## Q19 0.584  0.123 0.418
## Q20 0.552  0.128 0.382
## Q21 0.783 -0.202 0.518
## Q22 0.542  0.260 0.482
## Q23 0.173  0.124 0.064
## Q24 0.557 -0.038 0.294
## Q25 0.961 -0.240 0.784
## Q26 0.501  0.189 0.368
## Q27 0.484  0.533 0.741
## Q28 0.658  0.014 0.441
## Q29 0.666  0.272 0.674
## Q30 0.465  0.151 0.300
## Q31 0.623  0.209 0.544
## Q32 0.506  0.220 0.400
## Q33 0.675  0.164 0.578
## Q34 0.430  0.085 0.224
## 
## Rotated SS loadings:  10.463 3.87 
## 
## Factor correlations: 
## 
##       F1 F2
## F1 1.000   
## F2 0.429  1

• Notice we get a similar correlation between factors 0.4435671 is what we got from efa and 0.429 here

• Moderate suggestion for 2 factors. Lets further compare loadings

Aside

• positive factor loadings : Students with high factor loadings on this factor tend to get the item right

• negative factor loadings : Students with high factor (magnitude) loadings on this factor tend to get the item wrong

• Think of communality as how well the factors help explain response variation given those 2 factors

10.1.2 Information Criteria

• note that both rotations of 2 dim model have same information criteria, so this isnt how we can compare them. This is more of a dimensionality check

10.1.2.1 Does having discrimination help?

anova(
  Fall_models$Rasch,
  Fall_models$Dim1
)

• All four are information criteria that combine model fit with a penalty for additional parameters. We have some evidence for the descrimination parameter.

10.1.2.2 Does having dimensionality help?

anova(
  Fall_models$Dim1,
  Fall_models$Dim2
)

• Best-fitting under AIC, SABIC, HQ: 2D 2PL.

• Best-fitting under BIC: 1D 2PL.

• Therefore we have conflicting evidence and need to run further test(s).

10.1.2.3 absolute goodness-of-fit

• here we test if the model adequately reproduces the population response statistics tested.

# Limited-information model fit
M2(F_Dim1, type = "M2*")

M2(F_Dim2, type = "M2*")

• BIC favors 1D model

• Other measures favor 2D model

Lets assume the 2D model is favored, suppose we wish to further inspect this model itself.

11 Plotting 2D surface

11.1 Uncorrelated Structure

plot(F_Dim2, rotate = "varimax")

11.2 Correlated Structure

plot(F_Dim2, rotate = "oblimin")

11.3 Examine Coef

• here we assume uncorrelated structure.

coef(F_Dim2, simplify = TRUE)

## $items
##         a1     a2      d g u
## Q1  -0.671  0.272  2.043 0 1
## Q2  -1.925  0.156  3.120 0 1
## Q3  -1.856 -1.644  0.637 0 1
## Q4  -7.032 -6.365 18.699 0 1
## Q5  -1.541 -0.164  1.687 0 1
## Q6  -2.892 -0.252  6.190 0 1
## Q7  -2.014 -2.456  7.518 0 1
## Q8  -2.081 -0.603  3.959 0 1
## Q9  -0.859 -0.205  0.882 0 1
## Q10 -0.611  1.508  2.053 0 1
## Q11 -2.828 -0.501  3.701 0 1
## Q12 -0.165 -0.100  0.717 0 1
## Q13 -1.107 -0.656  0.535 0 1
## Q14 -1.077  0.070  0.457 0 1
## Q15 -1.214 -0.445  1.767 0 1
## Q16 -1.685 -0.114  3.601 0 1
## Q17 -1.448  0.047  3.265 0 1
## Q18 -1.825  0.536  3.879 0 1
## Q19 -1.443 -0.015  2.049 0 1
## Q20 -1.338 -0.034  1.324 0 1
## Q21 -1.612  0.718  4.070 0 1
## Q22 -1.615 -0.297  2.143 0 1
## Q23 -0.424 -0.130  0.875 0 1
## Q24 -1.070  0.247  1.724 0 1
## Q25 -2.971  1.294  4.709 0 1
## Q26 -1.290 -0.158  1.172 0 1
## Q27 -2.614 -1.202  4.867 0 1
## Q28 -1.497  0.217  0.444 0 1
## Q29 -2.422 -0.344  4.990 0 1
## Q30 -1.109 -0.099  1.806 0 1
## Q31 -1.850 -0.177  2.631 0 1
## Q32 -1.371 -0.216  1.411 0 1
## Q33 -1.991 -0.065  2.530 0 1
## Q34 -0.914  0.000 -0.040 0 1
## 
## $means
## F1 F2 
##  0  0 
## 
## $cov
##    F1 F2
## F1  1  0
## F2  0  1

11.4 Item Characteristic Curves (ICCs)

• here they are surfaces given we are in 2D

• Tells us how the prob of answering correctly changes at different abilities

plot(F_Dim2, type = "trace")

11.5 Information curves

• tell us the precision of measuring underlying train

• item wise : Shows how precisely an individual question measures proficiency at different ability levels. Identifies which students a question measures most effectively.

• exam wise : Shows how precisely the entire exam measures proficiency at different ability levels. Combines information across questions.

Generally, higher information means higher precision and vice versa

11.5.1 item 1

# Item information surface : Uncor()
itemplot(F_Dim2, item = 1, type = "info")

itemplot(F_Dim2, item = 1, type = "info", rotate="oblimin")

11.5.2 item 2

# Item information surface : Uncor()
itemplot(F_Dim2, item = 2, type = "info")

# correlated str 
itemplot(F_Dim2, item = 2, type = "info", rotate="oblimin")

11.5.3 Iterate across exam to generate plts :

make_item_plots <- function(model, n_items = 34) {
  list(
    Uncorrelated = lapply(seq_len(n_items), function(i) {
      local({
        item <- i
        function() itemplot(model, item = item, type = "info")
      })
    }),
    Correlated = lapply(seq_len(n_items), function(i) {
      local({
        item <- i
        function() itemplot(
          model, item = item,
          type = "info", rotate = "oblimin"
        )
      })
    })
  )
}

F_info_plots <- make_item_plots(F_Dim2)

# Test function
F_info_plots$Uncorrelated[[1]]()

F_info_plots$Correlated[[1]]()

11.6 Overall test information surface

# Independent factros
plot(F_Dim2, type = "info")

# corr() factros
plot(F_Dim2, type = "info", rotate = "oblimin")

11.7 Winter

12 Resources :

• Article :

  • Assessing Dimensionality of IRT Models Using Traditional and Revised Parallel Analyses
  • CAN WE DISTINGUISH STATISTICAL LITERACY AND STATISTICAL REASONING?
    • Used for

• YouTube :

  • Tetrachoric Correlation
  • StatQuest: Hierarchical Clustering
  • Flat and Hierarchical Clustering | The Dendrogram Explained

• Wiki :

  • Tetrachoric correlation
  • Dendrogram

• Forums / Blogs

  • Hierarchical clustering of correlation matrix
  • Exploratory Factor Analysis

• Textbooks

  • Biology 723: Statistical Computing for Biologists

13 Notes :

13.1 Assumption of Tetrachoric Corr() :

• Underlying variables come from a normal dist.

• there is a latent continuous scale

• errors are independent from raters and cases

• variance of errors is homogenous across levels of independent variable

In other words :

• although we only observe correct(1)/incorrect(0); we believe the constructs are normally distributed.

13.2 What does Tetrachoric Corr() measure?

• association between items

• does item_i and item_j measure the same construct?

13.3 Interpretation :

• when its 0, we suspect may measure separate underlying construct

• When positive, we believe they answer those correctly together

• When negative, we believe when they answer item_i correct, they get item_j wrong

• Absolute value of correlation denotes the strength of the relationship between items

13.4 Correlation clustering : hclust()

• here we use hclust , re-ordering according to clustering

• our notion of similarity comes from as.dist(1 - *_corr)

  • 1 - *_corr is our dissimilarity
    • So for ex. a correlation of .8 would be .2 dissimilarity
    • So a score of 0 would be similar

  print("correlation")

  ## [1] "correlation"

  S_corr[1:3, 1:3]

  ##           Q1        Q2        Q3
  ## Q1 1.0000000 0.4971214 0.1977192
  ## Q2 0.4971214 1.0000000 0.5150938
  ## Q3 0.1977192 0.5150938 1.0000000

  print("dissimilarity")

  ## [1] "dissimilarity"

  1 - S_corr[1:3, 1:3]

  ##           Q1        Q2        Q3
  ## Q1 0.0000000 0.5028786 0.8022808
  ## Q2 0.5028786 0.0000000 0.4849062
  ## Q3 0.8022808 0.4849062 0.0000000

  • as.dist() is just keeping off-diag [i,j entries for \(i\ne j\)]

as.dist(1 - S_corr[1:3, 1:3])

##           Q1        Q2
## Q2 0.5028786          
## Q3 0.8022808 0.4849062

13.4.1 Dendrogram

• Represent similarity and order of hiearchy

  • similarity is arbitrarily chosen

• height of branch defines order

• the ends of clusters is object broken down. We start from broad to specific groups

  • think animals –> mammals, non-mammals

13.5 IRT

• slope of sigmoid is discrimination (assumed constant 1 for Rasch mdl)

• intercept is difficulty