Check effect sizes are similar for continuous and categorical indicators

marklhc commented 1 year ago

Currently it seems numbers are smaller for categorical

YichiZhang2024 commented 1 year ago

The es_lavaan() function fails to compute the effect size for noninvariant loadings. See attached example:

library(tidyverse)
library(lavaan)
library(pinsearch)

## read in data
orig_dat <- haven::read_sav(file = "https://osf.io/download/wxjsg/")
dat <- orig_dat %>%
    filter_at(vars(class14, audit1:audit3), ~ !is.na(.)) %>%
    filter(eth %in% 1:4) %>%
    mutate_at(vars(class1:class15, audit1:audit10), as.numeric) %>%
    mutate(campuslive = factor(campuslive),
           eth = factor(eth, 1:4),
           group = case_when(
               # white male
               eth == 1 & gender == 1 ~ "1",
               # white female
               eth == 1 & gender == 2 ~ "2",
               # asian male
               eth == 2 & gender == 1 ~ "3",
               # asian female
               eth == 2 & gender == 2 ~ "4"
           )) %>%
    mutate(group = factor(group, 1:4)) %>%
    filter(group %in% 1:4) %>%
    arrange(group)
dat$class6 <- ifelse(dat$class6 > 3, 4, dat$class6)
dat[,9:23] <- sapply(dat %>% select("class1":"class15"), FUN = function(x) ifelse(x > 3, 4, x))

## specifying model

c_mod <- "
f1 =~ class1 + class2 + class3 + class4 + class5 + class6 +
      class7 + class8 + class9 + class10 + class11 + class12 +
      class13 + class14 + class15"

## continuous
res <- pinSearch(c_mod,
                 data = dat,
                 group = "group",
                 group.label = 1:4,
                 type = "residuals")
pinsearch:::es_lavaan(res[[1]])
#>       class15-f1 class1-f1 class14-f1 class9-f1 class10-f1 class12-f1
#> fmacs  0.1219231 0.1210571  0.1390922 0.1563828  0.1271502  0.1119919
#>        class4-f1
#> fmacs 0.09584454

## categorical

ps_thres <- pinSearch(c_mod,
                      data = dat,
                      group = "group",
                      ordered = TRUE,
                      group.label = 1:4,
                      type = "thresholds")
#> Unique variances are constrained to 1 for identification
ps_thres[[2]]
#>        lhs     rhs group       type
#> 1       f1  class3     2   loadings
#> 2       f1  class8     2   loadings
#> 3       f1  class6     1   loadings
#> 4       f1  class6     4   loadings
#> 5       f1 class14     1   loadings
#> 6  class14      t1     3 thresholds
#> 7   class3      t3     4 thresholds
#> 8   class9      t2     3 thresholds
#> 9  class15      t2     1 thresholds
#> 10 class15      t3     1 thresholds
#> 11 class14      t2     3 thresholds
#> 12 class14      t2     4 thresholds
#> 13 class14      t3     3 thresholds
#> 14  class8      t3     1 thresholds
pinsearch:::es_lavaan(ps_thres[[1]])
#>       class8-f1 class15-f1  class9-f1 class14-f1 class3-f1 class6-f1
#> fmacs 0.0878668 0.09903536 0.05507706  0.1611185 0.1285094 0.1528713

## using the same set of noninvariant parameters on continuous results

mod_psfit_ethc <- "
f1 =~ c(l1, l1, l1, l1) * class1 + c(l2, l2, l2, l2) * class2 +
c(l3, l3m, l3, l3) * class3 + c(l4, l4, l4, l4) * class4 +
c(l5, l5, l5, l5) * class5 + c(l6m, l6, l6, l6m) * class6 +
c(l7, l7, l7, l7) * class7 + c(l8, l8m, l8, l8) * class8 +
c(l9, l9, l9, l9) * class9 + c(l10, l10, l10, l10) * class10 +
c(l11, l11, l11, l11) * class11 + c(l12, l12, l12, l12) * class12 +
c(l13, l13, l13, l13) * class13 + c(l14m, l14, l14, l14) * class14 +
c(l15, l15, l15, l15) * class15

## common factor
f1 ~ c(0, NA, NA, NA) * 1
f1 ~~ c(1, NA, NA, NA) * f1

