SB life course

df$sex <- df$n622
df$sc_prenatal <- df$n492
df$smoking_preg <- df$n503
df$birth_weight <- df$n646

df$physical_activity_age7 <- df$n121
df$mon_read_child <- df$n179
df$dad_read_child <- df$n180
df$mon_read_book <- df$n186
df$mon_read_news <- df$n184
df$dad_read_book <- df$n187
df$dad_read_news <- df$n185
df$sc_age7 <- df$n190
df$breast_f <- df$n222
df$oral_ability_age7 <- df$n65
df$awareness_age7 <- df$n67
df$reading_ability_age7 <- df$n68
df$creativity_ability_age7 <- df$n69
df$numbers_ability_age7 <- df$n70
df$bsag_age7 <- df$n455
df$arithmetic_age7 <- df$n90
df$copying_age7 <- df$n457
df$draw_man_age7 <- df$n1840
df$reading_test_age7 <- df$n92

df$mon_walk_child_age11 <- df$n1145
df$dad_walk_child_age11 <- df$n1146
df$sc_age11 <- df$n1687
df$knowledge_age11 <- df$n876
df$numbers_ability_age11 <- df$n877
df$child_use_book_age11 <- df$n878
df$oral_ability_age11 <- df$n879
df$bsag_age11 <- df$n1008
df$verbal_ability_age11 <- df$n914
df$nonverbal_ability_age11 <- df$n917
df$general_ability_age11 <- df$n920
df$reading_ability_age11 <- df$n923
df$maths_ability_age11 <- df$n926
df$copying_age11 <- df$n929
df$child_bor_books_lib_age11 <- df$n931
df$child_read_books_age11 <- df$n935
df$child_read_news_age11 <- df$n936
df$physical_activity_age11 <- df$n941
df$child_writes_stories_age11 <- df$n942
df$child_draws_age11 <- df$n943
df$tv_age11 <- df$n949

df$sc_age16 <- df$n2384
df$iq_ag16 <- df$n1897
df$maths_ability_age16 <-df$n2244
df$english_ability_age16 <- df$n2245
df$language_age16 <- df$n2246
df$science_age16 <- df$n2247
df$practical_ability_age16 <- df$n2248
df$social_std_ability_age16 <- df$n2249
df$school_attendance_age16 <- df$n1722
df$child_read_age16 <- df$n2864
df$play_outdoor_age16 <- df$n2865
df$play_indoor_age16 <- df$n2867
df$swims_age16 <- df$n2866
df$tv_age16 <- df$n2868
df$smoking_age16 <- df$n2887
df$alcohol_age16 <- df$n2888
df$reading_test_age16 <- df$n2928
df$maths_test_age16 <- df$n2930
df$height_age16 <- df$dvht11
df$weight_age16 <- df$dvwt11

df$tv_age23 <- df$n5913_sweep4
df$read_book_age23 <- df$n5914_sweep4
df$read_news_age23 <- df$n5956_sweep4
df$physical_activity_age23 <- df$n5916_sweep4
df$alcohol_age23 <- df$n5920_sweep4
df$qual_age23 <- df$maxqual2_sweep4
df$smoking_age23 <- df$currentn_sweep4
df$income_age23 <- df$famnet_sweep4
df$job_time <- df$fullpart_sweep4
df$malaise_age23 <- df$malaise_sweep4

df$main_activity_age33 <- df$n500520_sweep5
df$sc_age33 <- df$n540056_sweep5
df$smoking_1_age33 <- df$n504262_sweep5
df$smoking_2_age33 <- df$n504265_sweep5
df$alcohol_age33 <- df$n504273_sweep5
df$occ_pa_age33 <- df$n504361_sweep5
df$physical_activity_1_age33 <- df$n504362_sweep5
df$physical_activity_2_age33 <- df$n504363_sweep5
df$qual_age33 <- df$hqual33_sweep5

df$occ_pa_age42 <- df$jdemand2_sweep6
df$smoking_age42 <- df$smoking_sweep6
df$alcohol_age42 <-df$drinks_sweep6
df$physical_activity_1_age42 <- df$exercise_sweep6
df$physical_activity_2_age42 <- df$breathls_sweep6
df$physical_activity_3_age42 <- df$sweat_sweep6
df$sc_age42 <- df$sc_sweep6

df$marital_age46 <- df$nd7ms_sweep7
df$income_age46 <- df$nd7iamtc_sweep7
df$main_activity_age46 <- df$nd7ecact_sweep7
df$qual_age46 <- df$nd7achq1_sweep7
df$computer_home_age46 <- df$n7pchome_sweep7
df$computer_home_freq_age46 <- df$n7hpcuse_sweep7
df$computer_work_age46 <- df$n7pcwork_sweep7
df$computer_work_freq_age46 <- df$n7wpcuse_sweep7
df$internet_leisure <- df$n7intacc_sweep7
df$read_news_age46 <- df$n7pread1_sweep7
df$read_book_age46 <- df$n7pread2_sweep7
df$diseases_age46 <- df$n7lsiany_sweep7
df$smoking_age46 <- df$nd7smoke_sweep7
df$alcohol_age46 <- df$n7drinks_sweep7
df$physical_activity_1_46 <- df$n7exers1_sweep7
df$physical_activity_2_46 <- df$n7breals_sweep7
df$physical_activity_3_46 <- df$n7sweat_sweep7
df$transport_age46 <- df$n7tranrt_sweep7

df$computer_home_age50 <- df$N8PCHOME_sweep8
df$computer_home_freq_age50 <- df$N8HPCUSE_sweep8
df$computer_work_age50 <- df$N8PCWORK_sweep8
df$computer_work_freq_age50 <- df$N8WPCUSE_sweep8
df$general_health_age50 <- df$N8HLTHGN_sweep8
df$diabetes_age50 <- df$N8KHPB03_sweep8
df$hypertension_age50 <- df$N8KHPB09_sweep8
df$words_immed_age50 <- df$N8CFLISN_sweep8
df$animals_named <- df$N8CFANI_sweep8
df$letter_canc_correct <- df$N8CFCOR_sweep8
df$letter_canc_accuracy <- df$N8CFMIS_sweep8
df$words_delayed <- df$N8CFLISD_sweep8
df$smoking_age50 <- df$N8SMOKIG_sweep8
df$alcohol_age50 <- df$N8DRINKS_sweep8
df$physical_activity_1_age50 <- df$N8EXERSE_sweep8
df$physical_activity_2_age50 <- df$N8BREALS_sweep8
df$physical_activity_3_age50 <- df$N8SWEAT_sweep8
df$occ_pa_age50 <- df$N8PHYSWK_sweep8
df$malaise_age50 <- df$ND8MALG_sweep8

df$physical_activity_age55 <- df$N9LEIS01_sweep9
df$internet_use_age55 <- df$N9INTRNT_sweep9
df$marital_age55 <- df$ND9MS_sweep9
df$main_activity_age55 <- df$ND9ECACT_sweep9
df$qual_age55 <- df$ND9HNVQ_sweep9
df$weight_age55 <- df$ND9WGHTK_sweep9
df$bmi_age55 <- df$ND9BMI_sweep9

df$marital_age62 <- df$n10hms_sweep10
df$main_activity_age62 <- df$n10econact2_sweep10
df$letter_can_accuracy <- df$n10cfmis_sweep10
df$animals_named_age62 <- df$n10cfani_sweep10
df$words_delayed_age62 <- df$n10cflisd_sweep10
df$general_health_age62 <- df$n10hlthgen_sweep10
df$diabetes_age62 <- df$n10khlprb04_sweep10
df$hypertension_age62 <- df$n10khlprb06_sweep10
df$dementia_age62 <- df$n10khlprb10_sweep10
df$physical_activity_1_age62 <- df$n10exercise_sweep10
df$physical_activity_2_age62 <- df$n10freqexer_sweep10
df$physical_activity_3_age62 <- df$n10sweat_sweep10
df$smoking_age62 <- df$n10smoking_sweep10
df$read_book_age62 <- df$n10ylq26_sweep10
df$tv_age62 <- df$n10ylq27_sweep10
df$bmi_age62 <- df$nd10bmi_rec_sweep10
df$malaise_age62 <- df$nd10mal_sweep10
df$qual_age62 <- df$nd10hnvq_sweep10

COGNITIVE FUNCTION AT AGE 7

library(broom)
library(dplyr)

df_complete <- df %>%
  filter(arithmetic_age7!=-1 &
           copying_age7!=-1  &
           draw_man_age7!=-1 &
           reading_test_age7!=-1)

# Age 7
# Problem Arithmetic test
df_complete$mean_arithmetic_age7 <- mean(df_complete$arithmetic_age7)
df_complete$sd_arithmetic_age7 <- sd(df_complete$arithmetic_age7)
df_complete$z_arithmetic_age7 <- (df_complete$arithmetic_age7-df_complete$mean_arithmetic_age7)/df_complete$sd_arithmetic_age7
df_complete$z_arithmetic_age7 = scale(df_complete$z_arithmetic_age7)[,1]

hist(df_complete$z_arithmetic_age7)

#Copying designs test
df_complete$mean_copying_age7 <- mean(df_complete$copying_age7)
df_complete$sd_copying_age7 <- sd(df_complete$copying_age7)
df_complete$z_copying_age7 <- (df_complete$copying_age7-df_complete$mean_copying_age7)/df_complete$sd_copying_age7
df_complete$z_copying_age7 = scale(df_complete$z_copying_age7)[,1]

hist(df_complete$z_copying_age7)

#Draw-a-man test
df_complete$mean_draw_man_age7 <- mean(df_complete$draw_man_age7)
df_complete$sd_draw_man_age7 <- sd(df_complete$draw_man_age7)
df_complete$z_draw_man_age7 <- (df_complete$draw_man_age7-df_complete$mean_draw_man_age7)/df_complete$sd_draw_man_age7
df_complete$z_draw_man_age7 = scale(df_complete$z_draw_man_age7)[,1]

hist(df_complete$z_draw_man_age7)

#Southgate Group Reading test
df_complete$mean_reading_test_age7 <- mean(df_complete$reading_test_age7)
df_complete$sd_reading_test_age7 <- sd(df_complete$reading_test_age7)
df_complete$z_reading_test_age7 <- (df_complete$reading_test_age7-df_complete$mean_reading_test_age7)/df_complete$sd_reading_test_age7
df_complete$z_reading_test_age7 = scale(df_complete$z_draw_man_age7)[,1]

hist(df_complete$z_reading_test_age7)

# averaging all z scores
df_complete$global_cog_age7 <- (df_complete$z_arithmetic_age7 + df_complete$z_copying_age7 + df_complete$z_draw_man_age7 + df_complete$z_reading_test_age7) / 4

# creating a z score for the averaged variable
df_complete$z_global_cog_age7 = scale(df_complete$global_cog_age7)[,1]

hist(df_complete$z_global_cog_age7)

# Reference levels
df_complete <- df_complete %>%
  mutate(
    dad_read_score = case_when(
      dad_read_book %in% c("Yes most weeks") ~ 2,
      dad_read_book %in% c("Yes occasionally") ~ 1,
      dad_read_book %in% c("Hardly ever") ~ 0,
      TRUE ~ NA_real_
    ),
    
    mom_read_score = case_when(
      mon_read_child %in% c("Yes,every week") ~ 2,
      mon_read_child %in% c("Yes,occasionally") ~ 1,
      mon_read_child %in% c("Hardly ever") ~ 0,
      TRUE ~ NA_real_
    )
  )

df_complete <- df_complete %>%
  mutate(
    parent_read = case_when(
      dad_read_score == 0 & mom_read_score == 0 ~ 0,
      dad_read_score > 0 | mom_read_score > 0 ~ 1,
      TRUE ~ NA_real_
    )
  )

df_complete <- df_complete %>%
  mutate(
    parent_read_level = pmax(dad_read_score, mom_read_score, na.rm = TRUE),
    parent_read_level = ifelse(
      is.infinite(parent_read_level), NA, parent_read_level
    )
  )


df_complete$parent_read_level <- factor(
  df_complete$parent_read_level,
  levels = c(0, 1, 2),
  labels = c("Hardly ever", "Occasionally", "Weekly"))

table(df_complete$parent_read_level)

## 
##  Hardly ever Occasionally       Weekly 
##         1022         3170         9459

global_age7 <- lm(z_global_cog_age7 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f, data = df_complete)

global_age7 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.401	0.501	-0.800	0.423	-1.382	0.581
sexFemale	0.006	0.017	0.337	0.736	-0.027	0.038
sc_prenatalUnemployed,sick	-0.593	0.484	-1.227	0.220	-1.541	0.355
sc_prenatalI	0.331	0.077	4.295	0.000	0.180	0.482
sc_prenatalII	0.231	0.064	3.620	0.000	0.106	0.356
sc_prenatalIII	0.172	0.058	2.955	0.003	0.058	0.287
sc_prenatalIV	0.089	0.062	1.422	0.155	-0.034	0.211
sc_prenatalV	-0.029	0.064	-0.453	0.650	-0.154	0.096
sc_prenatalStudents	0.254	0.274	0.927	0.354	-0.283	0.791
sc_prenatalDead or away	-0.671	0.683	-0.983	0.326	-2.010	0.668
sc_prenatalRetired	0.263	0.684	0.385	0.700	-1.076	1.603
sc_prenatalSingle,no husbnd	-0.095	0.079	-1.207	0.227	-0.250	0.059
smoking_preg5mths-no change	0.029	0.080	0.365	0.715	-0.128	0.186
smoking_pregNow non-smoker	0.092	0.085	1.074	0.283	-0.076	0.259
smoking_preg1 to 4 now	0.037	0.092	0.407	0.684	-0.143	0.218
smoking_preg5 to 9 now	-0.019	0.096	-0.202	0.840	-0.207	0.169
smoking_preg10 to 14 now	0.038	0.110	0.347	0.729	-0.177	0.254
smoking_preg15 to 19 now	-0.058	0.135	-0.429	0.668	-0.322	0.206
smoking_preg20 to 24 now	-0.233	0.155	-1.504	0.133	-0.537	0.071
smoking_preg25 to 29 now	0.029	0.240	0.122	0.903	-0.441	0.499
smoking_preg30 or more now	0.025	0.269	0.094	0.925	-0.501	0.552
smoking_pregVariable-5mths	0.010	0.087	0.116	0.908	-0.160	0.180
birth_weight	0.000	0.000	3.087	0.002	0.000	0.001
parent_read_levelOccasionally	0.122	0.035	3.461	0.001	0.053	0.191
parent_read_levelWeekly	0.199	0.033	6.040	0.000	0.134	0.263
sc_age7No male head	0.137	0.112	1.226	0.220	-0.082	0.357
sc_age7I	0.578	0.111	5.224	0.000	0.361	0.795
sc_age7II	0.460	0.104	4.438	0.000	0.257	0.664
sc_age7III non-manual	0.432	0.104	4.156	0.000	0.228	0.636
sc_age7III manual	0.187	0.101	1.854	0.064	-0.011	0.385
sc_age7IV non-manual	0.201	0.120	1.678	0.093	-0.034	0.435
sc_age7IV manual	0.097	0.102	0.949	0.342	-0.103	0.298
sc_age7V	-0.042	0.106	-0.396	0.692	-0.249	0.165
breast_fDont know	-0.444	0.493	-0.899	0.368	-1.411	0.523
breast_fNo	-0.323	0.480	-0.673	0.501	-1.265	0.618
breast_fUnder one month	-0.237	0.480	-0.494	0.621	-1.179	0.704
breast_fOver one month	-0.143	0.480	-0.297	0.767	-1.084	0.799

arith_age7 <- lm(z_arithmetic_age7 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f, data = df_complete)

arith_age7 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.416	0.506	-0.821	0.411	-1.407	0.576
sexFemale	-0.088	0.017	-5.197	0.000	-0.121	-0.055
sc_prenatalUnemployed,sick	0.384	0.489	0.786	0.432	-0.573	1.342
sc_prenatalI	0.313	0.078	4.024	0.000	0.161	0.466
sc_prenatalII	0.205	0.064	3.174	0.002	0.078	0.331
sc_prenatalIII	0.162	0.059	2.751	0.006	0.047	0.277
sc_prenatalIV	0.085	0.063	1.344	0.179	-0.039	0.208
sc_prenatalV	0.003	0.065	0.053	0.958	-0.123	0.130
sc_prenatalStudents	0.189	0.277	0.683	0.495	-0.353	0.732
sc_prenatalDead or away	-0.066	0.690	-0.096	0.923	-1.419	1.286
sc_prenatalRetired	-0.128	0.690	-0.186	0.853	-1.482	1.225
sc_prenatalSingle,no husbnd	-0.013	0.080	-0.160	0.872	-0.169	0.144
smoking_preg5mths-no change	-0.010	0.081	-0.126	0.900	-0.169	0.148
smoking_pregNow non-smoker	0.067	0.086	0.782	0.434	-0.101	0.236
smoking_preg1 to 4 now	-0.030	0.093	-0.321	0.748	-0.212	0.152
smoking_preg5 to 9 now	-0.005	0.097	-0.051	0.959	-0.195	0.185
smoking_preg10 to 14 now	-0.006	0.111	-0.056	0.955	-0.224	0.212
smoking_preg15 to 19 now	-0.148	0.136	-1.091	0.275	-0.415	0.118
smoking_preg20 to 24 now	-0.240	0.157	-1.534	0.125	-0.547	0.067
smoking_preg25 to 29 now	-0.533	0.242	-2.199	0.028	-1.008	-0.058
smoking_preg30 or more now	-0.131	0.271	-0.482	0.630	-0.663	0.401
smoking_pregVariable-5mths	-0.010	0.088	-0.117	0.907	-0.182	0.162
birth_weight	0.000	0.000	2.700	0.007	0.000	0.001
parent_read_levelOccasionally	0.094	0.036	2.625	0.009	0.024	0.163
parent_read_levelWeekly	0.155	0.033	4.672	0.000	0.090	0.221
sc_age7No male head	0.020	0.113	0.177	0.859	-0.202	0.242
sc_age7I	0.458	0.112	4.099	0.000	0.239	0.677
sc_age7II	0.332	0.105	3.171	0.002	0.127	0.538
sc_age7III non-manual	0.286	0.105	2.722	0.006	0.080	0.492
sc_age7III manual	0.043	0.102	0.425	0.671	-0.157	0.243
sc_age7IV non-manual	0.099	0.121	0.823	0.411	-0.137	0.336
sc_age7IV manual	-0.032	0.103	-0.310	0.757	-0.235	0.171
sc_age7V	-0.100	0.107	-0.935	0.350	-0.309	0.110
breast_fDont know	-0.271	0.498	-0.544	0.586	-1.248	0.705
breast_fNo	0.002	0.485	0.004	0.996	-0.949	0.953
breast_fUnder one month	0.041	0.485	0.085	0.932	-0.910	0.992
breast_fOver one month	0.100	0.485	0.206	0.837	-0.851	1.051

copying_age7 <- lm(z_copying_age7 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f, data = df_complete)

copying_age7 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.450	0.513	-0.877	0.381	-1.454	0.555
sexFemale	-0.048	0.017	-2.784	0.005	-0.081	-0.014
sc_prenatalUnemployed,sick	0.304	0.495	0.614	0.539	-0.666	1.274
sc_prenatalI	0.220	0.079	2.785	0.005	0.065	0.374
sc_prenatalII	0.103	0.065	1.575	0.115	-0.025	0.231
sc_prenatalIII	0.071	0.060	1.196	0.232	-0.046	0.188
sc_prenatalIV	0.007	0.064	0.117	0.907	-0.118	0.132
sc_prenatalV	-0.116	0.065	-1.771	0.077	-0.244	0.012
sc_prenatalStudents	0.267	0.281	0.952	0.341	-0.283	0.817
sc_prenatalDead or away	-0.503	0.699	-0.719	0.472	-1.874	0.868
sc_prenatalRetired	0.788	0.700	1.126	0.260	-0.584	2.160
sc_prenatalSingle,no husbnd	-0.154	0.081	-1.906	0.057	-0.312	0.004
smoking_preg5mths-no change	0.021	0.082	0.251	0.802	-0.140	0.181
smoking_pregNow non-smoker	0.024	0.087	0.274	0.784	-0.147	0.195
smoking_preg1 to 4 now	0.019	0.094	0.197	0.844	-0.166	0.203
smoking_preg5 to 9 now	-0.073	0.098	-0.745	0.456	-0.266	0.119
smoking_preg10 to 14 now	0.039	0.113	0.351	0.726	-0.181	0.260
smoking_preg15 to 19 now	0.051	0.138	0.372	0.710	-0.219	0.322
smoking_preg20 to 24 now	-0.114	0.159	-0.716	0.474	-0.425	0.198
smoking_preg25 to 29 now	0.069	0.246	0.279	0.780	-0.413	0.550
smoking_preg30 or more now	-0.224	0.275	-0.815	0.415	-0.763	0.315
smoking_pregVariable-5mths	-0.044	0.089	-0.490	0.624	-0.218	0.131
birth_weight	0.000	0.000	2.217	0.027	0.000	0.001
parent_read_levelOccasionally	0.053	0.036	1.459	0.144	-0.018	0.124
parent_read_levelWeekly	0.120	0.034	3.563	0.000	0.054	0.186
sc_age7No male head	0.042	0.115	0.370	0.711	-0.182	0.267
sc_age7I	0.350	0.113	3.091	0.002	0.128	0.572
sc_age7II	0.281	0.106	2.645	0.008	0.073	0.489
sc_age7III non-manual	0.257	0.106	2.410	0.016	0.048	0.465
sc_age7III manual	0.083	0.103	0.801	0.423	-0.120	0.285
sc_age7IV non-manual	0.087	0.122	0.713	0.476	-0.153	0.327
sc_age7IV manual	0.043	0.105	0.413	0.680	-0.162	0.249
sc_age7V	-0.078	0.108	-0.717	0.474	-0.290	0.135
breast_fDont know	0.200	0.505	0.395	0.692	-0.790	1.190
breast_fNo	0.076	0.492	0.155	0.877	-0.888	1.040
breast_fUnder one month	0.165	0.492	0.335	0.738	-0.799	1.129
breast_fOver one month	0.234	0.492	0.476	0.634	-0.730	1.198

draw_age7 <- lm(z_draw_man_age7 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f, data = df_complete)

draw_age7 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.183	0.510	-0.360	0.719	-1.182	0.816
sexFemale	0.076	0.017	4.488	0.000	0.043	0.110
sc_prenatalUnemployed,sick	-1.255	0.492	-2.551	0.011	-2.220	-0.291
sc_prenatalI	0.242	0.078	3.085	0.002	0.088	0.395
sc_prenatalII	0.201	0.065	3.097	0.002	0.074	0.329
sc_prenatalIII	0.148	0.059	2.494	0.013	0.032	0.264
sc_prenatalIV	0.090	0.063	1.420	0.156	-0.034	0.214
sc_prenatalV	0.012	0.065	0.180	0.857	-0.116	0.139
sc_prenatalStudents	0.162	0.279	0.582	0.561	-0.384	0.709
sc_prenatalDead or away	-0.747	0.695	-1.074	0.283	-2.109	0.616
sc_prenatalRetired	0.075	0.696	0.108	0.914	-1.288	1.438
sc_prenatalSingle,no husbnd	-0.063	0.080	-0.784	0.433	-0.220	0.094
smoking_preg5mths-no change	0.040	0.081	0.487	0.626	-0.120	0.199
smoking_pregNow non-smoker	0.095	0.087	1.095	0.274	-0.075	0.265
smoking_preg1 to 4 now	0.063	0.094	0.674	0.500	-0.120	0.246
smoking_preg5 to 9 now	0.009	0.098	0.096	0.924	-0.182	0.201
smoking_preg10 to 14 now	0.042	0.112	0.375	0.708	-0.177	0.261
smoking_preg15 to 19 now	-0.040	0.137	-0.294	0.769	-0.309	0.228
smoking_preg20 to 24 now	-0.181	0.158	-1.148	0.251	-0.490	0.128
smoking_preg25 to 29 now	0.277	0.244	1.136	0.256	-0.201	0.756
smoking_preg30 or more now	0.216	0.273	0.791	0.429	-0.320	0.752
smoking_pregVariable-5mths	0.042	0.088	0.480	0.631	-0.131	0.215
birth_weight	0.000	0.000	2.205	0.027	0.000	0.001
parent_read_levelOccasionally	0.115	0.036	3.188	0.001	0.044	0.185
parent_read_levelWeekly	0.168	0.033	5.006	0.000	0.102	0.233
sc_age7No male head	0.180	0.114	1.576	0.115	-0.044	0.403
sc_age7I	0.484	0.113	4.297	0.000	0.263	0.705
sc_age7II	0.401	0.106	3.795	0.000	0.194	0.607
sc_age7III non-manual	0.393	0.106	3.710	0.000	0.185	0.600
sc_age7III manual	0.225	0.103	2.185	0.029	0.023	0.426
sc_age7IV non-manual	0.215	0.122	1.766	0.077	-0.024	0.454
sc_age7IV manual	0.144	0.104	1.379	0.168	-0.060	0.348
sc_age7V	0.024	0.108	0.227	0.821	-0.187	0.235
breast_fDont know	-0.646	0.502	-1.286	0.198	-1.630	0.338
breast_fNo	-0.536	0.489	-1.096	0.273	-1.494	0.422
breast_fUnder one month	-0.467	0.489	-0.956	0.339	-1.426	0.491
breast_fOver one month	-0.386	0.489	-0.790	0.430	-1.344	0.572

reading_age7 <- lm(z_reading_test_age7 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f, data = df_complete)

