#================================================================== # Notes #================================================================== # Note that some of the cognitive test scores are scaled so that all test variances are # within a similar range. This is to aid model convergence. # Some of the test slopes' residual variances were near zero and estimated as negative. # To avoid these out of bounds estimates, residual variances that were estimated as # negative are fixed to zero in the model. # We use the marker variable method of scale setting. #### functions #### if (!require("pacman")) install.packages("pacman") pacman::p_load(dplyr, haven, knitr,lavaan, psych, qgraph, readxl, semPlot, semTools, skimr) #(.packages()) #Disambiguate functions select <- dplyr::select summarise <- dplyr::summarise filter <- dplyr::filter rescale <- psych::rescale #### red data #### data<- read.csv("", header=T, fileEncoding = "UTF-8-BOM")#add file name here #================================================================== # Select variables and recode missing #================================================================== data=data[,c("lbc36no", "matreas_w1", "matreas_w2", "matreas_w3", "matreas_w4", "blkdes_w1", "blkdes_w2", "blkdes_w3", "blkdes_w4", "spantot_w1", "spantot_w2", "spantot_w3", "spantot_w4", "vpatotal_w1", "vpatotal_w2", "vpatotal_w3", "vpatotal_w4", "lmtotal_w1", "lmtotal_w2", "lmtotal_w3", "lmtotal_w4", "digback_w1", "digback_w2", "digback_w3", "digback_w4", "nart_w1", "nart_w2", "nart_total_w3", "nart_total_w4", "wtar_w1", "wtar_w2", "wtar_total_w3", "wtar_total_w4", "vftot_w1", "vftot_w2", "vftot_w3" , "vftot_w4", "digsym_w1", "digsym_w2" , "digsym_w3", "digsym_w4", "symsear_w1", "symsear_w2", "symsear_w3", "symsear_w4", "crtmean_w1", "crtmean_w2", "crtmean_w3", "crtmean_w4", "ittotal_w1", "ittotal_w2" , "ittotal_w3", "ittotal_w4", "matreas_w5", "blkdes_w5", "spantot_w5", "vpa_total_w5", "lmtotal_w5", "digback_w5", "nart_total_w5", "wtar_total_w5", "vftot_w5", "digsym_w5", "symsear_w5", "crtmean_w5", "ittotal_w5")] #recode missing data[data == -999] <- NA data[data == -777] <- NA data[data == 999] <- NA #rescaled some of the cognitive test variables so that variances are within a similar #range see http://www.statmodel.com/discussion/messages/11/1615.html?1335376547 dset_mod <- mutate(data, blkdes_w1 = blkdes_w1/2, blkdes_w2 = blkdes_w2/2, blkdes_w3 = blkdes_w3/2, blkdes_w4 = blkdes_w4/2, blkdes_w5 = blkdes_w5/2, vftot_w1 = vftot_w1/2, vftot_w2 = vftot_w2/2, vftot_w3 = vftot_w3/2, vftot_w4 = vftot_w4/2, vftot_w5 = vftot_w5/2, lmtotal_w1 = lmtotal_w1/3, lmtotal_w2 = lmtotal_w2/3, lmtotal_w3 = lmtotal_w3/3, lmtotal_w4 = lmtotal_w4/3, lmtotal_w5 = lmtotal_w5/3, digback_w1 = 3*digback_w1, digback_w2 = 3*digback_w2, digback_w3 = 3*digback_w3, digback_w4 = 3*digback_w4, digback_w5 = 3*digback_w5, digsym_w1 = digsym_w1/2, digsym_w2 = digsym_w2/2, digsym_w3 = digsym_w3/2, digsym_w4 = digsym_w4/2, digsym_w5 = digsym_w5/2, ittotal_w1 = ittotal_w1/2, ittotal_w2 = ittotal_w2/2, ittotal_w3 = ittotal_w3/2, ittotal_w4 = ittotal_w4/2, ittotal_w5 = ittotal_w5/2, crtmean_w1 = -50 * crtmean_w1, crtmean_w2 = -50 * crtmean_w2, crtmean_w3 = -50 * crtmean_w3, crtmean_w4 = -50 * crtmean_w4, crtmean_w5 = -50 * crtmean_w5) #================================================================== # Hierarchical model for general cognitive ability #================================================================== #### growth curve models for individual cognitive tests #### pgmodel1 <- ' Imatreas =~ 1*matreas_w1 + 1*matreas_w2 + 1*matreas_w3 + 1*matreas_w4 + 1*matreas_w5 Smatreas =~ 0*matreas_w1 + 2.98*matreas_w2 + 6.75*matreas_w3 + 9.82*matreas_w4 + 12.54*matreas_w5 ' fit1 <- growth(pgmodel1, dset_mod, missing = "ml.x") summary(fit1, standardized = T) pgmodel2 <-' Iblkdes =~ 1*blkdes_w1 + 1*blkdes_w2 + 1*blkdes_w3 + 1*blkdes_w4 + 1*blkdes_w5 Sblkdes=~ 0*blkdes_w1 + 2.98*blkdes_w2 + 6.75*blkdes_w3 + 9.82*blkdes_w4 + 12.54*blkdes_w5 ' fit2 <- growth(pgmodel2, dset_mod, missing = "ml.x") summary(fit2, standardized = T) pgmodel3 <- ' Ispantot =~ 1*spantot_w1 + 1*spantot_w2 + 1*spantot_w3 + 1*spantot_w4 + 1*spantot_w5 Sspantot=~ 0*spantot_w1 + 2.98*spantot_w2 + 6.75*spantot_w3 + 9.82*spantot_w4 + 12.54*spantot_w5 ' fit3 <- growth(pgmodel3, dset_mod, missing = "ml.x") summary(fit3, standardized = T) pgmodel4 <- ' Inart =~ 1*nart_w1 + 1*nart_w2 + 1*nart_total_w3 + 1*nart_total_w4 + 1*nart_total_w5 Snart =~ 0*nart_w1 + 2.98*nart_w2 + 6.75*nart_total_w3 + 9.82*nart_total_w4 + 12.54*nart_total_w5 ' fit4 <- growth(pgmodel4, dset_mod, missing = "ml.x") summary(fit4, standardized = T) pgmodel5 <- ' Iwtar =~ 1*wtar_w1 + 1*wtar_w2 + 1*wtar_total_w3 + 1*wtar_total_w4 + 1*wtar_total_w5 Swtar =~ 0*wtar_w1 + 2.98*wtar_w2 + 6.75*wtar_total_w3 + 9.82*wtar_total_w4 + 12.54*wtar_total_w5 ' fit5 <- growth(pgmodel5, dset_mod, missing = "ml.x") summary(fit5, standardized = T) pgmodel6 <- ' Ivftot =~ 1*vftot_w1 + 1*vftot_w2 + 1*vftot_w3 + 1*vftot_w4 + 1*vftot_w5 Svftot =~ 0*vftot_w1 + 2.98*vftot_w2 + 6.75*vftot_w3 + 9.82*vftot_w4 + 12.54*vftot_w5 ' fit6 <- growth(pgmodel6, dset_mod, missing = "ml.x") summary(fit6, standardized = T) pgmodel7 <-' Ivpatotal =~ 1*vpatotal_w1 + 1*vpatotal_w2 + 1*vpatotal_w3 + 1*vpatotal_w4 + 1*vpa_total_w5 Svpatotal =~ 0*vpatotal_w1 + 2.98*vpatotal_w2 + 6.75*vpatotal_w3 + 9.82*vpatotal_w4 + 12.54*vpa_total_w5 ' fit7 <- growth(pgmodel7, dset_mod, missing = "ml.x") summary(fit7, standardized = T) pgmodel8 <- ' Ilmtotal =~ 1*lmtotal_w1 + 1*lmtotal_w2 + 1*lmtotal_w3 + 1*lmtotal_w4 + 1*lmtotal_w5 Slmtotal =~ 0*lmtotal_w1 + 2.98*lmtotal_w2 + 6.75*lmtotal_w3 + 9.82*lmtotal_w4 + 12.54*lmtotal_w5 ' fit8 <- growth(pgmodel8, dset_mod, missing = "ml.x") summary(fit8, standardized = T) pgmodel9 <- ' Idigback =~ 1*digback_w1 + 1*digback_w2 + 1*digback_w3 + 1*digback_w4 + 1*digback_w5 Sdigback =~ 0*digback_w1 + 2.98*digback_w2 + 6.75*digback_w3 + 9.82*digback_w4 + 12.54*digback_w5 ' fit9 <- growth(pgmodel9, dset_mod, missing = "ml.x") summary(fit9, standardized = T) pgmodel10 <- ' Isymsear =~ 1*symsear_w1 + 1*symsear_w2 + 1*symsear_w3 + 1*symsear_w4 + 1*symsear_w5 Ssymsear =~ 0*symsear_w1 + 2.98*symsear_w2 + 6.75*symsear_w3 + 9.82*symsear_w4 + 12.54*symsear_w5 ' fit10 <- growth(pgmodel10, dset_mod, missing = "ml.x") summary(fit10, standardized = T) pgmodel11 <- ' Idigsym =~ 1*digsym_w1 + 1*digsym_w2 + 1*digsym_w3 + 1*digsym_w4 + 1*digsym_w5 Sdigsym =~ 0*digsym_w1 + 2.98*digsym_w2 + 6.75*digsym_w3 + 9.82*digsym_w4 + 12.54*digsym_w5 ' fit11 <- growth(pgmodel11, dset_mod, missing = "ml.x") summary(fit11, standardized = T) pgmodel12 <- ' Iittotal =~ 1*ittotal_w1 + 1*ittotal_w2 + 1*ittotal_w3 + 1*ittotal_w4 + 1*ittotal_w5 Sittotal =~ 0*ittotal_w1 + 2.98*ittotal_w2 + 6.75*ittotal_w3 + 9.82*ittotal_w4 + 12.54*ittotal_w5 ' fit12 <- growth(pgmodel12, dset_mod, missing = "ml.x") summary(fit12, standardized = T) pgmodel13 <- ' Icrtmean =~ 1*crtmean_w1 + 1*crtmean_w2 + 1*crtmean_w3 + 1*crtmean_w4 + 1*crtmean_w5 Scrtmean =~ 0*crtmean_w1 + 2.98*crtmean_w2 + 6.75*crtmean_w3 + 9.82*crtmean_w4 + 12.54*crtmean_w5 ' fit13 <- growth(pgmodel13, dset_mod, missing = "ml.x") summary(fit13, standardized = T) #### cognitive domains #### #visuo dom_visL <- ' #domain parameters Ivis =~ Iblkdes + Imatreas + Ispantot Svis =~ Sblkdes + Smatreas + Sspantot #indicator as scaling reference: loading=1, int=0 Iblkdes ~ 0*1 Sblkdes ~ 0*1 # negative residual variances fixed at 0 Sspantot ~~ 0*Sspantot # estimate covariances between task parameters (not the fixed ones) Imatreas ~~ Smatreas Iblkdes ~~ Sblkdes # Ispantot ~~ Sspantot ' fitVisL <- growth(model = c(pgmodel1, pgmodel2, pgmodel3, #growth curves dom_visL), #domain parameters dset_mod, missing = "ml.x") fitmeasures(fitVisL, c("cfi", "tli", "RMSEA", "SRMR")) summary(fitVisL, standardized = T) #crystal dom_cryL <- ' #domain parameters Icryst =~ Iwtar + Inart + Ivftot Scryst =~ Swtar + Snart + Svftot #indicator as scaling reference: loading=1, int=0 Iwtar ~ 0*1 Swtar ~ 0*1 # negative residual variances fixed at 0 Snart ~~ 0*Snart Swtar ~~ 0*Swtar #estimate covariance within task parameters # Inart ~~ Snart # Iwtar ~~ Swtar Ivftot ~~ Svftot ' fitCryL <- growth(model = c(pgmodel4, pgmodel5, pgmodel6, #growth curves dom_cryL), #domain parameters dset_mod, missing = "ml.x") fitmeasures(fitCryL, c("cfi", "tli", "RMSEA", "SRMR")) summary(fitCryL, standardized =T) #memory dom_memL <-' #Domain parameters Imem =~ Ilmtotal + Ivpatotal + Idigback Smem =~ Slmtotal + Svpatotal + Sdigback #indicator as scaling reference: loading=1, int=0 Ilmtotal ~ 0*1 Slmtotal ~ 0*1 # negative residual variances fixed at 0 Slmtotal ~~ 0*Slmtotal #estimate covariance within task parameters Ivpatotal ~~ Svpatotal Idigback ~~ Sdigback # Ilmtotal ~~ Slmtotal ' fitMemL <- growth(model = c(pgmodel7, pgmodel8, pgmodel9, # growth curves dom_memL), #domain parameters dset_mod, missing = "ml.x") fitmeasures(fitMemL, c("cfi", "tli", "RMSEA", "SRMR")) summary(fitMemL, standardized = T) #speed dom_speL <- ' #Domain parameters Ispeed =~ Iittotal + Idigsym + Isymsear + Icrtmean Sspeed =~ Sittotal + Sdigsym + Ssymsear + Scrtmean #indicator as scaling reference: loading=1, int=0 Iittotal ~ 0*1 Sittotal ~ 0*1 # negative residual variances fixed at 0 Ssymsear~~ 0*Ssymsear #estimate covariance within task parameters #Isymsear ~~ Ssymsear Idigsym ~~ Sdigsym Iittotal ~~ Sittotal Icrtmean ~~ Scrtmean ' fitSpeL <- growth(model = c(pgmodel10, pgmodel11, pgmodel12, pgmodel13, # growth curves dom_speL), #domain parameters dset_mod, missing = "ml.x") fitmeasures(fitSpeL, c("cfi", "tli", "RMSEA", "SRMR")) summary(fitSpeL, standardized = T) #### General Cognitive Ability #### general_4p <- ' # General parameters Ig =~ Ivis + Icryst + Imem + Ispeed Sg =~ Svis + Scryst + Smem + Sspeed #indicator as scaling reference: loading=1, int=0 Ivis~ 0*1 Svis~ 0*1 #estimate covariance between domain parameters Ivis ~~ Svis Icryst ~~ Scryst Imem ~~ Smem Ispeed ~~ Sspeed ' options(max.print=10000) fitGen_4p <- growth(model = c(pgmodel1, pgmodel2, pgmodel3, pgmodel4, pgmodel5, pgmodel6, pgmodel7, pgmodel8, pgmodel9, pgmodel10, pgmodel11, pgmodel12, pgmodel13, dom_visL, dom_cryL, dom_memL, dom_speL, general_4p), dset_mod, missing = "ml.x") fitmeasures(fitGen_4p, c("cfi", "tli", "RMSEA", "SRMR")) summary(fitGen_4p, standardized = T) #to get fully standardized results including p-values use standardizedSolution(fitGen_4p)