## intercepts
class1 ~ c(I1, I1, I1, I1) * 1
class2 ~ c(I2, I2, I2, I2) * 1
class3 ~ c(I3, I3, I3, I3m) * 1
class4 ~ c(I4, I4, I4, I4) * 1
class5 ~ c(I5, I5, I5, I5) * 1
class6 ~ c(I6, I6, I6, I6) * 1
class7 ~ c(I7, I7, I7, I7) * 1
class8 ~ c(I8m, I8, I8, I8) * 1
class9 ~ c(I9, I9, I9m, I9) * 1
class10 ~ c(I10, I10, I10, I10) * 1
class11 ~ c(I11, I11, I11, I11) * 1
class12 ~ c(I12, I12, I12, I12) * 1
class13 ~ c(I13, I13, I13, I13) * 1
class14 ~ c(I14, I14, I14m, I14m) * 1
class15 ~ c(I15m, I15, I15, I15) * 1

## unique factor variances
class1 ~~ c(U1, U1, U1, U1) * class1
class2 ~~ c(U2, U2, U2, U2) * class2
class3 ~~ c(U3, U3, U3, U3) * class3
class4 ~~ c(U4, U4, U4, U4) * class4
class5 ~~ c(U5, U5, U5, U5) * class5
class6 ~~ c(U6, U6, U6, U6) * class6
class7 ~~ c(U7, U7, U7, U7) * class7
class8 ~~ c(U8, U8, U8, U8) * class8
class9 ~~ c(U9, U9, U9, U9) * class9
class10 ~~ c(U10, U10, U10, U10) * class10
class11 ~~ c(U11, U11, U11, U11) * class11
class12 ~~ c(U12, U12, U12, U12) * class12
class13 ~~ c(U13, U13, U13, U13) * class13
class14 ~~ c(U14, U14, U14, U14) * class14
class15 ~~ c(U15, U15, U15, U15) * class15
"

par_str_fitc <- cfa(model = mod_psfit_ethc,
                    data = dat,
                    group = "group",
                    group.label = 1:4)
## summary(par_str_fitc, fit.measures = TRUE)

## Only print out effect sizes for five items, missing item 6
pinsearch:::es_lavaan(par_str_fitc)
#>       class8-f1 class15-f1 class9-f1 class3-f1 class14-f1
#> fmacs 0.1180813   0.101678 0.1446231 0.1980019  0.1571596

^{Created on 2023-01-27 with reprex v2.0.2}

marklhc commented 1 year ago

Could you try again @YichiZhang2024 to see whether it has been fixed? Thank you.

YichiZhang2024 commented 1 year ago

Here are the results. The number of items matches now.

library(tidyverse)
library(lavaan)
#> This is lavaan 0.6-14
#> lavaan is FREE software! Please report any bugs.
library(reprex)
library(pinsearch)
devtools::install_github("marklhc/pinsearch", ref = "master")
#> Skipping install of 'pinsearch' from a github remote, the SHA1 (9db432ed) has not changed since last install.
#>   Use `force = TRUE` to force installation
## read in data
orig_dat <- haven::read_sav(file = "https://osf.io/download/wxjsg/")
dat <- orig_dat %>%
    filter_at(vars(class14, audit1:audit3), ~ !is.na(.)) %>%
    filter(eth %in% 1:4) %>%
    mutate_at(vars(class1:class15, audit1:audit10), as.numeric) %>%
    mutate(campuslive = factor(campuslive),
           eth = factor(eth, 1:4),
           group = case_when(
               # white male
               eth == 1 & gender == 1 ~ "1",
               # white female
               eth == 1 & gender == 2 ~ "2",
               # asian male
               eth == 2 & gender == 1 ~ "3",
               # asian female
               eth == 2 & gender == 2 ~ "4"
           )) %>%
    mutate(group = factor(group, 1:4)) %>%
    filter(group %in% 1:4) %>%
    arrange(group)
dat$class6 <- ifelse(dat$class6 > 3, 4, dat$class6)
dat[,9:23] <- sapply(dat %>% select("class1":"class15"), FUN = function(x) ifelse(x > 3, 4, x))

## specifying model

c_mod <- "
f1 =~ class1 + class2 + class3 + class4 + class5 + class6 +
      class7 + class8 + class9 + class10 + class11 + class12 +
      class13 + class14 + class15"

## continuous
res <- pinSearch(c_mod,
                 data = dat,
                 group = "group",
                 group.label = 1:4,
                 type = "residuals")
pinsearch:::es_lavaan(res[[1]])
#>       class15-f1 class1-f1 class14-f1 class9-f1 class10-f1 class12-f1
#> fmacs  0.1219231 0.1210571  0.1390922 0.1563828  0.1271502  0.1119919
#>        class4-f1
#> fmacs 0.09584454
#       class15-f1 class1-f1 class14-f1 class9-f1 class10-f1 class12-f1  class4-f1
# fmacs   0.121925 0.1210569  0.1390924  0.156384  0.1271512  0.1119923 0.09584571