reading_age7 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.183	0.510	-0.360	0.719	-1.182	0.816
sexFemale	0.076	0.017	4.488	0.000	0.043	0.110
sc_prenatalUnemployed,sick	-1.255	0.492	-2.551	0.011	-2.220	-0.291
sc_prenatalI	0.242	0.078	3.085	0.002	0.088	0.395
sc_prenatalII	0.201	0.065	3.097	0.002	0.074	0.329
sc_prenatalIII	0.148	0.059	2.494	0.013	0.032	0.264
sc_prenatalIV	0.090	0.063	1.420	0.156	-0.034	0.214
sc_prenatalV	0.012	0.065	0.180	0.857	-0.116	0.139
sc_prenatalStudents	0.162	0.279	0.582	0.561	-0.384	0.709
sc_prenatalDead or away	-0.747	0.695	-1.074	0.283	-2.109	0.616
sc_prenatalRetired	0.075	0.696	0.108	0.914	-1.288	1.438
sc_prenatalSingle,no husbnd	-0.063	0.080	-0.784	0.433	-0.220	0.094
smoking_preg5mths-no change	0.040	0.081	0.487	0.626	-0.120	0.199
smoking_pregNow non-smoker	0.095	0.087	1.095	0.274	-0.075	0.265
smoking_preg1 to 4 now	0.063	0.094	0.674	0.500	-0.120	0.246
smoking_preg5 to 9 now	0.009	0.098	0.096	0.924	-0.182	0.201
smoking_preg10 to 14 now	0.042	0.112	0.375	0.708	-0.177	0.261
smoking_preg15 to 19 now	-0.040	0.137	-0.294	0.769	-0.309	0.228
smoking_preg20 to 24 now	-0.181	0.158	-1.148	0.251	-0.490	0.128
smoking_preg25 to 29 now	0.277	0.244	1.136	0.256	-0.201	0.756
smoking_preg30 or more now	0.216	0.273	0.791	0.429	-0.320	0.752
smoking_pregVariable-5mths	0.042	0.088	0.480	0.631	-0.131	0.215
birth_weight	0.000	0.000	2.205	0.027	0.000	0.001
parent_read_levelOccasionally	0.115	0.036	3.188	0.001	0.044	0.185
parent_read_levelWeekly	0.168	0.033	5.006	0.000	0.102	0.233
sc_age7No male head	0.180	0.114	1.576	0.115	-0.044	0.403
sc_age7I	0.484	0.113	4.297	0.000	0.263	0.705
sc_age7II	0.401	0.106	3.795	0.000	0.194	0.607
sc_age7III non-manual	0.393	0.106	3.710	0.000	0.185	0.600
sc_age7III manual	0.225	0.103	2.185	0.029	0.023	0.426
sc_age7IV non-manual	0.215	0.122	1.766	0.077	-0.024	0.454
sc_age7IV manual	0.144	0.104	1.379	0.168	-0.060	0.348
sc_age7V	0.024	0.108	0.227	0.821	-0.187	0.235
breast_fDont know	-0.646	0.502	-1.286	0.198	-1.630	0.338
breast_fNo	-0.536	0.489	-1.096	0.273	-1.494	0.422
breast_fUnder one month	-0.467	0.489	-0.956	0.339	-1.426	0.491
breast_fOver one month	-0.386	0.489	-0.790	0.430	-1.344	0.572

extract_parent_read <- function(model, outcome_label) {
  
  coef_df <- tidy(model, conf.int = TRUE) %>%
    filter(term %in% c(
      "parent_read_levelOccasionally",
      "parent_read_levelWeekly"
    )) %>%
    mutate(
      exposure = case_when(
        term == "parent_read_levelOccasionally" ~ "Occasionally",
        term == "parent_read_levelWeekly" ~ "Weekly"
      ),
      outcome = outcome_label
    ) %>%
    dplyr::select(outcome, exposure, estimate, conf.low, conf.high)
  
  # referência (Hardly ever)
  ref_df <- tibble(
    outcome   = outcome_label,
    exposure  = "Hardly ever",
    estimate  = 0,
    conf.low  = 0,
    conf.high = 0
  )
  
  bind_rows(ref_df, coef_df)
}

plot_df <- bind_rows(
  extract_parent_read(global_age7,  "Global cognition"),
  extract_parent_read(arith_age7,   "Arithmetic"),
  extract_parent_read(copying_age7, "Copying"),
  extract_parent_read(draw_age7,    "Draw-a-man"),
  extract_parent_read(reading_age7, "Reading")
)

plot_df$exposure <- factor(
  plot_df$exposure,
  levels = c("Hardly ever", "Occasionally", "Weekly")
)


ggplot(plot_df, aes(x = estimate, y = outcome)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high, shape = exposure),
    size = 1, position = position_dodge2(width = 0.9)
  ) +
  
  scale_shape_manual(
    values = c(
      "Hardly ever" = 15,
      "Occasionally" = 16,
      "Weekly" = 17
    )
  ) +
  
  labs(
    x = "Standardized difference in Z score (95% CI)",
    y = NULL,
    shape = "Parental reading frequency"
  ) +
  
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

COGNITIVE FUNCTION AT AGE 11

df_complete_7_11 <- df_complete %>%
  filter(verbal_ability_age11!=-1 &
           nonverbal_ability_age11!=-1  &
           reading_ability_age11!=-1 &
           maths_ability_age11!=-1 &
           copying_age11!=-1)

# Age 11
# Verbal score – general ability test
df_complete_7_11$mean_verbal_ability_age11 <- mean(df_complete_7_11$verbal_ability_age11)
df_complete_7_11$sd_verbal_ability_age11 <- sd(df_complete_7_11$verbal_ability_age11)
df_complete_7_11$z_verbal_ability_age11 <- (df_complete_7_11$arithmetic_age7-df_complete_7_11$mean_verbal_ability_age11)/df_complete_7_11$sd_verbal_ability_age11
df_complete_7_11$z_verbal_ability_age11 = scale(df_complete_7_11$z_verbal_ability_age11)[,1]

hist(df_complete_7_11$z_verbal_ability_age11)

# Non-verbal score – general ability test
df_complete_7_11$mean_n_verbal_ability_age11 <- mean(df_complete_7_11$nonverbal_ability_age11)
df_complete_7_11$sd_n_verbal_ability_age11 <- sd(df_complete_7_11$nonverbal_ability_age11)
df_complete_7_11$z_n_verbal_ability_age11 <- (df_complete_7_11$nonverbal_ability_age11-df_complete_7_11$mean_n_verbal_ability_age11)/df_complete_7_11$sd_n_verbal_ability_age11
df_complete_7_11$z_n_verbal_ability_age11 = scale(df_complete_7_11$z_n_verbal_ability_age11)[,1]

hist(df_complete_7_11$z_n_verbal_ability_age11)

#Reading comprehension test
df_complete_7_11$mean_reading_ability_age11 <- mean(df_complete_7_11$reading_ability_age11)
df_complete_7_11$sd_reading_ability_age11 <- sd(df_complete_7_11$reading_ability_age11)
df_complete_7_11$z_reading_ability_age11 <- (df_complete_7_11$reading_ability_age11-df_complete_7_11$mean_reading_ability_age11)/df_complete_7_11$sd_reading_ability_age11
df_complete_7_11$z_reading_ability_age11 = scale(df_complete_7_11$z_reading_ability_age11)[,1]

hist(df_complete_7_11$z_reading_ability_age11)

#Mathematics test score
df_complete_7_11$mean_maths_ability_age11 <- mean(df_complete_7_11$maths_ability_age11)
df_complete_7_11$sd_maths_ability_age11 <- sd(df_complete_7_11$maths_ability_age11)
df_complete_7_11$z_maths_ability_age11 <- (df_complete_7_11$maths_ability_age11-df_complete_7_11$mean_maths_ability_age11)/df_complete_7_11$sd_maths_ability_age11
df_complete_7_11$z_maths_ability_age11 = scale(df_complete_7_11$z_maths_ability_age11)[,1]

hist(df_complete_7_11$z_maths_ability_age11)

#Copying designs test score
df_complete_7_11$mean_copying_age11 <- mean(df_complete_7_11$copying_age11)
df_complete_7_11$sd_copying_age11 <- sd(df_complete_7_11$copying_age11)
df_complete_7_11$z_copying_age11 <- (df_complete_7_11$copying_age11-df_complete_7_11$mean_copying_age11)/df_complete_7_11$sd_copying_age11
df_complete_7_11$z_copying_age11 = scale(df_complete_7_11$z_copying_age11)[,1]

hist(df_complete_7_11$z_copying_age11)

# averaging all z scores
df_complete_7_11$global_cog_age11 <- (df_complete_7_11$z_verbal_ability_age11 + df_complete_7_11$z_n_verbal_ability_age11 + df_complete_7_11$z_reading_ability_age11 + df_complete_7_11$z_maths_ability_age11 + df_complete_7_11$z_copying_age11) / 5

# creating a z score for the averaged variable
df_complete_7_11$z_global_cog_age11 = scale(df_complete_7_11$global_cog_age11)[,1]

hist(df_complete_7_11$z_global_cog_age11)

# Reference levels
library(dplyr)
table(df_complete_7_11$tv_age11)

## 
##          NA   most days   Sometimes Hardly ever 
##         438       10134        1430         299

table(df_complete_7_11$child_read_books_age11)

## 
##          NA   most days   Sometimes Hardly ever 
##         437        5223        5273        1368

table(df_complete_7_11$child_read_news_age11)

## 
##          NA   most days   Sometimes Hardly ever 
##         499        6156        4613        1033

table(df_complete_7_11$physical_activity_age11)

## 
##          NA   most days   Sometimes Hardly ever 
##         442        5476        5049        1334

table(df_complete_7_11$child_writes_stories_age11)

## 
##          NA   most days   Sometimes Hardly ever 
##         502        1082        4032        6685

table(df_complete_7_11$child_draws_age11)

## 
##          NA   most days   Sometimes Hardly ever 
##         473        2799        6007        3022

table(df_complete_7_11$child_writes_stories_age11)

## 
##          NA   most days   Sometimes Hardly ever 
##         502        1082        4032        6685

library(dplyr)

recode_freq <- function(x) {
  case_when(
    x == "Hardly ever" ~ 0,
    x == "Sometimes"  ~ 1,
    x == "most days"  ~ 2,
    TRUE ~ NA_real_
  )
}
df_complete_7_11 <- df_complete_7_11 %>%
  mutate(
    read_books_score = recode_freq(child_read_books_age11),
    read_news_score  = recode_freq(child_read_news_age11),
    
    child_reading_age11 = pmax(read_books_score, read_news_score, na.rm = TRUE),
    child_reading_age11 = ifelse(
      is.infinite(child_reading_age11), NA, child_reading_age11
    )
  )


df_complete_7_11$child_reading_age11 <- factor(
  df_complete_7_11$child_reading_age11,
  levels = c(0, 1, 2),
  labels = c("Hardly ever", "Sometimes", "Most days"))

table(df_complete_7_11$child_reading_age11, useNA = "ifany")

## 
## Hardly ever   Sometimes   Most days        <NA> 
##         233        3248        8536         284

df_complete_7_11 <- df_complete_7_11 %>%
  mutate(
    write_story_score = recode_freq(child_writes_stories_age11),
    draw_score        = recode_freq(child_draws_age11),
    
    child_creative_age11 = pmax(write_story_score, draw_score, na.rm = TRUE),
    child_creative_age11 = ifelse(
      is.infinite(child_creative_age11), NA, child_creative_age11
    )
  )

df_complete_7_11$child_creative_age11 <- factor(
  df_complete_7_11$child_creative_age11,
  levels = c(0, 1, 2),
  labels = c("Hardly ever", "Sometimes", "Most days"))

table(df_complete_7_11$child_creative_age11, useNA = "ifany")

## 
## Hardly ever   Sometimes   Most days        <NA> 
##        2074        6548        3391         288

df_complete_7_11$tv_age11 <- factor(
  df_complete_7_11$tv_age11,
  levels = c("Hardly ever", "Sometimes", "most days")
)

table(df_complete_7_11$tv_age11, useNA = "ifany")

## 
## Hardly ever   Sometimes   most days        <NA> 
##         299        1430       10134         438

df_complete_7_11$physical_activity_age11 <- factor(
  df_complete_7_11$physical_activity_age11,
  levels = c("Hardly ever", "Sometimes", "most days")
)

table(df_complete_7_11$physical_activity_age11, useNA = "ifany")

## 
## Hardly ever   Sometimes   most days        <NA> 
##        1334        5049        5476         442

global_age11 <- lm(z_global_cog_age11 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11, data = df_complete_7_11)

global_age11 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.758	0.465	-1.629	0.103	-1.670	0.154
sexFemale	-0.053	0.017	-3.182	0.001	-0.086	-0.020
sc_prenatalUnemployed,sick	0.251	0.873	0.288	0.774	-1.460	1.962
sc_prenatalI	0.418	0.077	5.429	0.000	0.267	0.569
sc_prenatalII	0.267	0.064	4.175	0.000	0.142	0.392
sc_prenatalIII	0.134	0.058	2.295	0.022	0.019	0.248
sc_prenatalIV	0.073	0.062	1.179	0.239	-0.049	0.195
sc_prenatalV	-0.090	0.064	-1.404	0.160	-0.216	0.036
sc_prenatalStudents	0.599	0.283	2.114	0.035	0.044	1.155
sc_prenatalDead or away	-1.144	0.875	-1.308	0.191	-2.859	0.570
sc_prenatalRetired	0.758	0.622	1.219	0.223	-0.461	1.978
sc_prenatalSingle,no husbnd	0.023	0.080	0.285	0.776	-0.134	0.179
smoking_preg5mths-no change	-0.033	0.082	-0.405	0.686	-0.194	0.128
smoking_pregNow non-smoker	0.005	0.087	0.057	0.955	-0.166	0.176
smoking_preg1 to 4 now	-0.111	0.094	-1.178	0.239	-0.295	0.074
smoking_preg5 to 9 now	-0.140	0.097	-1.440	0.150	-0.331	0.051
smoking_preg10 to 14 now	-0.152	0.111	-1.370	0.171	-0.370	0.066
smoking_preg15 to 19 now	-0.203	0.138	-1.468	0.142	-0.474	0.068
smoking_preg20 to 24 now	-0.373	0.156	-2.395	0.017	-0.679	-0.068
smoking_preg25 to 29 now	-0.178	0.264	-0.675	0.500	-0.697	0.340
smoking_preg30 or more now	-0.478	0.288	-1.663	0.096	-1.042	0.085
smoking_pregVariable-5mths	-0.124	0.089	-1.391	0.164	-0.298	0.051
birth_weight	0.000	0.000	1.759	0.079	0.000	0.000
parent_read_levelOccasionally	0.120	0.035	3.394	0.001	0.051	0.189
parent_read_levelWeekly	0.277	0.033	8.387	0.000	0.212	0.341
sc_age7No male head	0.007	0.118	0.061	0.952	-0.224	0.238
sc_age7I	0.444	0.117	3.803	0.000	0.215	0.672
sc_age7II	0.289	0.108	2.670	0.008	0.077	0.501
sc_age7III non-manual	0.282	0.109	2.589	0.010	0.068	0.496
sc_age7III manual	0.037	0.105	0.357	0.721	-0.168	0.242
sc_age7IV non-manual	0.048	0.122	0.394	0.694	-0.191	0.287
sc_age7IV manual	-0.047	0.106	-0.441	0.660	-0.255	0.161
sc_age7V	-0.201	0.110	-1.825	0.068	-0.416	0.015
breast_fDont know	-0.901	0.453	-1.990	0.047	-1.788	-0.013
breast_fNo	-0.749	0.436	-1.717	0.086	-1.604	0.106
breast_fUnder one month	-0.684	0.436	-1.567	0.117	-1.539	0.171
breast_fOver one month	-0.598	0.436	-1.371	0.171	-1.452	0.257
sc_age11I	0.314	0.057	5.555	0.000	0.203	0.425
sc_age11II	0.310	0.040	7.837	0.000	0.232	0.387
sc_age11III non manual	0.231	0.046	5.074	0.000	0.142	0.321
sc_age11III manual	0.086	0.032	2.640	0.008	0.022	0.149
sc_age11IV	0.006	0.037	0.163	0.871	-0.067	0.079
sc_age11V	-0.015	0.051	-0.293	0.769	-0.114	0.085
sc_age11No male head	0.012	0.055	0.225	0.822	-0.096	0.121
tv_age11Sometimes	0.173	0.058	2.958	0.003	0.058	0.287
tv_age11most days	0.282	0.054	5.227	0.000	0.176	0.388
child_reading_age11Sometimes	0.492	0.063	7.806	0.000	0.369	0.616
child_reading_age11Most days	0.782	0.062	12.642	0.000	0.661	0.904

arith_age11 <- lm(z_maths_ability_age11 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11, data = df_complete_7_11)

arith_age11 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.315	0.477	-0.660	0.509	-1.251	0.621
sexFemale	-0.056	0.017	-3.247	0.001	-0.089	-0.022
sc_prenatalUnemployed,sick	-0.213	0.896	-0.238	0.812	-1.970	1.543
sc_prenatalI	0.429	0.079	5.427	0.000	0.274	0.584
sc_prenatalII	0.316	0.066	4.809	0.000	0.187	0.444
sc_prenatalIII	0.120	0.060	2.003	0.045	0.003	0.237
sc_prenatalIV	0.070	0.064	1.094	0.274	-0.055	0.195
sc_prenatalV	-0.092	0.066	-1.401	0.161	-0.221	0.037
sc_prenatalStudents	0.784	0.291	2.694	0.007	0.213	1.354
sc_prenatalDead or away	-0.915	0.898	-1.019	0.308	-2.675	0.845
sc_prenatalRetired	0.799	0.638	1.252	0.211	-0.452	2.051
sc_prenatalSingle,no husbnd	0.004	0.082	0.054	0.957	-0.156	0.165
smoking_preg5mths-no change	0.032	0.084	0.378	0.705	-0.133	0.197
smoking_pregNow non-smoker	0.026	0.089	0.286	0.775	-0.150	0.201
smoking_preg1 to 4 now	-0.044	0.097	-0.457	0.647	-0.233	0.145
smoking_preg5 to 9 now	-0.059	0.100	-0.594	0.552	-0.255	0.136
smoking_preg10 to 14 now	-0.123	0.114	-1.076	0.282	-0.346	0.101
smoking_preg15 to 19 now	-0.150	0.142	-1.055	0.291	-0.428	0.128
smoking_preg20 to 24 now	-0.172	0.160	-1.078	0.281	-0.486	0.141
smoking_preg25 to 29 now	-0.112	0.271	-0.411	0.681	-0.643	0.420
smoking_preg30 or more now	-0.510	0.295	-1.730	0.084	-1.089	0.068
smoking_pregVariable-5mths	-0.072	0.091	-0.788	0.431	-0.251	0.107
birth_weight	0.000	0.000	2.098	0.036	0.000	0.001
parent_read_levelOccasionally	0.092	0.036	2.551	0.011	0.021	0.163
parent_read_levelWeekly	0.258	0.034	7.627	0.000	0.192	0.325
sc_age7No male head	-0.074	0.121	-0.616	0.538	-0.311	0.162
sc_age7I	0.381	0.120	3.184	0.001	0.147	0.616
sc_age7II	0.193	0.111	1.742	0.081	-0.024	0.411
sc_age7III non-manual	0.177	0.112	1.587	0.113	-0.042	0.397
sc_age7III manual	-0.064	0.107	-0.600	0.548	-0.275	0.146
sc_age7IV non-manual	-0.065	0.125	-0.521	0.602	-0.311	0.180
sc_age7IV manual	-0.139	0.109	-1.281	0.200	-0.353	0.074
sc_age7V	-0.258	0.113	-2.290	0.022	-0.479	-0.037
breast_fDont know	-1.297	0.465	-2.792	0.005	-2.208	-0.386
breast_fNo	-1.026	0.448	-2.292	0.022	-1.903	-0.149
breast_fUnder one month	-0.979	0.448	-2.187	0.029	-1.856	-0.101
breast_fOver one month	-0.901	0.448	-2.013	0.044	-1.778	-0.024
sc_age11I	0.255	0.058	4.385	0.000	0.141	0.369
sc_age11II	0.326	0.041	8.033	0.000	0.246	0.406
sc_age11III non manual	0.277	0.047	5.926	0.000	0.185	0.369
sc_age11III manual	0.112	0.033	3.357	0.001	0.046	0.177
sc_age11IV	0.014	0.038	0.367	0.714	-0.061	0.089
sc_age11V	0.010	0.052	0.188	0.851	-0.092	0.112
sc_age11No male head	-0.011	0.057	-0.197	0.843	-0.123	0.100
tv_age11Sometimes	0.173	0.060	2.890	0.004	0.056	0.290
tv_age11most days	0.261	0.055	4.719	0.000	0.153	0.370
child_reading_age11Sometimes	0.393	0.065	6.080	0.000	0.267	0.520
child_reading_age11Most days	0.676	0.064	10.646	0.000	0.552	0.801

copying_age11 <- lm(z_copying_age11 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11, data = df_complete_7_11)

copying_age11 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	0.094	0.504	0.186	0.852	-0.895	1.083
sexFemale	-0.049	0.018	-2.715	0.007	-0.085	-0.014
sc_prenatalUnemployed,sick	-0.115	0.947	-0.121	0.903	-1.971	1.741
sc_prenatalI	0.043	0.084	0.513	0.608	-0.121	0.207
sc_prenatalII	0.005	0.069	0.068	0.946	-0.131	0.141
sc_prenatalIII	-0.057	0.063	-0.905	0.365	-0.181	0.067
sc_prenatalIV	-0.047	0.067	-0.695	0.487	-0.179	0.085
sc_prenatalV	-0.190	0.070	-2.734	0.006	-0.326	-0.054
sc_prenatalStudents	0.055	0.307	0.178	0.859	-0.548	0.657
sc_prenatalDead or away	-0.972	0.949	-1.024	0.306	-2.832	0.888
sc_prenatalRetired	1.512	0.675	2.242	0.025	0.190	2.835
sc_prenatalSingle,no husbnd	-0.134	0.086	-1.555	0.120	-0.304	0.035
smoking_preg5mths-no change	-0.113	0.089	-1.269	0.205	-0.287	0.062
smoking_pregNow non-smoker	-0.122	0.095	-1.292	0.196	-0.308	0.063
smoking_preg1 to 4 now	-0.191	0.102	-1.869	0.062	-0.391	0.009
smoking_preg5 to 9 now	-0.206	0.106	-1.948	0.051	-0.413	0.001
smoking_preg10 to 14 now	-0.242	0.121	-2.011	0.044	-0.479	-0.006
smoking_preg15 to 19 now	-0.234	0.150	-1.561	0.119	-0.528	0.060
smoking_preg20 to 24 now	-0.610	0.169	-3.612	0.000	-0.942	-0.279
smoking_preg25 to 29 now	-0.425	0.287	-1.483	0.138	-0.987	0.137
smoking_preg30 or more now	-0.643	0.312	-2.063	0.039	-1.255	-0.032
smoking_pregVariable-5mths	-0.194	0.096	-2.010	0.044	-0.383	-0.005
birth_weight	0.000	0.000	0.326	0.745	0.000	0.000
parent_read_levelOccasionally	0.057	0.038	1.492	0.136	-0.018	0.132
parent_read_levelWeekly	0.096	0.036	2.685	0.007	0.026	0.166
sc_age7No male head	-0.090	0.128	-0.705	0.481	-0.340	0.160
sc_age7I	0.084	0.126	0.666	0.505	-0.164	0.332
sc_age7II	0.067	0.117	0.572	0.567	-0.163	0.297
sc_age7III non-manual	0.072	0.118	0.609	0.542	-0.160	0.304
sc_age7III manual	-0.002	0.113	-0.022	0.983	-0.225	0.220
sc_age7IV non-manual	-0.033	0.132	-0.248	0.804	-0.292	0.227
sc_age7IV manual	-0.073	0.115	-0.636	0.525	-0.299	0.152
sc_age7V	-0.142	0.119	-1.192	0.233	-0.376	0.091
breast_fDont know	-0.499	0.491	-1.018	0.309	-1.462	0.463
breast_fNo	-0.464	0.473	-0.980	0.327	-1.391	0.463
breast_fUnder one month	-0.402	0.473	-0.849	0.396	-1.329	0.526
breast_fOver one month	-0.367	0.473	-0.777	0.437	-1.294	0.560
sc_age11I	0.252	0.061	4.100	0.000	0.131	0.372
sc_age11II	0.155	0.043	3.622	0.000	0.071	0.239
sc_age11III non manual	0.080	0.049	1.615	0.106	-0.017	0.177
sc_age11III manual	0.041	0.035	1.165	0.244	-0.028	0.110
sc_age11IV	0.039	0.040	0.976	0.329	-0.040	0.118
sc_age11V	-0.010	0.055	-0.181	0.857	-0.118	0.098
sc_age11No male head	0.072	0.060	1.201	0.230	-0.046	0.190
tv_age11Sometimes	0.013	0.063	0.201	0.841	-0.111	0.137
tv_age11most days	0.088	0.058	1.510	0.131	-0.026	0.203
child_reading_age11Sometimes	0.295	0.068	4.317	0.000	0.161	0.429
child_reading_age11Most days	0.330	0.067	4.915	0.000	0.198	0.461

non_verbal_age11 <- lm(z_n_verbal_ability_age11 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11, data = df_complete_7_11)

non_verbal_age11 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.606	0.482	-1.258	0.208	-1.550	0.338
sexFemale	0.017	0.017	0.982	0.326	-0.017	0.051
sc_prenatalUnemployed,sick	0.198	0.904	0.219	0.827	-1.574	1.970
sc_prenatalI	0.320	0.080	4.012	0.000	0.164	0.476
sc_prenatalII	0.203	0.066	3.065	0.002	0.073	0.333
sc_prenatalIII	0.110	0.060	1.818	0.069	-0.009	0.228
sc_prenatalIV	0.033	0.064	0.519	0.604	-0.093	0.160
sc_prenatalV	-0.108	0.066	-1.622	0.105	-0.238	0.022
sc_prenatalStudents	0.370	0.293	1.262	0.207	-0.205	0.946
sc_prenatalDead or away	-0.693	0.906	-0.765	0.444	-2.469	1.083
sc_prenatalRetired	-0.236	0.644	-0.367	0.714	-1.499	1.027
sc_prenatalSingle,no husbnd	0.031	0.083	0.378	0.705	-0.131	0.193
smoking_preg5mths-no change	-0.004	0.085	-0.050	0.960	-0.171	0.162
smoking_pregNow non-smoker	0.021	0.090	0.237	0.813	-0.156	0.198
smoking_preg1 to 4 now	-0.086	0.097	-0.888	0.375	-0.277	0.104
smoking_preg5 to 9 now	-0.145	0.101	-1.440	0.150	-0.343	0.052
smoking_preg10 to 14 now	-0.077	0.115	-0.672	0.501	-0.303	0.148
smoking_preg15 to 19 now	-0.204	0.143	-1.427	0.154	-0.485	0.076
smoking_preg20 to 24 now	-0.107	0.161	-0.661	0.508	-0.423	0.210
smoking_preg25 to 29 now	0.176	0.274	0.643	0.520	-0.361	0.712
smoking_preg30 or more now	-0.249	0.298	-0.836	0.403	-0.833	0.335
smoking_pregVariable-5mths	-0.119	0.092	-1.288	0.198	-0.299	0.062
birth_weight	0.000	0.000	1.158	0.247	0.000	0.000
parent_read_levelOccasionally	0.160	0.037	4.386	0.000	0.089	0.232
parent_read_levelWeekly	0.292	0.034	8.544	0.000	0.225	0.359
sc_age7No male head	-0.039	0.122	-0.321	0.749	-0.278	0.200
sc_age7I	0.368	0.121	3.049	0.002	0.131	0.605
sc_age7II	0.221	0.112	1.968	0.049	0.001	0.440
sc_age7III non-manual	0.206	0.113	1.825	0.068	-0.015	0.427
sc_age7III manual	0.012	0.108	0.107	0.915	-0.201	0.224
sc_age7IV non-manual	0.009	0.126	0.068	0.946	-0.239	0.256
sc_age7IV manual	-0.058	0.110	-0.527	0.598	-0.273	0.157
sc_age7V	-0.245	0.114	-2.151	0.031	-0.468	-0.022
breast_fDont know	-0.788	0.469	-1.681	0.093	-1.706	0.131
breast_fNo	-0.844	0.452	-1.868	0.062	-1.729	0.042
breast_fUnder one month	-0.759	0.452	-1.681	0.093	-1.644	0.126
breast_fOver one month	-0.679	0.452	-1.504	0.133	-1.564	0.206
sc_age11I	0.304	0.059	5.182	0.000	0.189	0.419
sc_age11II	0.278	0.041	6.778	0.000	0.197	0.358
sc_age11III non manual	0.191	0.047	4.040	0.000	0.098	0.283
sc_age11III manual	0.079	0.034	2.346	0.019	0.013	0.145
sc_age11IV	0.024	0.038	0.619	0.536	-0.052	0.099
sc_age11V	-0.040	0.052	-0.754	0.451	-0.142	0.063
sc_age11No male head	0.013	0.057	0.230	0.818	-0.099	0.126
tv_age11Sometimes	0.207	0.060	3.425	0.001	0.089	0.325
tv_age11most days	0.308	0.056	5.523	0.000	0.199	0.418
child_reading_age11Sometimes	0.437	0.065	6.700	0.000	0.309	0.565
child_reading_age11Most days	0.677	0.064	10.570	0.000	0.552	0.803

