## White background for combined plotswhite_bg <-function(gg) {grobTree(rectGrob(gp =gpar(fill ='white', lwd =0)), gg) %>% invisible}## Pull out the legend of a plotg_legend <-function(p){ tmp <-ggplot_gtable(ggplot_build(p)) leg <-which(sapply(tmp$grobs, function(x) x$name) =="guide-box") tmp$grobs[[leg]]}
## Descriptive statistics of hyperactivity and inattentiondata %>%select(all_of(var_names) ) %>%rename_at(all_of(var_names),~ long_var_names[.] ) %>% descr_and_cor %>%kable2(register_output =TRUE,caption ='Descriptive statistics for hyperactivity and inattention (both twins)',footnote ='Pearson correlation with standard error is presented on the right',align ='c',digits =2 ) %>% kable_styling
Descriptive statistics for hyperactivity and inattention (both twins)
n
M
SD
median
range
1
2
3
1. Hyperactivity, early childhood
724
0.83
0.40
0.80
0-2
2. Inattention, early childhood
724
0.63
0.35
0.67
0-2
0.703 (0.000)
3. Hyperactivity, primary school
722
0.45
0.45
0.32
0-2
0.309 (0.000)
0.263 (0.000)
4. Inattention, primary school
722
0.77
0.55
0.67
0-2
0.244 (0.000)
0.259 (0.000)
0.723 (0.000)
Note:
Pearson correlation with standard error is presented on the right
## Distribution of untransformed ADHD scoresHs1 <-list()for (i in var_names) { Hs1[[i]] <- data %>%my_hist( i,title = long_var_names[i] )}grid.arrange(grobs = Hs1,ncol =2) %>% white_bg
Code
## Sqrt-transformation of ADHD scores and adjustment for sexfor (i in var_names) { data[, i] <-lm(sqrt(data[, i]) ~ data$gender, na.action = na.exclude ) %>% resid}
Code
## Distribution of adjusted and sqrt-transformed ADHD scoresHs2 <-list()for (i in var_names) { Hs2[[i]] <- data %>%my_hist( i,title = long_var_names[i] )}grid.arrange(grobs = Hs2,ncol =2) %>% white_bg
## Run the models# models <- models_noboot# # for (m in names(models)) {# print(paste('Processing', m))# M <- models[[m]]# M <- M %>%# mxBootstrap(# replications = 10000# )# # models[[m]] <- M# }# ## Save to cache# save(# models,# file = 'cache/boot_transmission_models.RData'# )## Load from cacheload('cache/boot_transmission_models.RData')
Parameter estimates
Code
## Extract the estimates with bootstrapped 95% CIs# par_table <- list()# # for (i in names(models)) {# print(paste('Processing', i))# par_table[[i]] <- cbind(# est = mxEval(# output,# models[[i]]# ),# mxBootstrapEval(# output,# models[[i]],# bq = c(0.025, 0.975)# )# ) %>%# `row.names<-`(# c(# 'g_tw',# 'g_m',# 'g_f',# 'R^2',# 'R^2_vt',# 'R^2_gn'# )# )# }# ## Save to cache# save(# par_table,# file = 'cache/par_table.RData'# )## Load from cacheload('cache/par_table.RData')
Code
## Function to extract and return the key parametersget_model_params <-function(model) { model %>% summary %>% .$CI %>%select( estimate,CI_low = lbound,CI_up = ubound, ) %>%mutate(SE =sapply(row.names(.), mxSE, model, silent =TRUE ),.after = estimate, ) %>%mutate(sig =ifelse(round(CI_up, 3) <0|round(CI_low, 3) >0, '*', '') ) %>%mutate(`95% CI`=sprintf('[%0.3f; %0.3f]', CI_low, CI_up ) ) %>%select( estimate, SE, `95% CI`, sig ) %>%`row.names<-`(c('g_tw','g_m','g_f','R^2','R^2_vt','R^2_gn' ) )}
Code
## The tables with parameter estimates## Table 3 and Table 4## Function that makes a tablepar_tables <- par_table %>%lapply(function(x) x %>% as.data.frame %>% rownames_to_column %>%mutate(sig =ifelse(round(`97.5%`, 3) <0|round(`2.5%`, 3) >0, '*', '') ) %>%mutate(`95% CI`=sprintf('[%0.3f; %0.3f]',`2.5%`,`97.5%` ) ) %>% column_to_rownames %>%select(est = V1, sig, `95% CI` ) )## ADHD-PGSlist('Early childhood'=mlth.data.frame(Hyperactivity = par_tables[['ADHD_PRS, Hyp1']],Inattention = par_tables[['ADHD_PRS, Ina1']] ),'Primary school'=mlth.data.frame(Hyperactivity = par_tables[['ADHD_PRS, Hyp2']],Inattention = par_tables[['ADHD_PRS, Ina2']] )) %>%kable2(register_output =TRUE,caption ='The estimates of direct ADHD PRS effects and explained variance',footnote ='est = estimate; SE = standard error; 95% CI = 95% confidence interval; sig = statistical significance based on 95% CI',align ='c',digits =3 ) %>% kable_styling
The estimates of direct ADHD PRS effects and explained variance
Hyperactivity
Inattention
est
sig
95% CI
est
sig
95% CI
Early childhood
g_tw
0.016
[-0.091; 0.116]
0.002
[-0.100; 0.099]
g_m
0.056
[-0.056; 0.161]
0.000
[-0.106; 0.105]
g_f
-0.038
[-0.164; 0.092]
-0.055
[-0.180; 0.068]
R^2
0.005
[0.000; 0.016]
0.003
[0.000; 0.010]
R^2_vt
0.000
[0.000; 0.003]
0.000
[0.000; 0.000]
R^2_gn
0.005
[-0.006; 0.026]
0.003
[-0.006; 0.019]
Primary school
g_tw
0.157
*
[0.055; 0.249]
0.120
*
[0.024; 0.217]
g_m
-0.089
[-0.187; 0.006]
-0.026
[-0.128; 0.072]
g_f
-0.015
[-0.114; 0.097]
-0.009
[-0.110; 0.099]
R^2
0.016
*
[0.003; 0.034]
0.011
*
[0.001; 0.023]
R^2_vt
0.025
*
[0.003; 0.062]
0.014
*
[0.001; 0.047]
R^2_gn
-0.008
[-0.034; 0.011]
-0.003
[-0.026; 0.011]
Note:
est = estimate; SE = standard error; 95% CI = 95% confidence interval; sig = statistical significance based on 95% CI
Code
## EA-PGSlist('Early childhood'=mlth.data.frame(Hyperactivity = par_tables[['Edu_PRS, Hyp1']],Inattention = par_tables[['Edu_PRS, Ina1']] ),'Primary school'=mlth.data.frame(Hyperactivity = par_tables[['Edu_PRS, Hyp2']],Inattention = par_tables[['Edu_PRS, Ina2']] )) %>%kable2(register_output =TRUE,caption ='The estimates of direct EA PRS effects and explained variance',footnote ='est = estimate; SE = standard error; 95% CI = 95% confidence interval; sig = statistical significance based on 95% CI',align ='c',digits =3 ) %>% kable_styling
The estimates of direct EA PRS effects and explained variance
Hyperactivity
Inattention
est
sig
95% CI
est
sig
95% CI
Early childhood
g_tw
-0.129
*
[-0.226; -0.030]
-0.150
*
[-0.251; -0.047]
g_m
-0.001
[-0.100; 0.090]
0.013
[-0.095; 0.116]
g_f
0.004
[-0.103; 0.108]
0.100
*
[0.001; 0.209]
R^2
0.016
*
[0.003; 0.030]
0.014
*
[0.002; 0.032]
R^2_vt
0.017
*
[0.001; 0.051]
0.022
*
[0.002; 0.063]
R^2_gn
0.000
[-0.026; 0.020]
-0.008
[-0.038; 0.011]
Primary school
g_tw
-0.100
*
[-0.196; -0.002]
-0.153
*
[-0.253; -0.051]
g_m
0.020
[-0.083; 0.122]
-0.070
[-0.165; 0.028]
g_f
-0.091
[-0.193; 0.016]
-0.060
[-0.173; 0.064]
R^2
0.026
*
[0.007; 0.049]
0.055
*
[0.028; 0.083]
R^2_vt
0.010
[0.000; 0.038]
0.023
*
[0.003; 0.064]
R^2_gn
0.016
[-0.013; 0.048]
0.032
[-0.014; 0.074]
Note:
est = estimate; SE = standard error; 95% CI = 95% confidence interval; sig = statistical significance based on 95% CI
## Percent of missing data in each transmission modelmodels_noboot %>%sapply( \(m) { NA_MZ <- m$MZ$data@observed %>% is.na %>% sum NA_DZ <- m$DZ$data@observed %>% is.na %>% sum NO_MZ <- m$MZ$data@observed %>% {nrow(.) *ncol(.)} NO_DZ <- m$DZ$data@observed %>% {nrow(.) *ncol(.)} (NA_MZ + NA_DZ) / (NO_MZ + NO_DZ) } ) %>% {100* .} %>% t %>% t %>%`colnames<-`('% of missing') %>%`row.names<-`(c('ADHD-PGS, hyperactivity in early childhood','EA-PGS, hyperactivity in early childhood','ADHD-PGS, inattention in early childhood','EA-PGS, inattention in early childhood','ADHD-PGS, hyperactivity in primary school','EA-PGS, hyperactivity in primary school','ADHD-PGS, inattention in primary school','EA-PGS, inattention in primary school' )) %>%kable2(register_output =TRUE,caption ='Percent of missing data in each transmission model',digits =1,align ='c' ) %>% kable_styling
Percent of missing data in each transmission model