## categorical

ps_thres <- pinSearch(c_mod,
                      data = dat,
                      group = "group",
                      ordered = TRUE,
                      group.label = 1:4,
                      type = "thresholds")
#> Unique variances are constrained to 1 for identification
ps_thres[[2]]
#>        lhs     rhs group       type
#> 1       f1  class3     2   loadings
#> 2       f1  class8     2   loadings
#> 3       f1  class6     1   loadings
#> 4       f1  class6     4   loadings
#> 5       f1 class14     1   loadings
#> 6  class14      t1     3 thresholds
#> 7   class3      t3     4 thresholds
#> 8   class9      t2     3 thresholds
#> 9  class15      t2     1 thresholds
#> 10 class15      t3     1 thresholds
#> 11 class14      t2     3 thresholds
#> 12 class14      t2     4 thresholds
#> 13 class14      t3     3 thresholds
#> 14  class8      t3     1 thresholds
pinsearch:::es_lavaan(ps_thres[[1]])
#>       class8-f1 class15-f1  class9-f1 class14-f1 class3-f1 class6-f1
#> fmacs 0.0878668 0.09903536 0.05507706  0.1611185 0.1285094 0.1528713
#        class8-f1 class15-f1  class9-f1 class14-f1 class3-f1 class6-f1
# fmacs 0.08786646 0.09903646 0.05507853  0.1611189 0.1285096 0.1528724

## using the same set of noninvariant parameters on continuous results
## Item class6 has noninvariant loadins, item class3, class8, class9 and class14
## have noninvariant thresholds/intercepts.

mod_psfit_ethc <- "
f1 =~ c(l1, l1, l1, l1) * class1 + c(l2, l2, l2, l2) * class2 +
c(l3, l3m, l3, l3) * class3 + c(l4, l4, l4, l4) * class4 +
c(l5, l5, l5, l5) * class5 + c(l6a, l6, l6, l6b) * class6 +
c(l7, l7, l7, l7) * class7 + c(l8, l8m, l8, l8) * class8 +
c(l9, l9, l9, l9) * class9 + c(l10, l10, l10, l10) * class10 +
c(l11, l11, l11, l11) * class11 + c(l12, l12, l12, l12) * class12 +
c(l13, l13, l13, l13) * class13 + c(l14m, l14, l14, l14) * class14 +
c(l15, l15, l15, l15) * class15

## common factor
f1 ~ c(0, NA, NA, NA) * 1
f1 ~~ c(1, NA, NA, NA) * f1

## intercepts
class1 ~ c(I1, I1, I1, I1) * 1
class2 ~ c(I2, I2, I2, I2) * 1
class3 ~ c(I3, I3, I3, I3m) * 1
class4 ~ c(I4, I4, I4, I4) * 1
class5 ~ c(I5, I5, I5, I5) * 1
class6 ~ c(I6, I6, I6, I6) * 1
class7 ~ c(I7, I7, I7, I7) * 1
class8 ~ c(I8m, I8, I8, I8) * 1
class9 ~ c(I9, I9, I9m, I9) * 1
class10 ~ c(I10, I10, I10, I10) * 1
class11 ~ c(I11, I11, I11, I11) * 1
class12 ~ c(I12, I12, I12, I12) * 1
class13 ~ c(I13, I13, I13, I13) * 1
class14 ~ c(I14, I14, I14a, I14b) * 1
class15 ~ c(I15m, I15, I15, I15) * 1

## unique factor variances
class1 ~~ c(U1, U1, U1, U1) * class1
class2 ~~ c(U2, U2, U2, U2) * class2
class3 ~~ c(U3, U3, U3, U3) * class3
class4 ~~ c(U4, U4, U4, U4) * class4
class5 ~~ c(U5, U5, U5, U5) * class5
class6 ~~ c(U6, U6, U6, U6) * class6
class7 ~~ c(U7, U7, U7, U7) * class7
class8 ~~ c(U8, U8, U8, U8) * class8
class9 ~~ c(U9, U9, U9, U9) * class9
class10 ~~ c(U10, U10, U10, U10) * class10
class11 ~~ c(U11, U11, U11, U11) * class11
class12 ~~ c(U12, U12, U12, U12) * class12
class13 ~~ c(U13, U13, U13, U13) * class13
class14 ~~ c(U14, U14, U14, U14) * class14
class15 ~~ c(U15, U15, U15, U15) * class15
"

par_str_fitc <- cfa(model = mod_psfit_ethc,
                    data = dat,
                    group = "group",
                    group.label = 1:4)