verbal_age11 <- lm(z_verbal_ability_age11 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11, data = df_complete_7_11)

verbal_age11 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-0.945	0.513	-1.842	0.066	-1.950	0.061
sexFemale	-0.104	0.018	-5.647	0.000	-0.140	-0.068
sc_prenatalUnemployed,sick	0.719	0.963	0.747	0.455	-1.167	2.606
sc_prenatalI	0.302	0.085	3.553	0.000	0.135	0.468
sc_prenatalII	0.225	0.071	3.190	0.001	0.087	0.363
sc_prenatalIII	0.200	0.064	3.113	0.002	0.074	0.325
sc_prenatalIV	0.149	0.069	2.177	0.030	0.015	0.283
sc_prenatalV	0.069	0.071	0.983	0.326	-0.069	0.208
sc_prenatalStudents	0.452	0.312	1.448	0.148	-0.160	1.065
sc_prenatalDead or away	-0.574	0.965	-0.595	0.552	-2.465	1.317
sc_prenatalRetired	-0.040	0.686	-0.058	0.954	-1.384	1.305
sc_prenatalSingle,no husbnd	0.103	0.088	1.170	0.242	-0.069	0.275
smoking_preg5mths-no change	-0.059	0.091	-0.651	0.515	-0.236	0.118
smoking_pregNow non-smoker	0.012	0.096	0.123	0.902	-0.177	0.200
smoking_preg1 to 4 now	-0.083	0.104	-0.799	0.424	-0.286	0.120
smoking_preg5 to 9 now	-0.081	0.107	-0.758	0.448	-0.292	0.129
smoking_preg10 to 14 now	-0.064	0.123	-0.521	0.602	-0.304	0.176
smoking_preg15 to 19 now	-0.205	0.152	-1.343	0.179	-0.503	0.094
smoking_preg20 to 24 now	-0.313	0.172	-1.823	0.068	-0.650	0.024
smoking_preg25 to 29 now	-0.500	0.291	-1.715	0.086	-1.071	0.072
smoking_preg30 or more now	-0.061	0.317	-0.192	0.848	-0.682	0.561
smoking_pregVariable-5mths	-0.084	0.098	-0.862	0.388	-0.277	0.108
birth_weight	0.000	0.000	1.171	0.242	0.000	0.000
parent_read_levelOccasionally	0.064	0.039	1.646	0.100	-0.012	0.140
parent_read_levelWeekly	0.125	0.036	3.437	0.001	0.054	0.196
sc_age7No male head	0.044	0.130	0.342	0.732	-0.210	0.299
sc_age7I	0.315	0.129	2.452	0.014	0.063	0.567
sc_age7II	0.205	0.119	1.722	0.085	-0.028	0.439
sc_age7III non-manual	0.176	0.120	1.464	0.143	-0.060	0.411
sc_age7III manual	0.009	0.115	0.081	0.935	-0.217	0.235
sc_age7IV non-manual	0.063	0.135	0.471	0.638	-0.200	0.327
sc_age7IV manual	-0.037	0.117	-0.316	0.752	-0.266	0.192
sc_age7V	-0.105	0.121	-0.869	0.385	-0.343	0.132
breast_fDont know	-0.286	0.499	-0.574	0.566	-1.265	0.692
breast_fNo	0.003	0.481	0.007	0.995	-0.939	0.946
breast_fUnder one month	0.014	0.481	0.030	0.976	-0.928	0.957
breast_fOver one month	0.073	0.481	0.151	0.880	-0.870	1.015
sc_age11I	0.199	0.062	3.187	0.001	0.077	0.321
sc_age11II	0.156	0.044	3.587	0.000	0.071	0.242
sc_age11III non manual	0.119	0.050	2.359	0.018	0.020	0.217
sc_age11III manual	0.053	0.036	1.480	0.139	-0.017	0.123
sc_age11IV	-0.012	0.041	-0.297	0.766	-0.092	0.068
sc_age11V	0.045	0.056	0.806	0.420	-0.064	0.155
sc_age11No male head	0.005	0.061	0.087	0.931	-0.114	0.125
tv_age11Sometimes	0.099	0.064	1.546	0.122	-0.027	0.226
tv_age11most days	0.174	0.059	2.920	0.004	0.057	0.290
child_reading_age11Sometimes	0.340	0.070	4.896	0.000	0.204	0.477
child_reading_age11Most days	0.517	0.068	7.576	0.000	0.383	0.651

reading_age11 <- lm(z_reading_ability_age11 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11, data = df_complete_7_11)

reading_age11 %>%
  broom::tidy(conf.int = TRUE) %>%
  knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	-1.141	0.469	-2.431	0.015	-2.061	-0.221
sexFemale	-0.012	0.017	-0.740	0.460	-0.046	0.021
sc_prenatalUnemployed,sick	0.376	0.881	0.427	0.669	-1.351	2.104
sc_prenatalI	0.514	0.078	6.613	0.000	0.362	0.667
sc_prenatalII	0.278	0.065	4.313	0.000	0.152	0.405
sc_prenatalIII	0.142	0.059	2.412	0.016	0.027	0.257
sc_prenatalIV	0.076	0.063	1.215	0.225	-0.047	0.199
sc_prenatalV	-0.025	0.065	-0.393	0.694	-0.152	0.101
sc_prenatalStudents	0.642	0.286	2.244	0.025	0.081	1.203
sc_prenatalDead or away	-1.246	0.883	-1.411	0.158	-2.977	0.485
sc_prenatalRetired	0.878	0.628	1.399	0.162	-0.352	2.109
sc_prenatalSingle,no husbnd	0.083	0.080	1.034	0.301	-0.075	0.241
smoking_preg5mths-no change	0.017	0.083	0.201	0.841	-0.146	0.179
smoking_pregNow non-smoker	0.082	0.088	0.937	0.349	-0.090	0.255
smoking_preg1 to 4 now	-0.022	0.095	-0.227	0.820	-0.208	0.165
smoking_preg5 to 9 now	-0.047	0.098	-0.481	0.631	-0.240	0.145
smoking_preg10 to 14 now	-0.079	0.112	-0.706	0.480	-0.299	0.141
smoking_preg15 to 19 now	0.012	0.140	0.088	0.930	-0.261	0.286
smoking_preg20 to 24 now	-0.232	0.157	-1.476	0.140	-0.541	0.076
smoking_preg25 to 29 now	0.175	0.267	0.655	0.513	-0.348	0.698
smoking_preg30 or more now	-0.375	0.290	-1.292	0.196	-0.944	0.194
smoking_pregVariable-5mths	-0.007	0.090	-0.073	0.942	-0.182	0.169
birth_weight	0.000	0.000	1.749	0.080	0.000	0.000
parent_read_levelOccasionally	0.087	0.036	2.433	0.015	0.017	0.156
parent_read_levelWeekly	0.292	0.033	8.786	0.000	0.227	0.358
sc_age7No male head	0.187	0.119	1.571	0.116	-0.046	0.419
sc_age7I	0.556	0.118	4.725	0.000	0.325	0.787
sc_age7II	0.424	0.109	3.884	0.000	0.210	0.638
sc_age7III non-manual	0.453	0.110	4.121	0.000	0.238	0.669
sc_age7III manual	0.190	0.106	1.796	0.073	-0.017	0.397
sc_age7IV non-manual	0.211	0.123	1.711	0.087	-0.031	0.452
sc_age7IV manual	0.128	0.107	1.194	0.233	-0.082	0.338
sc_age7V	-0.021	0.111	-0.186	0.852	-0.238	0.197
breast_fDont know	-0.592	0.457	-1.295	0.195	-1.487	0.304
breast_fNo	-0.548	0.440	-1.245	0.213	-1.411	0.315
breast_fUnder one month	-0.503	0.440	-1.142	0.253	-1.366	0.360
breast_fOver one month	-0.422	0.440	-0.960	0.337	-1.285	0.440
sc_age11I	0.200	0.057	3.499	0.000	0.088	0.312
sc_age11II	0.276	0.040	6.919	0.000	0.198	0.354
sc_age11III non manual	0.223	0.046	4.844	0.000	0.133	0.313
sc_age11III manual	0.045	0.033	1.369	0.171	-0.019	0.109
sc_age11IV	-0.042	0.037	-1.112	0.266	-0.115	0.032
sc_age11V	-0.062	0.051	-1.220	0.222	-0.163	0.038
sc_age11No male head	-0.032	0.056	-0.565	0.572	-0.141	0.078
tv_age11Sometimes	0.171	0.059	2.908	0.004	0.056	0.287
tv_age11most days	0.252	0.054	4.634	0.000	0.146	0.359
child_reading_age11Sometimes	0.426	0.064	6.691	0.000	0.301	0.551
child_reading_age11Most days	0.808	0.062	12.929	0.000	0.685	0.930

extract_age11_exposures <- function(model, outcome_label) {
  
  coef_df <- tidy(model, conf.int = TRUE) %>%
    filter(term %in% c(
      "tv_age11Sometimes",
      "tv_age11most days",
      "child_reading_age11Sometimes",
      "child_reading_age11Most days"
    )) %>%
    mutate(
      exposure = case_when(
        term == "tv_age11Sometimes" ~ "TV: Sometimes",
        term == "tv_age11most days" ~ "TV: Most days",
        term == "child_reading_age11Sometimes" ~ "Reading: Sometimes",
        term == "child_reading_age11Most days" ~ "Reading: Most days"
      ),
      domain = case_when(
        grepl("^tv_", term) ~ "Passive SB (TV)",
        grepl("^child_reading_", term) ~ "Cognitively engaged SB (Reading)"
      ),
      outcome = outcome_label
    ) %>%
    dplyr::select(outcome, domain, exposure, estimate, conf.low, conf.high)
  
  # referências (Hardly ever)
  ref_df <- tibble(
    outcome   = rep(outcome_label, 2),
    domain    = c("Passive SB (TV)", "Cognitively engaged SB (Reading)"),
    exposure  = c("TV: Hardly ever", "Reading: Hardly ever"),
    estimate  = 0,
    conf.low  = 0,
    conf.high = 0
  )
  
  bind_rows(ref_df, coef_df)
}

plot_df_age11 <- bind_rows(
  extract_age11_exposures(global_age11,     "Global cognition"),
  extract_age11_exposures(arith_age11,      "Arithmetic"),
  extract_age11_exposures(copying_age11,    "Copying"),
  extract_age11_exposures(non_verbal_age11, "Non-verbal ability"),
  extract_age11_exposures(verbal_age11,     "Verbal ability"),
  extract_age11_exposures(reading_age11,    "Reading ability")
)

plot_df_age11$exposure <- factor(
  plot_df_age11$exposure,
  levels = c(
    "TV: Hardly ever", "TV: Sometimes", "TV: Most days",
    "Reading: Hardly ever", "Reading: Sometimes", "Reading: Most days"
  )
)

plot_df_age11$dose <- NA

plot_df_age11$dose[plot_df_age11$exposure=="Reading: Hardly ever"] <- "Hardly ever"
plot_df_age11$dose[plot_df_age11$exposure=="TV: Hardly ever"] <- "Hardly ever"
plot_df_age11$dose[plot_df_age11$exposure=="Reading: Sometimes"] <- "Sometimes"
plot_df_age11$dose[plot_df_age11$exposure=="TV: Sometimes"] <- "Sometimes"
plot_df_age11$dose[plot_df_age11$exposure=="TV: Most days"] <- "Most days"
plot_df_age11$dose[plot_df_age11$exposure=="Reading: Most days"] <- "Most days"

library(ggplot2)
plot_df_age11$dose <- factor(plot_df_age11$dose, levels = c("Most days", 
                                                            "Sometimes",
                                                            "Hardly ever"), 
                             ordered = T)

ggplot(plot_df_age11,
       aes(x = estimate, y = outcome, colour = domain, fill=domain, shape = dose)) +
  
  scale_color_jama() +
  scale_fill_jama() +

  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high, shape = dose),
    size = 1, position = position_dodge2(width = 0.9)
  ) +
  
  facet_wrap(~ domain, ncol = 1, scales = "free_y") +
  
  labs(
    x = "Standardized difference in cognitive score (z)",
    y = NULL,
    shape = "Exposure category"
  ) +
  
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    strip.text = element_text(face = "bold")
  )

library(broom)
library(dplyr)

extract_age7_11_exposures <- function(model, outcome_label) {
  
  coef_df <- tidy(model, conf.int = TRUE) %>%
    filter(term %in% c(
      # Age 7
      "parent_read_levelOccasionally",
      "parent_read_levelWeekly",
      # Age 11
      "tv_age11Sometimes",
      "tv_age11most days",
      "child_reading_age11Sometimes",
      "child_reading_age11Most days"
    )) %>%
    mutate(
      exposure = case_when(
        term == "parent_read_levelOccasionally" ~ "Reading: Sometimes",
        term == "parent_read_levelWeekly"       ~ "Reading: Most days",
        term == "child_reading_age11Sometimes"  ~ "Reading: Sometimes",
        term == "child_reading_age11Most days"  ~ "Reading: Most days",
        term == "tv_age11Sometimes"              ~ "TV: Sometimes",
        term == "tv_age11most days"              ~ "TV: Most days"
      ),
      domain = case_when(
        grepl("^tv_", term) ~ "Passive SB (TV)",
        TRUE                ~ "Cognitively engaged SB (Reading)"
      ),
      age = case_when(
        grepl("^parent_", term)         ~ "Age 7",
        grepl("^tv_", term)             ~ "Age 11",
        grepl("^child_reading_", term)  ~ "Age 11"
      ),
      outcome = outcome_label
    ) %>%
    dplyr::select(outcome, domain, age, exposure, estimate, conf.low, conf.high)
  
  # referências (Hardly ever)
  ref_df <- tibble(
    outcome   = outcome_label,
    domain    = c(
      "Cognitively engaged SB (Reading)",
      "Cognitively engaged SB (Reading)",
      "Passive SB (TV)"
    ),
    age       = c("Age 7", "Age 11", "Age 11"),
    exposure  = c(
      "Reading: Hardly ever",
      "Reading: Hardly ever",
      "TV: Hardly ever"
    ),
    estimate  = 0,
    conf.low  = 0,
    conf.high = 0
  )
  
  bind_rows(ref_df, coef_df)
}


plot_df <- bind_rows(
  extract_age7_11_exposures(global_age7,  "Global cognition"),
  extract_age7_11_exposures(arith_age7,   "Arithmetic"),
  extract_age7_11_exposures(copying_age7, "Copying"),
  extract_age7_11_exposures(draw_age7,    "Draw-a-man"),
  extract_age7_11_exposures(reading_age7, "Reading"),
  
  extract_age7_11_exposures(global_age11,     "Global cognition"),
  extract_age7_11_exposures(arith_age11,      "Arithmetic"),
  extract_age7_11_exposures(copying_age11,    "Copying"),
  extract_age7_11_exposures(non_verbal_age11, "Non-verbal ability"),
  extract_age7_11_exposures(verbal_age11,     "Verbal ability"),
  extract_age7_11_exposures(reading_age11,    "Reading ability")
)

plot_df$dose <- factor(
  gsub(".*: ", "", plot_df$exposure),
  levels = c("Hardly ever", "Sometimes", "Most days"),
  ordered = TRUE
)


library(ggplot2)
library(ggsci)

ggplot(plot_df,
       aes(x = estimate, y = outcome,
           colour = age, shape = dose)) +
  
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high),
    position = position_dodge(width = 0.6),
    size = 0.9
  ) +
  
  facet_wrap(~ domain, ncol = 1, scales = "free_y") +
  
  scale_color_jama() +
  
  labs(
    x = "Standardized difference in cognitive score (z)",
    y = NULL,
    colour = "Age",
    shape  = "Exposure frequency"
  ) +
  
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    strip.text = element_text(face = "bold")
  )

COGNITIVE FUNCTION AT AGE 16

df_complete_7_16 <- df_complete_7_11 %>%
  filter(reading_test_age16!=-1 &
            maths_test_age16!=-1)

# Age 16
#Reading comprehension test
df_complete_7_16$mean_reading_ability_age16 <- mean(df_complete_7_16$reading_test_age16)
df_complete_7_16$sd_reading_ability_age16 <- sd(df_complete_7_16$reading_test_age16)
df_complete_7_16$z_reading_ability_age16 <- (df_complete_7_16$reading_test_age16-df_complete_7_16$mean_reading_ability_age16)/df_complete_7_16$sd_reading_ability_age16
df_complete_7_16$z_reading_ability_age16 = scale(df_complete_7_16$z_reading_ability_age16)[,1]

hist(df_complete_7_16$z_reading_ability_age16)

#Mathematics test score
df_complete_7_16$mean_maths_ability_age16 <- mean(df_complete_7_16$maths_test_age16)
df_complete_7_16$sd_maths_ability_age16 <- sd(df_complete_7_16$maths_test_age16)
df_complete_7_16$z_maths_ability_age16 <- (df_complete_7_16$maths_test_age16-df_complete_7_16$mean_maths_ability_age16)/df_complete_7_16$sd_maths_ability_age16
df_complete_7_16$z_maths_ability_age16 = scale(df_complete_7_16$z_maths_ability_age16)[,1]

hist(df_complete_7_16$z_maths_ability_age16)

# averaging all z scores
df_complete_7_16$global_cog_age16 <- (df_complete_7_16$z_reading_ability_age16 + df_complete_7_16$z_maths_ability_age16) / 2

# creating a z score for the averaged variable
df_complete_7_16$z_global_cog_age16 = scale(df_complete_7_16$global_cog_age16)[,1]

hist(df_complete_7_16$z_global_cog_age16)

# Reference levels
library(dplyr)
table(df_complete_7_16$tv_age16)

## 
##          NA       Often   Sometimes Hardly ever   No chance 
##         324        5736        2546         426          50

table(df_complete_7_16$child_read_age16)

## 
##          NA       Often   Sometimes Hardly ever   No chance 
##         592        2242        3838        2136         274

table(df_complete_7_16$play_indoor_age16)

## 
##          NA       Often   Sometimes Hardly ever   No chance 
##         764        2096        2655        2716         851

table(df_complete_7_16$play_outdoor_age16)

## 
##          NA       Often   Sometimes Hardly ever   No chance 
##         585        3176        3009        2048         264

df_complete_7_16$pa_age16 <- dplyr::case_when(
  
  is.na(df_complete_7_16$play_indoor_age16) &
  is.na(df_complete_7_16$play_outdoor_age16) ~ NA_character_,
  
  df_complete_7_16$play_indoor_age16 == "Most days" |
  df_complete_7_16$play_outdoor_age16 == "Most days" ~ "Most days",
  
  df_complete_7_16$play_indoor_age16 == "Sometimes" |
  df_complete_7_16$play_outdoor_age16 == "Sometimes" ~ "Sometimes",
  
  df_complete_7_16$play_indoor_age16 == "Hardly ever" &
  df_complete_7_16$play_outdoor_age16 == "Hardly ever" ~ "Hardly ever"
)

df_complete_7_16$pa_age16 <- factor(
  df_complete_7_16$pa_age16,
  levels = c("Hardly ever", "Sometimes", "Most days")
)

df_complete_7_16$tv_age16 <- factor(
  df_complete_7_16$tv_age16,
  levels = c("No chance", "Hardly ever", "Sometimes", "Often")
)

df_complete_7_16$tv_age16 <- fct_collapse(df_complete_7_16$tv_age16,
                                                "Hardly ever" = c("No chance", "Hardly ever"))


df_complete_7_16$child_read_age16 <- factor(
  df_complete_7_16$child_read_age16,
  levels = c("Hardly ever", "Sometimes", "Most days")
)


global_age16 <- lm(z_global_cog_age16 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11 +
                   sc_age16 + smoking_age16 + alcohol_age16 +
                   child_read_age16 + tv_age16 + school_attendance_age16, data = df_complete_7_16)

DescTools::Desc(df_complete_7_16$n2868)

## ────────────────────────────────────────────────────────────────────────────── 
## df_complete_7_16$n2868 (factor)
## 
##   length      n    NAs unique levels  dupes
##    9'082  9'082      0      5      5      y
##          100.0%   0.0%                     
## 
##          level   freq   perc  cumfreq  cumperc
## 1        Often  5'736  63.2%    5'736    63.2%
## 2    Sometimes  2'546  28.0%    8'282    91.2%
## 3  Hardly ever    426   4.7%    8'708    95.9%
## 4           NA    324   3.6%    9'032    99.4%
## 5    No chance     50   0.6%    9'082   100.0%

summary(global_age16)

## 
## Call:
## lm(formula = z_global_cog_age16 ~ sex + sc_prenatal + smoking_preg + 
##     birth_weight + parent_read_level + sc_age7 + breast_f + sc_age11 + 
##     tv_age11 + child_reading_age11 + sc_age16 + smoking_age16 + 
##     alcohol_age16 + child_read_age16 + tv_age16 + school_attendance_age16, 
##     data = df_complete_7_16)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.95817 -0.52850 -0.01963  0.56136  2.47544 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -2.1947999  0.6395462  -3.432 0.000604 ***
## sexFemale                     -0.2164616  0.0228955  -9.454  < 2e-16 ***
## sc_prenatalI                   0.3996873  0.1069935   3.736 0.000189 ***
## sc_prenatalII                  0.3583805  0.0859352   4.170 3.09e-05 ***
## sc_prenatalIII                 0.1854049  0.0781085   2.374 0.017647 *  
## sc_prenatalIV                  0.1261444  0.0830901   1.518 0.129033    
## sc_prenatalV                   0.0133601  0.0864026   0.155 0.877122    
## sc_prenatalStudents            0.5387014  0.3202801   1.682 0.092634 .  
## sc_prenatalRetired            -0.0560031  0.8251552  -0.068 0.945892    
## sc_prenatalSingle,no husbnd    0.2327157  0.1076485   2.162 0.030678 *  
## smoking_preg5mths-no change    0.1204601  0.1073726   1.122 0.261961    
## smoking_pregNow non-smoker     0.1331376  0.1146478   1.161 0.245583    
## smoking_preg1 to 4 now         0.0232709  0.1231596   0.189 0.850140    
## smoking_preg5 to 9 now         0.0748755  0.1290737   0.580 0.561873    
## smoking_preg10 to 14 now       0.0944198  0.1444866   0.653 0.513473    
## smoking_preg15 to 19 now       0.0378703  0.1937341   0.195 0.845028    
## smoking_preg20 to 24 now       0.2649211  0.2139438   1.238 0.215670    
## smoking_preg25 to 29 now       0.3798666  0.3471253   1.094 0.273864    
## smoking_preg30 or more now    -0.0032198  0.3480449  -0.009 0.992619    
## smoking_pregVariable-5mths     0.0788829  0.1159726   0.680 0.496417    
## birth_weight                   0.0003577  0.0001722   2.077 0.037855 *  
## parent_read_levelOccasionally  0.0397535  0.0473836   0.839 0.401523    
## parent_read_levelWeekly        0.2459046  0.0445804   5.516 3.63e-08 ***
## sc_age7No male head           -0.0091009  0.1575697  -0.058 0.953944    
## sc_age7I                       0.4244963  0.1564308   2.714 0.006677 ** 
## sc_age7II                      0.2760123  0.1421373   1.942 0.052206 .  
## sc_age7III non-manual          0.2622384  0.1430707   1.833 0.066870 .  
## sc_age7III manual              0.0946882  0.1367200   0.693 0.488610    
## sc_age7IV non-manual           0.0713224  0.1617748   0.441 0.659322    
## sc_age7IV manual              -0.0025879  0.1387269  -0.019 0.985118    
## sc_age7V                      -0.0285052  0.1450572  -0.197 0.844219    
## breast_fDont know             -0.7598063  0.5934819  -1.280 0.200513    
## breast_fNo                    -0.6884009  0.5733045  -1.201 0.229899    
## breast_fUnder one month       -0.6399981  0.5734778  -1.116 0.264476    
## breast_fOver one month        -0.5426867  0.5733388  -0.947 0.343918    
## sc_age11I                      0.2308798  0.0818674   2.820 0.004818 ** 
## sc_age11II                     0.2518359  0.0566068   4.449 8.81e-06 ***
## sc_age11III non manual         0.2216753  0.0642269   3.451 0.000562 ***
## sc_age11III manual             0.1379524  0.0455862   3.026 0.002489 ** 
## sc_age11IV                     0.0359138  0.0516552   0.695 0.486923    
## sc_age11V                      0.0526049  0.0709497   0.741 0.458460    
## sc_age11No male head           0.1483636  0.0774030   1.917 0.055322 .  
## tv_age11Sometimes              0.1529831  0.0811524   1.885 0.059467 .  
## tv_age11most days              0.2482052  0.0753031   3.296 0.000987 ***
## child_reading_age11Sometimes   0.3471052  0.0813162   4.269 2.00e-05 ***
## child_reading_age11Most days   0.5417820  0.0801444   6.760 1.53e-11 ***
## sc_age16I                      0.2665387  0.0802698   3.321 0.000905 ***
## sc_age16II                     0.1558685  0.0426578   3.654 0.000261 ***
## sc_age16III non-manual         0.1587386  0.0536288   2.960 0.003091 ** 
## sc_age16III manual             0.0681656  0.0321079   2.123 0.033799 *  
## sc_age16IV non-manual          0.0926711  0.1166108   0.795 0.426822    
## sc_age16IV manual              0.0387056  0.0461611   0.838 0.401794    
## sc_age16V                     -0.0540261  0.0667078  -0.810 0.418040    
## sc_age16Unclear                0.0520467  0.1285402   0.405 0.685563    
## smoking_age16None,dont smoke   0.6768602  0.1279508   5.290 1.27e-07 ***
## smoking_age16Less 1 a week     0.5496281  0.1429525   3.845 0.000122 ***
## smoking_age161 to 9 a week     0.4151491  0.1347945   3.080 0.002082 ** 
## smoking_age1610 to 19 a week   0.3929308  0.1369732   2.869 0.004139 ** 
## smoking_age1620 to 29 a week   0.3217505  0.1368990   2.350 0.018796 *  
## smoking_age1630 to 39 a week   0.3217600  0.1370252   2.348 0.018902 *  
## smoking_age1640 to 49 a week   0.2909721  0.1384094   2.102 0.035578 *  
## smoking_age1650 to 59 a week   0.3650651  0.1399947   2.608 0.009141 ** 
## smoking_age1660 or more a wk   0.3214966  0.1353845   2.375 0.017599 *  
## alcohol_age16Less than 1 week  0.5857616  0.1778508   3.294 0.000996 ***
## alcohol_age162 to 4 weeks      0.5158579  0.1789393   2.883 0.003957 ** 
## alcohol_age165 to 8 weeks      0.4593757  0.1823387   2.519 0.011787 *  
## alcohol_age169 to 12 weeks     0.3875317  0.1878067   2.063 0.039118 *  
## alcohol_age16Over 12 weeks     0.3637565  0.1828431   1.989 0.046704 *  
## alcohol_age16Do not remember   0.3697804  0.1795718   2.059 0.039521 *  
## alcohol_age16Never had one     0.1569093  0.1838743   0.853 0.393504    
## child_read_age16Sometimes      0.2077748  0.0237090   8.764  < 2e-16 ***
## tv_age16Sometimes             -0.0040151  0.0528784  -0.076 0.939476    
## tv_age16Often                 -0.0241693  0.0504498  -0.479 0.631904    
## school_attendance_age16        0.0016305  0.0003254   5.011 5.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8047 on 5235 degrees of freedom
##   (3773 observations deleted due to missingness)
## Multiple R-squared:  0.2707, Adjusted R-squared:  0.2605 
## F-statistic: 26.61 on 73 and 5235 DF,  p-value: < 2.2e-16

