library(psych) 
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.0     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.1     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x ggplot2::%+%()   masks psych::%+%()
## x ggplot2::alpha() masks psych::alpha()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
library(ggplot2)
library(corrplot)
## corrplot 0.84 loaded
library(lavaan)
## This is lavaan 0.6-8
## lavaan is FREE software! Please report any bugs.
## 
## Attaching package: 'lavaan'
## The following object is masked from 'package:psych':
## 
##     cor2cov
library(semPlot)
library(here)
## here() starts at /Users/caoanjie/Desktop/projects/psych262/assignment_3

Reading in data

d <- read_csv(here("data/new_rt.csv"))
## Parsed with column specification:
## cols(
##   subject = col_character(),
##   block_number = col_double(),
##   block_type = col_character(),
##   trial_number = col_double(),
##   item_type = col_character(),
##   trial_type = col_character(),
##   trial_complexity = col_character(),
##   item_id = col_character(),
##   rt = col_double(),
##   exposure_type = col_character(),
##   half = col_character(),
##   block_deviant_number = col_double(),
##   trial_type_index = col_character(),
##   first_dev_position = col_double(),
##   second_dev_position = col_double()
## )

data preprocessing

Here the data processing separates the original dataframe into two halves. The four factors of interests are the type of trials: background simple trial, background complex trial, deviant simple trial, and deviant complex trial. I separated the looking time data by half to examine whether the factorial invariances hold throughout the

d_first_half <- d %>% 
  filter(block_number  < 6) %>% 
  group_by(subject, trial_complexity, item_type) %>% 
  summarise(rt = mean(rt)) %>% 
  unite(trial_type, trial_complexity, item_type) %>% 
  pivot_wider(names_from = trial_type, 
              values_from = rt, 
              names_prefix = "1_")
## `summarise()` has grouped output by 'subject', 'trial_complexity'. You can override using the `.groups` argument.
d_second_half <- d %>% 
  filter(block_number > 6 | block_number == 6) %>% 
  group_by(subject, trial_complexity, item_type) %>% 
  summarise(rt = mean(rt)) %>% 
  unite(trial_type, trial_complexity, item_type) %>% 
  pivot_wider(names_from = trial_type, 
              values_from = rt, 
               names_prefix = "2_")
## `summarise()` has grouped output by 'subject', 'trial_complexity'. You can override using the `.groups` argument.
d_processed <- left_join(d_first_half, 
                         d_second_half, 
                         by = "subject") %>% 
  drop_na() %>% 
  janitor::clean_names()
describe(d_processed[, -1])
##                               vars   n    mean      sd  median trimmed     mad
## x1_complex_background            1 158 4347.57 2985.79 3849.77 4005.00 2688.26
## x1_complex_dissimilar_deviant    2 158 5001.80 3553.25 4288.57 4610.72 3177.34
## x1_simple_background             3 158 3916.27 2916.20 2986.15 3537.82 2429.04
## x1_simple_dissimilar_deviant     4 158 4321.52 3942.97 3342.61 3673.98 2781.77
## x2_complex_background            5 158 2895.70 1920.90 2544.34 2711.47 2014.10
## x2_complex_dissimilar_deviant    6 158 3480.01 3050.52 2876.37 3031.52 2491.88
## x2_simple_background             7 158 2629.62 2149.83 2032.67 2311.78 1674.62
## x2_simple_dissimilar_deviant     8 158 2820.68 2579.64 2232.56 2446.80 2346.30
##                                  min      max    range skew kurtosis     se
## x1_complex_background         184.19 15799.46 15615.27 1.31     2.33 237.54
## x1_complex_dissimilar_deviant  44.94 16776.66 16731.71 1.01     0.73 282.68
## x1_simple_background          225.62 16662.87 16437.25 1.44     2.58 232.00
## x1_simple_dissimilar_deviant   93.99 24488.67 24394.68 2.12     6.00 313.69
## x2_complex_background         183.11  9676.58  9493.48 0.92     0.67 152.82
## x2_complex_dissimilar_deviant 112.91 17792.69 17679.78 1.86     4.92 242.69
## x2_simple_background          163.92 12435.71 12271.78 1.87     4.78 171.03
## x2_simple_dissimilar_deviant   78.57 13313.68 13235.11 1.39     2.02 205.23
corrplot(cor(d_processed[,-1]), order = "original", tl.col='black', tl.cex=.75) 