summary(par_str_fitc, fit.measures = TRUE)
#> lavaan 0.6.14 ended normally after 53 iterations
#> 
#>   Estimator                                         ML
#>   Optimization method                           NLMINB
#>   Number of model parameters                       182
#>   Number of equality constraints                   121
#> 
#>   Number of observations per group:                   
#>     1                                              133
#>     2                                              259
#>     3                                               78
#>     4                                              106
#> 
#> Model Test User Model:
#>                                                       
#>   Test statistic                               954.760
#>   Degrees of freedom                               479
#>   P-value (Chi-square)                           0.000
#>   Test statistic for each group:
#>     1                                          307.203
#>     2                                          223.359
#>     3                                          228.104
#>     4                                          196.094
#> 
#> Model Test Baseline Model:
#> 
#>   Test statistic                              3642.119
#>   Degrees of freedom                               420
#>   P-value                                        0.000
#> 
#> User Model versus Baseline Model:
#> 
#>   Comparative Fit Index (CFI)                    0.852
#>   Tucker-Lewis Index (TLI)                       0.871
#> 
#> Loglikelihood and Information Criteria:
#> 
#>   Loglikelihood user model (H0)             -11506.839
#>   Loglikelihood unrestricted model (H1)     -11029.458
#>                                                       
#>   Akaike (AIC)                               23135.677
#>   Bayesian (BIC)                             23401.400
#>   Sample-size adjusted Bayesian (SABIC)      23207.750
#> 
#> Root Mean Square Error of Approximation:
#> 
#>   RMSEA                                          0.083
#>   90 Percent confidence interval - lower         0.075
#>   90 Percent confidence interval - upper         0.091
#>   P-value H_0: RMSEA <= 0.050                    0.000
#>   P-value H_0: RMSEA >= 0.080                    0.748
#> 
#> Standardized Root Mean Square Residual:
#> 
#>   SRMR                                           0.155
#> 
#> Parameter Estimates:
#> 
#>   Standard errors                             Standard
#>   Information                                 Expected
#>   Information saturated (h1) model          Structured
#> 
#> 
#> Group 1 [1]:
#> 
#> Latent Variables:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>   f1 =~                                               
#>     class1    (l1)    1.000                           
#>     class2    (l2)    1.141    0.054   21.146    0.000
#>     class3    (l3)   -0.393    0.057   -6.853    0.000
#>     class4    (l4)    1.078    0.053   20.430    0.000
#>     class5    (l5)    0.847    0.053   16.050    0.000
#>     class6   (l6a)    0.670    0.071    9.403    0.000
#>     class7    (l7)    1.141    0.054   21.037    0.000
#>     class8    (l8)   -0.323    0.062   -5.227    0.000
#>     class9    (l9)    0.835    0.053   15.628    0.000
#>     class10  (l10)    0.823    0.055   14.956    0.000
#>     class11  (l11)    0.927    0.048   19.362    0.000
#>     class12  (l12)    0.935    0.051   18.345    0.000
#>     class13  (l13)    0.843    0.053   16.048    0.000
#>     class14 (l14m)    0.904    0.077   11.774    0.000
#>     class15  (l15)    0.904    0.057   15.838    0.000
#> 
#> Intercepts:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     f1                0.000                           
#>    .class1    (I1)    3.719    0.090   41.424    0.000
#>    .class2    (I2)    3.089    0.101   30.451    0.000
#>    .class3    (I3)    2.645    0.064   41.176    0.000
#>    .class4    (I4)    2.779    0.097   28.730    0.000
#>    .class5    (I5)    2.471    0.082   30.186    0.000
#>    .class6    (I6)    2.016    0.067   30.060    0.000
#>    .class7    (I7)    3.403    0.102   33.492    0.000
#>    .class8   (I8m)    2.467    0.092   26.883    0.000
#>    .class9    (I9)    3.158    0.082   38.398    0.000
#>    .class10  (I10)    2.834    0.082   34.727    0.000
#>    .class11  (I11)    3.469    0.084   41.111    0.000
#>    .class12  (I12)    3.241    0.086   37.497    0.000
#>    .class13  (I13)    2.409    0.081   29.587    0.000
#>    .class14  (I14)    3.585    0.083   43.097    0.000
#>    .class15 (I15m)    2.472    0.110   22.426    0.000
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     f1                1.000                           
#>    .class1    (U1)    0.560    0.037   15.283    0.000
#>    .class2    (U2)    0.609    0.041   14.915    0.000
#>    .class3    (U3)    0.808    0.049   16.604    0.000
#>    .class4    (U4)    0.619    0.041   15.166    0.000
#>    .class5    (U5)    0.813    0.050   16.124    0.000
#>    .class6    (U6)    0.570    0.035   16.203    0.000
#>    .class7    (U7)    0.622    0.042   14.956    0.000
#>    .class8    (U8)    1.027    0.061   16.820    0.000
#>    .class9    (U9)    0.838    0.052   16.173    0.000
#>    .class10  (U10)    0.927    0.057   16.271    0.000
#>    .class11  (U11)    0.552    0.036   15.475    0.000
#>    .class12  (U12)    0.669    0.043   15.716    0.000
#>    .class13  (U13)    0.804    0.050   16.124    0.000
#>    .class14  (U14)    0.617    0.039   15.804    0.000
#>    .class15  (U15)    0.889    0.055   16.085    0.000
#> 
#> 
#> Group 2 [2]:
#> 
#> Latent Variables:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>   f1 =~                                               
#>     class1    (l1)    1.000                           
#>     class2    (l2)    1.141    0.054   21.146    0.000
#>     class3   (l3m)   -0.711    0.078   -9.103    0.000
#>     class4    (l4)    1.078    0.053   20.430    0.000
#>     class5    (l5)    0.847    0.053   16.050    0.000
#>     class6    (l6)    0.715    0.059   12.192    0.000
#>     class7    (l7)    1.141    0.054   21.037    0.000
#>     class8   (l8m)   -0.480    0.085   -5.621    0.000
#>     class9    (l9)    0.835    0.053   15.628    0.000
#>     class10  (l10)    0.823    0.055   14.956    0.000
#>     class11  (l11)    0.927    0.048   19.362    0.000
#>     class12  (l12)    0.935    0.051   18.345    0.000
#>     class13  (l13)    0.843    0.053   16.048    0.000
#>     class14  (l14)    0.850    0.058   14.747    0.000
#>     class15  (l15)    0.904    0.057   15.838    0.000
#> 
#> Intercepts:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     f1               -0.415    0.097   -4.263    0.000
#>    .class1    (I1)    3.719    0.090   41.424    0.000
#>    .class2    (I2)    3.089    0.101   30.451    0.000
#>    .class3    (I3)    2.645    0.064   41.176    0.000
#>    .class4    (I4)    2.779    0.097   28.730    0.000
#>    .class5    (I5)    2.471    0.082   30.186    0.000
#>    .class6    (I6)    2.016    0.067   30.060    0.000
#>    .class7    (I7)    3.403    0.102   33.492    0.000
#>    .class8    (I8)    2.708    0.066   40.952    0.000
#>    .class9    (I9)    3.158    0.082   38.398    0.000
#>    .class10  (I10)    2.834    0.082   34.727    0.000
#>    .class11  (I11)    3.469    0.084   41.111    0.000
#>    .class12  (I12)    3.241    0.086   37.497    0.000
#>    .class13  (I13)    2.409    0.081   29.587    0.000
#>    .class14  (I14)    3.585    0.083   43.097    0.000
#>    .class15  (I15)    2.747    0.092   29.750    0.000
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     f1                0.562    0.065    8.638    0.000
#>    .class1    (U1)    0.560    0.037   15.283    0.000
#>    .class2    (U2)    0.609    0.041   14.915    0.000
#>    .class3    (U3)    0.808    0.049   16.604    0.000
#>    .class4    (U4)    0.619    0.041   15.166    0.000
#>    .class5    (U5)    0.813    0.050   16.124    0.000
#>    .class6    (U6)    0.570    0.035   16.203    0.000
#>    .class7    (U7)    0.622    0.042   14.956    0.000
#>    .class8    (U8)    1.027    0.061   16.820    0.000
#>    .class9    (U9)    0.838    0.052   16.173    0.000
#>    .class10  (U10)    0.927    0.057   16.271    0.000
#>    .class11  (U11)    0.552    0.036   15.475    0.000
#>    .class12  (U12)    0.669    0.043   15.716    0.000
#>    .class13  (U13)    0.804    0.050   16.124    0.000
#>    .class14  (U14)    0.617    0.039   15.804    0.000
#>    .class15  (U15)    0.889    0.055   16.085    0.000
#> 
#> 
#> Group 3 [3]:
#> 
#> Latent Variables:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>   f1 =~                                               
#>     class1    (l1)    1.000                           
#>     class2    (l2)    1.141    0.054   21.146    0.000
#>     class3    (l3)   -0.393    0.057   -6.853    0.000
#>     class4    (l4)    1.078    0.053   20.430    0.000
#>     class5    (l5)    0.847    0.053   16.050    0.000
#>     class6    (l6)    0.715    0.059   12.192    0.000
#>     class7    (l7)    1.141    0.054   21.037    0.000
#>     class8    (l8)   -0.323    0.062   -5.227    0.000
#>     class9    (l9)    0.835    0.053   15.628    0.000
#>     class10  (l10)    0.823    0.055   14.956    0.000
#>     class11  (l11)    0.927    0.048   19.362    0.000
#>     class12  (l12)    0.935    0.051   18.345    0.000
#>     class13  (l13)    0.843    0.053   16.048    0.000
#>     class14  (l14)    0.850    0.058   14.747    0.000
#>     class15  (l15)    0.904    0.057   15.838    0.000
#> 
#> Intercepts:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     f1               -0.640    0.128   -5.002    0.000
#>    .class1    (I1)    3.719    0.090   41.424    0.000
#>    .class2    (I2)    3.089    0.101   30.451    0.000
#>    .class3    (I3)    2.645    0.064   41.176    0.000
#>    .class4    (I4)    2.779    0.097   28.730    0.000
#>    .class5    (I5)    2.471    0.082   30.186    0.000
#>    .class6    (I6)    2.016    0.067   30.060    0.000
#>    .class7    (I7)    3.403    0.102   33.492    0.000
#>    .class8    (I8)    2.708    0.066   40.952    0.000
#>    .class9   (I9m)    3.521    0.131   26.894    0.000
#>    .class10  (I10)    2.834    0.082   34.727    0.000
#>    .class11  (I11)    3.469    0.084   41.111    0.000
#>    .class12  (I12)    3.241    0.086   37.497    0.000
#>    .class13  (I13)    2.409    0.081   29.587    0.000
#>    .class14 (I14a)    3.159    0.121   26.045    0.000
#>    .class15  (I15)    2.747    0.092   29.750    0.000
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     f1                0.625    0.116    5.377    0.000
#>    .class1    (U1)    0.560    0.037   15.283    0.000
#>    .class2    (U2)    0.609    0.041   14.915    0.000
#>    .class3    (U3)    0.808    0.049   16.604    0.000
#>    .class4    (U4)    0.619    0.041   15.166    0.000
#>    .class5    (U5)    0.813    0.050   16.124    0.000
#>    .class6    (U6)    0.570    0.035   16.203    0.000
#>    .class7    (U7)    0.622    0.042   14.956    0.000
#>    .class8    (U8)    1.027    0.061   16.820    0.000
#>    .class9    (U9)    0.838    0.052   16.173    0.000
#>    .class10  (U10)    0.927    0.057   16.271    0.000
#>    .class11  (U11)    0.552    0.036   15.475    0.000
#>    .class12  (U12)    0.669    0.043   15.716    0.000
#>    .class13  (U13)    0.804    0.050   16.124    0.000
#>    .class14  (U14)    0.617    0.039   15.804    0.000
#>    .class15  (U15)    0.889    0.055   16.085    0.000
#> 
#> 
#> Group 4 [4]:
#> 
#> Latent Variables:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>   f1 =~                                               
#>     class1    (l1)    1.000                           
#>     class2    (l2)    1.141    0.054   21.146    0.000
#>     class3    (l3)   -0.393    0.057   -6.853    0.000
#>     class4    (l4)    1.078    0.053   20.430    0.000
#>     class5    (l5)    0.847    0.053   16.050    0.000
#>     class6   (l6b)    0.542    0.077    7.089    0.000
#>     class7    (l7)    1.141    0.054   21.037    0.000
#>     class8    (l8)   -0.323    0.062   -5.227    0.000
#>     class9    (l9)    0.835    0.053   15.628    0.000
#>     class10  (l10)    0.823    0.055   14.956    0.000
#>     class11  (l11)    0.927    0.048   19.362    0.000
#>     class12  (l12)    0.935    0.051   18.345    0.000
#>     class13  (l13)    0.843    0.053   16.048    0.000
#>     class14  (l14)    0.850    0.058   14.747    0.000
#>     class15  (l15)    0.904    0.057   15.838    0.000
#> 
#> Intercepts:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     f1               -0.727    0.120   -6.057    0.000
#>    .class1    (I1)    3.719    0.090   41.424    0.000
#>    .class2    (I2)    3.089    0.101   30.451    0.000
#>    .class3   (I3m)    2.960    0.103   28.753    0.000
#>    .class4    (I4)    2.779    0.097   28.730    0.000
#>    .class5    (I5)    2.471    0.082   30.186    0.000
#>    .class6    (I6)    2.016    0.067   30.060    0.000
#>    .class7    (I7)    3.403    0.102   33.492    0.000
#>    .class8    (I8)    2.708    0.066   40.952    0.000
#>    .class9    (I9)    3.158    0.082   38.398    0.000
#>    .class10  (I10)    2.834    0.082   34.727    0.000
#>    .class11  (I11)    3.469    0.084   41.111    0.000
#>    .class12  (I12)    3.241    0.086   37.497    0.000
#>    .class13  (I13)    2.409    0.081   29.587    0.000
#>    .class14 (I14b)    3.344    0.113   29.573    0.000
#>    .class15  (I15)    2.747    0.092   29.750    0.000
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     f1                0.652    0.106    6.150    0.000
#>    .class1    (U1)    0.560    0.037   15.283    0.000
#>    .class2    (U2)    0.609    0.041   14.915    0.000
#>    .class3    (U3)    0.808    0.049   16.604    0.000
#>    .class4    (U4)    0.619    0.041   15.166    0.000
#>    .class5    (U5)    0.813    0.050   16.124    0.000
#>    .class6    (U6)    0.570    0.035   16.203    0.000
#>    .class7    (U7)    0.622    0.042   14.956    0.000
#>    .class8    (U8)    1.027    0.061   16.820    0.000
#>    .class9    (U9)    0.838    0.052   16.173    0.000
#>    .class10  (U10)    0.927    0.057   16.271    0.000
#>    .class11  (U11)    0.552    0.036   15.475    0.000
#>    .class12  (U12)    0.669    0.043   15.716    0.000
#>    .class13  (U13)    0.804    0.050   16.124    0.000
#>    .class14  (U14)    0.617    0.039   15.804    0.000
#>    .class15  (U15)    0.889    0.055   16.085    0.000