tab_model(global_age16)

	z global cog age 16
Predictors	Estimates	CI	p
(Intercept)	-2.19	-3.45 – -0.94	0.001
sex [Female]	-0.22	-0.26 – -0.17	<0.001
sc prenatal [I]	0.40	0.19 – 0.61	<0.001
sc prenatal [II]	0.36	0.19 – 0.53	<0.001
sc prenatal [III]	0.19	0.03 – 0.34	0.018
sc prenatal [IV]	0.13	-0.04 – 0.29	0.129
sc prenatal [V]	0.01	-0.16 – 0.18	0.877
sc prenatal [Students]	0.54	-0.09 – 1.17	0.093
sc prenatal [Retired]	-0.06	-1.67 – 1.56	0.946
sc prenatal [Single,no husbnd]	0.23	0.02 – 0.44	0.031
smoking preg [5mths-no change]	0.12	-0.09 – 0.33	0.262
smoking preg [Now non-smoker]	0.13	-0.09 – 0.36	0.246
smoking preg [1 to 4 now]	0.02	-0.22 – 0.26	0.850
smoking preg [5 to 9 now]	0.07	-0.18 – 0.33	0.562
smoking preg [10 to 14 now]	0.09	-0.19 – 0.38	0.513
smoking preg [15 to 19 now]	0.04	-0.34 – 0.42	0.845
smoking preg [20 to 24 now]	0.26	-0.15 – 0.68	0.216
smoking preg [25 to 29 now]	0.38	-0.30 – 1.06	0.274
smoking preg [30 or more now]	-0.00	-0.69 – 0.68	0.993
smoking_pregVariable-5mths	0.08	-0.15 – 0.31	0.496
birth weight	0.00	0.00 – 0.00	0.038
parent read level [Occasionally]	0.04	-0.05 – 0.13	0.402
parent read level [Weekly]	0.25	0.16 – 0.33	<0.001
sc age7 [No male head]	-0.01	-0.32 – 0.30	0.954
sc age7 [I]	0.42	0.12 – 0.73	0.007
sc age7 [II]	0.28	-0.00 – 0.55	0.052
sc age7 [III non-manual]	0.26	-0.02 – 0.54	0.067
sc age7 [III manual]	0.09	-0.17 – 0.36	0.489
sc age7 [IV non-manual]	0.07	-0.25 – 0.39	0.659
sc age7 [IV manual]	-0.00	-0.27 – 0.27	0.985
sc age7 [V]	-0.03	-0.31 – 0.26	0.844
breast f [Dont know]	-0.76	-1.92 – 0.40	0.201
breast f [No]	-0.69	-1.81 – 0.44	0.230
breast f [Under one month]	-0.64	-1.76 – 0.48	0.264
breast f [Over one month]	-0.54	-1.67 – 0.58	0.344
sc age11 [I]	0.23	0.07 – 0.39	0.005
sc age11 [II]	0.25	0.14 – 0.36	<0.001
sc age11 [III non manual]	0.22	0.10 – 0.35	0.001
sc age11 [III manual]	0.14	0.05 – 0.23	0.002
sc age11 [IV]	0.04	-0.07 – 0.14	0.487
sc age11 [V]	0.05	-0.09 – 0.19	0.458
sc age11 [No male head]	0.15	-0.00 – 0.30	0.055
tv age11 [Sometimes]	0.15	-0.01 – 0.31	0.059
tv age11 [most days]	0.25	0.10 – 0.40	0.001
child reading age11 [Sometimes]	0.35	0.19 – 0.51	<0.001
child reading age11 [Most days]	0.54	0.38 – 0.70	<0.001
sc age16 [I]	0.27	0.11 – 0.42	0.001
sc age16 [II]	0.16	0.07 – 0.24	<0.001
sc age16 [III non-manual]	0.16	0.05 – 0.26	0.003
sc age16 [III manual]	0.07	0.01 – 0.13	0.034
sc age16 [IV non-manual]	0.09	-0.14 – 0.32	0.427
sc age16 [IV manual]	0.04	-0.05 – 0.13	0.402
sc age16 [V]	-0.05	-0.18 – 0.08	0.418
sc age16 [Unclear]	0.05	-0.20 – 0.30	0.686
smoking age16 [None,dont smoke]	0.68	0.43 – 0.93	<0.001
smoking age16 [Less 1 a week]	0.55	0.27 – 0.83	<0.001
smoking age16 [1 to 9 a week]	0.42	0.15 – 0.68	0.002
smoking age16 [10 to 19 a week]	0.39	0.12 – 0.66	0.004
smoking age16 [20 to 29 a week]	0.32	0.05 – 0.59	0.019
smoking age16 [30 to 39 a week]	0.32	0.05 – 0.59	0.019
smoking age16 [40 to 49 a week]	0.29	0.02 – 0.56	0.036
smoking age16 [50 to 59 a week]	0.37	0.09 – 0.64	0.009
smoking age16 [60 or more a wk]	0.32	0.06 – 0.59	0.018
alcohol age16 [Less than 1 week]	0.59	0.24 – 0.93	0.001
alcohol age16 [2 to 4 weeks]	0.52	0.17 – 0.87	0.004
alcohol age16 [5 to 8 weeks]	0.46	0.10 – 0.82	0.012
alcohol age16 [9 to 12 weeks]	0.39	0.02 – 0.76	0.039
alcohol age16 [Over 12 weeks]	0.36	0.01 – 0.72	0.047
alcohol age16 [Do not remember]	0.37	0.02 – 0.72	0.040
alcohol age16 [Never had one]	0.16	-0.20 – 0.52	0.394
child read age16 [Sometimes]	0.21	0.16 – 0.25	<0.001
tv age16 [Sometimes]	-0.00	-0.11 – 0.10	0.939
tv age16 [Often]	-0.02	-0.12 – 0.07	0.632
school attendance age16	0.00	0.00 – 0.00	<0.001
Observations	5309
R2 / R2 adjusted	0.271 / 0.260

arith_age16 <- lm(z_maths_ability_age16 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11 +
                   sc_age16 + smoking_age16 + alcohol_age16 +
                   child_read_age16 + tv_age16 + school_attendance_age16, data = df_complete_7_16)

summary(arith_age16)

## 
## Call:
## lm(formula = z_maths_ability_age16 ~ sex + sc_prenatal + smoking_preg + 
##     birth_weight + parent_read_level + sc_age7 + breast_f + sc_age11 + 
##     tv_age11 + child_reading_age11 + sc_age16 + smoking_age16 + 
##     alcohol_age16 + child_read_age16 + tv_age16 + school_attendance_age16, 
##     data = df_complete_7_16)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4953 -0.6086 -0.1169  0.5697  2.6204 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -1.8867619  0.6733895  -2.802 0.005099 ** 
## sexFemale                     -0.2708661  0.0241071 -11.236  < 2e-16 ***
## sc_prenatalI                   0.3082450  0.1126553   2.736 0.006237 ** 
## sc_prenatalII                  0.2881823  0.0904827   3.185 0.001456 ** 
## sc_prenatalIII                 0.0701856  0.0822418   0.853 0.393473    
## sc_prenatalIV                  0.0336025  0.0874870   0.384 0.700931    
## sc_prenatalV                  -0.0735968  0.0909749  -0.809 0.418563    
## sc_prenatalStudents            0.3801157  0.3372286   1.127 0.259720    
## sc_prenatalRetired            -0.5481982  0.8688204  -0.631 0.528089    
## sc_prenatalSingle,no husbnd    0.1507726  0.1133450   1.330 0.183507    
## smoking_preg5mths-no change    0.1779728  0.1130545   1.574 0.115497    
## smoking_pregNow non-smoker     0.1342059  0.1207147   1.112 0.266292    
## smoking_preg1 to 4 now         0.0864359  0.1296769   0.667 0.505090    
## smoking_preg5 to 9 now         0.1377830  0.1359040   1.014 0.310713    
## smoking_preg10 to 14 now       0.0930547  0.1521325   0.612 0.540784    
## smoking_preg15 to 19 now       0.0947958  0.2039860   0.465 0.642153    
## smoking_preg20 to 24 now       0.2500078  0.2252652   1.110 0.267120    
## smoking_preg25 to 29 now       0.4003578  0.3654943   1.095 0.273397    
## smoking_preg30 or more now     0.0154052  0.3664626   0.042 0.966470    
## smoking_pregVariable-5mths     0.1032489  0.1221096   0.846 0.397846    
## birth_weight                   0.0003136  0.0001813   1.730 0.083772 .  
## parent_read_levelOccasionally -0.0022678  0.0498910  -0.045 0.963746    
## parent_read_levelWeekly        0.1798358  0.0469395   3.831 0.000129 ***
## sc_age7No male head            0.0063936  0.1659079   0.039 0.969261    
## sc_age7I                       0.3479025  0.1647087   2.112 0.034714 *  
## sc_age7II                      0.1945533  0.1496589   1.300 0.193666    
## sc_age7III non-manual          0.1555995  0.1506416   1.033 0.301693    
## sc_age7III manual              0.0381168  0.1439549   0.265 0.791187    
## sc_age7IV non-manual          -0.0295637  0.1703355  -0.174 0.862217    
## sc_age7IV manual              -0.0426324  0.1460680  -0.292 0.770400    
## sc_age7V                      -0.0414111  0.1527333  -0.271 0.786299    
## breast_fDont know             -0.5553548  0.6248875  -0.889 0.374190    
## breast_fNo                    -0.2911293  0.6036424  -0.482 0.629622    
## breast_fUnder one month       -0.2668352  0.6038248  -0.442 0.658574    
## breast_fOver one month        -0.1632619  0.6036785  -0.270 0.786828    
## sc_age11I                      0.2219384  0.0861996   2.575 0.010060 *  
## sc_age11II                     0.2646350  0.0596022   4.440 9.18e-06 ***
## sc_age11III non manual         0.2604881  0.0676256   3.852 0.000119 ***
## sc_age11III manual             0.1467802  0.0479985   3.058 0.002239 ** 
## sc_age11IV                     0.0679723  0.0543886   1.250 0.211446    
## sc_age11V                      0.1008294  0.0747042   1.350 0.177166    
## sc_age11No male head           0.0974909  0.0814990   1.196 0.231664    
## tv_age11Sometimes              0.1049854  0.0854468   1.229 0.219253    
## tv_age11most days              0.1665095  0.0792880   2.100 0.035771 *  
## child_reading_age11Sometimes   0.2395832  0.0856193   2.798 0.005157 ** 
## child_reading_age11Most days   0.3772555  0.0843855   4.471 7.96e-06 ***
## sc_age16I                      0.3553822  0.0845175   4.205 2.66e-05 ***
## sc_age16II                     0.1575510  0.0449152   3.508 0.000456 ***
## sc_age16III non-manual         0.1607104  0.0564667   2.846 0.004443 ** 
## sc_age16III manual             0.0772373  0.0338069   2.285 0.022373 *  
## sc_age16IV non-manual          0.1221131  0.1227815   0.995 0.319998    
## sc_age16IV manual              0.0651165  0.0486038   1.340 0.180387    
## sc_age16V                     -0.0147199  0.0702378  -0.210 0.834010    
## sc_age16Unclear                0.0734136  0.1353422   0.542 0.587546    
## smoking_age16None,dont smoke   0.5436481  0.1347216   4.035 5.53e-05 ***
## smoking_age16Less 1 a week     0.4067117  0.1505172   2.702 0.006913 ** 
## smoking_age161 to 9 a week     0.2982309  0.1419274   2.101 0.035663 *  
## smoking_age1610 to 19 a week   0.2852239  0.1442215   1.978 0.048017 *  
## smoking_age1620 to 29 a week   0.1957739  0.1441433   1.358 0.174462    
## smoking_age1630 to 39 a week   0.2083018  0.1442762   1.444 0.148863    
## smoking_age1640 to 49 a week   0.0916798  0.1457337   0.629 0.529317    
## smoking_age1650 to 59 a week   0.1900891  0.1474028   1.290 0.197250    
## smoking_age1660 or more a wk   0.1219100  0.1425487   0.855 0.392470    
## alcohol_age16Less than 1 week  0.4416303  0.1872623   2.358 0.018393 *  
## alcohol_age162 to 4 weeks      0.3823805  0.1884083   2.030 0.042455 *  
## alcohol_age165 to 8 weeks      0.3853068  0.1919877   2.007 0.044808 *  
## alcohol_age169 to 12 weeks     0.3208912  0.1977449   1.623 0.104702    
## alcohol_age16Over 12 weeks     0.2969705  0.1925187   1.543 0.122999    
## alcohol_age16Do not remember   0.3102729  0.1890743   1.641 0.100855    
## alcohol_age16Never had one     0.2309384  0.1936045   1.193 0.232988    
## child_read_age16Sometimes      0.1386275  0.0249636   5.553 2.94e-08 ***
## tv_age16Sometimes             -0.0223446  0.0556766  -0.401 0.688195    
## tv_age16Often                 -0.0501302  0.0531194  -0.944 0.345354    
## school_attendance_age16        0.0019117  0.0003426   5.580 2.52e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8473 on 5235 degrees of freedom
##   (3773 observations deleted due to missingness)
## Multiple R-squared:  0.2211, Adjusted R-squared:  0.2102 
## F-statistic: 20.35 on 73 and 5235 DF,  p-value: < 2.2e-16

reading_age16 <- lm(z_reading_ability_age16 ~ sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 + tv_age11 + child_reading_age11 +
                   sc_age16 + smoking_age16 + alcohol_age16 +
                   child_read_age16 + tv_age16 + school_attendance_age16, data = df_complete_7_16)

summary(reading_age16)

## 
## Call:
## lm(formula = z_reading_ability_age16 ~ sex + sc_prenatal + smoking_preg + 
##     birth_weight + parent_read_level + sc_age7 + breast_f + sc_age11 + 
##     tv_age11 + child_reading_age11 + sc_age16 + smoking_age16 + 
##     alcohol_age16 + child_read_age16 + tv_age16 + school_attendance_age16, 
##     data = df_complete_7_16)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4561 -0.4848  0.1047  0.5853  2.1646 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -2.0981228  0.6528153  -3.214 0.001317 ** 
## sexFemale                     -0.1221422  0.0233705  -5.226 1.80e-07 ***
## sc_prenatalI                   0.4174284  0.1092133   3.822 0.000134 ***
## sc_prenatalII                  0.3624943  0.0877182   4.132 3.64e-05 ***
## sc_prenatalIII                 0.2664360  0.0797290   3.342 0.000838 ***
## sc_prenatalIV                  0.1954256  0.0848140   2.304 0.021252 *  
## sc_prenatalV                   0.0978534  0.0881953   1.110 0.267262    
## sc_prenatalStudents            0.5979520  0.3269252   1.829 0.067454 .  
## sc_prenatalRetired             0.4465187  0.8422752   0.530 0.596042    
## sc_prenatalSingle,no husbnd    0.2717467  0.1098819   2.473 0.013427 *  
## smoking_preg5mths-no change    0.0407349  0.1096003   0.372 0.710155    
## smoking_pregNow non-smoker     0.1075191  0.1170265   0.919 0.358264    
## smoking_preg1 to 4 now        -0.0441851  0.1257149  -0.351 0.725249    
## smoking_preg5 to 9 now        -0.0018389  0.1317517  -0.014 0.988865    
## smoking_preg10 to 14 now       0.0783742  0.1474844   0.531 0.595160    
## smoking_preg15 to 19 now      -0.0260384  0.1977536  -0.132 0.895250    
## smoking_preg20 to 24 now       0.2309837  0.2183827   1.058 0.290240    
## smoking_preg25 to 29 now       0.2893292  0.3543273   0.817 0.414218    
## smoking_preg30 or more now    -0.0212510  0.3552660  -0.060 0.952303    
## smoking_pregVariable-5mths     0.0399711  0.1183787   0.338 0.735637    
## birth_weight                   0.0003358  0.0001758   1.910 0.056158 .  
## parent_read_levelOccasionally  0.0744444  0.0483667   1.539 0.123824    
## parent_read_levelWeekly        0.2666293  0.0455053   5.859 4.93e-09 ***
## sc_age7No male head           -0.0229171  0.1608389  -0.142 0.886702    
## sc_age7I                       0.4228142  0.1596764   2.648 0.008123 ** 
## sc_age7II                      0.3065754  0.1450863   2.113 0.034643 *  
## sc_age7III non-manual          0.3205212  0.1460391   2.195 0.028224 *  
## sc_age7III manual              0.1337993  0.1395566   0.959 0.337731    
## sc_age7IV non-manual           0.1590568  0.1651312   0.963 0.335484    
## sc_age7IV manual               0.0379339  0.1416052   0.268 0.788799    
## sc_age7V                      -0.0103431  0.1480668  -0.070 0.944313    
## breast_fDont know             -0.8241516  0.6057952  -1.360 0.173747    
## breast_fNo                    -0.9587332  0.5851992  -1.638 0.101419    
## breast_fUnder one month       -0.8951471  0.5853761  -1.529 0.126280    
## breast_fOver one month        -0.8220416  0.5852342  -1.405 0.160189    
## sc_age11I                      0.1972476  0.0835660   2.360 0.018293 *  
## sc_age11II                     0.1925990  0.0577812   3.333 0.000864 ***
## sc_age11III non manual         0.1419861  0.0655594   2.166 0.030374 *  
## sc_age11III manual             0.1036866  0.0465320   2.228 0.025904 *  
## sc_age11IV                    -0.0027671  0.0527269  -0.052 0.958148    
## sc_age11V                     -0.0053198  0.0724217  -0.073 0.941446    
## sc_age11No male head           0.1718785  0.0790090   2.175 0.029642 *  
## tv_age11Sometimes              0.1727712  0.0828361   2.086 0.037054 *  
## tv_age11most days              0.2841326  0.0768655   3.696 0.000221 ***
## child_reading_age11Sometimes   0.3906219  0.0830033   4.706 2.59e-06 ***
## child_reading_age11Most days   0.6064054  0.0818072   7.413 1.44e-13 ***
## sc_age16I                      0.1285462  0.0819352   1.569 0.116737    
## sc_age16II                     0.1254443  0.0435429   2.881 0.003981 ** 
## sc_age16III non-manual         0.1274959  0.0547415   2.329 0.019894 *  
## sc_age16III manual             0.0465245  0.0327740   1.420 0.155798    
## sc_age16IV non-manual          0.0461408  0.1190302   0.388 0.698299    
## sc_age16IV manual              0.0051575  0.0471188   0.109 0.912845    
## sc_age16V                     -0.0833700  0.0680919  -1.224 0.220866    
## sc_age16Unclear                0.0210825  0.1312071   0.161 0.872351    
## smoking_age16None,dont smoke   0.6852611  0.1306055   5.247 1.61e-07 ***
## smoking_age16Less 1 a week     0.5911945  0.1459184   4.052 5.16e-05 ***
## smoking_age161 to 9 a week     0.4555148  0.1375911   3.311 0.000937 ***
## smoking_age1610 to 19 a week   0.4281823  0.1398151   3.062 0.002206 ** 
## smoking_age1620 to 29 a week   0.3883972  0.1397393   2.779 0.005465 ** 
## smoking_age1630 to 39 a week   0.3758865  0.1398681   2.687 0.007223 ** 
## smoking_age1640 to 49 a week   0.4366099  0.1412811   3.090 0.002010 ** 
## smoking_age1650 to 59 a week   0.4727241  0.1428992   3.308 0.000946 ***
## smoking_age1660 or more a wk   0.4618000  0.1381934   3.342 0.000839 ***
## alcohol_age16Less than 1 week  0.6218800  0.1815408   3.426 0.000618 ***
## alcohol_age162 to 4 weeks      0.5542126  0.1826519   3.034 0.002423 ** 
## alcohol_age165 to 8 weeks      0.4487371  0.1861218   2.411 0.015944 *  
## alcohol_age169 to 12 weeks     0.3827124  0.1917032   1.996 0.045944 *  
## alcohol_age16Over 12 weeks     0.3634668  0.1866367   1.947 0.051533 .  
## alcohol_age16Do not remember   0.3611014  0.1832975   1.970 0.048888 *  
## alcohol_age16Never had one     0.0539465  0.1876893   0.287 0.773799    
## child_read_age16Sometimes      0.2386090  0.0242009   9.860  < 2e-16 ***
## tv_age16Sometimes              0.0150546  0.0539755   0.279 0.780320    
## tv_age16Often                  0.0062483  0.0514965   0.121 0.903431    
## school_attendance_age16        0.0010486  0.0003321   3.157 0.001601 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8214 on 5235 degrees of freedom
##   (3773 observations deleted due to missingness)
## Multiple R-squared:  0.2374, Adjusted R-squared:  0.2268 
## F-statistic: 22.33 on 73 and 5235 DF,  p-value: < 2.2e-16

library(broom)
library(dplyr)

extract_age7_11_16_exposures <- function(model, outcome_label) {

  coef_df <- tidy(model, conf.int = TRUE) %>%
    filter(term %in% c(
      # Age 7
      "parent_read_levelOccasionally",
      "parent_read_levelWeekly",
      
      # Age 11
      "child_reading_age11Sometimes",
      "child_reading_age11Most days",
      "tv_age11Sometimes",
      "tv_age11most days",
      
      # Age 16
      "child_read_age16Sometimes",
      "tv_age16Sometimes"
    )) %>%
    mutate(
      exposure = case_when(
        term == "parent_read_levelOccasionally" ~ "Reading: Sometimes",
        term == "parent_read_levelWeekly"       ~ "Reading: Most days",
        
        term == "child_reading_age11Sometimes"  ~ "Reading: Sometimes",
        term == "child_reading_age11Most days"  ~ "Reading: Most days",
        
        term == "child_read_age16Sometimes"     ~ "Reading: Sometimes",
        
        term == "tv_age11Sometimes"             ~ "TV: Sometimes",
        term == "tv_age11most days"             ~ "TV: Most days",
        term == "tv_age16Sometimes"             ~ "TV: Sometimes"
      ),
      domain = case_when(
        grepl("^tv_", term) ~ "Passive SB (TV)",
        TRUE                ~ "Cognitively engaged SB (Reading)"
      ),
      age = case_when(
        grepl("^parent_", term)         ~ "Age 7",
        grepl("age11", term)            ~ "Age 11",
        grepl("age16", term)            ~ "Age 16"
      ),
      outcome = outcome_label
    ) %>%
    dplyr::select(outcome, domain, age, exposure, estimate, conf.low, conf.high)

  # referências
  ref_df <- tibble(
    outcome  = outcome_label,
    domain   = c(
      "Cognitively engaged SB (Reading)",
      "Cognitively engaged SB (Reading)",
      "Cognitively engaged SB (Reading)",
      "Passive SB (TV)",
      "Passive SB (TV)"
    ),
    age      = c("Age 7", "Age 11", "Age 16", "Age 11", "Age 16"),
    exposure = c(
      "Reading: Hardly ever",
      "Reading: Hardly ever",
      "Reading: Hardly ever",
      "TV: Hardly ever",
      "TV: Hardly ever"
    ),
    estimate = 0,
    conf.low = 0,
    conf.high = 0
  )

  bind_rows(ref_df, coef_df)
}

plot_df <- bind_rows(
  extract_age7_11_16_exposures(global_age16,  "Global cognition"),
  extract_age7_11_16_exposures(arith_age16,   "Maths ability"),
  extract_age7_11_16_exposures(reading_age16, "Reading ability")
)


plot_df$dose <- factor(
  gsub(".*: ", "", plot_df$exposure),
  levels = c("Most days", "Sometimes", "Hardly ever"),
  ordered = TRUE
)

plot_df$age <- factor(plot_df$age, 
  levels = c("Age 7", "Age 11", "Age 16"),
  ordered = TRUE
)

library(ggplot2)
library(ggsci)
png("Figure_read_tv_age16.tiff", units = "in", width = 10, height = 14, res=300)
plot_age16 <- ggplot(plot_df,
       aes(x = estimate,
           y = outcome,
           colour = age,
           shape = dose)) +

  geom_vline(xintercept = 0,
             linetype = "dashed",
             color = "grey50") +

  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high),
    position = position_dodge(width = 0.6),
    size = 0.9
  ) +

  facet_wrap(~ domain, ncol = 2, scales = "free_y") +

  scale_color_jama() +

  labs(
    x = "Standardized difference in cognitive score (z)",
    y = NULL,
    colour = "Age",
    shape  = "Exposure frequency"
  ) +

  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    strip.text = element_text(face = "bold")
  )

