Load libraries

library(plyr)
library(dplyr)
library(psych)
library(lmtest)  # for bptest (testing for heteroscedasticity)
library(jtools)  # summ.lm (summary of regression); effectplot (scatterplot)
library(lme4)    # to use lmer for mixed models (testing between and within effects)
library(nortest) # to run anderson-darling normality of residuals test (ad.test)
library(car)     #for calculating VIF

Load data

data <- read.csv("/Users/emurbanw/OneDrive-MichiganMedicine/OneDriveDocuments/EFDC/Incubator/Wei_Zhao_PSID/Merged_2014_2019_CDS_genomic_data.3.01.23.csv", header=TRUE)

Subpopulations & mean-centering age

# RQ2_BPI: N = 540; MAge: 10.4, SD = 3.61, range: 4-17
RQ2_BPI <- data[ which(data$RQ2_BPI == 1),]
describe(RQ2_BPI$BPI_Age)
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 540 10.4 3.61     10   10.29 4.45   4  17    13 0.21     -1.1 0.16
MAge <- mean(RQ2_BPI$BPI_Age)
RQ2_BPI$MAge_RQ2_BPI <- RQ2_BPI$BPI_Age - MAge
describe(RQ2_BPI$MAge_RQ2_BPI)
##    vars   n mean   sd median trimmed  mad  min max range skew kurtosis   se
## X1    1 540    0 3.61   -0.4    -0.1 4.45 -6.4 6.6    13 0.21     -1.1 0.16
RQ2_BPI$MAge_RQ2_BPI_sq <- RQ2_BPI$MAge_RQ2_BPI ^ 2

# RQ2_CDI: N = 278; MAge: 14.28, SD = 1.58, range: 12-17
RQ2_CDI <- data[ which(data$RQ2_CDI == 1),]
describe(RQ2_CDI$CDI_Age)
##    vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 278 14.28 1.58     14   14.23 1.48  12  17     5 0.15    -1.11 0.09
MAge <- mean(RQ2_CDI$CDI_Age)
RQ2_CDI$MAge_RQ2_CDI <- RQ2_CDI$CDI_Age - MAge
describe(RQ2_CDI$MAge_RQ2_CDI)
##    vars   n mean   sd median trimmed  mad   min  max range skew kurtosis   se
## X1    1 278    0 1.58  -0.28   -0.05 1.48 -2.28 2.72     5 0.15    -1.11 0.09
RQ2_CDI$MAge_RQ2_CDI_sq <- RQ2_CDI$MAge_RQ2_CDI ^ 2

Analyses

BPI

# Looking at bivariate relationship
cor.test(RQ2_BPI$z_new_PRS, RQ2_BPI$BPI)
## 
##  Pearson's product-moment correlation
## 
## data:  RQ2_BPI$z_new_PRS and RQ2_BPI$BPI
## t = 2.6337, df = 538, p-value = 0.008689
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.02871752 0.19533975
## sample estimates:
##       cor 
## 0.1128217
# running the models
child_BPI <- lmer(BPI ~ z_new_PRS + MAge_RQ2_BPI + MAge_RQ2_BPI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_BPI + (1|Fam_ID_14), data = RQ2_BPI)

parent_BPI <- lmer(BPI ~ bioparent_z_new_PRS + MAge_RQ2_BPI + MAge_RQ2_BPI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_BPI + (1|Fam_ID_14), data = RQ2_BPI)

pair_BPI <- lmer(BPI ~ z_new_PRS + bioparent_z_new_PRS + MAge_RQ2_BPI + MAge_RQ2_BPI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_BPI + (1|Fam_ID_14), data = RQ2_BPI)

