data("bfi")
dta <- bfi

Load data

str(dta)
## 'data.frame':    2800 obs. of  28 variables:
##  $ A1       : int  2 2 5 4 2 6 2 4 4 2 ...
##  $ A2       : int  4 4 4 4 3 6 5 3 3 5 ...
##  $ A3       : int  3 5 5 6 3 5 5 1 6 6 ...
##  $ A4       : int  4 2 4 5 4 6 3 5 3 6 ...
##  $ A5       : int  4 5 4 5 5 5 5 1 3 5 ...
##  $ C1       : int  2 5 4 4 4 6 5 3 6 6 ...
##  $ C2       : int  3 4 5 4 4 6 4 2 6 5 ...
##  $ C3       : int  3 4 4 3 5 6 4 4 3 6 ...
##  $ C4       : int  4 3 2 5 3 1 2 2 4 2 ...
##  $ C5       : int  4 4 5 5 2 3 3 4 5 1 ...
##  $ E1       : int  3 1 2 5 2 2 4 3 5 2 ...
##  $ E2       : int  3 1 4 3 2 1 3 6 3 2 ...
##  $ E3       : int  3 6 4 4 5 6 4 4 NA 4 ...
##  $ E4       : int  4 4 4 4 4 5 5 2 4 5 ...
##  $ E5       : int  4 3 5 4 5 6 5 1 3 5 ...
##  $ N1       : int  3 3 4 2 2 3 1 6 5 5 ...
##  $ N2       : int  4 3 5 5 3 5 2 3 5 5 ...
##  $ N3       : int  2 3 4 2 4 2 2 2 2 5 ...
##  $ N4       : int  2 5 2 4 4 2 1 6 3 2 ...
##  $ N5       : int  3 5 3 1 3 3 1 4 3 4 ...
##  $ O1       : int  3 4 4 3 3 4 5 3 6 5 ...
##  $ O2       : int  6 2 2 3 3 3 2 2 6 1 ...
##  $ O3       : int  3 4 5 4 4 5 5 4 6 5 ...
##  $ O4       : int  4 3 5 3 3 6 6 5 6 5 ...
##  $ O5       : int  3 3 2 5 3 1 1 3 1 2 ...
##  $ gender   : int  1 2 2 2 1 2 1 1 1 2 ...
##  $ education: int  NA NA NA NA NA 3 NA 2 1 NA ...
##  $ age      : int  16 18 17 17 17 21 18 19 19 17 ...
summary(dta)
##        A1              A2              A3              A4            A5      
##  Min.   :1.000   Min.   :1.000   Min.   :1.000   Min.   :1.0   Min.   :1.00  
##  1st Qu.:1.000   1st Qu.:4.000   1st Qu.:4.000   1st Qu.:4.0   1st Qu.:4.00  
##  Median :2.000   Median :5.000   Median :5.000   Median :5.0   Median :5.00  
##  Mean   :2.413   Mean   :4.802   Mean   :4.604   Mean   :4.7   Mean   :4.56  
##  3rd Qu.:3.000   3rd Qu.:6.000   3rd Qu.:6.000   3rd Qu.:6.0   3rd Qu.:5.00  
##  Max.   :6.000   Max.   :6.000   Max.   :6.000   Max.   :6.0   Max.   :6.00  
##  NA's   :16      NA's   :27      NA's   :26      NA's   :19    NA's   :16    
##        C1              C2             C3              C4              C5       
##  Min.   :1.000   Min.   :1.00   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:4.000   1st Qu.:4.00   1st Qu.:4.000   1st Qu.:1.000   1st Qu.:2.000  
##  Median :5.000   Median :5.00   Median :5.000   Median :2.000   Median :3.000  
##  Mean   :4.502   Mean   :4.37   Mean   :4.304   Mean   :2.553   Mean   :3.297  
##  3rd Qu.:5.000   3rd Qu.:5.00   3rd Qu.:5.000   3rd Qu.:4.000   3rd Qu.:5.000  
##  Max.   :6.000   Max.   :6.00   Max.   :6.000   Max.   :6.000   Max.   :6.000  
##  NA's   :21      NA's   :24     NA's   :20      NA's   :26      NA's   :16     
##        E1              E2              E3              E4       
##  Min.   :1.000   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:2.000   1st Qu.:2.000   1st Qu.:3.000   1st Qu.:4.000  
##  Median :3.000   Median :3.000   Median :4.000   Median :5.000  
##  Mean   :2.974   Mean   :3.142   Mean   :4.001   Mean   :4.422  
##  3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:5.000   3rd Qu.:6.000  
##  Max.   :6.000   Max.   :6.000   Max.   :6.000   Max.   :6.000  
##  NA's   :23      NA's   :16      NA's   :25      NA's   :9      
##        E5              N1              N2              N3       
##  Min.   :1.000   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:4.000   1st Qu.:2.000   1st Qu.:2.000   1st Qu.:2.000  
##  Median :5.000   Median :3.000   Median :4.000   Median :3.000  
##  Mean   :4.416   Mean   :2.929   Mean   :3.508   Mean   :3.217  
##  3rd Qu.:5.000   3rd Qu.:4.000   3rd Qu.:5.000   3rd Qu.:4.000  
##  Max.   :6.000   Max.   :6.000   Max.   :6.000   Max.   :6.000  
##  NA's   :21      NA's   :22      NA's   :21      NA's   :11     
##        N4              N5             O1              O2              O3       
##  Min.   :1.000   Min.   :1.00   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:2.000   1st Qu.:2.00   1st Qu.:4.000   1st Qu.:1.000   1st Qu.:4.000  
##  Median :3.000   Median :3.00   Median :5.000   Median :2.000   Median :5.000  
##  Mean   :3.186   Mean   :2.97   Mean   :4.816   Mean   :2.713   Mean   :4.438  
##  3rd Qu.:4.000   3rd Qu.:4.00   3rd Qu.:6.000   3rd Qu.:4.000   3rd Qu.:5.000  
##  Max.   :6.000   Max.   :6.00   Max.   :6.000   Max.   :6.000   Max.   :6.000  
##  NA's   :36      NA's   :29     NA's   :22                      NA's   :28     
##        O4              O5           gender        education         age       
##  Min.   :1.000   Min.   :1.00   Min.   :1.000   Min.   :1.00   Min.   : 3.00  
##  1st Qu.:4.000   1st Qu.:1.00   1st Qu.:1.000   1st Qu.:3.00   1st Qu.:20.00  
##  Median :5.000   Median :2.00   Median :2.000   Median :3.00   Median :26.00  
##  Mean   :4.892   Mean   :2.49   Mean   :1.672   Mean   :3.19   Mean   :28.78  
##  3rd Qu.:6.000   3rd Qu.:3.00   3rd Qu.:2.000   3rd Qu.:4.00   3rd Qu.:35.00  
##  Max.   :6.000   Max.   :6.00   Max.   :2.000   Max.   :5.00   Max.   :86.00  
##  NA's   :14      NA's   :20                     NA's   :223
head(dta)
##       A1 A2 A3 A4 A5 C1 C2 C3 C4 C5 E1 E2 E3 E4 E5 N1 N2 N3 N4 N5 O1 O2 O3 O4
## 61617  2  4  3  4  4  2  3  3  4  4  3  3  3  4  4  3  4  2  2  3  3  6  3  4
## 61618  2  4  5  2  5  5  4  4  3  4  1  1  6  4  3  3  3  3  5  5  4  2  4  3
## 61620  5  4  5  4  4  4  5  4  2  5  2  4  4  4  5  4  5  4  2  3  4  2  5  5
## 61621  4  4  6  5  5  4  4  3  5  5  5  3  4  4  4  2  5  2  4  1  3  3  4  3
## 61622  2  3  3  4  5  4  4  5  3  2  2  2  5  4  5  2  3  4  4  3  3  3  4  3
## 61623  6  6  5  6  5  6  6  6  1  3  2  1  6  5  6  3  5  2  2  3  4  3  5  6
##       O5 gender education age
## 61617  3      1        NA  16
## 61618  3      2        NA  18
## 61620  2      2        NA  17
## 61621  5      2        NA  17
## 61622  3      1        NA  17
## 61623  1      2         3  21
dta <- dta[,1:25]