plot_age16
dev.off()

## quartz_off_screen 
##                 2

plot_age16

DescTools::Desc(df_complete_7_16$z_global_cog_age16)

## ────────────────────────────────────────────────────────────────────────────── 
## df_complete_7_16$z_global_cog_age16 (numeric)
## 
##        length            n              NAs      unique          0s'
##         9'082        9'082                0         667           0
##                     100.0%             0.0%                    0.0%
##                                                                    
##           .05          .10              .25      median         .75
##   -1.65856251  -1.32768872      -0.69028072  0.02376084  0.75550390
##                                                                    
##         range           sd            vcoef         mad         IQR
##    5.28129608   1.00000000  -3.36053075e+14  1.07832120  1.44578462
##                                                                    
##              mean       meanCI
##   -2.97572043e-15  -0.02056910
##                     0.02056910
##                               
##               .90          .95
##        1.36887466   1.54617901
##                               
##              skew         kurt
##       -0.16693366  -0.53942019
##                               
## lowest : -3.10213444 (4), -3.02107551 (3), -2.94001658 (6), -2.86117033 (3), -2.85895765 (6)
## highest: 2.01925647 (17), 2.02146916 (8), 2.09810271, 2.1003154 (4), 2.17916164 (2)
## 
## ' 95%-CI (classic)

COGNITIVE DECLINE

library(dplyr)
library(tidyr)
library(lme4)
library(emmeans)
library(ggplot2)
library(ggsci)

df_long <- df_complete_7_16 %>%
  dplyr::select(
    ncdsid,
    z_global_cog_age7,
    z_global_cog_age11,
    z_global_cog_age16,
    parent_read_level,
    sex,
    sc_prenatal,
    smoking_preg,
    birth_weight,
    sc_age7,
    breast_f,
    sc_age11,
    sc_age16,
    smoking_age16,
    alcohol_age16,
    school_attendance_age16
  ) %>%
  pivot_longer(
    cols = starts_with("z_global_cog"),
    names_to = "wave",
    values_to = "z_global_cog"
  ) %>%
  mutate(
    age = case_when(
      wave == "z_global_cog_age7"  ~ 7,
      wave == "z_global_cog_age11" ~ 11,
      wave == "z_global_cog_age16" ~ 16
    ),
    id_num = as.numeric(as.factor(ncdsid)),
    parent_read_level = factor(
      parent_read_level,
      levels = c("Hardly ever", "Occasionally", "Weekly")
    )
  ) %>%
  filter(!is.na(z_global_cog), !is.na(parent_read_level))

mixed <- lmer(
  z_global_cog ~ age * parent_read_level +
    sex + sc_prenatal + smoking_preg +
    birth_weight + sc_age7 + breast_f +
    sc_age11 + sc_age16 +
    smoking_age16 + alcohol_age16 +
    school_attendance_age16 +
    (1 | id_num),
  data = df_long,
  REML = TRUE
)

summary(mixed)

## Linear mixed model fit by REML ['lmerMod']
## Formula: z_global_cog ~ age * parent_read_level + sex + sc_prenatal +  
##     smoking_preg + birth_weight + sc_age7 + breast_f + sc_age11 +  
##     sc_age16 + smoking_age16 + alcohol_age16 + school_attendance_age16 +  
##     (1 | id_num)
##    Data: df_long
## 
## REML criterion at convergence: 57248.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1707 -0.5580  0.0087  0.5556  3.7615 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id_num   (Intercept) 0.4748   0.6891  
##  Residual             0.3104   0.5571  
## Number of obs: 25419, groups:  id_num, 8473
## 
## Fixed effects:
##                                     Estimate Std. Error t value
## (Intercept)                       -0.4771652  0.4726989  -1.009
## age                               -0.0276125  0.0035877  -7.696
## parent_read_levelOccasionally      0.0588150  0.0588333   1.000
## parent_read_levelWeekly           -0.1367652  0.0543543  -2.516
## sexFemale                         -0.0138493  0.0169086  -0.819
## sc_prenatalI                       0.3405575  0.0759160   4.486
## sc_prenatalII                      0.2687267  0.0632984   4.245
## sc_prenatalIII                     0.1299338  0.0574254   2.263
## sc_prenatalIV                      0.0674120  0.0613240   1.099
## sc_prenatalV                      -0.0431959  0.0636922  -0.678
## sc_prenatalStudents                0.4447746  0.2627576   1.693
## sc_prenatalRetired                 0.4148887  0.5448411   0.761
## sc_prenatalSingle,no husbnd       -0.0254288  0.0789177  -0.322
## smoking_preg5mths-no change       -0.0119955  0.0829882  -0.145
## smoking_pregNow non-smoker         0.0243965  0.0880821   0.277
## smoking_preg1 to 4 now            -0.0559987  0.0944689  -0.593
## smoking_preg5 to 9 now            -0.0618972  0.0979754  -0.632
## smoking_preg10 to 14 now          -0.0335994  0.1115498  -0.301
## smoking_preg15 to 19 now          -0.1950853  0.1375448  -1.418
## smoking_preg20 to 24 now          -0.3710597  0.1638393  -2.265
## smoking_preg25 to 29 now           0.3684535  0.2676413   1.377
## smoking_preg30 or more now        -0.4230639  0.2822802  -1.499
## smoking_pregVariable-5mths        -0.0720824  0.0891213  -0.809
## birth_weight                       0.0004528  0.0001247   3.632
## sc_age7No male head               -0.0412713  0.1159071  -0.356
## sc_age7I                           0.2906540  0.1131384   2.569
## sc_age7II                          0.1834090  0.1042092   1.760
## sc_age7III non-manual              0.2010335  0.1050740   1.913
## sc_age7III manual                 -0.0121083  0.1004322  -0.121
## sc_age7IV non-manual               0.0346083  0.1186395   0.292
## sc_age7IV manual                  -0.0584296  0.1020653  -0.572
## sc_age7V                          -0.1804525  0.1066851  -1.691
## breast_fDont know                 -0.6295774  0.4553967  -1.382
## breast_fNo                        -0.5320111  0.4403901  -1.208
## breast_fUnder one month           -0.4635105  0.4404788  -1.052
## breast_fOver one month            -0.3665665  0.4403660  -0.832
## sc_age11I                          0.2734182  0.0596603   4.583
## sc_age11II                         0.2189058  0.0415800   5.265
## sc_age11III non manual             0.1930274  0.0470315   4.104
## sc_age11III manual                 0.0665467  0.0337184   1.974
## sc_age11IV                         0.0022970  0.0383053   0.060
## sc_age11V                         -0.0152743  0.0532941  -0.287
## sc_age11No male head               0.0262475  0.0575054   0.456
## sc_age16I                          0.1850237  0.0563972   3.281
## sc_age16II                         0.0976992  0.0318229   3.070
## sc_age16III non-manual             0.0210368  0.0398973   0.527
## sc_age16III manual                 0.0422884  0.0241738   1.749
## sc_age16IV non-manual              0.0453397  0.0837153   0.542
## sc_age16IV manual                  0.0113622  0.0341105   0.333
## sc_age16V                         -0.1004395  0.0509113  -1.973
## sc_age16Unclear                    0.2001783  0.0913158   2.192
## smoking_age16None,dont smoke       0.4046931  0.0902181   4.486
## smoking_age16Less 1 a week         0.3069917  0.1009775   3.040
## smoking_age161 to 9 a week         0.1966095  0.0956604   2.055
## smoking_age1610 to 19 a week       0.1349575  0.0973387   1.386
## smoking_age1620 to 29 a week       0.1157112  0.0972128   1.190
## smoking_age1630 to 39 a week       0.1650466  0.0975903   1.691
## smoking_age1640 to 49 a week       0.1027948  0.0987487   1.041
## smoking_age1650 to 59 a week       0.2151257  0.1009549   2.131
## smoking_age1660 or more a wk       0.2029623  0.0960843   2.112
## alcohol_age16Less than 1 week      0.4775481  0.1012484   4.717
## alcohol_age162 to 4 weeks          0.4084158  0.1024215   3.988
## alcohol_age165 to 8 weeks          0.3499986  0.1052116   3.327
## alcohol_age169 to 12 weeks         0.2919475  0.1100385   2.653
## alcohol_age16Over 12 weeks         0.2438006  0.1055246   2.310
## alcohol_age16Do not remember       0.1955057  0.1029448   1.899
## alcohol_age16Never had one        -0.1264345  0.1067089  -1.185
## school_attendance_age16            0.0012446  0.0002428   5.127
## age:parent_read_levelOccasionally  0.0006549  0.0041040   0.160
## age:parent_read_levelWeekly        0.0309402  0.0037620   8.224

tab_model(mixed, digits = 3, digits.p = 3)

	z global cog
Predictors	Estimates	CI	p
(Intercept)	-0.477	-1.404 – 0.449	0.313
age	-0.028	-0.035 – -0.021	<0.001
parent read level [Occasionally]	0.059	-0.057 – 0.174	0.317
parent read level [Weekly]	-0.137	-0.243 – -0.030	0.012
sex [Female]	-0.014	-0.047 – 0.019	0.413
sc prenatal [I]	0.341	0.192 – 0.489	<0.001
sc prenatal [II]	0.269	0.145 – 0.393	<0.001
sc prenatal [III]	0.130	0.017 – 0.242	0.024
sc prenatal [IV]	0.067	-0.053 – 0.188	0.272
sc prenatal [V]	-0.043	-0.168 – 0.082	0.498
sc prenatal [Students]	0.445	-0.070 – 0.960	0.091
sc prenatal [Retired]	0.415	-0.653 – 1.483	0.446
sc prenatal [Single,no husbnd]	-0.025	-0.180 – 0.129	0.747
smoking preg [5mths-no change]	-0.012	-0.175 – 0.151	0.885
smoking preg [Now non-smoker]	0.024	-0.148 – 0.197	0.782
smoking preg [1 to 4 now]	-0.056	-0.241 – 0.129	0.553
smoking preg [5 to 9 now]	-0.062	-0.254 – 0.130	0.528
smoking preg [10 to 14 now]	-0.034	-0.252 – 0.185	0.763
smoking preg [15 to 19 now]	-0.195	-0.465 – 0.075	0.156
smoking preg [20 to 24 now]	-0.371	-0.692 – -0.050	0.024
smoking preg [25 to 29 now]	0.368	-0.156 – 0.893	0.169
smoking preg [30 or more now]	-0.423	-0.976 – 0.130	0.134
smoking_pregVariable-5mths	-0.072	-0.247 – 0.103	0.419
birth weight	0.000	0.000 – 0.001	<0.001
sc age7 [No male head]	-0.041	-0.268 – 0.186	0.722
sc age7 [I]	0.291	0.069 – 0.512	0.010
sc age7 [II]	0.183	-0.021 – 0.388	0.078
sc age7 [III non-manual]	0.201	-0.005 – 0.407	0.056
sc age7 [III manual]	-0.012	-0.209 – 0.185	0.904
sc age7 [IV non-manual]	0.035	-0.198 – 0.267	0.771
sc age7 [IV manual]	-0.058	-0.258 – 0.142	0.567
sc age7 [V]	-0.180	-0.390 – 0.029	0.091
breast f [Dont know]	-0.630	-1.522 – 0.263	0.167
breast f [No]	-0.532	-1.395 – 0.331	0.227
breast f [Under one month]	-0.464	-1.327 – 0.400	0.293
breast f [Over one month]	-0.367	-1.230 – 0.497	0.405
sc age11 [I]	0.273	0.156 – 0.390	<0.001
sc age11 [II]	0.219	0.137 – 0.300	<0.001
sc age11 [III non manual]	0.193	0.101 – 0.285	<0.001
sc age11 [III manual]	0.067	0.000 – 0.133	0.048
sc age11 [IV]	0.002	-0.073 – 0.077	0.952
sc age11 [V]	-0.015	-0.120 – 0.089	0.774
sc age11 [No male head]	0.026	-0.086 – 0.139	0.648
sc age16 [I]	0.185	0.074 – 0.296	0.001
sc age16 [II]	0.098	0.035 – 0.160	0.002
sc age16 [III non-manual]	0.021	-0.057 – 0.099	0.598
sc age16 [III manual]	0.042	-0.005 – 0.090	0.080
sc age16 [IV non-manual]	0.045	-0.119 – 0.209	0.588
sc age16 [IV manual]	0.011	-0.055 – 0.078	0.739
sc age16 [V]	-0.100	-0.200 – -0.001	0.049
sc age16 [Unclear]	0.200	0.021 – 0.379	0.028
smoking age16 [None,dont smoke]	0.405	0.228 – 0.582	<0.001
smoking age16 [Less 1 a week]	0.307	0.109 – 0.505	0.002
smoking age16 [1 to 9 a week]	0.197	0.009 – 0.384	0.040
smoking age16 [10 to 19 a week]	0.135	-0.056 – 0.326	0.166
smoking age16 [20 to 29 a week]	0.116	-0.075 – 0.306	0.234
smoking age16 [30 to 39 a week]	0.165	-0.026 – 0.356	0.091
smoking age16 [40 to 49 a week]	0.103	-0.091 – 0.296	0.298
smoking age16 [50 to 59 a week]	0.215	0.017 – 0.413	0.033
smoking age16 [60 or more a wk]	0.203	0.015 – 0.391	0.035
alcohol age16 [Less than 1 week]	0.478	0.279 – 0.676	<0.001
alcohol age16 [2 to 4 weeks]	0.408	0.208 – 0.609	<0.001
alcohol age16 [5 to 8 weeks]	0.350	0.144 – 0.556	0.001
alcohol age16 [9 to 12 weeks]	0.292	0.076 – 0.508	0.008
alcohol age16 [Over 12 weeks]	0.244	0.037 – 0.451	0.021
alcohol age16 [Do not remember]	0.196	-0.006 – 0.397	0.058
alcohol age16 [Never had one]	-0.126	-0.336 – 0.083	0.236
school attendance age16	0.001	0.001 – 0.002	<0.001
age × parent read level [Occasionally]	0.001	-0.007 – 0.009	0.873
age × parent read level [Weekly]	0.031	0.024 – 0.038	<0.001
Random Effects
σ2	0.31
τ00 id_num	0.47
ICC	0.60
N id_num	8473
Observations	25419
Marginal R2 / Conditional R2	0.189 / 0.680

emm <- emmeans(
  mixed,
  ~ parent_read_level | age,
  at = list(age = c(7, 11, 16)),
  nuisance = c(
    "sex", "sc_prenatal", "smoking_preg", "sc_age7", "sc_age11", "sc_age16",
    "smoking_age16", "alcohol_age16", "school_attendance_age16", "breast_f"
  )
)

emm_df <- as.data.frame(emm)

plot_trajectory_7_16 <- ggplot(
  emm_df,
  aes(
    x = age,
    y = emmean,
    colour = parent_read_level,
    fill   = parent_read_level
  )
) +
  geom_line(size = 1.1) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
  scale_x_continuous(breaks = c(7, 11, 16)) +
  labs(
    x = "Age (years)",
    y = "Predicted global cognitive score (z)",
    colour = "Parental reading frequency",
    fill   = "Parental reading frequency"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  ) +
  scale_color_bmj() +
  scale_fill_bmj()

plot_trajectory_7_16

GROUP-BASED TRAJECTORY MODELLING

library(lcmm)

df_read_long <- df_complete_7_16 %>%
  mutate(
    reading_age7  = parent_read_level,
    reading_age11 = child_reading_age11,
    reading_age16 = child_read_age16
  ) %>%
  mutate(
    read7 = as.numeric(
      factor(
        reading_age7,
        levels = c("Hardly ever", "Occasionally", "Weekly")
      )
    ) - 1,
    read11 = as.numeric(
      factor(
        reading_age11,
        levels = c("Hardly ever", "Sometimes", "Most days")
      )
    ) - 1,
    read16 = as.numeric(
      factor(
        reading_age16,
        levels = c("Hardly ever", "Sometimes", "Most days")
      )
    ) - 1
  ) %>%
  dplyr::select(
    ncdsid,
    read7,
    read11,
    read16
  ) %>%
  pivot_longer(
    cols = starts_with("read"),
    names_to = "wave",
    values_to = "reading_num"
  ) %>%
  mutate(
    age = case_when(
      wave == "read7"  ~ 7,
      wave == "read11" ~ 11,
      wave == "read16" ~ 16
    ),
    id_num = as.numeric(as.factor(ncdsid))
  ) %>%
  filter(!is.na(reading_num))

m1 <- hlme(
  fixed   = reading_num ~ age,
  random  = ~ 1,
  subject = "id_num",
  ng      = 1,
  data    = df_read_long
)

summary(m1)

## Heterogenous linear mixed model 
##      fitted by maximum likelihood method 
##  
## hlme(fixed = reading_num ~ age, random = ~1, subject = "id_num", 
##     ng = 1, data = df_read_long)
##  
## Statistical Model: 
##      Dataset: df_read_long 
##      Number of subjects: 9080 
##      Number of observations: 23579 
##      Number of latent classes: 1 
##      Number of parameters: 4  
##  
## Iteration process: 
##      Convergence criteria satisfied 
##      Number of iterations:  34 
##      Convergence criteria: parameters= 5e-13 
##                          : likelihood= 4.9e-07 
##                          : second derivatives= 7.4e-06 
##  
## Goodness-of-fit statistics: 
##      maximum log-likelihood: -21067.55  
##      AIC: 42143.11  
##      BIC: 42171.56  
##  
##  
## Maximum Likelihood Estimates: 
##  
## Fixed effects in the longitudinal model:
## 
##               coef      Se    Wald p-value
## intercept  2.55734 0.01250 204.535 0.00000
## age       -0.10674 0.00110 -96.814 0.00000
## 
## 
## Variance-covariance matrix of the random-effects:
##           intercept
## intercept         0
## 
##                              coef      Se
## Residual standard error:  0.59129 0.00272

m2 <- hlme(
  fixed   = reading_num ~ age,
  mixture = ~ age,
  random  = ~ 1,
  subject = "id_num",
  ng      = 2,
  data    = df_read_long,
  nwg     = TRUE,
  B       = m1
)

summary(m2)

## Heterogenous linear mixed model 
##      fitted by maximum likelihood method 
##  
## hlme(fixed = reading_num ~ age, mixture = ~age, random = ~1, 
##     subject = "id_num", ng = 2, nwg = TRUE, data = df_read_long)
##  
## Statistical Model: 
##      Dataset: df_read_long 
##      Number of subjects: 9080 
##      Number of observations: 23579 
##      Number of latent classes: 2 
##      Number of parameters: 8  
##  
## Iteration process: 
##      Convergence criteria satisfied 
##      Number of iterations:  206 
##      Convergence criteria: parameters= 4.8e-10 
##                          : likelihood= 1.6e-06 
##                          : second derivatives= 9e-07 
##  
## Goodness-of-fit statistics: 
##      maximum log-likelihood: -20627.31  
##      AIC: 41270.62  
##      BIC: 41327.53  
##  
##  
## Maximum Likelihood Estimates: 
##  
## Fixed effects in the class-membership model:
## (the class of reference is the last class) 
## 
##                      coef       Se     Wald p-value
## intercept class1 -2.21667  0.05951  -37.249 0.00000
## 
## Fixed effects in the longitudinal model:
## 
##                      coef       Se     Wald p-value
## intercept class1  0.55147  0.05724    9.635 0.00000
## intercept class2  2.77983  0.01391  199.814 0.00000
## age class1        0.02617  0.00484    5.404 0.00000
## age class2       -0.12158  0.00117 -104.255 0.00000
## 
## 
## Variance-covariance matrix of the random-effects:
##           intercept
## intercept         0
## 
##                                     coef       Se
## Proportional coefficient class1  0.95921 11.73219
## Residual standard error:         0.53844  0.00279

postprob(m2)

##  
## Posterior classification: 
##   class1  class2
## N 660.00 8420.00
## %   7.27   92.73
##  
## Posterior classification table: 
##      --> mean of posterior probabilities in each class 
##         prob1  prob2
## class1 0.8537 0.1463
## class2 0.0390 0.9610
##  
## Posterior probabilities above a threshold (%): 
##          class1 class2
## prob>0.7  98.79  95.63
## prob>0.8  83.64  92.66
## prob>0.9  33.03  91.69
## 

age_grid <- data.frame(
  age = seq(7, 16, by = 0.1)
)

pred <- predictY(
  m2,
  newdata = age_grid,
  var.time = "age",
  draws = TRUE,
  ndraws = 1000
)

str(pred)

## List of 2
##  $ pred : num [1:91, 1:6] 0.735 0.737 0.74 0.743 0.745 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:6] "Ypred_class1" "Ypred_class2" "lower.Ypred_class1" "lower.Ypred_class2" ...
##  $ times:'data.frame':   91 obs. of  1 variable:
##   ..$ age: num [1:91] 7 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 ...
##  - attr(*, "class")= chr "predictY"

traj_df <- bind_rows(

  data.frame(
    age   = pred$times$age,
    mean  = pred$pred[, "Ypred_class1"],
    lower = pred$pred[, "lower.Ypred_class1"],
    upper = pred$pred[, "upper.Ypred_class1"],
    class = "Class 1"
  ),

  data.frame(
    age   = pred$times$age,
    mean  = pred$pred[, "Ypred_class2"],
    lower = pred$pred[, "lower.Ypred_class2"],
    upper = pred$pred[, "upper.Ypred_class2"],
    class = "Class 2"
  )
)

traj_df$class <- factor(
  traj_df$class,
  levels = c("Class 1", "Class 2"),
  labels = c(
    "High initial / rapidly decreasing",
    "Moderate initial / slowly increasing"
  )
)

fig_traj_reading <- ggplot(traj_df, aes(x = age, y = mean, colour = class, fill = class)) +

  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.20,
    colour = NA
  ) +

  geom_line(size = 1.2) +

  scale_x_continuous(
    breaks = c(7, 11, 16)
  ) +

  scale_y_continuous(
    breaks = c(0, 1, 2),
    labels = c("Low", "Moderate", "High")
  ) +

  labs(
    x = "Age (years)",
    y = "Reading frequency",
    colour = NULL,
    fill   = NULL
  ) +

  theme_classic(base_size = 14) +
  theme(
    legend.position = "top",
    axis.title = element_text(face = "bold")
  ) +
  coord_cartesian(ylim = c(0,2)) +
  scale_color_bmj() +
  scale_fill_bmj()

fig_traj_reading

GBTM + COGNITIVE DECLINE

library(lme4)
library(ggeffects)
library(ggplot2)
library(ggsci)

pp <- postprob(m2)

##  
## Posterior classification: 
##   class1  class2
## N 660.00 8420.00
## %   7.27   92.73
##  
## Posterior classification table: 
##      --> mean of posterior probabilities in each class 
##         prob1  prob2
## class1 0.8537 0.1463
## class2 0.0390 0.9610
##  
## Posterior probabilities above a threshold (%): 
##          class1 class2
## prob>0.7  98.79  95.63
## prob>0.8  83.64  92.66
## prob>0.9  33.03  91.69
## 

str(pp)

## List of 3
##  $ : num [1:2, 1:2] 660 7.27 8420 92.73
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "N" "%"
##   .. ..$ : chr [1:2] "class1" "class2"
##  $ : num [1:2, 1:2] 0.854 0.039 0.146 0.961
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "class1" "class2"
##   .. ..$ : chr [1:2] "prob1" "prob2"
##  $ : num [1:3, 1:2] 98.8 83.6 33 95.6 92.7 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "prob>0.7" "prob>0.8" "prob>0.9"
##   .. ..$ : chr [1:2] "class1" "class2"

class_df <- m2$pprob %>%
  mutate(
    read_class = factor(
      class,
      labels = c(
        "High initial / rapidly decreasing",
        "Moderate initial / slowly increasing"
      )
    )
  ) %>%
  dplyr::select(id_num, read_class)

table(class_df$read_class)

## 
##    High initial / rapidly decreasing Moderate initial / slowly increasing 
##                                  660                                 8420

df_cog <- df_long %>%
  left_join(class_df, by = "id_num")

mixed_cog <- lmer(
  z_global_cog ~ age * read_class + 
                sex + sc_prenatal  + smoking_preg + 
                 birth_weight +  parent_read_level  + sc_age7  + breast_f +
                   sc_age11 +sc_age16 + smoking_age16 + alcohol_age16 +
                   school_attendance_age16 + (1 | ncdsid), data = df_cog
)
summary(mixed_cog)