# examining results
summ(child_BPI, digits=3) # child: p = .01
## MODEL INFO:
## Observations: 540
## Dependent Variable: BPI
## Type: Mixed effects linear regression 
## 
## MODEL FIT:
## AIC = 2762.327, BIC = 2813.826
## Pseudo-R² (fixed effects) = 0.037
## Pseudo-R² (total) = 0.406 
## 
## FIXED EFFECTS:
## -----------------------------------------------------------------
##                           Est.    S.E.   t val.      d.f.       p
## --------------------- -------- ------- -------- --------- -------
## (Intercept)              3.567   0.245   14.543   512.506   0.000
## z_new_PRS                0.369   0.145    2.546   521.646   0.011
## MAge_RQ2_BPI             0.049   0.039    1.266   529.793   0.206
## MAge_RQ2_BPI_sq         -0.035   0.011   -3.232   494.432   0.001
## Sex2                    -0.278   0.263   -1.058   487.944   0.291
## PC1                     -6.058   5.522   -1.097   407.870   0.273
## PC2                      0.304   3.959    0.077   366.457   0.939
## PC3                     -3.516   4.277   -0.822   424.117   0.411
## PC4                      3.963   4.356    0.910   421.553   0.363
## Which_BPI                1.769   1.963    0.901   513.555   0.368
## -----------------------------------------------------------------
## 
## p values calculated using Satterthwaite d.f.
## 
## RANDOM EFFECTS:
## -------------------------------------
##    Group      Parameter    Std. Dev. 
## ----------- ------------- -----------
##  Fam_ID_14   (Intercept)     1.986   
##  Residual                    2.518   
## -------------------------------------
## 
## Grouping variables:
## ------------------------------
##    Group     # groups    ICC  
## ----------- ---------- -------
##  Fam_ID_14     341      0.384 
## ------------------------------
summ(parent_BPI, digits=3) # parent: p = .003
## MODEL INFO:
## Observations: 540
## Dependent Variable: BPI
## Type: Mixed effects linear regression 
## 
## MODEL FIT:
## AIC = 2759.735, BIC = 2811.234
## Pseudo-R² (fixed effects) = 0.045
## Pseudo-R² (total) = 0.421 
## 
## FIXED EFFECTS:
## ---------------------------------------------------------------------
##                               Est.    S.E.   t val.      d.f.       p
## ------------------------- -------- ------- -------- --------- -------
## (Intercept)                  3.613   0.245   14.746   514.469   0.000
## bioparent_z_new_PRS          0.461   0.154    2.989   342.581   0.003
## MAge_RQ2_BPI                 0.051   0.039    1.323   529.707   0.186
## MAge_RQ2_BPI_sq             -0.036   0.011   -3.313   493.814   0.001
## Sex2                        -0.277   0.262   -1.056   486.688   0.291
## PC1                         -5.777   5.517   -1.047   410.629   0.296
## PC2                          0.174   3.966    0.044   374.051   0.965
## PC3                         -3.057   4.280   -0.714   425.779   0.475
## PC4                          4.043   4.358    0.928   423.733   0.354
## Which_BPI                    1.190   1.970    0.604   518.682   0.546
## ---------------------------------------------------------------------
## 
## p values calculated using Satterthwaite d.f.
## 
## RANDOM EFFECTS:
## -------------------------------------
##    Group      Parameter    Std. Dev. 
## ----------- ------------- -----------
##  Fam_ID_14   (Intercept)     2.013   
##  Residual                    2.495   
## -------------------------------------
## 
## Grouping variables:
## ------------------------------
##    Group     # groups    ICC  
## ----------- ---------- -------
##  Fam_ID_14     341      0.394 
## ------------------------------
summ(pair_BPI, digits=3) # child: p = .19; parent: p = .04
## MODEL INFO:
## Observations: 540
## Dependent Variable: BPI
## Type: Mixed effects linear regression 
## 
## MODEL FIT:
## AIC = 2761.811, BIC = 2817.601
## Pseudo-R² (fixed effects) = 0.047
## Pseudo-R² (total) = 0.418 
## 
## FIXED EFFECTS:
## ---------------------------------------------------------------------
##                               Est.    S.E.   t val.      d.f.       p
## ------------------------- -------- ------- -------- --------- -------
## (Intercept)                  3.596   0.245   14.676   512.572   0.000
## z_new_PRS                    0.214   0.163    1.313   528.881   0.190
## bioparent_z_new_PRS          0.356   0.173    2.050   388.081   0.041
## MAge_RQ2_BPI                 0.049   0.039    1.250   528.744   0.212
## MAge_RQ2_BPI_sq             -0.036   0.011   -3.335   493.367   0.001
## Sex2                        -0.272   0.262   -1.039   486.757   0.299
## PC1                         -6.252   5.515   -1.134   409.114   0.258
## PC2                          0.014   3.958    0.004   370.587   0.997
## PC3                         -3.227   4.272   -0.755   423.906   0.450
## PC4                          4.225   4.352    0.971   422.424   0.332
## Which_BPI                    1.320   1.970    0.670   516.783   0.503
## ---------------------------------------------------------------------
## 
## p values calculated using Satterthwaite d.f.
## 
## RANDOM EFFECTS:
## -------------------------------------
##    Group      Parameter    Std. Dev. 
## ----------- ------------- -----------
##  Fam_ID_14   (Intercept)     1.996   
##  Residual                    2.502   
## -------------------------------------
## 
## Grouping variables:
## ------------------------------
##    Group     # groups    ICC  
## ----------- ---------- -------
##  Fam_ID_14     341      0.389 
## ------------------------------
#Find change between models that don't vs. do include PRS_parent
anova(child_BPI, pair_BPI) # p = .04
## refitting model(s) with ML (instead of REML)
## Data: RQ2_BPI
## Models:
## child_BPI: BPI ~ z_new_PRS + MAge_RQ2_BPI + MAge_RQ2_BPI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_BPI + (1 | Fam_ID_14)
## pair_BPI: BPI ~ z_new_PRS + bioparent_z_new_PRS + MAge_RQ2_BPI + MAge_RQ2_BPI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_BPI + (1 | Fam_ID_14)
##           npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## child_BPI   12 2768.1 2819.6 -1372.0   2744.1                       
## pair_BPI    13 2765.8 2821.6 -1369.9   2739.8 4.2494  1    0.03927 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(parent_BPI, pair_BPI) # p = .18
## refitting model(s) with ML (instead of REML)
## Data: RQ2_BPI
## Models:
## parent_BPI: BPI ~ bioparent_z_new_PRS + MAge_RQ2_BPI + MAge_RQ2_BPI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_BPI + (1 | Fam_ID_14)
## pair_BPI: BPI ~ z_new_PRS + bioparent_z_new_PRS + MAge_RQ2_BPI + MAge_RQ2_BPI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_BPI + (1 | Fam_ID_14)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## parent_BPI   12 2765.6 2817.1 -1370.8   2741.6                     
## pair_BPI     13 2765.8 2821.6 -1369.9   2739.8 1.7673  1     0.1837
# Scatter plot
effect_plot(child_BPI, pred = z_new_PRS, interval = FALSE, plot.points = TRUE, x.label = "Child PRS score", y.label = "BPI Internalizing Score", main.title = "Unadjusted Assocation Between Child PRS and Child BPI")