pinsearch:::es_lavaan(par_str_fitc)
#>       class8-f1 class15-f1 class9-f1 class14-f1 class3-f1  class6-f1
#> fmacs 0.1180371  0.1018944 0.1410199   0.174897   0.19675 0.07738711

^{Created on 2023-02-27 with reprex v2.0.2}

marklhc commented 1 year ago

Great. Looks like the results are somewhat different between categorical and continuous, but there is no systematic pattern of one larger than the other or vice versa.

marklhc commented 1 year ago

Bug reported:

library(tidyverse)
library(lavaan)
#> This is lavaan 0.6-14
#> lavaan is FREE software! Please report any bugs.
library(pinsearch)
# devtools::install_github("marklhc/pinsearch", ref = "contrast")
## read in data
orig_dat <- haven::read_sav(file = "https://osf.io/download/wxjsg/")
dat <- orig_dat %>%
    filter_at(vars(class14, audit1:audit3), ~ !is.na(.)) %>%
    filter(eth %in% 1:4) %>%
    mutate_at(vars(class1:class15, audit1:audit10), as.numeric) %>%
    mutate(campuslive = factor(campuslive),
           eth = factor(eth, 1:4),
           group = case_when(
               # white male
               eth == 1 & gender == 1 ~ "1",
               # white female
               eth == 1 & gender == 2 ~ "2",
               # asian male
               eth == 2 & gender == 1 ~ "3",
               # asian female
               eth == 2 & gender == 2 ~ "4"
           )) %>%
    mutate(group = factor(group, 1:4)) %>%
    filter(group %in% 1:4) %>%
    arrange(group)