Remove useless data

dta$A1 <- 7 - dta$A1
dta$C4 <- 7 - dta$C4
dta$C5 <- 7 - dta$C5
dta$E1 <- 7 - dta$E1
dta$E2 <- 7 - dta$E2
dta$O2 <- 7 - dta$O2
dta$O5 <- 7 - dta$O5

reverse certain questions

my_summary <- function(x) {
 require(moments)
 funs <- c(mean, sd, skewness, kurtosis)
 sapply(funs, function(f) f(x, na.rm = T))
 }
dta_desc <- apply(dta, 2, my_summary)
head(dta_desc)
##              A1        A2        A3        A4        A5         C1         C2
## [1,]  4.5865661  4.802380  4.603821  4.699748  4.560345  4.5023390  4.3699568
## [2,]  1.4077372  1.172020  1.301834  1.479633  1.258512  1.2413465  1.3183465
## [3,] -0.8254883 -1.124894 -0.998997 -1.031499 -0.847690 -0.8551631 -0.7422207
## [4,]  2.6942957  4.057765  3.444524  3.042640  3.161176  3.3068088  2.8656243
##              C3         C4          C5         E1         E2         E3
## [1,]  4.3039568  4.4466474  3.70330460  4.0255672  3.8581178  4.0007207
## [2,]  1.2885518  1.3751181  1.62854187  1.6315055  1.6052103  1.3527188
## [3,] -0.6918287 -0.5964955 -0.06620282 -0.3736569 -0.2209396 -0.4706335
## [4,]  2.8697332  2.3802970  1.78461246  1.9090390  1.8526925  2.5367154
##              E4        E5        N1          N2        N3        N4        N5
## [1,]  4.4224292  4.416337 2.9290857  3.50773660 3.2165651 3.1856006 2.9696860
## [2,]  1.4575174  1.334768 1.5709175  1.52594359 1.6029021 1.5696851 1.6186474
## [3,] -0.8241831 -0.777486 0.3714298 -0.07698521 0.1506797 0.1969966 0.3744599
## [4,]  2.6977079  2.908401 1.9885722  1.95035250 1.8227046 1.9090309 1.9401121
##              O1        O2         O3        O4         O5
## [1,]  4.8160547  4.286786  4.4383117  4.892319  4.5104317
## [2,]  1.1295303  1.565152  1.2209011  1.221250  1.3279590
## [3,] -0.8973669 -0.585679 -0.7730516 -1.218247 -0.7384818
## [4,]  3.4277033  2.188889  3.3043641  4.082686  2.7630094
rownames(dta_desc) <- c("mean", "sd", "skewness", "kurtosis")
rslt1 <- as.data.frame(t(dta_desc))
rslt1 |> knitr::kable()
mean sd skewness kurtosis
A1 4.586566 1.407737 -0.8254883 2.694296
A2 4.802380 1.172020 -1.1248938 4.057765
A3 4.603821 1.301834 -0.9989970 3.444524
A4 4.699748 1.479633 -1.0314991 3.042640
A5 4.560345 1.258512 -0.8476900 3.161176
C1 4.502339 1.241346 -0.8551631 3.306809
C2 4.369957 1.318347 -0.7422207 2.865624
C3 4.303957 1.288552 -0.6918287 2.869733
C4 4.446647 1.375118 -0.5964955 2.380297
C5 3.703305 1.628542 -0.0662028 1.784612
E1 4.025567 1.631506 -0.3736569 1.909039
E2 3.858118 1.605210 -0.2209396 1.852693
E3 4.000721 1.352719 -0.4706335 2.536715
E4 4.422429 1.457517 -0.8241831 2.697708
E5 4.416337 1.334768 -0.7774860 2.908401
N1 2.929086 1.570917 0.3714298 1.988572
N2 3.507737 1.525944 -0.0769852 1.950352
N3 3.216565 1.602902 0.1506797 1.822705
N4 3.185601 1.569685 0.1969966 1.909031
N5 2.969686 1.618647 0.3744599 1.940112
O1 4.816055 1.129530 -0.8973669 3.427703
O2 4.286786 1.565152 -0.5856790 2.188889
O3 4.438312 1.220901 -0.7730516 3.304364
O4 4.892319 1.221250 -1.2182471 4.082686
O5 4.510432 1.327959 -0.7384818 2.763009 
dtal_desc <- melt(dta_desc)
names(dtal_desc)[1:2] <- c("moments", "items")
head(dtal_desc)
##    moments items      value
## 1     mean    A1  4.5865661
## 2       sd    A1  1.4077372
## 3 skewness    A1 -0.8254883
## 4 kurtosis    A1  2.6942957
## 5     mean    A2  4.8023801
## 6       sd    A2  1.1720199
ggplot(data = subset(dtal_desc, moments == "mean"),
 aes(x = reorder(items, value, max), y = value, group = moments)) +
 geom_point(size = 3) +
 geom_hline(yintercept = mean(t(dta_desc["mean",])) +
 c(-1.5, 0, 1.5) * sd(t(dta_desc["mean", ])), linetype = "dashed") +
 coord_flip() +
 labs(x = "items",  y = "mean") +
  theme_bw()