# configural invariance model

configural <- '
# Define the latent factors.
first_half_rt =~ NA*x1_complex_background    + lambda1*x1_complex_background     + x1_complex_dissimilar_deviant + x1_simple_background  + x1_simple_dissimilar_deviant 
second_half_rt =~NA*x2_complex_background    + lambda1*x2_complex_background     + x2_complex_dissimilar_deviant + x2_simple_background  + x2_simple_dissimilar_deviant

# Intercepts
x1_complex_background ~ i1*1
x1_complex_dissimilar_deviant ~ 1
x1_simple_background ~ 1
x1_simple_dissimilar_deviant ~ 1
x2_complex_background ~ i1*1
x2_complex_dissimilar_deviant ~ 1
x2_simple_background ~ 1
x2_simple_dissimilar_deviant ~ 1

# Unique Variances
x1_complex_background ~~ x1_complex_background
x1_complex_dissimilar_deviant ~~ x1_complex_dissimilar_deviant
x1_simple_background ~~ x1_simple_background
x1_simple_dissimilar_deviant ~~ x1_simple_dissimilar_deviant
x2_complex_background ~~ x2_complex_background
x2_complex_dissimilar_deviant ~~ x2_complex_dissimilar_deviant
x2_simple_background ~~ x2_simple_background
x2_simple_dissimilar_deviant ~~ x2_simple_dissimilar_deviant

# Latent Variable Means
first_half_rt ~ 0*1
second_half_rt ~ 1

# Latent Variable Variances and Covariance
first_half_rt ~~ 1*first_half_rt
second_half_rt ~~ second_half_rt
first_half_rt ~~ second_half_rt
'
fit_configural <- cfa(configural, data = d_processed, mimic = "mplus")
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some observed variances are larger than 1000000
##   lavaan NOTE: use varTable(fit) to investigate
## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
##     The variance-covariance matrix of the estimated parameters (vcov)
##     does not appear to be positive definite! The smallest eigenvalue
##     (= -1.604676e-09) is smaller than zero. This may be a symptom that
##     the model is not identified.
summary(fit_configural, fit.measures = TRUE)
## lavaan 0.6-8 ended normally after 226 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        27
##   Number of equality constraints                     2
##                                                       
##   Number of observations                           158
##   Number of missing patterns                         1
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                92.043
##   Degrees of freedom                                19
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               711.134
##   Degrees of freedom                                28
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.893
##   Tucker-Lewis Index (TLI)                       0.842
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)             -11520.427
##   Loglikelihood unrestricted model (H1)     -11474.406
##                                                       
##   Akaike (AIC)                               23090.855
##   Bayesian (BIC)                             23167.420
##   Sample-size adjusted Bayesian (BIC)        23088.283
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.156
##   90 Percent confidence interval - lower         0.125
##   90 Percent confidence interval - upper         0.189
##   P-value RMSEA <= 0.05                          0.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.054
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                     Estimate     Std.Err  z-value       P(>|z|)
##   first_half_rt =~                                             
##     x1_cmp_ (lmb1)     2252.806  207.933        10.834    0.000
##     x1_cm__            2729.991  245.118        11.137    0.000
##     x1_smp_            2265.608  199.989        11.329    0.000
##     x1_sm__            2917.764  277.109        10.529    0.000
##   second_half_rt =~                                            
##     x2_cmp_ (lmb1)     2252.806  207.933        10.834    0.000
##     x2_cm__            3199.752  436.219         7.335    0.000
##     x2_smp_            2399.933  316.163         7.591    0.000
##     x2_sm__            2831.382  378.332         7.484    0.000
## 
## Covariances:
##                    Estimate     Std.Err  z-value       P(>|z|)
##   first_half_rt ~~                                            
##     second_half_rt       0.562    0.067         8.341    0.000
## 
## Intercepts:
##                    Estimate     Std.Err  z-value       P(>|z|)
##    .x1_cmplx_ (i1)    4347.569  236.784        18.361    0.000
##    .x1_cmpl__         5001.802  281.786        17.750    0.000
##    .x1_smpl_b         3916.266  231.265        16.934    0.000
##    .x1_smpl__         4321.519  312.692        13.820    0.000
##    .x2_cmplx_ (i1)    4347.569  236.784        18.361    0.000
##    .x2_cmpl__         5542.159  447.105        12.396    0.000
##    .x2_smpl_b         4176.310  325.091        12.847    0.000
##    .x2_smpl__         4645.422  387.602        11.985    0.000
##     frst_hlf_            0.000                                
##     scnd_hlf_           -0.644    0.110        -5.874    0.000
## 
## Variances:
##                    Estimate     Std.Err  z-value       P(>|z|)
##    .x1_cmplx_bckgr 3783408.245    0.123  30677043.728    0.000
##    .x1_cmplx_dssm_ 5092843.263    0.029 176555677.856    0.000
##    .x1_smpl_bckgrn 3317433.381    0.047  70933293.867    0.000
##    .x1_smpl_dssml_ 6935285.200    0.015 450673522.922    0.000
##    .x2_cmplx_bckgr 1255962.988    0.459   2733515.769    0.000
##    .x2_cmplx_dssm_ 4383793.132    0.047  92437556.001    0.000
##    .x2_smpl_bckgrn 1856821.298    0.158  11772348.015    0.000
##    .x2_smpl_dssml_ 2804703.972    0.121  23114397.414    0.000
##     first_half_rt        1.000                                
##     second_half_rt       0.475    0.100         4.738    0.000
semPaths(fit_configural, what="est", 
         sizeLat = 7, sizeMan = 7, edge.label.cex = .75)