dat$class6 <- ifelse(dat$class6 > 3, 4, dat$class6)
dat[,9:23] <- sapply(dat %>% select("class1":"class15"), FUN = function(x) ifelse(x > 3, 4, x))

## specifying model

c_mod <- "
f1 =~ class1 + class2 + class3 + class4 + class5 + class6 +
      class7 + class8 + class9 + class10 + class11 + class12 +
      class13 + class14 + class15"

## categorical

ps_thres <- pinSearch(c_mod,
                      data = dat,
                      group = "group",
                      ordered = TRUE,
                      group.label = 1:4,
                      type = "thresholds")
#> Unique variances are constrained to 1 for identification
ps_thres[[2]]
#>        lhs     rhs group       type
#> 1       f1  class3     2   loadings
#> 2       f1  class8     2   loadings
#> 3       f1  class6     1   loadings
#> 4       f1  class6     4   loadings
#> 5       f1 class14     1   loadings
#> 6  class14      t1     3 thresholds
#> 7   class3      t3     4 thresholds
#> 8   class9      t2     3 thresholds
#> 9  class15      t2     1 thresholds
#> 10 class15      t3     1 thresholds
#> 11 class14      t2     3 thresholds
#> 12 class14      t2     4 thresholds
#> 13 class14      t3     3 thresholds
#> 14  class8      t3     1 thresholds
pinsearch:::es_lavaan(ps_thres[[1]])
#> Error in dimnames(x) <- dn: length of 'dimnames' [2] not equal to array extent