Mean of N3, N4, N5, N1 is lower than total mean

ggplot(data = subset(dtal_desc, moments == "sd"),
 aes(x = reorder(items, value, max), y = value, group = moments)) +
 geom_point(size = 3) +
 geom_hline(yintercept = mean(t(dta_desc["sd",])) +
 c(-1.5, 0, 1.5) * sd(t(dta_desc["sd", ])), linetype = "dashed") +
 coord_flip() +
 labs(x = "items",  y = "sd") +
  theme_bw()

O1 has way lower sd than total sd

ggplot(data = subset(dtal_desc, moments == "skewness"),
 aes(x = reorder(items, value, max), y = value, group = moments)) +
 geom_point(size = 3) +
 geom_hline(yintercept = mean(t(dta_desc["skewness",])) +
 c(-1.5, 0, 1.5) * sd(t(dta_desc["skewness", ])), linetype = "dashed") +
 coord_flip() +
 labs(x = "items",  y = "skewness") +
  theme_bw()

Skewness of N5, N1, N4 is higher than skewness of all

ggplot(data = subset(dtal_desc, moments == "kurtosis"),
 aes(x = reorder(items, value, max), y = value, group = moments)) +
 geom_point(size = 3) +
 geom_hline(yintercept = mean(t(dta_desc["kurtosis",])) +
 c(-1.5, 0, 1.5) * sd(t(dta_desc["kurtosis", ])), linetype = "dashed") +
 coord_flip() +
 labs(x = "items",  y = "kurtosis") +
  theme_bw()