configural invariance model with correlated uniqueness

corrunique <- '
# Define the latent factors.
first_half_rt =~ NA*x1_complex_background    + lambda1*x1_complex_background     + x1_complex_dissimilar_deviant + x1_simple_background  + x1_simple_dissimilar_deviant 
second_half_rt =~NA*x2_complex_background    + lambda1*x2_complex_background     + x2_complex_dissimilar_deviant + x2_simple_background  + x2_simple_dissimilar_deviant

# Intercepts
x1_complex_background ~ i1*1
x1_complex_dissimilar_deviant ~ 1
x1_simple_background ~ 1
x1_simple_dissimilar_deviant ~ 1
x2_complex_background ~ i1*1
x2_complex_dissimilar_deviant ~ 1
x2_simple_background ~ 1
x2_simple_dissimilar_deviant ~ 1

# Unique Variances and Covariances
x1_complex_background ~~ x1_complex_background
x1_complex_dissimilar_deviant ~~ x1_complex_dissimilar_deviant
x1_simple_background ~~ x1_simple_background
x1_simple_dissimilar_deviant ~~ x1_simple_dissimilar_deviant
x2_complex_background ~~ x2_complex_background
x2_complex_dissimilar_deviant ~~ x2_complex_dissimilar_deviant
x2_simple_background ~~ x2_simple_background
x2_simple_dissimilar_deviant ~~ x2_simple_dissimilar_deviant

x1_complex_background ~~ x2_complex_background
x1_complex_dissimilar_deviant ~~ x2_complex_dissimilar_deviant
x1_simple_background ~~ x2_simple_background
x1_simple_dissimilar_deviant ~~ x2_simple_dissimilar_deviant


# Latent Variable Means
first_half_rt ~ 0*1
second_half_rt ~ 1