thre_mat <- rbind(lavInspect(ps_thres[[1]], what = "est")$`1`$tau[c(40:42, 7:9, 25:27, 43:45, 22:24, 16:18)],
                  lavInspect(ps_thres[[1]], what = "est")$`2`$tau[c(40:42, 7:9, 25:27, 43:45, 22:24, 16:18)],
                  lavInspect(ps_thres[[1]], what = "est")$`3`$tau[c(40:42, 7:9, 25:27, 43:45, 22:24, 16:18)],
                  lavInspect(ps_thres[[1]], what = "est")$`4`$tau[c(40:42, 7:9, 25:27, 43:45, 22:24, 16:18)])
colnames(thre_mat) <- rep(c(14, 3, 9, 15, 8, 6), each  = 3)
rownames(thre_mat) <- NULL

lambda_mat <- rbind(lavInspect(ps_thres[[1]], what = "est")$`1`$lambda[c(14, 3, 9, 15, 8, 6)],
                    lavInspect(ps_thres[[1]], what = "est")$`2`$lambda[c(14, 3, 9, 15, 8, 6)],
                    lavInspect(ps_thres[[1]], what = "est")$`3`$lambda[c(14, 3, 9, 15, 8, 6)],
                    lavInspect(ps_thres[[1]], what = "est")$`4`$lambda[c(14, 3, 9, 15, 8, 6)])