Kurtosis of O4 and A2 is way higher than Kurtosis of all

dta <- na.omit(dta)
dta$tot <- apply(dta, 1, sum)
dta$grp <- NA
dta$grp[rank(dta$tot) < 2800*.33] <- "L"
dta$grp[rank(dta$tot) > 2800*.66] <- "H"
dta$grp <- factor(dta$grp)
head(dta)
##       A1 A2 A3 A4 A5 C1 C2 C3 C4 C5 E1 E2 E3 E4 E5 N1 N2 N3 N4 N5 O1 O2 O3 O4
## 61617  5  4  3  4  4  2  3  3  3  3  4  4  3  4  4  3  4  2  2  3  3  1  3  4
## 61618  5  4  5  2  5  5  4  4  4  3  6  6  6  4  3  3  3  3  5  5  4  5  4  3
## 61620  2  4  5  4  4  4  5  4  5  2  5  3  4  4  5  4  5  4  2  3  4  5  5  5
## 61621  3  4  6  5  5  4  4  3  2  2  2  4  4  4  4  2  5  2  4  1  3  4  4  3
## 61622  5  3  3  4  5  4  4  5  4  5  5  5  5  4  5  2  3  4  4  3  3  4  4  3
## 61623  1  6  5  6  5  6  6  6  6  4  5  6  6  5  6  3  5  2  2  3  4  4  5  6
##       O5 tot  grp
## 61617  4  82    L
## 61618  4 105 <NA>
## 61620  5 102 <NA>
## 61621  2  86    L
## 61622  4 100    L
## 61623  6 119    H

Assign items to different group according to total score