effect_plot(parent_BPI, pred = bioparent_z_new_PRS, interval = FALSE, plot.points = TRUE, x.label = "Parent PRS score", y.label = "BPI Internalizing Score", main.title = "Unadjusted Assocation Between Parent PRS and Child BPI")

effect_plot(pair_BPI, pred = z_new_PRS, interval = FALSE, plot.points = TRUE, x.label = "Child PRS score", y.label = "BPI Internalizing Score", main.title = "Adjusted Assocation Between Child PRS and Child BPI")

effect_plot(pair_BPI, pred = bioparent_z_new_PRS, interval = FALSE, plot.points = TRUE, x.label = "Parent PRS score", y.label = "BPI Internalizing Score", main.title = "Adjusted Assocation Between Parent PRS and Child BPI")

CDI

# Looking at bivariate relationship
cor.test(RQ2_CDI$z_new_PRS, RQ2_CDI$CDI)
## 
##  Pearson's product-moment correlation
## 
## data:  RQ2_CDI$z_new_PRS and RQ2_CDI$CDI
## t = 1.5876, df = 276, p-value = 0.1135
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.02276771  0.21041822
## sample estimates:
##        cor 
## 0.09513009
# running the models
child_CDI <- lmer(CDI ~ z_new_PRS + MAge_RQ2_CDI + MAge_RQ2_CDI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_CDI + (1|Fam_ID_14), data = RQ2_CDI)