rownames(lambda_mat) <- NULL
colnames(lambda_mat) <- c(14, 3, 9, 15, 8, 6)
compute_pvar <- function(x, g, na.rm = TRUE) {
  mean((x - ave(x, g, FUN = \(xx) mean(xx, na.rm = na.rm)))^2,
       na.rm = na.rm)
}
pvar <- apply(
  dat[,9:23],
  MARGIN = 2,
  FUN = compute_pvar,
  g = dat$group
)
# fmacs (white vs asian)
fmacs_eth <-pinsearch::fmacs_ordered(
  thresholds = thre_mat,
  loadings = lambda_mat,
  pooled_item_sd = sqrt(pvar)[c(14, 3, 9, 15, 8, 6)],
  latent_mean = 0,
  latent_sd = 1,
  group_factor = c(1,1,2,2),
  num_obs = summary(ps_thres[[1]])$data$nobs
)
# fmacs (male vs female)
fmacs_gen <- pinsearch::fmacs_ordered(
  thresholds = thre_mat,
  loadings = lambda_mat,
  pooled_item_sd = sqrt(pvar)[c(14, 3, 9, 15, 8, 6)],
  latent_mean = 0,
  latent_sd = 1,
  group_factor = c(1,2,1,2),
  num_obs = summary(ps_thres[[1]])$data$nobs
)
# overall 
fmacs <- pinsearch::fmacs_ordered(
  thresholds = thre_mat,
  loadings = lambda_mat,
  pooled_item_sd = sqrt(pvar)[c(14, 3, 9, 15, 8, 6)],
  latent_mean = 0,
  latent_sd = 1,
  num_obs = summary(ps_thres[[1]])$data$nobs
)
# interaction
fmacs_int <- pinsearch::fmacs_ordered(
  thresholds = thre_mat,
  loadings = lambda_mat,
  pooled_item_sd = sqrt(pvar)[c(14, 3, 9, 15, 8, 6)],
  latent_mean = 0,
  latent_sd = 1,
  group_factor = c(1,2,2,1),
  num_obs = summary(ps_thres[[1]])$data$nobs
)
print("ethnicity")
#> [1] "ethnicity"
round(fmacs_eth, 3)
#>         14     3     9    15     8     6
#> fmacs 0.04 0.077 0.014 0.024 0.071 0.055
print("gender")
#> [1] "gender"
round(fmacs_gen, 3)
#>          14     3     9    15     8     6
#> fmacs 0.059 0.085 0.026 0.022 0.042 0.113
print("general")
#> [1] "general"
round(fmacs, 3)
#>          14     3     9    15     8     6
#> fmacs 0.095 0.124 0.037 0.042 0.105 0.199
print("interaction")
#> [1] "interaction"
round(fmacs_int, 3)
#>          14     3     9    15     8    6
#> fmacs 0.068 0.094 0.024 0.014 0.047 0.17

^{Created on 2023-03-09 with reprex v2.0.2}

marklhc commented 4 months ago

Results seem reasonable so far. Will close for now.

marklhc / pinsearch

Check effect sizes are similar for continuous and categorical indicators #1