dtam <- aggregate(dta[, 1:25], by=list(dta$grp), mean)
print(dtam)
##   Group.1       A1       A2       A3       A4       A5       C1       C2
## 1       H 5.082601 5.504394 5.397188 5.349736 5.281195 5.177504 5.159930
## 2       L 4.265487 4.189159 3.899336 4.078540 3.953540 3.987832 3.731195
##         C3       C4       C5       E1       E2       E3       E4       E5
## 1 4.866432 5.117750 4.365554 4.794376 4.678383 4.945518 5.138840 5.261863
## 2 3.794248 3.884956 3.092920 3.245575 3.091814 3.234513 3.658186 3.608407
##         N1       N2       N3       N4       N5       O1       O2       O3
## 1 3.337434 3.959578 3.708260 3.307557 3.302285 5.395431 4.829525 5.203866
## 2 2.690265 3.238938 2.892699 3.196903 2.762168 4.357301 3.976770 3.845133
##         O4       O5
## 1 5.351494 5.045694
## 2 4.662611 4.168142

Mean of different groups

dtam <- t(dtam[, -1])

t-test

item_t <- sapply(dta[, 1:25], function(x) t.test(x ~ dta$grp)$statistic)
print(item_t)
##      A1.t      A2.t      A3.t      A4.t      A5.t      C1.t      C2.t      C3.t 
## 11.489868 25.195291 25.832339 18.156276 22.154179 19.903818 23.719711 16.384197 
##      C4.t      C5.t      E1.t      E2.t      E3.t      E4.t      E5.t      N1.t 
## 18.247550 15.204577 19.539964 20.135917 28.895154 20.871644 27.790835  7.575812 
##      N2.t      N3.t      N4.t      N5.t      O1.t      O2.t      O3.t      O4.t 
##  8.847026  9.616210  1.283397  5.988082 19.630090 10.562105 24.535055 11.736758 
##      O5.t 
## 13.186029
rslt2 <- data.frame(Item = rownames(dtam), low.mean.score = dtam[, 2],high.mean.score = dtam[, 1], mean.dif = dtam[, 1]-dtam[,2], t.value = item_t)
rslt2 |> knitr::kable()
Item low.mean.score high.mean.score mean.dif t.value
A1 A1 4.265487 5.082601 0.8171143 11.489869
A2 A2 4.189159 5.504394 1.3152344 25.195291
A3 A3 3.899336 5.397188 1.4978518 25.832339
A4 A4 4.078540 5.349736 1.2711966 18.156276
A5 A5 3.953540 5.281195 1.3276553 22.154179
C1 C1 3.987832 5.177504 1.1896725 19.903818
C2 C2 3.731195 5.159930 1.4287350 23.719711
C3 C3 3.794248 4.866432 1.0721845 16.384197
C4 C4 3.884956 5.117750 1.2327947 18.247550
C5 C5 3.092920 4.365554 1.2726332 15.204577
E1 E1 3.245575 4.794376 1.5488009 19.539964
E2 E2 3.091814 4.678383 1.5865690 20.135917
E3 E3 3.234513 4.945519 1.7110052 28.895154
E4 E4 3.658186 5.138840 1.4806542 20.871644
E5 E5 3.608407 5.261863 1.6534558 27.790834
N1 N1 2.690266 3.337434 0.6471686 7.575812
N2 N2 3.238938 3.959578 0.7206402 8.847026
N3 N3 2.892699 3.708260 0.8155610 9.616210
N4 N4 3.196903 3.307557 0.1106545 1.283397
N5 N5 2.762168 3.302285 0.5401166 5.988082
O1 O1 4.357301 5.395431 1.0381297 19.630090
O2 O2 3.976770 4.829525 0.8527556 10.562105
O3 O3 3.845133 5.203866 1.3587337 24.535055
O4 O4 4.662611 5.351494 0.6888832 11.736758
O5 O5 4.168142 5.045694 0.8775526 13.186030 
ggplot(data = rslt2, aes(x = reorder(Item, t.value, max), y = t.value)) +
 geom_point() +
 geom_hline(yintercept = 2, linetype = "dashed") +
 coord_flip() +
 labs(x = "Items", y = "t-value") +
 theme_bw()

