library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(sandwich)
library(clubSandwich)
## Registered S3 method overwritten by 'clubSandwich':
##   method    from    
##   bread.mlm sandwich
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(interactions)
library(sjPlot)

data0 <- read.delim("C:/Users/Lenovo/Dropbox/STATISTICS/R DATA/KNCAST_PDC_RECODE_R2.dat", row.names=1)

vars <- c("BLYCOGC", "PDC24", "NCIIA_24", 
          "BBSEX", "RAN", "INC_IMP_R", 
          "EDU_B", "MSTATUS", "SEFKOR")
data0$complete_flag <- as.integer(complete.cases(data0[, vars]))
table(data0$complete_flag)
## 
##   0   1 
## 110 722
data<- data0[data0$complete_flag == 1, ]


#######################Centering
data$PDC24_c   <- data$PDC24   - mean(data$PDC24,   na.rm = TRUE)
data$PDC12_c   <- data$PDC12   - mean(data$PDC12,   na.rm = TRUE)
data$NCIIA24_c <- data$NCIIA_24 - mean(data$NCIIA_24, na.rm = TRUE)
data$NCIIA12_c <- data$NCIIA_12 - mean(data$NCIIA_12, na.rm = TRUE)
data$INC_IMP_RC <- data$INC_IMP_R - mean(data$INC_IMP_R, na.rm = TRUE)

########################FACTOR
data <- data %>%
  mutate(
    BBSEX   = factor(BBSEX),     # 아동 성별
    EDU_B   = factor(EDU_B),     # 엄마 교육수준
    MSTATUS = factor(MSTATUS),   # 혼인상태
    SEFKOR  = factor(SEFKOR))     # 국적/한국인 여부

data <- data %>%mutate(BBSEX = factor(BBSEX,levels = c(0, 1),labels = c("Female", "Male")),
                       EDU_B = factor(EDU_B,levels = c(0, 1), labels = c("High school or less", "College or higher")),
                       MSTATUS = factor(MSTATUS, levels = c(0, 1),labels = c("Not married", "Married")),
                       RAN = factor(RAN, levels = c(0, 1),labels = c("Control", "Intervention")),
                       SEFKOR = factor(SEFKOR, levels = c(0, 1),labels = c("Non-Korean", "Korean")))


COG_B<-lmer(BLYCOGC~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)
LAN_B<-lmer(BLYLANC~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)
LAN_R_B<-lmer(BLYLREB~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
              data = data,REML = TRUE)
LAN_E_B<-lmer(BLYLEXB~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
              data = data,REML = TRUE)
MOT_B<-lmer(BLYMOTC~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)
MOT_F_B<-lmer(BLYMOTFB~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
              data = data,REML = TRUE)
MOT_G_B<-lmer(BLYMOTGB~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
              data = data,REML = TRUE)
SOC_B<-lmer(BLYSEC~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)
ADP_B<-lmer(BLYADC~PDC24_c+NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)