## Linear mixed model fit by REML ['lmerMod']
## Formula: z_global_cog ~ age * read_class + sex + sc_prenatal + smoking_preg +  
##     birth_weight + parent_read_level + sc_age7 + breast_f + sc_age11 +  
##     sc_age16 + smoking_age16 + alcohol_age16 + school_attendance_age16 +  
##     (1 | ncdsid)
##    Data: df_cog
## 
## REML criterion at convergence: 57396.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0857 -0.5576  0.0064  0.5539  3.7452 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ncdsid   (Intercept) 0.4729   0.6877  
##  Residual             0.3135   0.5599  
## Number of obs: 25419, groups:  ncdsid, 8473
## 
## Fixed effects:
##                                                      Estimate Std. Error
## (Intercept)                                        -0.2841875  0.4749314
## age                                                -0.0276449  0.0035011
## read_classModerate initial / slowly increasing      0.2255078  0.1359681
## sexFemale                                          -0.0140386  0.0168950
## sc_prenatalI                                        0.3427848  0.0758572
## sc_prenatalII                                       0.2702546  0.0632488
## sc_prenatalIII                                      0.1323723  0.0573828
## sc_prenatalIV                                       0.0685472  0.0612754
## sc_prenatalV                                       -0.0374185  0.0636590
## sc_prenatalStudents                                 0.4469663  0.2625468
## sc_prenatalRetired                                  0.4132950  0.5444029
## sc_prenatalSingle,no husbnd                        -0.0244175  0.0788547
## smoking_preg5mths-no change                        -0.0210636  0.0829555
## smoking_pregNow non-smoker                          0.0147442  0.0880476
## smoking_preg1 to 4 now                             -0.0674254  0.0944404
## smoking_preg5 to 9 now                             -0.0710731  0.0979262
## smoking_preg10 to 14 now                           -0.0452350  0.1115018
## smoking_preg15 to 19 now                           -0.1886392  0.1374445
## smoking_preg20 to 24 now                           -0.3840253  0.1637427
## smoking_preg25 to 29 now                            0.3548898  0.2674496
## smoking_preg30 or more now                         -0.4344657  0.2820689
## smoking_pregVariable-5mths                         -0.0818783  0.0890867
## birth_weight                                        0.0004520  0.0001246
## parent_read_levelOccasionally                      -0.4186086  0.1321426
## parent_read_levelWeekly                            -0.2797175  0.1337540
## sc_age7No male head                                -0.0407403  0.1158139
## sc_age7I                                            0.2942855  0.1130514
## sc_age7II                                           0.1885984  0.1041343
## sc_age7III non-manual                               0.2042509  0.1049929
## sc_age7III manual                                  -0.0074551  0.1003588
## sc_age7IV non-manual                                0.0346207  0.1185440
## sc_age7IV manual                                   -0.0541220  0.1019894
## sc_age7V                                           -0.1745230  0.1066106
## breast_fDont know                                  -0.7807967  0.4567551
## breast_fNo                                         -0.6943870  0.4420912
## breast_fUnder one month                            -0.6259690  0.4421815
## breast_fOver one month                             -0.5290790  0.4420707
## sc_age11I                                           0.2685328  0.0596261
## sc_age11II                                          0.2127385  0.0415780
## sc_age11III non manual                              0.1874149  0.0470167
## sc_age11III manual                                  0.0602206  0.0337321
## sc_age11IV                                         -0.0021447  0.0382922
## sc_age11V                                          -0.0240777  0.0533012
## sc_age11No male head                                0.0217706  0.0574711
## sc_age16I                                           0.1849779  0.0563518
## sc_age16II                                          0.0974785  0.0317973
## sc_age16III non-manual                              0.0214950  0.0398654
## sc_age16III manual                                  0.0428163  0.0241547
## sc_age16IV non-manual                               0.0532748  0.0836738
## sc_age16IV manual                                   0.0105365  0.0340838
## sc_age16V                                          -0.1005486  0.0508703
## sc_age16Unclear                                     0.2038535  0.0912474
## smoking_age16None,dont smoke                        0.3977581  0.0901639
## smoking_age16Less 1 a week                          0.3002868  0.1009115
## smoking_age161 to 9 a week                          0.1899048  0.0955996
## smoking_age1610 to 19 a week                        0.1263498  0.0972866
## smoking_age1620 to 29 a week                        0.1073334  0.0971594
## smoking_age1630 to 39 a week                        0.1608836  0.0975179
## smoking_age1640 to 49 a week                        0.0937298  0.0986979
## smoking_age1650 to 59 a week                        0.2108023  0.1008800
## smoking_age1660 or more a wk                        0.1947475  0.0960312
## alcohol_age16Less than 1 week                       0.4613152  0.1012565
## alcohol_age162 to 4 weeks                           0.3928616  0.1024203
## alcohol_age165 to 8 weeks                           0.3325654  0.1052263
## alcohol_age169 to 12 weeks                          0.2739161  0.1100516
## alcohol_age16Over 12 weeks                          0.2267445  0.1055346
## alcohol_age16Do not remember                        0.1809893  0.1029324
## alcohol_age16Never had one                         -0.1438544  0.1067209
## school_attendance_age16                             0.0012460  0.0002426
## age:read_classModerate initial / slowly increasing  0.0236965  0.0036387
##                                                    t value
## (Intercept)                                         -0.598
## age                                                 -7.896
## read_classModerate initial / slowly increasing       1.659
## sexFemale                                           -0.831
## sc_prenatalI                                         4.519
## sc_prenatalII                                        4.273
## sc_prenatalIII                                       2.307
## sc_prenatalIV                                        1.119
## sc_prenatalV                                        -0.588
## sc_prenatalStudents                                  1.702
## sc_prenatalRetired                                   0.759
## sc_prenatalSingle,no husbnd                         -0.310
## smoking_preg5mths-no change                         -0.254
## smoking_pregNow non-smoker                           0.167
## smoking_preg1 to 4 now                              -0.714
## smoking_preg5 to 9 now                              -0.726
## smoking_preg10 to 14 now                            -0.406
## smoking_preg15 to 19 now                            -1.372
## smoking_preg20 to 24 now                            -2.345
## smoking_preg25 to 29 now                             1.327
## smoking_preg30 or more now                          -1.540
## smoking_pregVariable-5mths                          -0.919
## birth_weight                                         3.628
## parent_read_levelOccasionally                       -3.168
## parent_read_levelWeekly                             -2.091
## sc_age7No male head                                 -0.352
## sc_age7I                                             2.603
## sc_age7II                                            1.811
## sc_age7III non-manual                                1.945
## sc_age7III manual                                   -0.074
## sc_age7IV non-manual                                 0.292
## sc_age7IV manual                                    -0.531
## sc_age7V                                            -1.637
## breast_fDont know                                   -1.709
## breast_fNo                                          -1.571
## breast_fUnder one month                             -1.416
## breast_fOver one month                              -1.197
## sc_age11I                                            4.504
## sc_age11II                                           5.117
## sc_age11III non manual                               3.986
## sc_age11III manual                                   1.785
## sc_age11IV                                          -0.056
## sc_age11V                                           -0.452
## sc_age11No male head                                 0.379
## sc_age16I                                            3.283
## sc_age16II                                           3.066
## sc_age16III non-manual                               0.539
## sc_age16III manual                                   1.773
## sc_age16IV non-manual                                0.637
## sc_age16IV manual                                    0.309
## sc_age16V                                           -1.977
## sc_age16Unclear                                      2.234
## smoking_age16None,dont smoke                         4.412
## smoking_age16Less 1 a week                           2.976
## smoking_age161 to 9 a week                           1.986
## smoking_age1610 to 19 a week                         1.299
## smoking_age1620 to 29 a week                         1.105
## smoking_age1630 to 39 a week                         1.650
## smoking_age1640 to 49 a week                         0.950
## smoking_age1650 to 59 a week                         2.090
## smoking_age1660 or more a wk                         2.028
## alcohol_age16Less than 1 week                        4.556
## alcohol_age162 to 4 weeks                            3.836
## alcohol_age165 to 8 weeks                            3.160
## alcohol_age169 to 12 weeks                           2.489
## alcohol_age16Over 12 weeks                           2.149
## alcohol_age16Do not remember                         1.758
## alcohol_age16Never had one                          -1.348
## school_attendance_age16                              5.137
## age:read_classModerate initial / slowly increasing   6.512

tab_model(mixed_cog)

	z global cog
Predictors	Estimates	CI	p
(Intercept)	-0.28	-1.22 – 0.65	0.550
age	-0.03	-0.03 – -0.02	<0.001
read class [Moderate initial / slowly increasing]	0.23	-0.04 – 0.49	0.097
sex [Female]	-0.01	-0.05 – 0.02	0.406
sc prenatal [I]	0.34	0.19 – 0.49	<0.001
sc prenatal [II]	0.27	0.15 – 0.39	<0.001
sc prenatal [III]	0.13	0.02 – 0.24	0.021
sc prenatal [IV]	0.07	-0.05 – 0.19	0.263
sc prenatal [V]	-0.04	-0.16 – 0.09	0.557
sc prenatal [Students]	0.45	-0.07 – 0.96	0.089
sc prenatal [Retired]	0.41	-0.65 – 1.48	0.448
sc prenatal [Single,no husbnd]	-0.02	-0.18 – 0.13	0.757
smoking preg [5mths-no change]	-0.02	-0.18 – 0.14	0.800
smoking preg [Now non-smoker]	0.01	-0.16 – 0.19	0.867
smoking preg [1 to 4 now]	-0.07	-0.25 – 0.12	0.475
smoking preg [5 to 9 now]	-0.07	-0.26 – 0.12	0.468
smoking preg [10 to 14 now]	-0.05	-0.26 – 0.17	0.685
smoking preg [15 to 19 now]	-0.19	-0.46 – 0.08	0.170
smoking preg [20 to 24 now]	-0.38	-0.70 – -0.06	0.019
smoking preg [25 to 29 now]	0.35	-0.17 – 0.88	0.185
smoking preg [30 or more now]	-0.43	-0.99 – 0.12	0.124
smoking_pregVariable-5mths	-0.08	-0.26 – 0.09	0.358
birth weight	0.00	0.00 – 0.00	<0.001
parent read level [Occasionally]	-0.42	-0.68 – -0.16	0.002
parent read level [Weekly]	-0.28	-0.54 – -0.02	0.037
sc age7 [No male head]	-0.04	-0.27 – 0.19	0.725
sc age7 [I]	0.29	0.07 – 0.52	0.009
sc age7 [II]	0.19	-0.02 – 0.39	0.070
sc age7 [III non-manual]	0.20	-0.00 – 0.41	0.052
sc age7 [III manual]	-0.01	-0.20 – 0.19	0.941
sc age7 [IV non-manual]	0.03	-0.20 – 0.27	0.770
sc age7 [IV manual]	-0.05	-0.25 – 0.15	0.596
sc age7 [V]	-0.17	-0.38 – 0.03	0.102
breast f [Dont know]	-0.78	-1.68 – 0.11	0.087
breast f [No]	-0.69	-1.56 – 0.17	0.116
breast f [Under one month]	-0.63	-1.49 – 0.24	0.157
breast f [Over one month]	-0.53	-1.40 – 0.34	0.231
sc age11 [I]	0.27	0.15 – 0.39	<0.001
sc age11 [II]	0.21	0.13 – 0.29	<0.001
sc age11 [III non manual]	0.19	0.10 – 0.28	<0.001
sc age11 [III manual]	0.06	-0.01 – 0.13	0.074
sc age11 [IV]	-0.00	-0.08 – 0.07	0.955
sc age11 [V]	-0.02	-0.13 – 0.08	0.651
sc age11 [No male head]	0.02	-0.09 – 0.13	0.705
sc age16 [I]	0.18	0.07 – 0.30	0.001
sc age16 [II]	0.10	0.04 – 0.16	0.002
sc age16 [III non-manual]	0.02	-0.06 – 0.10	0.590
sc age16 [III manual]	0.04	-0.00 – 0.09	0.076
sc age16 [IV non-manual]	0.05	-0.11 – 0.22	0.524
sc age16 [IV manual]	0.01	-0.06 – 0.08	0.757
sc age16 [V]	-0.10	-0.20 – -0.00	0.048
sc age16 [Unclear]	0.20	0.03 – 0.38	0.025
smoking age16 [None,dont smoke]	0.40	0.22 – 0.57	<0.001
smoking age16 [Less 1 a week]	0.30	0.10 – 0.50	0.003
smoking age16 [1 to 9 a week]	0.19	0.00 – 0.38	0.047
smoking age16 [10 to 19 a week]	0.13	-0.06 – 0.32	0.194
smoking age16 [20 to 29 a week]	0.11	-0.08 – 0.30	0.269
smoking age16 [30 to 39 a week]	0.16	-0.03 – 0.35	0.099
smoking age16 [40 to 49 a week]	0.09	-0.10 – 0.29	0.342
smoking age16 [50 to 59 a week]	0.21	0.01 – 0.41	0.037
smoking age16 [60 or more a wk]	0.19	0.01 – 0.38	0.043
alcohol age16 [Less than 1 week]	0.46	0.26 – 0.66	<0.001
alcohol age16 [2 to 4 weeks]	0.39	0.19 – 0.59	<0.001
alcohol age16 [5 to 8 weeks]	0.33	0.13 – 0.54	0.002
alcohol age16 [9 to 12 weeks]	0.27	0.06 – 0.49	0.013
alcohol age16 [Over 12 weeks]	0.23	0.02 – 0.43	0.032
alcohol age16 [Do not remember]	0.18	-0.02 – 0.38	0.079
alcohol age16 [Never had one]	-0.14	-0.35 – 0.07	0.178
school attendance age16	0.00	0.00 – 0.00	<0.001
age × read class [Moderate initial / slowly increasing]	0.02	0.02 – 0.03	<0.001
Random Effects
σ2	0.31
τ00 ncdsid	0.47
ICC	0.60
N ncdsid	8473
Observations	25419
Marginal R2 / Conditional R2	0.188 / 0.676

pred_cog <- ggpredict(
  mixed_cog,
  terms = c("age [7:16]", "read_class")
)

fig2 <- ggplot(
  pred_cog,
  aes(
    x = x,
    y = predicted,
    colour = group,
    fill = group
  )
)  +
  geom_line(size = 1.2) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
  scale_color_jama() +
  scale_fill_jama() +
  scale_x_continuous(breaks = c(7, 11, 16)) +
  labs(
    x = "Age (years)",
    y = "Global cognitive function (z-score)",
    colour = "Reading trajectory",
    fill   = "Reading trajectory"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank()
  )

fig2

SB ACROSS ADULTHOOD

table(df_complete_7_16$tv_age11)

## 
## Hardly ever   Sometimes   most days 
##         208        1043        7536

table(df_complete_7_16$tv_age16)

## 
## Hardly ever   Sometimes       Often 
##         476        2546        5736

table(df_complete_7_16$tv_age23)

## 
## 5 TIMES A WEEK +  3-4TIMES A WEEK  1-2TIMES A WEEK 2-3TIMES LST4WKS 
##             4976             1288              728              236 
## ONCE IN LST 4WKS NOT DONE LST4WKS 
##               59              106

table(df_complete_7_16$read_book_age23)

## 
## 5 TIMES A WEEK +  3-4TIMES A WEEK  1-2TIMES A WEEK 2-3TIMES LST4WKS 
##             1564              816              941              579 
## ONCE IN LST 4WKS NOT DONE LST4WKS 
##              706             2780

table(df_complete_7_16$read_book_age46)

## 
##                           Refused                        Don't know 
##                                 0                                 0 
##                     Other missing                    Not applicable 
##                                 1                                 0 
##                         Every day Most days (4 or 5 days each week) 
##                              1407                               757 
##              At least once a week             At least once a month 
##                               772                               684 
##              Less often than that                             Never 
##                              1412                               701

table(df_complete_7_16$read_news_age46)

## 
##                           Refused                        Don't know 
##                                 0                                 0 
##                     Other missing                    Not applicable 
##                                 1                                 0 
##                         Every day Most days (4 or 5 days each week) 
##                              2713                              1030 
##              At least once a week             At least once a month 
##                              1431                               280 
##              Less often than that                             Never 
##                               146                               133

table(df_complete_7_16$computer_home_age46)

## 
##        Refused     Don't know  Other missing Not applicable             No 
##              0              0              1              0            763 
##            yes 
##           4970

table(df_complete_7_16$computer_home_freq_age46)

## 
##               Refused            Don't know         Other missing 
##                     0                     2                     0 
##        Not applicable                 Daily      2-4 times a week 
##                   764                  1731                  1243 
##           Once a week Less than once a week                 Never 
##                   603                   750                   641

table(df_complete_7_16$computer_work_age46)

## 
##        Refused     Don't know  Other missing Not applicable             No 
##              0              0              1            690           1417 
##            yes 
##           3626

table(df_complete_7_16$computer_work_freq_age46)

## 
##               Refused            Don't know         Other missing 
##                     0                     2                     0 
##        Not applicable                 Daily      2-4 times a week 
##                  2108                  3165                   246 
##           Once a week Less than once a week 
##                   123                    90

table(df_complete_7_16$computer_home_age50)

## 
##                 Refusal              Don't Know Schedule not applicable 
##                       0                       0                       0 
##     Item not applicable                     Yes                      No 
##                       8                    5288                     517

table(df_complete_7_16$computer_home_freq_age46)

## 
##               Refused            Don't know         Other missing 
##                     0                     2                     0 
##        Not applicable                 Daily      2-4 times a week 
##                   764                  1731                  1243 
##           Once a week Less than once a week                 Never 
##                   603                   750                   641

table(df_complete_7_16$computer_work_age50)

## 
##                 Refusal              Don't Know Schedule not applicable 
##                       0                       0                       0 
##     Item not applicable                     Yes                      No 
##                     818                    3751                    1244

table(df_complete_7_16$computer_work_freq_age46)

## 
##               Refused            Don't know         Other missing 
##                     0                     2                     0 
##        Not applicable                 Daily      2-4 times a week 
##                  2108                  3165                   246 
##           Once a week Less than once a week 
##                   123                    90

table(df_complete_7_16$internet_use_age55)

## 
##                     Refused                  Don't know 
##                           3                           2 
##              Not applicable Everyday or almost everyday 
##                          68                        3962 
##        Several times a week        Once or twice a week 
##                         478                         373 
##       At least once a month         Less often or never 
##                         164                         453

table(df_complete_7_16$read_book_age62)

## 
##                  Not answered                   Multi-coded 
##                            10                             2 
##                Not applicable Every day or almost every day 
##                           745                          1469 
##          Several times a week          Once or twice a week 
##                           525                           392 
##         At least once a month              Every few months 
##                           315                           536 
##          At least once a year           Less often or never 
##                           314                           722

table(df_complete_7_16$tv_age62)

## 
##            Not answered             Multi-coded          Not applicable 
##                      12                       2                     745 
##                    None  Less than 1 hour a day      1 to 2 hours a day 
##                      76                     289                     781 
##      2 to 3 hours a day      3 to 4 hours a day More than 4 hours a day 
##                    1071                     967                    1087

df_complete_7_16$read3_age11 <- factor(
  df_complete_7_16$child_read_books_age11,
  levels = c("Hardly ever", "Sometimes", "most days"),
  labels = c("Low", "Moderate", "High")
)

df_complete_7_16$read3_age16 <- factor(
  df_complete_7_16$child_read_age16,
  levels = c("Hardly ever", "Sometimes", "Most days"),
  labels = c("Low", "Moderate", "High")
)

df_complete_7_16$tv3_age11 <- factor(
  df_complete_7_16$tv_age11,
  levels = c("Hardly ever", "Sometimes", "most days"),
  labels = c("Low", "Moderate", "High")
)

df_complete_7_16$tv3_age16 <- factor(
  df_complete_7_16$tv_age16,
  levels = c("Hardly ever", "Sometimes", "Often"),
  labels = c("Low", "Moderate", "High")
)

df_complete_7_16$tv3_age23 <- case_when(
  df_complete_7_16$tv_age23 %in% c("NOT DONE LST4WKS", "ONCE IN LST 4WKS") ~ "Low",
  df_complete_7_16$tv_age23 %in% c("2-3TIMES LST4WKS", "1-2TIMES A WEEK") ~ "Moderate",
  df_complete_7_16$tv_age23 %in% c("3-4TIMES A WEEK", "5 TIMES A WEEK +") ~ "High",
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "Moderate", "High"))

df_complete_7_16$tv3_age62 <- case_when(
  df_complete_7_16$tv_age62 %in% c("None", "Less than 1 hour a day") ~ "Low",
  df_complete_7_16$tv_age62 %in% c("1 to 2 hours a day", "2 to 3 hours a day") ~ "Moderate",
  df_complete_7_16$tv_age62 %in% c("3 to 4 hours a day", "More than 4 hours a day") ~ "High",
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "Moderate", "High"))


df_complete_7_16$read3_age23 <- case_when(
  df_complete_7_16$read_book_age23 %in% c("NOT DONE LST4WKS") ~ "Low",
  df_complete_7_16$read_book_age23 %in% c("ONCE IN LST 4WKS", 
                                          "1-2TIMES A WEEK") ~ "Moderate",
  df_complete_7_16$read_book_age23 %in% c("2-3TIMES LST4WKS",
                                          "3-4TIMES A WEEK", 
                                          "5 TIMES A WEEK +") ~ "High",
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "Moderate", "High"))


df_complete_7_16$read3_age46 <- case_when(
  df_complete_7_16$read_book_age46 %in% c("Never", "Less often than that") &
    df_complete_7_16$read_news_age46 %in% c("Never", "Less often than that") ~ "Low",
  
  df_complete_7_16$read_book_age46 %in% c("At least once a week", "At least once a month") |
    df_complete_7_16$read_news_age46 %in% c("At least once a week", "At least once a month") ~ "Moderate",
  
  df_complete_7_16$read_book_age46 %in% c("Every day", "Most days (4 or 5 days each week)") |
    df_complete_7_16$read_news_age46 %in% c("Every day", "Most days (4 or 5 days each week)") ~ "High",
  
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "Moderate", "High"))


df_complete_7_16$read3_age62 <- case_when(
  df_complete_7_16$read_book_age62 %in% c("Less often or never", "Every few months" ,"At least once a year") ~ "Low",
  df_complete_7_16$read_book_age62 %in% c("At least once a month", "At least once a year") ~ "Moderate",
  df_complete_7_16$read_book_age62 %in% c("Several times a week", "Once or twice a week", "Every day or almost every day") ~ "High",
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "Moderate", "High"))


df_complete_7_16$computer3_home_age46 <- case_when(
  df_complete_7_16$computer_home_age46 == "No" ~ "Low",
  
  df_complete_7_16$computer_home_freq_age46 %in%
    c("Never", "Less than once a week") ~ "Low",
  
  df_complete_7_16$computer_home_freq_age46 %in%
    c("Once a week", "2-4 times a week") ~ "Moderate",
  
  df_complete_7_16$computer_home_freq_age46 %in%
    c("Daily") ~ "High",
  
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "Moderate", "High"))


df_complete_7_16$computer3_work_age46 <- case_when(
  df_complete_7_16$computer_work_age46 == "No" ~ "Low",
  
  df_complete_7_16$computer_work_freq_age46 %in%
    c("Less than once a week") ~ "Low",
  
  df_complete_7_16$computer_work_freq_age46 %in%
    c("Once a week") ~ "Moderate",
  
  df_complete_7_16$computer_work_freq_age46 %in%
    c("2-4 times a week", "Daily") ~ "High",
  
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "Moderate", "High"))


df_complete_7_16$computer3_age46 <- pmax(
  as.numeric(df_complete_7_16$computer3_home_age46),
  as.numeric(df_complete_7_16$computer3_work_age46),
  na.rm = TRUE
)

df_complete_7_16$computer3_age46 <- factor(
  df_complete_7_16$computer3_age46,
  levels = 1:3,
  labels = c("Low", "Moderate", "High")
)


df_complete_7_16$computer3_age50 <- case_when(
  df_complete_7_16$computer_home_age50 == "No" ~ "Low",
  df_complete_7_16$computer_home_age50 == "Yes" ~ "High",
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "High"))


df_complete_7_16$internet3_age55 <- case_when(
  df_complete_7_16$internet_use_age55 %in%
    c("Less often or never") ~ "Low",
  
  df_complete_7_16$internet_use_age55 %in%
    c("Once or twice a week", "Several times a week") ~ "Moderate",
  
  df_complete_7_16$internet_use_age55 %in%
    c("Everyday or almost everyday") ~ "High",
  
  TRUE ~ NA_character_
) |> factor(levels = c("Low", "Moderate", "High"))

GROUP-BASED TRAJECTORY MODELLING - ADULTHOOD - TV

library(tidyr)
library(lcmm)
library(dplyr)

df_tv <- df_complete_7_16 %>%
  dplyr::select(
    ncdsid,
    tv3_age11,
    tv3_age16,
    tv3_age23,
    tv3_age62
  ) %>%
  mutate(
    tv11 = as.numeric(tv3_age11) - 1,
    tv16 = as.numeric(tv3_age16) - 1,
    tv23 = as.numeric(tv3_age23) - 1,
    tv62 = as.numeric(tv3_age62) - 1
  )

df_tv_long <- df_tv %>%
  pivot_longer(
    cols = c(tv11, tv16, tv23, tv62),
    names_to  = "wave",
    values_to = "tv_num"
  ) %>%
  mutate(
    age = case_when(
      wave == "tv11" ~ 11,
      wave == "tv16" ~ 16,
      wave == "tv23" ~ 23,
      wave == "tv62" ~ 62
    )
  ) %>%
  filter(!is.na(tv_num))

df_tv_long <- df_tv_long %>%
  mutate(
    id_num = as.numeric(as.factor(ncdsid))
  )

m1 <- hlme(
  fixed   = tv_num ~ age,
  random  = ~ 1,
  subject = "id_num",
  ng      = 1,
  data    = df_tv_long
)

summary(m1)

## Heterogenous linear mixed model 
##      fitted by maximum likelihood method 
##  
## hlme(fixed = tv_num ~ age, random = ~1, subject = "id_num", ng = 1, 
##     data = df_tv_long)
##  
## Statistical Model: 
##      Dataset: df_tv_long 
##      Number of subjects: 9079 
##      Number of observations: 29209 
##      Number of latent classes: 1 
##      Number of parameters: 4  
##  
## Iteration process: 
##      Convergence criteria satisfied 
##      Number of iterations:  11 
##      Convergence criteria: parameters= 2.4e-13 
##                          : likelihood= 1.1e-10 
##                          : second derivatives= 4.6e-14 
##  
## Goodness-of-fit statistics: 
##      maximum log-likelihood: -22691.85  
##      AIC: 45391.7  
##      BIC: 45420.16  
##  
##  
## Maximum Likelihood Estimates: 
##  
## Fixed effects in the longitudinal model:
## 
##               coef      Se     Wald p-value
## intercept  1.86267 0.00531  350.680 0.00000
## age       -0.00720 0.00018  -40.286 0.00000
## 
## 
## Variance-covariance matrix of the random-effects:
##           intercept
## intercept   0.02617
## 
##                              coef      Se
## Residual standard error:  0.50339 0.00248

m2 <- hlme(
  fixed   = tv_num ~ age,
  mixture = ~ age,
  random  = ~ 1,
  subject = "id_num",
  ng      = 2,
  nwg     = TRUE,
  data    = df_tv_long,
  B       = random(m1)
)

summary(m2)

## Heterogenous linear mixed model 
##      fitted by maximum likelihood method 
##  
## hlme(fixed = tv_num ~ age, mixture = ~age, random = ~1, subject = "id_num", 
##     ng = 2, nwg = TRUE, data = df_tv_long)
##  
## Statistical Model: 
##      Dataset: df_tv_long 
##      Number of subjects: 9079 
##      Number of observations: 29209 
##      Number of latent classes: 2 
##      Number of parameters: 8  
##  
## Iteration process: 
##      Maximum number of iteration reached without convergence 
##      Number of iterations:  500 
##      Convergence criteria: parameters= 0.00021 
##                          : likelihood= 8.1e-06 
##                          : second derivatives= 1 
##  
## Goodness-of-fit statistics: 
##      maximum log-likelihood: -22019.76  
##      AIC: 44055.51  
##      BIC: 44112.42  
##  
##  
## Maximum Likelihood Estimates: 
##  
## Fixed effects in the class-membership model:
## (the class of reference is the last class) 
## 
##                       coef Se Wald p-value
## intercept class1  -1.16542                
## 
## Fixed effects in the longitudinal model:
## 
##                       coef Se Wald p-value
## intercept class1   1.78087                
## intercept class2   1.88621                
## age class1        -0.01716                
## age class2        -0.00402                
## 
## 
## Variance-covariance matrix of the random-effects:
##           intercept
## intercept         0
## 
##                                      coef Se
## Proportional coefficient class1 141.44814   
## Residual standard error:          0.47405