Sadly N4 couldn’t show the differences between high score group and low score group

itotr <- psych::alpha(dta[, 1:25])$item.stats[, "r.drop"]
## Warning in psych::alpha(dta[, 1:25]): Some items were negatively correlated with the total scale and probably 
## should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( N1 N2 N3 N4 N5 O4 ) were negatively correlated with the total scale and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
 ldta <- list(A = dta[, 1:5], C  = dta[, 6:10], E = dta[ ,11:15], N = dta[ ,16:20], O = dta[ ,21:25])
isubalpha <- lapply(ldta, psych::alpha)
isubr <- c(isubalpha$A$item.stats[, "r.drop"],
           isubalpha$C$item.stats[, "r.drop"],
           isubalpha$E$item.stats[, "r.drop"],
           isubalpha$N$item.stats[, "r.drop"],
           isubalpha$O$item.stats[, "r.drop"])
rslt3 <- as.data.frame(t(rbind(itotr, isubr)))
names(rslt3) <- c("Item-total correlation without item", "Item-subscale correlation without item")
rslt3 |> knitr::kable()
Item-total correlation without item Item-subscale correlation without item
0.1383165 0.3190962
0.4398713 0.5759233
0.4523314 0.6035693
0.2952496 0.4145254
0.3953830 0.5004352
0.3326445 0.4654163
0.3656616 0.5128535
0.2582766 0.4769297
0.2770146 0.5731250
0.2141296 0.4860793
0.3033276 0.5153693
0.3233658 0.6142087
0.4724368 0.5049821
0.3756075 0.5827738
0.4717161 0.4634332
0.0514879 0.6778437
0.0671842 0.6548330
0.0937758 0.6781411
-0.1057078 0.5485366
0.0126812 0.4874632
0.3168253 0.3981233
0.1056166 0.3509392
0.4191609 0.4546552
0.1236295 0.2167170
0.1866414 0.4197456

A1 and A4 has good correlation to subscale A
E2 has not so good correlation to subscale E
N1, N2, N3 has not so good correlation to subscale N
O1, O2, O4 has good correlation to subscale O

itotalpha <- psych::alpha(dta[, 1:25])$alpha.drop[, "raw_alpha"]
## Warning in psych::alpha(dta[, 1:25]): Some items were negatively correlated with the total scale and probably 
## should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( N1 N2 N3 N4 N5 O4 ) were negatively correlated with the total scale and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
ialphad <- c(isubalpha$A$alpha.drop[, "raw_alpha"],
             isubalpha$C$alpha.drop[, "raw_alpha"],
             isubalpha$E$alpha.drop[, "raw_alpha"],
             isubalpha$N$alpha.drop[, "raw_alpha"],
             isubalpha$C$alpha.drop[, "raw_alpha"])
rslt4 <- as.data.frame(t(rbind(itotalpha, ialphad)))
names(rslt4) <- c("Main Reliability(item drop)", "Sub Reliability (item drop)")
rslt4 |> knitr::kable()
Main Reliability(item drop) Sub Reliability (item drop)
0.6979605 0.7314607
0.6762664 0.6332003
0.6734577 0.6150842
0.6850740 0.6963136
0.6783576 0.6582421
0.6833182 0.7044913
0.6801730 0.6869874
0.6885173 0.7000899
0.6868631 0.6630853
0.6925096 0.7031818
0.6840630 0.7312733
0.6821853 0.6924954
0.6712580 0.7329196
0.6781666 0.7056423
0.6714425 0.7457366
0.7066986 0.7598167
0.7048568 0.7673732
0.7032029 0.7594660
0.7199257 0.7982138
0.7107374 0.8167652
0.6852134 0.7044913
0.7017610 0.6869874
0.6773863 0.7000899
0.6979447 0.6630853
0.6938900 0.7031818

dropping A3 will cause severe reliability to subscale.
dropping A1, E1, E3, E4, E5, N1, N2, N3, N4, N5 will have relatively small effet on reliability to subscale.

So to sum things up, I will pick
A1, A3, A4 for subscale A
C1, C3, C5 for subscale C
E1, E3, E5 for subscale E
N2, N4, N5 for subscale N
O1, O2, O4 for subscale O