round(summary(COG_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              94.190     10.485 681.270   8.983    0.000
## PDC24_c                  -1.075      0.371 706.826  -2.895    0.004
## NCIIA24_c                -0.380      0.486 702.093  -0.783    0.434
## BBSEXMale                -3.874      1.006 674.232  -3.851    0.000
## RANIntervention          -0.664      1.035 667.292  -0.642    0.521
## INC_IMP_RC                0.798      0.368 674.053   2.167    0.031
## EDU_BCollege or higher    3.250      1.924 660.320   1.690    0.092
## MSTATUSMarried           -1.070      9.694 687.265  -0.110    0.912
## SEFKORKorean              9.220      3.527 632.420   2.614    0.009
round(summary(LAN_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              67.832     13.156 679.345   5.156    0.000
## PDC24_c                  -1.056      0.456 656.358  -2.314    0.021
## NCIIA24_c                 0.886      0.595 639.799   1.488    0.137
## BBSEXMale                -7.277      1.224 568.597  -5.946    0.000
## RANIntervention          -0.611      1.301 671.693  -0.470    0.639
## INC_IMP_RC                1.901      0.462 674.679   4.111    0.000
## EDU_BCollege or higher    6.676      2.420 668.340   2.759    0.006
## MSTATUSMarried           14.222     12.154 681.941   1.170    0.242
## SEFKORKorean             16.209      4.450 655.019   3.642    0.000
round(summary(LAN_R_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)               4.684      2.349 679.334   1.994    0.047
## PDC24_c                  -0.175      0.082 684.287  -2.125    0.034
## NCIIA24_c                 0.196      0.108 673.182   1.824    0.069
## BBSEXMale                -1.014      0.222 620.091  -4.575    0.000
## RANIntervention          -0.222      0.232 669.613  -0.958    0.338
## INC_IMP_RC                0.329      0.083 673.852   3.987    0.000
## EDU_BCollege or higher    1.103      0.432 665.078   2.554    0.011
## MSTATUSMarried            1.609      2.171 683.154   0.741    0.459
## SEFKORKorean              3.143      0.793 647.126   3.963    0.000
round(summary(LAN_E_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)               4.438      2.623 680.723   1.692    0.091
## PDC24_c                  -0.182      0.091 662.555  -1.998    0.046
## NCIIA24_c                 0.088      0.119 647.284   0.736    0.462
## BBSEXMale                -1.412      0.245 580.430  -5.775    0.000
## RANIntervention           0.047      0.259 673.162   0.180    0.857
## INC_IMP_RC                0.317      0.092 676.161   3.436    0.001
## EDU_BCollege or higher    1.163      0.482 669.817   2.410    0.016
## MSTATUSMarried            3.212      2.424 683.347   1.325    0.186
## SEFKORKorean              2.255      0.887 656.529   2.541    0.011
round(summary(MOT_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              96.718     11.365 686.571   8.510    0.000
## PDC24_c                  -0.697      0.403 707.690  -1.731    0.084
## NCIIA24_c                 1.136      0.526 703.687   2.159    0.031
## BBSEXMale                -3.304      1.090 680.060  -3.032    0.003
## RANIntervention          -0.366      1.122 674.901  -0.326    0.744
## INC_IMP_RC                0.765      0.399 680.547   1.917    0.056
## EDU_BCollege or higher    1.883      2.085 669.055   0.903    0.367
## MSTATUSMarried           -4.830     10.507 691.549  -0.460    0.646
## SEFKORKorean             10.944      3.824 645.473   2.862    0.004
round(summary(MOT_F_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              11.175      1.957 684.640   5.712    0.000
## PDC24_c                  -0.036      0.069 704.534  -0.523    0.601
## NCIIA24_c                 0.172      0.090 699.279   1.906    0.057
## BBSEXMale                -1.083      0.187 669.992  -5.792    0.000
## RANIntervention          -0.099      0.193 673.324  -0.513    0.608
## INC_IMP_RC                0.190      0.069 678.698   2.764    0.006
## EDU_BCollege or higher    0.003      0.359 667.736   0.009    0.993
## MSTATUSMarried           -1.089      1.809 689.433  -0.602    0.547
## SEFKORKorean              1.392      0.659 645.337   2.113    0.035
round(summary(MOT_G_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)               7.741      2.627 691.472   2.947    0.003
## PDC24_c                  -0.185      0.094 712.994  -1.976    0.049
## NCIIA24_c                 0.201      0.123 712.704   1.634    0.103
## BBSEXMale                -0.059      0.255 705.832  -0.230    0.818
## RANIntervention          -0.018      0.259 674.373  -0.070    0.944
## INC_IMP_RC                0.064      0.092 683.238   0.691    0.490
## EDU_BCollege or higher    0.651      0.481 665.277   1.355    0.176
## MSTATUSMarried           -0.523      2.431 698.466  -0.215    0.830
## SEFKORKorean              2.257      0.879 627.246   2.569    0.010
round(summary(SOC_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)             104.670     12.153 685.060   8.613    0.000
## PDC24_c                   0.350      0.422 667.835   0.831    0.406
## NCIIA24_c                 0.947      0.551 654.145   1.721    0.086
## BBSEXMale                -2.646      1.132 593.508  -2.338    0.020
## RANIntervention          -0.350      1.202 678.535  -0.291    0.771
## INC_IMP_RC                1.053      0.427 681.110   2.467    0.014
## EDU_BCollege or higher    4.395      2.235 675.652   1.967    0.050
## MSTATUSMarried            0.724     11.227 687.300   0.065    0.949
## SEFKORKorean              3.747      4.110 664.157   0.912    0.362
round(summary(ADP_B)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              98.189      8.999 680.951  10.911    0.000
## PDC24_c                  -0.470      0.310 643.484  -1.514    0.131
## NCIIA24_c                 0.094      0.405 625.047   0.232    0.817
## BBSEXMale                -4.162      0.831 548.783  -5.008    0.000
## RANIntervention          -0.032      0.890 674.292  -0.036    0.972
## INC_IMP_RC                1.269      0.316 676.732   4.011    0.000
## EDU_BCollege or higher   -0.128      1.656 671.466  -0.077    0.939
## MSTATUSMarried           -4.779      8.313 683.004  -0.575    0.566
## SEFKORKorean              7.513      3.046 660.181   2.466    0.014
COG<-lmer(BLYCOGC~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
          data = data,REML = TRUE)
LAN<-lmer(BLYLANC~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
          data = data,REML = TRUE)
LAN_R<-lmer(BLYLREB~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)
LAN_E<-lmer(BLYLEXB~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)
MOT<-lmer(BLYMOTC~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
          data = data,REML = TRUE)
MOT_F<-lmer(BLYMOTFB~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)
MOT_G<-lmer(BLYMOTGB~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
            data = data,REML = TRUE)
SOC<-lmer(BLYSEC~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
          data = data,REML = TRUE)
ADP<-lmer(BLYADC~PDC24_c+NCIIA24_c+PDC24_c:NCIIA24_c+BBSEX+RAN+INC_IMP_RC+EDU_B+MSTATUS+SEFKOR+(1|SUBJNO),
          data = data,REML = TRUE)

round(summary(COG)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              94.264     10.496 680.429   8.981    0.000
## PDC24_c                  -1.095      0.380 699.914  -2.881    0.004
## NCIIA24_c                -0.378      0.486 701.073  -0.778    0.437
## BBSEXMale                -3.858      1.009 675.546  -3.824    0.000
## RANIntervention          -0.664      1.036 666.417  -0.641    0.522
## INC_IMP_RC                0.800      0.368 673.370   2.171    0.030
## EDU_BCollege or higher    3.259      1.925 659.574   1.693    0.091
## MSTATUSMarried           -1.068      9.700 686.385  -0.110    0.912
## SEFKORKorean              9.198      3.531 631.112   2.605    0.009
## PDC24_c:NCIIA24_c         0.099      0.399 675.375   0.248    0.804
round(summary(LAN)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              68.597     13.159 677.563   5.213    0.000
## PDC24_c                  -1.270      0.462 608.459  -2.749    0.006
## NCIIA24_c                 0.959      0.591 612.353   1.621    0.106
## BBSEXMale                -7.024      1.216 536.721  -5.774    0.000
## RANIntervention          -0.619      1.301 670.220  -0.475    0.635
## INC_IMP_RC                1.925      0.462 672.985   4.163    0.000
## EDU_BCollege or higher    6.754      2.421 667.358   2.790    0.005
## MSTATUSMarried           14.324     12.151 679.434   1.179    0.239
## SEFKORKorean             15.919      4.455 655.382   3.574    0.000
## PDC24_c:NCIIA24_c         1.077      0.481 536.152   2.242    0.025
round(summary(LAN_R)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)               4.808      2.349 678.075   2.047    0.041
## PDC24_c                  -0.207      0.084 658.255  -2.476    0.014
## NCIIA24_c                 0.205      0.107 661.428   1.909    0.057
## BBSEXMale                -0.979      0.221 606.226  -4.425    0.000
## RANIntervention          -0.224      0.232 668.813  -0.963    0.336
## INC_IMP_RC                0.333      0.083 672.842   4.034    0.000
## EDU_BCollege or higher    1.117      0.432 664.774   2.588    0.010
## MSTATUSMarried            1.621      2.170 681.323   0.747    0.455
## SEFKORKorean              3.100      0.794 648.096   3.906    0.000
## PDC24_c:NCIIA24_c         0.173      0.087 605.938   1.978    0.048
round(summary(LAN_E)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)               4.565      2.625 679.305   1.739    0.082
## PDC24_c                  -0.218      0.093 627.206  -2.354    0.019
## NCIIA24_c                 0.096      0.118 630.872   0.807    0.420
## BBSEXMale                -1.376      0.244 562.648  -5.643    0.000
## RANIntervention           0.046      0.260 671.888   0.178    0.859
## INC_IMP_RC                0.321      0.092 674.806   3.477    0.001
## EDU_BCollege or higher    1.175      0.483 668.899   2.434    0.015
## MSTATUSMarried            3.223      2.424 681.404   1.330    0.184
## SEFKORKorean              2.209      0.888 656.440   2.487    0.013
## PDC24_c:NCIIA24_c         0.169      0.096 562.183   1.759    0.079
round(summary(MOT)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              97.914     11.307 683.741   8.659    0.000
## PDC24_c                  -1.027      0.407 690.815  -2.523    0.012
## NCIIA24_c                 1.200      0.521 692.520   2.304    0.022
## BBSEXMale                -2.991      1.079 660.633  -2.773    0.006
## RANIntervention          -0.380      1.116 673.534  -0.340    0.734
## INC_IMP_RC                0.807      0.397 678.388   2.033    0.042
## EDU_BCollege or higher    1.996      2.076 668.725   0.962    0.337
## MSTATUSMarried           -4.750     10.446 687.864  -0.455    0.649
## SEFKORKorean             10.554      3.811 648.734   2.769    0.006
## PDC24_c:NCIIA24_c         1.587      0.426 660.471   3.725    0.000
round(summary(MOT_F)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              11.244      1.956 683.839   5.747    0.000
## PDC24_c                  -0.054      0.071 697.322  -0.770    0.442
## NCIIA24_c                 0.174      0.090 698.631   1.926    0.055
## BBSEXMale                -1.068      0.187 672.144  -5.701    0.000
## RANIntervention          -0.099      0.193 672.383  -0.514    0.607
## INC_IMP_RC                0.192      0.069 677.976   2.794    0.005
## EDU_BCollege or higher    0.012      0.359 666.858   0.032    0.974
## MSTATUSMarried           -1.086      1.808 688.615  -0.601    0.548
## SEFKORKorean              1.371      0.659 643.814   2.082    0.038
## PDC24_c:NCIIA24_c         0.091      0.074 671.991   1.233    0.218
round(summary(MOT_G)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)               8.055      2.609 686.992   3.087    0.002
## PDC24_c                  -0.277      0.095 708.275  -2.914    0.004
## NCIIA24_c                 0.217      0.121 708.700   1.785    0.075
## BBSEXMale                 0.021      0.252 695.084   0.081    0.935
## RANIntervention          -0.022      0.257 672.716  -0.084    0.933
## INC_IMP_RC                0.075      0.092 680.058   0.822    0.411
## EDU_BCollege or higher    0.680      0.478 665.465   1.423    0.155
## MSTATUSMarried           -0.504      2.412 693.107  -0.209    0.834
## SEFKORKorean              2.160      0.876 634.693   2.467    0.014
## PDC24_c:NCIIA24_c         0.416      0.100 694.938   4.170    0.000
round(summary(SOC)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)             104.484     12.159 684.215   8.593    0.000
## PDC24_c                   0.399      0.431 652.289   0.925    0.355
## NCIIA24_c                 0.947      0.551 655.308   1.718    0.086
## BBSEXMale                -2.696      1.136 600.591  -2.373    0.018
## RANIntervention          -0.348      1.202 677.358  -0.290    0.772
## INC_IMP_RC                1.048      0.427 680.186   2.453    0.014
## EDU_BCollege or higher    4.373      2.235 674.486   1.956    0.051
## MSTATUSMarried            0.714     11.228 686.362   0.064    0.949
## SEFKORKorean              3.809      4.111 662.524   0.926    0.355
## PDC24_c:NCIIA24_c        -0.249      0.449 600.256  -0.555    0.579
round(summary(ADP)$coefficients, 3)
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              98.413      9.010 679.767  10.922    0.000
## PDC24_c                  -0.531      0.316 611.053  -1.680    0.093
## NCIIA24_c                 0.104      0.405 614.762   0.257    0.798
## BBSEXMale                -4.090      0.832 541.243  -4.917    0.000
## RANIntervention          -0.034      0.891 672.957  -0.038    0.970
## INC_IMP_RC                1.276      0.317 675.496   4.029    0.000
## EDU_BCollege or higher   -0.104      1.658 670.322  -0.063    0.950
## MSTATUSMarried           -4.758      8.320 681.455  -0.572    0.568
## SEFKORKorean              7.431      3.050 659.263   2.436    0.015
## PDC24_c:NCIIA24_c         0.307      0.329 540.673   0.935    0.350
#MODEL COMPARSION
anova(COG_B,COG)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## COG_B: BLYCOGC ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## COG: BLYCOGC ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##       npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
## COG_B   11 5805.7 5856.1 -2891.8    5783.7                     
## COG     12 5807.6 5862.6 -2891.8    5783.6 0.0624  1     0.8028
anova(LAN_B,LAN)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## LAN_B: BLYLANC ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## LAN: BLYLANC ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##       npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)  
## LAN_B   11 6116.7 6167.1 -3047.3    6094.7                       
## LAN     12 6113.8 6168.8 -3044.9    6089.8 4.8805  1    0.02716 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(LAN_R_B,LAN_R)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## LAN_R_B: BLYLREB ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## LAN_R: BLYLREB ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##         npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)  
## LAN_R_B   11 3635.8 3686.2 -1806.9    3613.8                       
## LAN_R     12 3634.0 3688.9 -1805.0    3610.0 3.8765  1    0.04897 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(LAN_E_B,LAN_E)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## LAN_E_B: BLYLEXB ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## LAN_E: BLYLEXB ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##         npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)  
## LAN_E_B   11 3789.2 3839.6 -1883.6    3767.2                       
## LAN_E     12 3788.1 3843.1 -1882.1    3764.1 3.0448  1    0.08099 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(MOT_B,MOT)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## MOT_B: BLYMOTC ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## MOT: BLYMOTC ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##       npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## MOT_B   11 5921.8 5972.2 -2949.9    5899.8                         
## MOT     12 5910.2 5965.2 -2943.1    5886.2 13.518  1  0.0002363 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(MOT_F_B,MOT_F)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## MOT_F_B: BLYMOTFB ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## MOT_F: BLYMOTFB ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##         npar    AIC    BIC  logLik -2*log(L) Chisq Df Pr(>Chisq)
## MOT_F_B   11 3379.2 3429.6 -1678.6    3357.2                    
## MOT_F     12 3379.7 3434.7 -1677.8    3355.7  1.54  1     0.2146
anova(MOT_G_B,MOT_G)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## MOT_G_B: BLYMOTGB ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## MOT_G: BLYMOTGB ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##         npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## MOT_G_B   11 3816.9 3867.3 -1897.4    3794.9                         
## MOT_G     12 3802.1 3857.1 -1889.0    3778.1 16.769  1  4.222e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(SOC_B,SOC)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## SOC_B: BLYSEC ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## SOC: BLYSEC ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##       npar    AIC    BIC  logLik -2*log(L) Chisq Df Pr(>Chisq)
## SOC_B   11 6002.7 6053.1 -2990.3    5980.7                    
## SOC     12 6004.4 6059.4 -2990.2    5980.4 0.307  1     0.5795
anova(ADP_B,ADP)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## ADP_B: BLYADC ~ PDC24_c + NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
## ADP: BLYADC ~ PDC24_c + NCIIA24_c + PDC24_c:NCIIA24_c + BBSEX + RAN + INC_IMP_RC + EDU_B + MSTATUS + SEFKOR + (1 | SUBJNO)
##       npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
## ADP_B   11 5565.2 5615.6 -2771.6    5543.2                     
## ADP     12 5566.4 5621.3 -2771.2    5542.4 0.8673  1     0.3517
#PICK A POINT PLOT
sd_NCIIA24_c <- sd(data$NCIIA24_c, na.rm = TRUE)
sd_NCIIA12_c <- sd(data$NCIIA12_c, na.rm = TRUE)

mean_NCIIA24_c <- mean(data$NCIIA24_c, na.rm = TRUE)  # ≈ 0
sd_NCIIA24_c   <- sd(data$NCIIA24_c,   na.rm = TRUE)  # 1.14

PLAN <- plot_model(LAN,type  = "pred",terms=c("PDC24_c",paste0("NCIIA24_c [",round(-sd_NCIIA24_c, 2), ", 0, ", round(sd_NCIIA24_c, 2), "]")))
PLAN_R <- plot_model(LAN_R,type  = "pred",terms=c("PDC24_c",paste0("NCIIA24_c [",round(-sd_NCIIA24_c, 2), ", 0, ", round(sd_NCIIA24_c, 2), "]")))
PLAN_E <- plot_model(LAN_E,type  = "pred",terms=c("PDC24_c",paste0("NCIIA24_c [",round(-sd_NCIIA24_c, 2), ", 0, ", round(sd_NCIIA24_c, 2), "]")))

PMOT <- plot_model(MOT,type  = "pred",terms=c("PDC24_c",paste0("NCIIA24_c [",round(-sd_NCIIA24_c, 2), ", 0, ", round(sd_NCIIA24_c, 2), "]")))
PMOT_F <- plot_model(MOT_F,type  = "pred",terms=c("PDC24_c",paste0("NCIIA24_c [",round(-sd_NCIIA24_c, 2), ", 0, ", round(sd_NCIIA24_c, 2), "]")))
PMOT_G <- plot_model(MOT_G,type  = "pred",terms=c("PDC24_c",paste0("NCIIA24_c [",round(-sd_NCIIA24_c, 2), ", 0, ", round(sd_NCIIA24_c, 2), "]")))

library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:sjPlot':
## 
##     set_theme
#DETAIL
PLAN2<-PLAN+ theme_classic(base_size = 12,base_family = "Serif")+labs(x = "Number of PDC", y = "Bayley Language Development",
                                                                      color = "Caregiver RTC",  title = " ") +scale_color_discrete(labels = c("−1 SD","Mean", "+1 SD"))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
PLAN_R2<-PLAN_R+ theme_classic(base_size = 12,base_family = "Serif")+labs(x = "Number of PDC", y = "Bayley Receptive Language Development",
                                                                          color = "Caregiver RTC",  title = " ") +scale_color_discrete(labels = c("−1 SD","Mean", "+1 SD"))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
PLAN_E2<-PLAN_E+ theme_classic(base_size = 12,base_family = "Serif")+labs(x = "Number of PDC", y = "Bayley Expresive Language Development",
                                                                          color = "Caregiver RTC",  title = " ") +scale_color_discrete(labels = c("−1 SD","Mean", "+1 SD"))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
PMOT2<-PMOT+ theme_classic(base_size = 12,base_family = "Serif")+labs(x = "Number of PDC", y = "Bayley Motor Development",
                                                                      color = "Caregiver RTC",  title = " ") +scale_color_discrete(labels = c("−1 SD","Mean", "+1 SD"))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
PMOT_F2<-PMOT_F+ theme_classic(base_size = 12,base_family = "Serif")+labs(x = "Number of PDC", y = "Bayley Fine Motor Development",
                                                                          color = "Caregiver RTC",  title = " ") +scale_color_discrete(labels = c("−1 SD","Mean", "+1 SD"))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
PMOT_G2<-PMOT_G+ theme_classic(base_size = 12,base_family = "Serif")+labs(x = "Number of PDC", y = "Bayley Gross Motor Development", 
                                                                          color = "Caregiver RTC",  title = " ") +scale_color_discrete(labels = c("−1 SD","Mean", "+1 SD"))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
##JN
MOT_JP<-johnson_neyman(MOT, pred=PDC24_c, modx = NCIIA24_c, vmat = NULL, alpha = 0.05,
                       plot = TRUE, control.fdr = FALSE, line.thickness = 0.5,
                       df = "residual", digits = getOption("jtools-digits", 2),
                       critical.t = NULL, sig.color = "#00BFC4", insig.color = "#F8766D",
                       mod.range = NULL, title = "")

MOT_JP_F<-MOT_JP$plot + xlab("Responsiveness to Distress") + ylab("PDC effect on Child Motor Development")

MOT_G_JP<-johnson_neyman(MOT_G, pred=PDC24_c, modx = NCIIA24_c, vmat = NULL, alpha = 0.05,
                         plot = TRUE, control.fdr = FALSE, line.thickness = 0.5,
                         df = "residual", digits = getOption("jtools-digits", 2),
                         critical.t = NULL, sig.color = "#00BFC4", insig.color = "#F8766D",
                         mod.range = NULL, title = "")

MOT_G_JP_F<-MOT_G_JP$plot + xlab("Responsiveness to Distress") + ylab("PDC effect on Gross Motor Development")

LAN_R_JP<-johnson_neyman(LAN_R, pred=PDC24_c, modx = NCIIA24_c, vmat = NULL, alpha = 0.05,
                         plot = TRUE, control.fdr = FALSE, line.thickness = 0.5,
                         df = "residual", digits = getOption("jtools-digits", 2),
                         critical.t = NULL, sig.color = "#00BFC4", insig.color = "#F8766D",
                         mod.range = NULL, title = "")

LAN_R_JP_F<-LAN_R_JP$plot + xlab("Responsiveness to Distress") + ylab("PDC effect on Receptive Language Development")