parent_CDI <- lmer(CDI ~ bioparent_z_new_PRS + MAge_RQ2_CDI + MAge_RQ2_CDI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_CDI + (1|Fam_ID_14), data = RQ2_CDI)
## boundary (singular) fit: see help('isSingular')
pair_CDI <- lmer(CDI ~ z_new_PRS + bioparent_z_new_PRS + MAge_RQ2_CDI + MAge_RQ2_CDI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_CDI + (1|Fam_ID_14), data = RQ2_CDI)
## boundary (singular) fit: see help('isSingular')
# examining results
summ(child_CDI, digits=3) # child: p = .10
## MODEL INFO:
## Observations: 278
## Dependent Variable: CDI
## Type: Mixed effects linear regression 
## 
## MODEL FIT:
## AIC = 1298.961, BIC = 1342.492
## Pseudo-R² (fixed effects) = 0.093
## Pseudo-R² (total) = 0.095 
## 
## FIXED EFFECTS:
## -----------------------------------------------------------------
##                           Est.    S.E.   t val.      d.f.       p
## --------------------- -------- ------- -------- --------- -------
## (Intercept)              3.463   0.301   11.500   263.402   0.000
## z_new_PRS                0.254   0.155    1.641   242.682   0.102
## MAge_RQ2_CDI             0.186   0.097    1.910   266.380   0.057
## MAge_RQ2_CDI_sq         -0.040   0.065   -0.614   261.347   0.539
## Sex2                    -1.145   0.312   -3.671   267.902   0.000
## PC1                      3.797   5.248    0.724   224.131   0.470
## PC2                      5.390   4.521    1.192   190.957   0.235
## PC3                     -5.370   4.938   -1.087   216.411   0.278
## PC4                      0.552   5.605    0.099   209.405   0.922
## Which_CDI                0.732   0.306    2.393   260.527   0.017
## -----------------------------------------------------------------
## 
## p values calculated using Satterthwaite d.f.
## 
## RANDOM EFFECTS:
## -------------------------------------
##    Group      Parameter    Std. Dev. 
## ----------- ------------- -----------
##  Fam_ID_14   (Intercept)     0.104   
##  Residual                    2.477   
## -------------------------------------
## 
## Grouping variables:
## ------------------------------
##    Group     # groups    ICC  
## ----------- ---------- -------
##  Fam_ID_14     211      0.002 
## ------------------------------
summ(parent_CDI, digits=3) # parent: p = .25
## MODEL INFO:
## Observations: 278
## Dependent Variable: CDI
## Type: Mixed effects linear regression 
## 
## MODEL FIT:
## AIC = 1300.396, BIC = 1343.928
## Pseudo-R² (fixed effects) = 0.089
## Pseudo-R² (total) = 0.089 
## 
## FIXED EFFECTS:
## ---------------------------------------------------------------------
##                               Est.    S.E.   t val.      d.f.       p
## ------------------------- -------- ------- -------- --------- -------
## (Intercept)                  3.507   0.300   11.693   268.000   0.000
## bioparent_z_new_PRS          0.171   0.148    1.156   268.000   0.249
## MAge_RQ2_CDI                 0.194   0.097    1.999   268.000   0.047
## MAge_RQ2_CDI_sq             -0.039   0.065   -0.606   268.000   0.545
## Sex2                        -1.141   0.313   -3.645   268.000   0.000
## PC1                          4.109   5.253    0.782   268.000   0.435
## PC2                          5.320   4.538    1.172   268.000   0.242
## PC3                         -5.326   4.950   -1.076   268.000   0.283
## PC4                          0.378   5.619    0.067   268.000   0.946
## Which_CDI                    0.684   0.305    2.245   268.000   0.026
## ---------------------------------------------------------------------
## 
## p values calculated using Satterthwaite d.f.
## 
## RANDOM EFFECTS:
## -------------------------------------
##    Group      Parameter    Std. Dev. 
## ----------- ------------- -----------
##  Fam_ID_14   (Intercept)     0.000   
##  Residual                    2.485   
## -------------------------------------
## 
## Grouping variables:
## ------------------------------
##    Group     # groups    ICC  
## ----------- ---------- -------
##  Fam_ID_14     211      0.000 
## ------------------------------
summ(pair_CDI, digits=3) # child: p = .22; parent: p = .66
## MODEL INFO:
## Observations: 278
## Dependent Variable: CDI
## Type: Mixed effects linear regression 
## 
## MODEL FIT:
## AIC = 1302.501, BIC = 1349.660
## Pseudo-R² (fixed effects) = 0.094
## Pseudo-R² (total) = 0.094 
## 
## FIXED EFFECTS:
## ---------------------------------------------------------------------
##                               Est.    S.E.   t val.      d.f.       p
## ------------------------- -------- ------- -------- --------- -------
## (Intercept)                  3.464   0.302   11.486   267.000   0.000
## z_new_PRS                    0.218   0.176    1.239   267.000   0.216
## bioparent_z_new_PRS          0.073   0.168    0.437   267.000   0.663
## MAge_RQ2_CDI                 0.185   0.098    1.895   267.000   0.059
## MAge_RQ2_CDI_sq             -0.039   0.065   -0.597   267.000   0.551
## Sex2                        -1.137   0.313   -3.638   267.000   0.000
## PC1                          3.727   5.257    0.709   267.000   0.479
## PC2                          5.282   4.533    1.165   267.000   0.245
## PC3                         -5.427   4.945   -1.097   267.000   0.273
## PC4                          0.697   5.619    0.124   267.000   0.901
## Which_CDI                    0.728   0.307    2.376   267.000   0.018
## ---------------------------------------------------------------------
## 
## p values calculated using Satterthwaite d.f.
## 
## RANDOM EFFECTS:
## -------------------------------------
##    Group      Parameter    Std. Dev. 
## ----------- ------------- -----------
##  Fam_ID_14   (Intercept)     0.000   
##  Residual                    2.483   
## -------------------------------------
## 
## Grouping variables:
## ------------------------------
##    Group     # groups    ICC  
## ----------- ---------- -------
##  Fam_ID_14     211      0.000 
## ------------------------------
#Find change between models that don't vs. do include PRS_parent
anova(child_CDI, pair_CDI) # p = .66
## refitting model(s) with ML (instead of REML)
## Data: RQ2_CDI
## Models:
## child_CDI: CDI ~ z_new_PRS + MAge_RQ2_CDI + MAge_RQ2_CDI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_CDI + (1 | Fam_ID_14)
## pair_CDI: CDI ~ z_new_PRS + bioparent_z_new_PRS + MAge_RQ2_CDI + MAge_RQ2_CDI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_CDI + (1 | Fam_ID_14)
##           npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## child_CDI   12 1307.5 1351.0 -641.74   1283.5                     
## pair_CDI    13 1309.3 1356.4 -641.64   1283.3 0.1985  1     0.6559
anova(parent_CDI, pair_CDI) # p = .21
## refitting model(s) with ML (instead of REML)
## Data: RQ2_CDI
## Models:
## parent_CDI: CDI ~ bioparent_z_new_PRS + MAge_RQ2_CDI + MAge_RQ2_CDI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_CDI + (1 | Fam_ID_14)
## pair_CDI: CDI ~ z_new_PRS + bioparent_z_new_PRS + MAge_RQ2_CDI + MAge_RQ2_CDI_sq + Sex2 + PC1 + PC2 + PC3 + PC4 + Which_CDI + (1 | Fam_ID_14)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## parent_CDI   12 1308.9 1352.4 -642.44   1284.9                     
## pair_CDI     13 1309.3 1356.4 -641.64   1283.3 1.5935  1     0.2068
# Scatter plot
effect_plot(child_CDI, pred = z_new_PRS, interval = FALSE, plot.points = TRUE, x.label = "Child PRS score", y.label = "CDI Internalizing Score", main.title = "Unadjusted Assocation Between Child PRS and Child CDI")

effect_plot(parent_CDI, pred = bioparent_z_new_PRS, interval = FALSE, plot.points = TRUE, x.label = "Parent PRS score", y.label = "CDI Internalizing Score", main.title = "Unadjusted Assocation Between Parent PRS and Child CDI")

effect_plot(pair_CDI, pred = z_new_PRS, interval = FALSE, plot.points = TRUE, x.label = "Child PRS score", y.label = "CDI Internalizing Score", main.title = "Adjusted Assocation Between Child PRS and Child CDI")

effect_plot(pair_CDI, pred = bioparent_z_new_PRS, interval = FALSE, plot.points = TRUE, x.label = "Parent PRS score", y.label = "CDI Internalizing Score", main.title = "Adjusted Assocation Between Parent PRS and Child CDI")