# Latent Variable Variances and Covariance
first_half_rt ~~ 1*first_half_rt
second_half_rt ~~ second_half_rt
first_half_rt ~~ second_half_rt
'
fit_corrunique <- cfa(corrunique, data = d_processed, mimic = "mplus")
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some observed variances are larger than 1000000
##   lavaan NOTE: use varTable(fit) to investigate
## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
##     The variance-covariance matrix of the estimated parameters (vcov)
##     does not appear to be positive definite! The smallest eigenvalue
##     (= -8.183856e-09) is smaller than zero. This may be a symptom that
##     the model is not identified.
summary(fit_corrunique, fit.measures = TRUE)
## lavaan 0.6-8 ended normally after 227 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        31
##   Number of equality constraints                     2
##                                                       
##   Number of observations                           158
##   Number of missing patterns                         1
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                92.043
##   Degrees of freedom                                15
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               711.134
##   Degrees of freedom                                28
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.887
##   Tucker-Lewis Index (TLI)                       0.789
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)             -11520.427
##   Loglikelihood unrestricted model (H1)     -11474.406
##                                                       
##   Akaike (AIC)                               23098.855
##   Bayesian (BIC)                             23187.670
##   Sample-size adjusted Bayesian (BIC)        23095.871
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.180
##   90 Percent confidence interval - lower         0.146
##   90 Percent confidence interval - upper         0.217
##   P-value RMSEA <= 0.05                          0.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.054
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                     Estimate     Std.Err  z-value       P(>|z|)
##   first_half_rt =~                                             
##     x1_cmp_ (lmb1)     2252.807  207.933        10.834    0.000
##     x1_cm__            2729.992  245.118        11.137    0.000
##     x1_smp_            2265.609  199.989        11.329    0.000
##     x1_sm__            2917.765  277.109        10.529    0.000
##   second_half_rt =~                                            
##     x2_cmp_ (lmb1)     2252.807  207.933        10.834    0.000
##     x2_cm__            3199.756  436.220         7.335    0.000
##     x2_smp_            2399.934  316.164         7.591    0.000
##     x2_sm__            2831.385  378.332         7.484    0.000
## 
## Covariances:
##                                    Estimate     Std.Err  z-value       P(>|z|)
##  .x1_complex_background ~~                                                    
##    .x2_cmplx_bckgr                       0.743    0.188         3.952    0.000
##  .x1_complex_dissimilar_deviant ~~                                            
##    .x2_cmplx_dssm_                       0.493    0.018        27.265    0.000
##  .x1_simple_background ~~                                                     
##    .x2_smpl_bckgrn                       0.756    0.068        11.100    0.000
##  .x1_simple_dissimilar_deviant ~~                                             
##    .x2_smpl_dssml_                       0.489    0.017        29.080    0.000
##   first_half_rt ~~                                                            
##     second_half_rt                       0.562    0.067         8.341    0.000
## 
## Intercepts:
##                    Estimate     Std.Err  z-value       P(>|z|)
##    .x1_cmplx_ (i1)    4347.569  236.784        18.361    0.000
##    .x1_cmpl__         5001.803  281.786        17.750    0.000
##    .x1_smpl_b         3916.266  231.265        16.934    0.000
##    .x1_smpl__         4321.520  312.692        13.820    0.000
##    .x2_cmplx_ (i1)    4347.569  236.784        18.361    0.000
##    .x2_cmpl__         5542.160  447.105        12.396    0.000
##    .x2_smpl_b         4176.311  325.091        12.847    0.000
##    .x2_smpl__         4645.423  387.602        11.985    0.000
##     frst_hlf_            0.000                                
##     scnd_hlf_           -0.644    0.110        -5.874    0.000
## 
## Variances:
##                    Estimate     Std.Err  z-value       P(>|z|)
##    .x1_cmplx_bckgr 3783407.800    0.123  30747903.856    0.000
##    .x1_cmplx_dssm_ 5092843.877    0.028 179657753.169    0.000
##    .x1_smpl_bckgrn 3317434.225    0.046  71606383.547    0.000
##    .x1_smpl_dssml_ 6935280.360    0.015 450377019.404    0.000
##    .x2_cmplx_bckgr 1255963.078    0.461   2726648.720    0.000
##    .x2_cmplx_dssm_ 4383794.462    0.048  92225942.188    0.000
##    .x2_smpl_bckgrn 1856820.366    0.158  11769475.236    0.000
##    .x2_smpl_dssml_ 2804703.557    0.121  23172034.788    0.000
##     first_half_rt        1.000                                
##     second_half_rt       0.475    0.100         4.738    0.000
semPaths(fit_corrunique, what="est",
         sizeLat = 7, sizeMan = 7, edge.label.cex = .75)