age_grid <- data.frame(
  age = c(11, 16, 23, 62)
)

pred_tv <- predictY(
  m2,
  newdata = age_grid,
  var.time = "age",
  draws = TRUE
)

str(pred_tv)

## List of 2
##  $ pred : num [1:4, 1:6] 1.592 1.506 1.386 0.717 1.842 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:6] "Ypred_class1" "Ypred_class2" "lower.Ypred_class1" "lower.Ypred_class2" ...
##  $ times:'data.frame':   4 obs. of  1 variable:
##   ..$ age: num [1:4] 11 16 23 62
##  - attr(*, "class")= chr "predictY"

pred_df <- data.frame(
  age = pred_tv$times$age,
  pred_tv$pred
)

traj_df <- pred_df %>%
  pivot_longer(
    cols = -age,
    names_to = "term",
    values_to = "value"
  ) %>%
  separate(term, into = c("type", "class"), sep = "\\.") %>%
  pivot_wider(
    names_from  = type,
    values_from = value
  ) %>%
  mutate(
  class = dplyr::recode(
    class,
    "class1" = "Low / slowly increasing TV",
    "class2" = "High / rapidly increasing TV"
  )
)

# Grid de idades observado
age_grid <- data.frame(
  age = seq(
    min(df_tv_long$age, na.rm = TRUE),
    max(df_tv_long$age, na.rm = TRUE),
    by = 0.5
  )
)

# Predições com IC95%
pred_tv <- predictY(
  m2,
  newdata = age_grid,
  var.time = "age",
  draws = TRUE
)

traj_df <- data.frame(
  age = pred_tv$times$age,

  mean_class1  = pred_tv$pred[, "Ypred_class1"],
  lower_class1 = pred_tv$pred[, "lower.Ypred_class1"],
  upper_class1 = pred_tv$pred[, "upper.Ypred_class1"],

  mean_class2  = pred_tv$pred[, "Ypred_class2"],
  lower_class2 = pred_tv$pred[, "lower.Ypred_class2"],
  upper_class2 = pred_tv$pred[, "upper.Ypred_class2"]
) %>%
  tidyr::pivot_longer(
    cols = -age,
    names_to = c("stat", "class"),
    names_sep = "_class",
    values_to = "value"
  ) %>%
  tidyr::pivot_wider(
    names_from  = stat,
    values_from = value
  ) %>%
  mutate(
    class = factor(
      class,
      levels = c("1", "2"),
      labels = c(
        "High / slowly declining",
        "High / rapidly declining"
      )
    )
  )

fig_tv <- ggplot(traj_df, aes(x = age, y = mean, color = class, fill = class)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.25,
    color = NA
  ) +
  geom_line(size = 1.2) +
  scale_color_bmj() +
  scale_fill_bmj() +
  labs(
    x = "Age (years)",
    y = "TV watching (ordinal score)",
    color = "Trajectory TV watching",
    fill  = "Trajectory TV watching"
  ) +
  theme_classic(base_size = 14) +
  theme(
    legend.position = "top",
    legend.title = element_text(size = 12),
    legend.text  = element_text(size = 11),
    axis.title   = element_text(size = 13),
    axis.text    = element_text(size = 11)
  ) +

  scale_x_continuous(
    breaks = c(11, 16, 23, 62)
  ) +
  scale_y_continuous(
    breaks = c(0, 1, 2),
    labels = c("Low", "Mid", "High")
  )  +
  coord_cartesian(ylim = c(0,2))
fig_tv

COGNITIVE FUNCION - ADULTHOOD

df_complete_7_62 <- df_complete_7_16 %>%
  filter(words_immed_age50>-1 &
            words_delayed>-1 & 
            letter_canc_correct>-1 &
            letter_canc_accuracy>-1 &
            animals_named>-1 & 
            
            words_delayed_age62>-1 &
            letter_can_accuracy>-1 &
            animals_named_age62>-1)

df_complete_7_50 <- df_complete_7_16 %>%
  filter(words_immed_age50>-1 &
            words_delayed>-1 & 
            letter_canc_correct>-1 &
            letter_canc_accuracy>-1 &
            animals_named>-1)

            #### #####

# Age 50
#Word recall test - immediate recall
df_complete_7_50$mean_words_immed_age50 <- mean(df_complete_7_50$words_immed_age50)
df_complete_7_50$sd_words_immed_age50 <- sd(df_complete_7_50$words_immed_age50)
df_complete_7_50$z_words_immed_age50 <- (df_complete_7_50$words_immed_age50-df_complete_7_50$mean_words_immed_age50)/df_complete_7_50$sd_words_immed_age50
df_complete_7_50$z_words_immed_age50 = scale(df_complete_7_50$z_words_immed_age50)[,1]

hist(df_complete_7_50$z_words_immed_age50)

#Word recall test - Delayed recall
df_complete_7_50$mean_words_delayed <- mean(df_complete_7_50$words_delayed)
df_complete_7_50$sd_words_delayed <- sd(df_complete_7_50$words_delayed)
df_complete_7_50$z_words_delayed <- (df_complete_7_50$words_delayed-df_complete_7_50$mean_words_delayed)/df_complete_7_50$sd_words_delayed
df_complete_7_50$z_words_delayed = scale(df_complete_7_50$z_words_delayed)[,1]

hist(df_complete_7_50$z_words_delayed)

#Animal naming - Number of animals named
df_complete_7_50$mean_animals_named <- mean(df_complete_7_50$animals_named)
df_complete_7_50$sd_animals_named <- sd(df_complete_7_50$animals_named)
df_complete_7_50$z_animals_named <- (df_complete_7_50$animals_named-df_complete_7_50$mean_animals_named)/df_complete_7_50$sd_animals_named
df_complete_7_50$z_animals_named = scale(df_complete_7_50$z_animals_named)[,1]

hist(df_complete_7_50$z_animals_named)

#Letter cancelation - accuracy
df_complete_7_50$mean_letter_canc_accuracy <- mean(df_complete_7_50$letter_canc_accuracy)
df_complete_7_50$sd_letter_canc_accuracy <- sd(df_complete_7_50$letter_canc_accuracy)
df_complete_7_50$z_letter_canc_accuracy <- (df_complete_7_50$letter_canc_accuracy-df_complete_7_50$mean_letter_canc_accuracy)/df_complete_7_50$sd_letter_canc_accuracy
df_complete_7_50$z_letter_canc_accuracy = scale(df_complete_7_50$z_letter_canc_accuracy)[,1]

df_complete_7_50$z_letter_canc_accuracy <- df_complete_7_50$z_letter_canc_accuracy*-1
hist(df_complete_7_50$z_letter_canc_accuracy)

# averaging all z scores
df_complete_7_50$global_cog_age50 <- (df_complete_7_50$z_words_immed_age50 + df_complete_7_50$z_words_delayed  + df_complete_7_50$z_animals_named  + df_complete_7_50$z_letter_canc_accuracy) / 4

# creating a z score for the averaged variable
df_complete_7_50$global_cog_age50 = scale(df_complete_7_50$global_cog_age50)[,1]

hist(df_complete_7_50$global_cog_age50)

# Age 50
#Word recall test - immediate recall
df_complete_7_62$mean_words_immed_age50 <- mean(df_complete_7_62$words_immed_age50)
df_complete_7_62$sd_words_immed_age50 <- sd(df_complete_7_62$words_immed_age50)
df_complete_7_62$z_words_immed_age50 <- (df_complete_7_62$words_immed_age50-df_complete_7_62$mean_words_immed_age50)/df_complete_7_62$sd_words_immed_age50
df_complete_7_62$z_words_immed_age50 = scale(df_complete_7_62$z_words_immed_age50)[,1]

hist(df_complete_7_62$z_words_immed_age50)

#Word recall test - Delayed recall
df_complete_7_62$mean_words_delayed <- mean(df_complete_7_62$words_delayed)
df_complete_7_62$sd_words_delayed <- sd(df_complete_7_62$words_delayed)
df_complete_7_62$z_words_delayed <- (df_complete_7_62$words_delayed-df_complete_7_62$mean_words_delayed)/df_complete_7_62$sd_words_delayed
df_complete_7_62$z_words_delayed = scale(df_complete_7_62$z_words_delayed)[,1]

hist(df_complete_7_62$z_words_delayed)

#Animal naming - Number of animals named
df_complete_7_62$mean_animals_named <- mean(df_complete_7_62$animals_named)
df_complete_7_62$sd_animals_named <- sd(df_complete_7_62$animals_named)
df_complete_7_62$z_animals_named <- (df_complete_7_62$animals_named-df_complete_7_62$mean_animals_named)/df_complete_7_62$sd_animals_named
df_complete_7_62$z_animals_named = scale(df_complete_7_62$z_animals_named)[,1]

hist(df_complete_7_62$z_animals_named)

#Letter cancelation - accuracy
df_complete_7_62$mean_letter_canc_accuracy <- mean(df_complete_7_62$letter_canc_accuracy)
df_complete_7_62$sd_letter_canc_accuracy <- sd(df_complete_7_62$letter_canc_accuracy)
df_complete_7_62$z_letter_canc_accuracy <- (df_complete_7_62$letter_canc_accuracy-df_complete_7_62$mean_letter_canc_accuracy)/df_complete_7_62$sd_letter_canc_accuracy
df_complete_7_62$z_letter_canc_accuracy = scale(df_complete_7_62$z_letter_canc_accuracy)[,1]

df_complete_7_62$z_letter_canc_accuracy <- df_complete_7_62$z_letter_canc_accuracy*-1
hist(df_complete_7_62$z_letter_canc_accuracy)

# averaging all z scores
df_complete_7_62$global_cog_age50 <- (df_complete_7_62$z_words_immed_age50 + df_complete_7_62$z_words_delayed  + df_complete_7_62$z_animals_named  + df_complete_7_62$z_letter_canc_accuracy) / 4

# creating a z score for the averaged variable
df_complete_7_62$global_cog_age50 = scale(df_complete_7_62$global_cog_age50)[,1]

hist(df_complete_7_62$global_cog_age50)

# Age 62
#Word recall test - Delayed recall
df_complete_7_62$mean_words_delayed_age62 <- mean(df_complete_7_62$words_delayed_age62)
df_complete_7_62$sd_words_delayed_age62 <- sd(df_complete_7_62$words_delayed_age62)
df_complete_7_62$z_words_delayed_age62 <- (df_complete_7_62$words_delayed_age62-df_complete_7_62$mean_words_delayed_age62)/df_complete_7_62$sd_words_delayed_age62
df_complete_7_62$z_words_delayed_age62 = scale(df_complete_7_62$z_words_delayed_age62)[,1]

hist(df_complete_7_62$z_words_delayed_age62)

#Animal naming - Number of animals named
df_complete_7_62$mean_animals_named_age62 <- mean(df_complete_7_62$animals_named_age62)
df_complete_7_62$sd_animals_named_age62 <- sd(df_complete_7_62$animals_named_age62)
df_complete_7_62$z_animals_named_age62 <- (df_complete_7_62$animals_named_age62-df_complete_7_62$mean_animals_named_age62)/df_complete_7_62$sd_animals_named_age62
df_complete_7_62$z_animals_named_age62 = scale(df_complete_7_62$z_animals_named_age62)[,1]

hist(df_complete_7_62$z_animals_named_age62)

#Letter cancelation - accuracy
df_complete_7_62$mean_letter_can_accuracy <- mean(df_complete_7_62$letter_can_accuracy)
df_complete_7_62$sd_letter_can_accuracy <- sd(df_complete_7_62$letter_can_accuracy)
df_complete_7_62$z_letter_can_accuracy <- (df_complete_7_62$letter_canc_accuracy-df_complete_7_62$mean_letter_can_accuracy)/df_complete_7_62$sd_letter_can_accuracy
df_complete_7_62$z_letter_can_accuracy = scale(df_complete_7_62$z_letter_can_accuracy)[,1]

df_complete_7_62$z_letter_can_accuracy <- df_complete_7_62$z_letter_can_accuracy*-1
hist(df_complete_7_62$z_letter_can_accuracy)

# averaging all z scores
df_complete_7_62$global_cog_age62 <- (df_complete_7_62$z_words_delayed_age62 + df_complete_7_62$z_animals_named_age62  + df_complete_7_62$z_letter_can_accuracy) / 3

# creating a z score for the averaged variable
df_complete_7_62$global_cog_age62 = scale(df_complete_7_62$global_cog_age62)[,1]

hist(df_complete_7_62$global_cog_age62)

COGNITIVE DECLINE - ADULTHOOD

library(dplyr)
library(tidyr)
library(ggeffects)
library(lme4)
library(broom.mixed)

tv_class_df <- m2$pprob %>%
  mutate(
    tv_class = factor(
      class,
      levels = c(1, 2),
      labels = c(
        "High / slowly declining TV",
        "High / rapidly declining TV"
      )
    )
  ) %>%
  dplyr::select(id_num, tv_class)

df_cog_long <- df_complete_7_62 %>%
  dplyr::select(
    ncdsid,
    z_global_cog_age7,
    z_global_cog_age11,
    z_global_cog_age16,
    global_cog_age50,
    global_cog_age62
  ) %>%
  pivot_longer(
    cols = -ncdsid,
    names_to  = "wave",
    values_to = "z_cog"
  ) %>%
  mutate(
    age = case_when(
      wave == "z_global_cog_age7"   ~ 7,
      wave == "z_global_cog_age11"  ~ 11,
      wave == "z_global_cog_age16"  ~ 16,
      wave == "global_cog_age50"    ~ 50,
      wave == "global_cog_age62"    ~ 62
    )
  ) %>%
  filter(!is.na(z_cog)) %>%
  mutate(
    id_num = as.numeric(as.factor(ncdsid))
  )

df_cov <- df_complete_7_62 %>%
  dplyr::select(
    ncdsid,

    # Demográficas
    sex,

    # Gestação / início da vida
    sc_prenatal,
    smoking_preg,
    birth_weight,
    breast_f,

    # Ambiente cognitivo precoce
    parent_read_level,

    # SEP ao longo da infância/adolescência
    sc_age7,
    sc_age11,
    sc_age16,

    # Comportamentos na adolescência
    smoking_age16,
    alcohol_age16,

    # Escolaridade
    school_attendance_age16
  ) %>%
  mutate(
    id_num = as.numeric(as.factor(ncdsid))
  )

df_tv_class <- m2$pprob %>%
  mutate(
    tv_class = factor(
      class,
      labels = c(
        "High / slowly decline TV",
        "High / rapidly decline TV"
      )
    )
  ) %>%
  dplyr::select(id_num, tv_class)

df_cog_tv <- df_cog_long %>%
  left_join(df_tv_class, by = "id_num") %>%
  left_join(df_cov,       by = "id_num")

cog_mixed_adj <- lmer(
  z_cog ~ age * tv_class +

    sex +
    sc_prenatal +
    smoking_preg +
    birth_weight +
    breast_f +
    parent_read_level +

    sc_age7 +
    sc_age11 +
    sc_age16 +

    smoking_age16 +
    alcohol_age16 +
    school_attendance_age16 +

    (1 | id_num),

  data = df_cog_tv,
  REML = TRUE
)

summary(cog_mixed_adj)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## z_cog ~ age * tv_class + sex + sc_prenatal + smoking_preg + birth_weight +  
##     breast_f + parent_read_level + sc_age7 + sc_age11 + sc_age16 +  
##     smoking_age16 + alcohol_age16 + school_attendance_age16 +      (1 | id_num)
##    Data: df_cog_tv
## 
## REML criterion at convergence: 40425.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1433 -0.5896  0.0236  0.6150  4.2980 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id_num   (Intercept) 0.3282   0.5729  
##  Residual             0.4906   0.7004  
## Number of obs: 16615, groups:  id_num, 3323
## 
## Fixed effects:
##                                         Estimate Std. Error t value
## (Intercept)                           -0.3912235  0.5383939  -0.727
## age                                   -0.0066663  0.0006609 -10.086
## tv_classHigh / rapidly decline TV     -0.0576405  0.0394470  -1.461
## sexFemale                              0.0219067  0.0233323   0.939
## sc_prenatalI                           0.3160972  0.0991735   3.187
## sc_prenatalII                          0.2517210  0.0856966   2.937
## sc_prenatalIII                         0.1238566  0.0796812   1.554
## sc_prenatalIV                          0.0769719  0.0855353   0.900
## sc_prenatalV                           0.0452001  0.0919379   0.492
## sc_prenatalStudents                    0.2099551  0.3408081   0.616
## sc_prenatalSingle,no husbnd            0.0433839  0.1146246   0.378
## smoking_preg5mths-no change           -0.0794388  0.1249160  -0.636
## smoking_pregNow non-smoker            -0.0306702  0.1312993  -0.234
## smoking_preg1 to 4 now                -0.1781563  0.1400026  -1.273
## smoking_preg5 to 9 now                -0.0568105  0.1457027  -0.390
## smoking_preg10 to 14 now              -0.1306281  0.1639510  -0.797
## smoking_preg15 to 19 now              -0.0939355  0.2163372  -0.434
## smoking_preg20 to 24 now              -0.2864433  0.2271215  -1.261
## smoking_preg25 to 29 now              -0.3097697  0.3993037  -0.776
## smoking_preg30 or more now            -0.1555445  0.4021066  -0.387
## smoking_pregVariable-5mths            -0.1147373  0.1332362  -0.861
## birth_weight                           0.0002363  0.0001840   1.284
## breast_fDont know                     -0.0043236  0.5001469  -0.009
## breast_fNo                            -0.1082451  0.4666688  -0.232
## breast_fUnder one month               -0.0570621  0.4668196  -0.122
## breast_fOver one month                 0.0314924  0.4666089   0.067
## parent_read_levelOccasionally          0.0242604  0.0563269   0.431
## parent_read_levelWeekly                0.1618518  0.0526955   3.071
## sc_age7No male head                    0.2602114  0.1878907   1.385
## sc_age7I                               0.1946954  0.1738999   1.120
## sc_age7II                              0.1838100  0.1657158   1.109
## sc_age7III non-manual                  0.2020526  0.1673749   1.207
## sc_age7III manual                      0.0239362  0.1625361   0.147
## sc_age7IV non-manual                   0.2785051  0.1874317   1.486
## sc_age7IV manual                       0.0290632  0.1654240   0.176
## sc_age7V                              -0.0279753  0.1728439  -0.162
## sc_age11I                              0.1496386  0.0785989   1.904
## sc_age11II                             0.0231250  0.0583108   0.397
## sc_age11III non manual                 0.0754331  0.0652675   1.156
## sc_age11III manual                    -0.0076475  0.0520950  -0.147
## sc_age11IV                            -0.1071075  0.0589509  -1.817
## sc_age11V                             -0.1007663  0.0831150  -1.212
## sc_age11No male head                  -0.1313120  0.0885925  -1.482
## sc_age16I                              0.1855246  0.0691549   2.683
## sc_age16II                             0.1364621  0.0412928   3.305
## sc_age16III non-manual                -0.0182039  0.0541276  -0.336
## sc_age16III manual                     0.0441456  0.0353300   1.250
## sc_age16IV non-manual                  0.2170316  0.1096862   1.979
## sc_age16IV manual                     -0.0137819  0.0516702  -0.267
## sc_age16V                             -0.0838931  0.0800080  -1.049
## sc_age16Unclear                        0.2143485  0.1225825   1.749
## smoking_age16None,dont smoke           0.3878464  0.1530211   2.535
## smoking_age16Less 1 a week             0.3408968  0.1627244   2.095
## smoking_age161 to 9 a week             0.2338096  0.1599430   1.462
## smoking_age1610 to 19 a week           0.2276830  0.1636879   1.391
## smoking_age1620 to 29 a week           0.1841900  0.1617031   1.139
## smoking_age1630 to 39 a week           0.1613804  0.1620844   0.996
## smoking_age1640 to 49 a week           0.1927983  0.1644404   1.172
## smoking_age1650 to 59 a week           0.2704211  0.1697456   1.593
## smoking_age1660 or more a wk           0.2738217  0.1619680   1.691
## alcohol_age16Less than 1 week          0.0841146  0.1662780   0.506
## alcohol_age162 to 4 weeks              0.0612540  0.1673520   0.366
## alcohol_age165 to 8 weeks             -0.0171581  0.1710058  -0.100
## alcohol_age169 to 12 weeks            -0.0039322  0.1767265  -0.022
## alcohol_age16Over 12 weeks            -0.0290028  0.1704998  -0.170
## alcohol_age16Do not remember          -0.1199315  0.1678537  -0.715
## alcohol_age16Never had one            -0.2527040  0.1750479  -1.444
## school_attendance_age16                0.0007009  0.0003374   2.077
## age:tv_classHigh / rapidly decline TV  0.0015164  0.0007106   2.134

tab_model(cog_mixed_adj)

	z cog
Predictors	Estimates	CI	p
(Intercept)	-0.39	-1.45 – 0.66	0.467
age	-0.01	-0.01 – -0.01	<0.001
tv class [High / rapidly decline TV]	-0.06	-0.13 – 0.02	0.144
sex [Female]	0.02	-0.02 – 0.07	0.348
sc prenatal [I]	0.32	0.12 – 0.51	0.001
sc prenatal [II]	0.25	0.08 – 0.42	0.003
sc prenatal [III]	0.12	-0.03 – 0.28	0.120
sc prenatal [IV]	0.08	-0.09 – 0.24	0.368
sc prenatal [V]	0.05	-0.14 – 0.23	0.623
sc prenatal [Students]	0.21	-0.46 – 0.88	0.538
sc prenatal [Single,no husbnd]	0.04	-0.18 – 0.27	0.705
smoking preg [5mths-no change]	-0.08	-0.32 – 0.17	0.525
smoking preg [Now non-smoker]	-0.03	-0.29 – 0.23	0.815
smoking preg [1 to 4 now]	-0.18	-0.45 – 0.10	0.203
smoking preg [5 to 9 now]	-0.06	-0.34 – 0.23	0.697
smoking preg [10 to 14 now]	-0.13	-0.45 – 0.19	0.426
smoking preg [15 to 19 now]	-0.09	-0.52 – 0.33	0.664
smoking preg [20 to 24 now]	-0.29	-0.73 – 0.16	0.207
smoking preg [25 to 29 now]	-0.31	-1.09 – 0.47	0.438
smoking preg [30 or more now]	-0.16	-0.94 – 0.63	0.699
smoking_pregVariable-5mths	-0.11	-0.38 – 0.15	0.389
birth weight	0.00	-0.00 – 0.00	0.199
breast f [Dont know]	-0.00	-0.98 – 0.98	0.993
breast f [No]	-0.11	-1.02 – 0.81	0.817
breast f [Under one month]	-0.06	-0.97 – 0.86	0.903
breast f [Over one month]	0.03	-0.88 – 0.95	0.946
parent read level [Occasionally]	0.02	-0.09 – 0.13	0.667
parent read level [Weekly]	0.16	0.06 – 0.27	0.002
sc age7 [No male head]	0.26	-0.11 – 0.63	0.166
sc age7 [I]	0.19	-0.15 – 0.54	0.263
sc age7 [II]	0.18	-0.14 – 0.51	0.267
sc age7 [III non-manual]	0.20	-0.13 – 0.53	0.227
sc age7 [III manual]	0.02	-0.29 – 0.34	0.883
sc age7 [IV non-manual]	0.28	-0.09 – 0.65	0.137
sc age7 [IV manual]	0.03	-0.30 – 0.35	0.861
sc age7 [V]	-0.03	-0.37 – 0.31	0.871
sc age11 [I]	0.15	-0.00 – 0.30	0.057
sc age11 [II]	0.02	-0.09 – 0.14	0.692
sc age11 [III non manual]	0.08	-0.05 – 0.20	0.248
sc age11 [III manual]	-0.01	-0.11 – 0.09	0.883
sc age11 [IV]	-0.11	-0.22 – 0.01	0.069
sc age11 [V]	-0.10	-0.26 – 0.06	0.225
sc age11 [No male head]	-0.13	-0.30 – 0.04	0.138
sc age16 [I]	0.19	0.05 – 0.32	0.007
sc age16 [II]	0.14	0.06 – 0.22	0.001
sc age16 [III non-manual]	-0.02	-0.12 – 0.09	0.737
sc age16 [III manual]	0.04	-0.03 – 0.11	0.211
sc age16 [IV non-manual]	0.22	0.00 – 0.43	0.048
sc age16 [IV manual]	-0.01	-0.12 – 0.09	0.790
sc age16 [V]	-0.08	-0.24 – 0.07	0.294
sc age16 [Unclear]	0.21	-0.03 – 0.45	0.080
smoking age16 [None,dont smoke]	0.39	0.09 – 0.69	0.011
smoking age16 [Less 1 a week]	0.34	0.02 – 0.66	0.036
smoking age16 [1 to 9 a week]	0.23	-0.08 – 0.55	0.144
smoking age16 [10 to 19 a week]	0.23	-0.09 – 0.55	0.164
smoking age16 [20 to 29 a week]	0.18	-0.13 – 0.50	0.255
smoking age16 [30 to 39 a week]	0.16	-0.16 – 0.48	0.319
smoking age16 [40 to 49 a week]	0.19	-0.13 – 0.52	0.241
smoking age16 [50 to 59 a week]	0.27	-0.06 – 0.60	0.111
smoking age16 [60 or more a wk]	0.27	-0.04 – 0.59	0.091
alcohol age16 [Less than 1 week]	0.08	-0.24 – 0.41	0.613
alcohol age16 [2 to 4 weeks]	0.06	-0.27 – 0.39	0.714
alcohol age16 [5 to 8 weeks]	-0.02	-0.35 – 0.32	0.920
alcohol age16 [9 to 12 weeks]	-0.00	-0.35 – 0.34	0.982
alcohol age16 [Over 12 weeks]	-0.03	-0.36 – 0.31	0.865
alcohol age16 [Do not remember]	-0.12	-0.45 – 0.21	0.475
alcohol age16 [Never had one]	-0.25	-0.60 – 0.09	0.149
school attendance age16	0.00	0.00 – 0.00	0.038
age × tv class [High / rapidly decline TV]	0.00	0.00 – 0.00	0.033
Random Effects
σ2	0.49
τ00 id_num	0.33
ICC	0.40
N id_num	3323
Observations	16615
Marginal R2 / Conditional R2	0.115 / 0.470

df_cog_tv$tv_class <- factor(
  df_cog_tv$tv_class,
  levels = c(
    "High / slowly declining TV",
    "High / rapidly declining TV"
  )
)

pred_cog <- ggpredict(
  cog_mixed_adj,
  terms = c("age [7,11,16,50,62]", "tv_class")
)

pred_df <- as.data.frame(pred_cog)

fig_cog <- ggplot(
  pred_df,
  aes(
    x = x,
    y = predicted,
    color = group,
    fill  = group
  )
) +
  geom_line(size = 1.3) +

  scale_color_manual(
    values = c("#1B9E77", "#D95F02")
  ) +
  scale_fill_manual(
    values = c("#1B9E77", "#D95F02")
  ) +

  scale_x_continuous(
    breaks = c(7, 11, 16, 50, 62)
  ) +

  labs(
    x = "Age (years)",
    y = "Global cognitive function (z-score)",
    color = "TV viewing trajectory",
    fill  = "TV viewing trajectory"
  ) +

  theme_classic(base_size = 14) +
  theme(
    legend.position = "top",
    legend.title = element_text(size = 12),
    legend.text  = element_text(size = 11),
    axis.title   = element_text(size = 13),
    axis.text    = element_text(size = 11)
  )

fig_cog

int_term <- tidy(
  cog_mixed_adj,
  effects = "fixed",
  conf.int = TRUE
) %>%
  dplyr::filter(grepl("^age:tv_class", term)) %>%
  mutate(
    p.value = 2 * (1 - pnorm(abs(statistic)))
  )

int_term

## # A tibble: 1 × 8
##   effect term            estimate std.error statistic conf.low conf.high p.value
##   <chr>  <chr>              <dbl>     <dbl>     <dbl>    <dbl>     <dbl>   <dbl>
## 1 fixed  age:tv_classHi…  0.00152  0.000711      2.13 0.000124   0.00291  0.0328

int_label <- with(int_term, {
  paste0(
    "Interaction (age × TV trajectory):\n",
    "β = ", sprintf("%.4f", estimate),
    " (95% CI ", sprintf("%.4f", conf.low),
    ", ", sprintf("%.4f", conf.high), "), ",
    "p = ", ifelse(p.value < 0.001, "<0.001", sprintf("%.3f", p.value))
  )
})

int_label

## [1] "Interaction (age × TV trajectory):\nβ = 0.0015 (95% CI 0.0001, 0.0029), p = 0.033"

fig_cog <- fig_cog +
  annotate(
    "text",
    x = 8,
    y = min(pred_df$predicted) + 0.05,
    label = int_label,
    hjust = 0,
    size = 4.2
  )

fig_cog

GROUP-BASED TRAJECTORY MODELLING - READING

library(lcmm)

# df_complete_7_62$read3_age11 <- factor(
#   df_complete_7_62$child_read_books_age11,
#   levels = c("Hardly ever", "Sometimes", "most days"),
#   labels = c("Low", "Moderate", "High")
# )
# 
# df_complete_7_62$read3_age16 <- factor(
#   df_complete_7_62$child_read_age16,
#   levels = c("Hardly ever", "Sometimes", "Most days"),
#   labels = c("Low", "Moderate", "High")
# )


df_read <- df_complete_7_62 %>%
  dplyr::select(
    ncdsid,
    read3_age11,
    read3_age16,
    read3_age23,
    read3_age46,
    read3_age62
  ) %>%
  mutate(
    rdbk11 = as.numeric(read3_age11) - 1,
    rdbk16 = as.numeric(read3_age16) - 1,
    rdbk23 = as.numeric(read3_age23) - 1,
    rdbk46 = as.numeric(read3_age46) - 1,
    rdbk62 = as.numeric(read3_age62) - 1
  )

df_read_long <- df_read %>%
  pivot_longer(
    cols = c(rdbk11, rdbk16, rdbk23, rdbk46, rdbk62),
    names_to  = "wave",
    values_to = "rdbk_num"
  ) %>%
  mutate(
    age = case_when(
      wave == "rdbk11" ~ 11,
      wave == "rdbk16" ~ 16,
      wave == "rdbk23" ~ 23,
      wave == "rdbk46" ~ 46,
      wave == "rdbk62" ~ 62
    )
  ) %>%
  filter(!is.na(rdbk_num))


df_read_long <- df_read_long %>%
  mutate(
    id_num = as.numeric(as.factor(ncdsid))
  )

m1 <- hlme(
  fixed   = rdbk_num ~ age,
  random  = ~ 1,
  subject = "id_num",
  ng      = 1,
  data    = df_read_long
)

summary(m1)

## Heterogenous linear mixed model 
##      fitted by maximum likelihood method 
##  
## hlme(fixed = rdbk_num ~ age, random = ~1, subject = "id_num", 
##     ng = 1, data = df_read_long)
##  
## Statistical Model: 
##      Dataset: df_read_long 
##      Number of subjects: 3548 
##      Number of observations: 15512 
##      Number of latent classes: 1 
##      Number of parameters: 4  
##  
## Iteration process: 
##      Convergence criteria satisfied 
##      Number of iterations:  10 
##      Convergence criteria: parameters= 6.6e-12 
##                          : likelihood= 4.3e-08 
##                          : second derivatives= 2.4e-16 
##  
## Goodness-of-fit statistics: 
##      maximum log-likelihood: -17712.83  
##      AIC: 35433.65  
##      BIC: 35458.35  
##  
##  
## Maximum Likelihood Estimates: 
##  
## Fixed effects in the longitudinal model:
## 
##              coef      Se    Wald p-value
## intercept 1.11983 0.01224  91.507 0.00000
## age       0.00314 0.00029  10.867 0.00000
## 
## 
## Variance-covariance matrix of the random-effects:
##           intercept
## intercept   0.09587
## 
##                             coef      Se
## Residual standard error: 0.70727 0.00458

m2 <- hlme(
  fixed   = rdbk_num ~ age,
  mixture = ~ age,
  random  = ~ 1,
  subject = "id_num",
  ng      = 2,
  nwg     = TRUE,
  data    = df_read_long,
  B       = random(m1)
)

summary(m2)

## Heterogenous linear mixed model 
##      fitted by maximum likelihood method 
##  
## hlme(fixed = rdbk_num ~ age, mixture = ~age, random = ~1, subject = "id_num", 
##     ng = 2, nwg = TRUE, data = df_read_long)
##  
## Statistical Model: 
##      Dataset: df_read_long 
##      Number of subjects: 3548 
##      Number of observations: 15512 
##      Number of latent classes: 2 
##      Number of parameters: 8  
##  
## Iteration process: 
##      Convergence criteria satisfied 
##      Number of iterations:  19 
##      Convergence criteria: parameters= 3.6e-07 
##                          : likelihood= 7.2e-05 
##                          : second derivatives= 2.5e-05 
##  
## Goodness-of-fit statistics: 
##      maximum log-likelihood: -17536.53  
##      AIC: 35089.07  
##      BIC: 35138.46  
##  
##  
## Maximum Likelihood Estimates: 
##  
## Fixed effects in the class-membership model:
## (the class of reference is the last class) 
## 
##                      coef       Se    Wald p-value
## intercept class1  0.55493  0.05741   9.666 0.00000
## 
## Fixed effects in the longitudinal model:
## 
##                      coef       Se    Wald p-value
## intercept class1  1.19188  0.01695  70.337 0.00000
## intercept class2  0.96950  0.02273  42.645 0.00000
## age class1        0.00833  0.00045  18.701 0.00000
## age class2       -0.00548  0.00066  -8.308 0.00000
## 
## 
## Variance-covariance matrix of the random-effects:
##           intercept
## intercept         0
## 
##                                     coef       Se
## Proportional coefficient class1  1.15835 10.28575
## Residual standard error:         0.68754  0.00430

age_grid_rdbk <- data.frame(
  age = c(11, 16, 23, 46, 62)
)

pred_read <- predictY(
  m2,
  newdata = age_grid_rdbk,
  var.time = "age",
  draws = TRUE
)

str(pred_read)

## List of 2
##  $ pred : num [1:5, 1:6] 1.28 1.33 1.38 1.57 1.71 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:6] "Ypred_class1" "Ypred_class2" "lower.Ypred_class1" "lower.Ypred_class2" ...
##  $ times:'data.frame':   5 obs. of  1 variable:
##   ..$ age: num [1:5] 11 16 23 46 62
##  - attr(*, "class")= chr "predictY"

pred_df <- data.frame(
  age = pred_read$times$age,
  pred_read$pred
)

traj_df <- pred_df %>%
  pivot_longer(
    cols = -age,
    names_to = "term",
    values_to = "value"
  ) %>%
  separate(term, into = c("type", "class"), sep = "\\.") %>%
  pivot_wider(
    names_from  = type,
    values_from = value
  ) %>%
  mutate(
  class = dplyr::recode(
    class,
    "class1" = "Mid / slowly increasing",
    "class2" = "Mid / slowly declining"
  )
)

# Grid de idades observado
age_grid_rdbk <- data.frame(
  age = seq(
    min(df_read_long$age, na.rm = TRUE),
    max(df_read_long$age, na.rm = TRUE),
    by = 0.5
  )
)

# Predições com IC95%
pred_read <- predictY(
  m2,
  newdata = age_grid_rdbk,
  var.time = "age",
  draws = TRUE
)

traj_df <- data.frame(
  age = pred_read$times$age,

  mean_class1  = pred_read$pred[, "Ypred_class1"],
  lower_class1 = pred_read$pred[, "lower.Ypred_class1"],
  upper_class1 = pred_read$pred[, "upper.Ypred_class1"],

  mean_class2  = pred_read$pred[, "Ypred_class2"],
  lower_class2 = pred_read$pred[, "lower.Ypred_class2"],
  upper_class2 = pred_read$pred[, "upper.Ypred_class2"]
) %>%
  tidyr::pivot_longer(
    cols = -age,
    names_to = c("stat", "class"),
    names_sep = "_class",
    values_to = "value"
  ) %>%
  tidyr::pivot_wider(
    names_from  = stat,
    values_from = value
  ) %>%
  mutate(
    class = factor(
      class,
      levels = c("1", "2"),
      labels = c(
        "Mid / slowly increasing",
        "Mid / slowly declining"
      )
    )
  )

fig_read <- ggplot(traj_df, aes(x = age, y = mean, color = class, fill = class)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.25,
    color = NA
  ) +
  geom_line(size = 1.2) +
  scale_color_bmj() +
  scale_fill_bmj() +
  labs(
    x = "Age (years)",
    y = "Reading books (ordinal score)",
    color = "Trajectory",
    fill  = "Trajectory"
  ) +
  theme_classic(base_size = 14) +
  theme(
    legend.position = "top",
    legend.title = element_text(size = 12),
    legend.text  = element_text(size = 11),
    axis.title   = element_text(size = 13),
    axis.text    = element_text(size = 11)
  ) +

  scale_x_continuous(
    breaks = c(11, 16, 23, 46, 62)
  ) +
  scale_y_continuous(
    breaks = c(0, 1, 2),
    labels = c("Low", "Mid", "High")
  )  +
  coord_cartesian(ylim = c(0,2))
fig_read

COGNITIVE FUNCION - READING ACROSS LIFE COURSE

read_class_df <- m2$pprob %>%
  mutate(
    class = factor(
      class,
      levels = c(1, 2),
      labels = c("Mid / slowly increasing",
                 "Mid / slowly declining"
      )
    )
  ) %>%
  dplyr::select(id_num, class)

df_cog_long <- df_complete_7_62 %>%
  dplyr::select(
    ncdsid,
    z_global_cog_age7,
    z_global_cog_age11,
    z_global_cog_age16,
    global_cog_age50,
    global_cog_age62
  ) %>%
  pivot_longer(
    cols = -ncdsid,
    names_to  = "wave",
    values_to = "z_cog"
  ) %>%
  mutate(
    age = case_when(
      wave == "z_global_cog_age7"   ~ 7,
      wave == "z_global_cog_age11"  ~ 11,
      wave == "z_global_cog_age16"  ~ 16,
      wave == "global_cog_age50"    ~ 50,
      wave == "global_cog_age62"    ~ 62
    )
  ) %>%
  filter(!is.na(z_cog)) %>%
  mutate(
    id_num = as.numeric(as.factor(ncdsid))
  )

df_cov <- df_complete_7_62 %>%
  dplyr::select(
    ncdsid,

    # Demográficas
    sex,

    # Gestação / início da vida
    sc_prenatal,
    smoking_preg,
    birth_weight,
    breast_f,

    # Ambiente cognitivo precoce
    parent_read_level,

    # SEP ao longo da infância/adolescência
    sc_age7,
    sc_age11,
    sc_age16,

    # Comportamentos na adolescência
    smoking_age16,
    alcohol_age16,

    # Escolaridade
    school_attendance_age16
  ) %>%
  mutate(
    id_num = as.numeric(as.factor(ncdsid))
  )

df_read_class <- m2$pprob %>%
  mutate(
    class = factor(
      class,
      labels = c("Mid / slowly increasing",
                 "Mid / slowly declining"
      )
    )
  ) %>%
  dplyr::select(id_num, class)

df_cog_read <- df_cog_long %>%
  left_join(df_read_class, by = "id_num") %>%
  left_join(df_cov,       by = "id_num")


library(lme4)

cog_mixed_adj_read <- lmer(
  z_cog ~ age * class +

    sex +
    sc_prenatal +
    smoking_preg +
    birth_weight +
    breast_f +
    parent_read_level +

    sc_age7 +
    sc_age11 +
    sc_age16 +

    smoking_age16 +
    alcohol_age16 +
    school_attendance_age16 +

    (1 | id_num),

  data = df_cog_read,
  REML = TRUE
)

summary(cog_mixed_adj_read)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## z_cog ~ age * class + sex + sc_prenatal + smoking_preg + birth_weight +  
##     breast_f + parent_read_level + sc_age7 + sc_age11 + sc_age16 +  
##     smoking_age16 + alcohol_age16 + school_attendance_age16 +      (1 | id_num)
##    Data: df_cog_read
## 
## REML criterion at convergence: 40257.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1896 -0.5894  0.0254  0.6156  4.3684 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id_num   (Intercept) 0.3061   0.5533  
##  Residual             0.4907   0.7005  
## Number of obs: 16615, groups:  id_num, 3323
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                     -0.3481967  0.5235888  -0.665
## age                             -0.0052528  0.0003027 -17.351
## classMid / slowly declining     -0.3121482  0.0281932 -11.072
## sexFemale                       -0.0359514  0.0231250  -1.555
## sc_prenatalI                     0.2904172  0.0965589   3.008
## sc_prenatalII                    0.2468640  0.0834457   2.958
## sc_prenatalIII                   0.1126232  0.0775782   1.452
## sc_prenatalIV                    0.0662626  0.0832904   0.796
## sc_prenatalV                     0.0362757  0.0895236   0.405
## sc_prenatalStudents              0.2229251  0.3317859   0.672
## sc_prenatalSingle,no husbnd      0.0403501  0.1116042   0.362
## smoking_preg5mths-no change     -0.0861584  0.1216342  -0.708
## smoking_pregNow non-smoker      -0.0324485  0.1278469  -0.254
## smoking_preg1 to 4 now          -0.1998548  0.1363114  -1.466
## smoking_preg5 to 9 now          -0.0948206  0.1418933  -0.668
## smoking_preg10 to 14 now        -0.1609347  0.1596467  -1.008
## smoking_preg15 to 19 now        -0.0623483  0.2105964  -0.296
## smoking_preg20 to 24 now        -0.2606441  0.2210911  -1.179
## smoking_preg25 to 29 now        -0.1828599  0.3889067  -0.470
## smoking_preg30 or more now      -0.0652964  0.3915678  -0.167
## smoking_pregVariable-5mths      -0.1172919  0.1297335  -0.904
## birth_weight                     0.0002057  0.0001791   1.149
## breast_fDont know                0.0795266  0.4870190   0.163
## breast_fNo                       0.0084575  0.4544764   0.019
## breast_fUnder one month          0.0452303  0.4546001   0.099
## breast_fOver one month           0.1410480  0.4544063   0.310
## parent_read_levelOccasionally    0.0074227  0.0548468   0.135
## parent_read_levelWeekly          0.1203771  0.0513865   2.343
## sc_age7No male head              0.2814421  0.1828456   1.539
## sc_age7I                         0.1753646  0.1692044   1.036
## sc_age7II                        0.1662552  0.1611662   1.032
## sc_age7III non-manual            0.1788856  0.1628196   1.099
## sc_age7III manual                0.0185931  0.1581256   0.118
## sc_age7IV non-manual             0.2469120  0.1823950   1.354
## sc_age7IV manual                 0.0297135  0.1608884   0.185
## sc_age7V                        -0.0323077  0.1681291  -0.192
## sc_age11I                        0.1438080  0.0765241   1.879
## sc_age11II                       0.0293913  0.0567625   0.518
## sc_age11III non manual           0.0856989  0.0635541   1.348
## sc_age11III manual               0.0118617  0.0507378   0.234
## sc_age11IV                      -0.0716551  0.0574624  -1.247
## sc_age11V                       -0.0755440  0.0809533  -0.933
## sc_age11No male head            -0.1193322  0.0862423  -1.384
## sc_age16I                        0.1605842  0.0673443   2.385
## sc_age16II                       0.1424864  0.0401989   3.545
## sc_age16III non-manual          -0.0303708  0.0527018  -0.576
## sc_age16III manual               0.0446027  0.0344016   1.297
## sc_age16IV non-manual            0.1930562  0.1068169   1.807
## sc_age16IV manual               -0.0237786  0.0503138  -0.473
## sc_age16V                       -0.0891071  0.0779066  -1.144
## sc_age16Unclear                  0.1860697  0.1193716   1.559
## smoking_age16None,dont smoke     0.4028653  0.1488658   2.706
## smoking_age16Less 1 a week       0.3669589  0.1583311   2.318
## smoking_age161 to 9 a week       0.2496321  0.1556662   1.604
## smoking_age1610 to 19 a week     0.2654166  0.1592950   1.666
## smoking_age1620 to 29 a week     0.2214319  0.1573228   1.408
## smoking_age1630 to 39 a week     0.1847409  0.1576866   1.172
## smoking_age1640 to 49 a week     0.2346005  0.1600647   1.466
## smoking_age1650 to 59 a week     0.3128703  0.1652177   1.894
## smoking_age1660 or more a wk     0.3167267  0.1575958   2.010
## alcohol_age16Less than 1 week    0.0595445  0.1618000   0.368
## alcohol_age162 to 4 weeks        0.0262022  0.1628272   0.161
## alcohol_age165 to 8 weeks       -0.0421912  0.1664045  -0.254
## alcohol_age169 to 12 weeks      -0.0194674  0.1718986  -0.113
## alcohol_age16Over 12 weeks      -0.0708415  0.1658964  -0.427
## alcohol_age16Do not remember    -0.1223890  0.1632964  -0.749
## alcohol_age16Never had one      -0.2705589  0.1703399  -1.588
## school_attendance_age16          0.0006436  0.0003284   1.960
## age:classMid / slowly declining -0.0002845  0.0005065  -0.562

tab_model(cog_mixed_adj_read)

	z cog
Predictors	Estimates	CI	p
(Intercept)	-0.35	-1.37 – 0.68	0.506
age	-0.01	-0.01 – -0.00	<0.001
class [Mid / slowly declining]	-0.31	-0.37 – -0.26	<0.001
sex [Female]	-0.04	-0.08 – 0.01	0.120
sc prenatal [I]	0.29	0.10 – 0.48	0.003
sc prenatal [II]	0.25	0.08 – 0.41	0.003
sc prenatal [III]	0.11	-0.04 – 0.26	0.147
sc prenatal [IV]	0.07	-0.10 – 0.23	0.426
sc prenatal [V]	0.04	-0.14 – 0.21	0.685
sc prenatal [Students]	0.22	-0.43 – 0.87	0.502
sc prenatal [Single,no husbnd]	0.04	-0.18 – 0.26	0.718
smoking preg [5mths-no change]	-0.09	-0.32 – 0.15	0.479
smoking preg [Now non-smoker]	-0.03	-0.28 – 0.22	0.800
smoking preg [1 to 4 now]	-0.20	-0.47 – 0.07	0.143
smoking preg [5 to 9 now]	-0.09	-0.37 – 0.18	0.504
smoking preg [10 to 14 now]	-0.16	-0.47 – 0.15	0.313
smoking preg [15 to 19 now]	-0.06	-0.48 – 0.35	0.767
smoking preg [20 to 24 now]	-0.26	-0.69 – 0.17	0.238
smoking preg [25 to 29 now]	-0.18	-0.95 – 0.58	0.638
smoking preg [30 or more now]	-0.07	-0.83 – 0.70	0.868
smoking_pregVariable-5mths	-0.12	-0.37 – 0.14	0.366
birth weight	0.00	-0.00 – 0.00	0.251
breast f [Dont know]	0.08	-0.88 – 1.03	0.870
breast f [No]	0.01	-0.88 – 0.90	0.985
breast f [Under one month]	0.05	-0.85 – 0.94	0.921
breast f [Over one month]	0.14	-0.75 – 1.03	0.756
parent read level [Occasionally]	0.01	-0.10 – 0.11	0.892
parent read level [Weekly]	0.12	0.02 – 0.22	0.019
sc age7 [No male head]	0.28	-0.08 – 0.64	0.124
sc age7 [I]	0.18	-0.16 – 0.51	0.300
sc age7 [II]	0.17	-0.15 – 0.48	0.302
sc age7 [III non-manual]	0.18	-0.14 – 0.50	0.272
sc age7 [III manual]	0.02	-0.29 – 0.33	0.906
sc age7 [IV non-manual]	0.25	-0.11 – 0.60	0.176
sc age7 [IV manual]	0.03	-0.29 – 0.35	0.853
sc age7 [V]	-0.03	-0.36 – 0.30	0.848
sc age11 [I]	0.14	-0.01 – 0.29	0.060
sc age11 [II]	0.03	-0.08 – 0.14	0.605
sc age11 [III non manual]	0.09	-0.04 – 0.21	0.178
sc age11 [III manual]	0.01	-0.09 – 0.11	0.815
sc age11 [IV]	-0.07	-0.18 – 0.04	0.212
sc age11 [V]	-0.08	-0.23 – 0.08	0.351
sc age11 [No male head]	-0.12	-0.29 – 0.05	0.166
sc age16 [I]	0.16	0.03 – 0.29	0.017
sc age16 [II]	0.14	0.06 – 0.22	<0.001
sc age16 [III non-manual]	-0.03	-0.13 – 0.07	0.564
sc age16 [III manual]	0.04	-0.02 – 0.11	0.195
sc age16 [IV non-manual]	0.19	-0.02 – 0.40	0.071
sc age16 [IV manual]	-0.02	-0.12 – 0.07	0.637
sc age16 [V]	-0.09	-0.24 – 0.06	0.253
sc age16 [Unclear]	0.19	-0.05 – 0.42	0.119
smoking age16 [None,dont smoke]	0.40	0.11 – 0.69	0.007
smoking age16 [Less 1 a week]	0.37	0.06 – 0.68	0.020
smoking age16 [1 to 9 a week]	0.25	-0.06 – 0.55	0.109
smoking age16 [10 to 19 a week]	0.27	-0.05 – 0.58	0.096
smoking age16 [20 to 29 a week]	0.22	-0.09 – 0.53	0.159
smoking age16 [30 to 39 a week]	0.18	-0.12 – 0.49	0.241
smoking age16 [40 to 49 a week]	0.23	-0.08 – 0.55	0.143
smoking age16 [50 to 59 a week]	0.31	-0.01 – 0.64	0.058
smoking age16 [60 or more a wk]	0.32	0.01 – 0.63	0.044
alcohol age16 [Less than 1 week]	0.06	-0.26 – 0.38	0.713
alcohol age16 [2 to 4 weeks]	0.03	-0.29 – 0.35	0.872
alcohol age16 [5 to 8 weeks]	-0.04	-0.37 – 0.28	0.800
alcohol age16 [9 to 12 weeks]	-0.02	-0.36 – 0.32	0.910
alcohol age16 [Over 12 weeks]	-0.07	-0.40 – 0.25	0.669
alcohol age16 [Do not remember]	-0.12	-0.44 – 0.20	0.454
alcohol age16 [Never had one]	-0.27	-0.60 – 0.06	0.112
school attendance age16	0.00	-0.00 – 0.00	0.050
age × class [Mid / slowly declining]	-0.00	-0.00 – 0.00	0.574
Random Effects
σ2	0.49
τ00 id_num	0.31
ICC	0.38
N id_num	3323
Observations	16615
Marginal R2 / Conditional R2	0.138 / 0.469

df_cog_read$class <- factor(
  df_cog_read$class,
  levels = c(
    "High / stable",
    "Low / slow increase"
  )
)

pred_cog <- ggpredict(
  cog_mixed_adj_read,
  terms = c("age [11,16,23,46,50, 62]", "class")
)

pred_df <- as.data.frame(pred_cog)

fig_cog_read <- ggplot(
  pred_df,
  aes(
    x = x,
    y = predicted,
    color = group,
    fill  = group
  )
) +
  geom_line(size = 1.3) +

  scale_color_bmj() +
  scale_fill_bmj() +

  scale_x_continuous(
    breaks = c(11, 16, 23, 46, 50, 62)
  ) +

  labs(
    x = "Age (years)",
    y = "Global cognitive function (z-score)",
    color = "Reading books trajectory",
    fill  = "Reading books trajectory"
  ) +

  theme_classic(base_size = 14) +
  theme(
    legend.position = "top",
    legend.title = element_text(size = 12),
    legend.text  = element_text(size = 11),
    axis.title   = element_text(size = 13),
    axis.text    = element_text(size = 11)
  )

fig_cog_read

int_term <- tidy(
  cog_mixed_adj_read,
  effects = "fixed",
  conf.int = TRUE
) %>%
  dplyr::filter(grepl("^age:class", term)) %>%
  mutate(
    p.value = 2 * (1 - pnorm(abs(statistic)))
  )

int_term

## # A tibble: 1 × 8
##   effect term            estimate std.error statistic conf.low conf.high p.value
##   <chr>  <chr>              <dbl>     <dbl>     <dbl>    <dbl>     <dbl>   <dbl>
## 1 fixed  age:classMid /… -2.85e-4  0.000507    -0.562 -0.00128  0.000708   0.574

int_label <- with(int_term, {
  paste0(
    "Interaction (age × reading books trajectory):\n",
    "β = ", sprintf("%.4f", estimate),
    " (95% CI ", sprintf("%.4f", conf.low),
    ", ", sprintf("%.4f", conf.high), "), ",
    "p = ", ifelse(p.value < 0.001, "<0.001", sprintf("%.3f", p.value))
  )
})

int_label

## [1] "Interaction (age × reading books trajectory):\nβ = -0.0003 (95% CI -0.0013, 0.0007), p = 0.574"

fig_cog_read <- fig_cog_read +
  annotate(
    "text",
    x = 8,
    y = min(pred_df$predicted) + 0.05,
    label = int_label,
    hjust = 0,
    size = 4.2
  )

fig_cog_read