Amniotic Testosterone & Autistic Traits in Adolescence
Niamh Dooley
Introduction
This R-code was written by Niamh Dooley (niamhdooley@rcsi.com) and analyzes data from the Longitudinal Foetal Testosterone Study in Cambridge (PI: Prof. Simon Baron-Cohen)
Data Imported: Two waves of the FT longitudinal study were used in this analysis, referred to as the 2013 & 2016 waves. Imported here are those that met inclusion criteria (valid pT value, at least one AQ questionnaire completed, no genetic abnormalities).
Load libraries & data
library(tidyverse)
library(tidyr)
library(broom)
library(knitr)
library(lme4)
library(ggplot2); theme_set(theme_minimal())
library(readxl)
library(RColorBrewer)
library(dvmisc)
library(ggeffects)
library(stargazer)
library(ggpubr)
library(dplyr)
library(qwraps2)
library(emmeans)
library(gridExtra)
library(MASS)
library(effects)
library(sjmisc)
library(lsr)load("ftdata.Rda")
head(smallft[,c(2,3,5,6,7,8, 15, 16, 17)]) # show a selection of non-identifiable info## GENDER AGE_YRS fT_NUM gen_puberty_score AQ_chrep AQ_parrep chAQ_SBC1
## 1 boy 13 0.65 2.4 13 3 0
## 2 boy 16 0.79 3.0 19 17 4
## 3 boy 15 1.25 2.8 15 5 1
## 4 girl 15 0.40 3.6 14 6 0
## 5 boy 15 1.10 2.8 20 7 5
## 6 girl 16 0.35 3.6 23 10 2
## chAQ_SBC2 chAQ_SBC3
## 1 3 0
## 2 4 5
## 3 4 2
## 4 4 1
## 5 5 3
## 6 7 5# fT_NUM: amniotic testosterone concentration (nmol/L)
# gen_puberty_score: pubertal stage as determined by the Pubertal Development Scale (Peterson, 1988)
# AQ_chrep: child-report AQ total score (from the Adult AQ; Baron-Cohen et al.,2002)
# AQ_parrep: parent-report AQ total score (from the Adolescent AQ; Baron-Cohen et al.,2006)
# SBC1-5 refer to the subscales social skills, attention switching, communication, imagination & attention to detailExplore & Clean data
Age
One individual was missing age. Their age were mean imputed.
# find subjects with missing ages
miss_age = sum(is.na(smallft$AGE_MTHS)) # n=1
miss_age_ids = which(is.na(smallft$AGE_MTHS))
# Mean Impute...
smallft$AGE_MTHS.imp <- smallft$AGE_MTHS
for(i in 1:length(miss_age_ids)) {
smallft$AGE_MTHS.imp[miss_age_ids[i]] = mean(smallft$AGE_MTHS, na.rm = TRUE)
}
smallft <- move_columns(smallft, AGE_MTHS.imp, .after = "AGE_MTHS")
# use this mean-imputed age (in months) variable (AGE_MTHS.imp) for analysis
smallft$agemz <- scale(smallft$AGE_MTHS.imp)Prenatal Testosterone
ggplot(data = smallft, aes(x = fT_NUM) )+
geom_histogram(aes(x=fT_NUM), alpha=0.5, colour="black",
fill = "lightblue") +
geom_density(alpha = 0.5, fill = "blue") +
labs(y = "Count", x = "Prenatal Testosterone (nmol/L)") +
scale_fill_brewer(palette = "Set1") +
theme_minimal()+
theme(legend.title = element_blank())+
facet_grid(~GENDER)
Histograms above show a female participant with a pT value way outside the typical female range. We decided to exclude her from the analysis.
rows <- which(smallft$fT_NUM==2.3)
smallft <- smallft[-c(rows), ]
# re-create histograms again
ggplot(data = smallft, aes(x = fT_NUM) )+
geom_histogram(aes(x=fT_NUM), alpha=0.5, colour="black",
fill = "lightblue") +
geom_density(alpha = 0.5, fill = "blue") +
labs(y = "Count", x = "Prenatal Testosterone (nmol/L)") +
scale_fill_brewer(palette = "Set1") +
theme_minimal()+
theme(legend.title = element_blank())+
facet_grid(~GENDER)
Pubertal stage
From PDS questionnaire. Note male mean PDS score is the ave of items on (1) height (2) body hair (3)skin (4)voice-drop (5)facial hair
While female mean PDS score reflects items on: (1)height (2)body hair (3)skin (4)breasts (5)menstruation (y/n)
# distribution of PDS stage for males & females separately
ggplot(data = smallft, aes(x = gen_puberty_score, fill=GENDER)) +
geom_histogram(color="black", binwidth=0.10, alpha=0.5) +
labs(y = "N", x = " ", title= "Pubertal Stage") +
facet_wrap(~GENDER) +
theme_minimal() +
scale_fill_brewer(palette = "Set1") +
theme(legend.position="none")
# what is the r'ship between pubertal stage and age (for boys & girls)?
ggplot(smallft, aes(AGE_MTHS.imp , gen_puberty_score, colour=GENDER)) +
geom_point(aes(colour=GENDER), alpha = 0.8) +
theme_classic()+
geom_smooth(method = 'lm', formula= y~x) 
# within-sex correlations
boys <- smallft[smallft$GENDER=="boy",]
girls <- smallft[smallft$GENDER=="girl",]
cor(x=boys$gen_puberty_score, y=boys$AGE_MTHS.imp, use="complete.obs") # male r=.64## [1] 0.6376187cor(x=girls$gen_puberty_score, y=girls$AGE_MTHS.imp, use="complete.obs") # female r=.51## [1] 0.5074126Pubertal Timing
Create pubertal timing variables as per https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2756189/
# steps:
# create sex-specific dataset for puberty analysis
# remove all NA cases of puberty score
# compute age (agemz) ~ PDS stage (gen_puberty_score), extract residuals lm[["residuals"]] and add to df as a variable
# Males
boys <- smallft[smallft$GENDER=="boy",]
rows <-which(is.na(boys$gen_puberty_score))
if (!is_empty(rows)) {
boys <- boys[-c(rows), ]}
mod <- lm(gen_puberty_score ~ agemz, data=boys)
boys$pub.timing <- mod$residuals
boys <- move_columns(boys, pub.timing, .after = "gen_puberty_score")
# Females
girls <- smallft[smallft$GENDER=="girl",]
rows <-which(is.na(girls$gen_puberty_score))
if (!is_empty(rows)) {
girls <- girls[-c(rows), ]}
mod <- lm(gen_puberty_score ~ agemz, data=girls)
girls$pub.timing <- mod$resid
girls <- move_columns(girls, pub.timing, .after = "gen_puberty_score")
# put male and female df's together into 1 df
ft4pub <- rbind(boys, girls)
# using a loop, copy pub.timing data back into the main df "smallft" (issues using direct merge due to missing puberty data)
smallft$pub.timing <- NA
for (i in 1:length(smallft$ID)) {
row <- which(ft4pub$ID == smallft$ID[i])
if (!is_empty(row)) {
smallft$pub.timing[i] <-ft4pub$pub.timing[row]}
}
smallft <- move_columns(smallft, pub.timing, .before = "gen_puberty_score")Standardize/centre variables
# pubertal stage (note pubertal timing is already a z-score)
smallft$pub.stage <- smallft$gen_puberty_score
smallft$pub.stage.z <- scale(smallft$gen_puberty_score, center = T, scale=F)
smallft$mat.age.z <- scale(smallft$mat.age, center = T, scale=F)
smallft$GA.z <- scale(smallft$GA, center = T, scale=F)
smallft$ftz <- scale(smallft$fT_NUM)Descriptives
| Sample Summary | Full Sample | Girls | Boys |
|---|---|---|---|
| Age (years) | Â Â | Â Â | Â Â |
| Â Â Min | 13 | 13 | 13 |
| Â Â Max | 21 | 21 | 21 |
| Â Â Mean (SD) | 95; 15.60 (1.78) | 15.80 (1.80) | 55; 15.46 (1.77) |
| Maternal Age (years) | Â Â | Â Â | Â Â |
| Â Â Min | 23 | 25 | 23 |
| Â Â Max | 45 | 45 | 45 |
| Â Â Mean (SD) | 73; 34.96 (4.75) | 30; 35.87 (4.36) | 43; 34.33 (4.96) |
| GA @ amnio | Â Â | Â Â | Â Â |
| Â Â Min | 14 | 15 | 14 |
| Â Â Max | 21 | 18 | 21 |
| Â Â Mean (SD) | 50; 16.64 (1.48) | 15; 16.47 (1.06) | 35; 16.71 (1.64) |
| Stage of Puberty | Â Â | Â Â | Â Â |
| Â Â Min | 1.8 | 3.0 | 1.8 |
| Â Â Max | 4 | 4.0 | 3.8 |
| Â Â Mean (SD) | 95; 3.30 (0.45) | 40; 3.62 (0.31) | 55; 3.07 (0.39) |
| Puberty Timing | Â Â | Â Â | Â Â |
| Â Â Min | -1.04864 | -0.4506405 | -1.0486402 |
| Â Â Max | 0.6444229 | 0.4978616 | 0.6444229 |
| Â Â Mean (SD) | 95; 0.00 (0.29) | 40; 0.00 (0.27) | 55; 0.00 (0.30) |
| Foetal Testosterone levels (nmol/L) | Â Â | Â Â | Â Â |
| Â Â Min | 0.1 | 0.15 | 0.10 |
| Â Â Max | 2.05 | 0.80 | 2.05 |
| Â Â Mean (SD) | 96; 0.66 (0.45) | 40; 0.29 (0.13) | 56; 0.92 (0.41) |
| Autism Quotient (self-report) | Â Â | Â Â | Â Â |
| Â Â Min | 4 | 6 | 4 |
| Â Â Max | 39 | 27 | 39 |
| Â Â Mean (SD) | 86; 16.23 (6.71) | 39; 14.67 (5.54) | 47; 17.53 (7.35) |
| Autism Quotient (parent-report) | Â Â | Â Â | Â Â |
| Â Â Min | 3 | 4 | 3 |
| Â Â Max | 40 | 15 | 40 |
| Â Â Mean (SD) | 85; 12.98 (8.11) | 35; 9.06 (2.97) | 50; 15.72 (9.38) |
| SOCIAL (self-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 0 | 0 |
| Â Â Max | 10 | 7 | 10 |
| Â Â Mean (SD) | 2.34 (2.35) | 1.87 (1.85) | 2.72 (2.65) |
| SOCIAL (parent-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 0 | 0 |
| Â Â Max | 10 | 3 | 10 |
| Â Â Mean (SD) | 2.01 (2.40) | 0.91 (0.98) | 2.78 (2.78) |
| ATT-SWITCH (self-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 0 | 1 |
| Â Â Max | 10 | 7 | 10 |
| Â Â Mean (SD) | 4.08 (2.00) | 3.67 (1.75) | 4.43 (2.14) |
| ATT-SWITCH (parent-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 0 | 0 |
| Â Â Max | 9 | 5 | 9 |
| Â Â Mean (SD) | 2.86 (2.09) | 2.09 (1.48) | 3.40 (2.29) |
| COMMUNIC (self-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 0 | 0 |
| Â Â Max | 8 | 6 | 8 |
| Â Â Mean (SD) | 2.23 (2.00) | 1.79 (1.92) | 2.60 (2.01) |
| COMMUNIC (parent-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 0 | 0 |
| Â Â Max | 10 | 6 | 10 |
| Â Â Mean (SD) | 2.01 (2.20) | 1.17 (1.22) | 2.60 (2.53) |
| IMAGIN (self-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 0 | 0 |
| Â Â Max | 7 | 6 | 7 |
| Â Â Mean (SD) | 2.50 (1.73) | 2.05 (1.43) | 2.87 (1.87) |
| IMAGIN (parent-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 0 | 0 |
| Â Â Max | 9 | 4 | 9 |
| Â Â Mean (SD) | 2.08 (1.96) | 0.94 (0.94) | 2.88 (2.10) |
| DETAIL (self-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 1 | 0 |
| Â Â Max | 9 | 9 | 9 |
| Â Â Mean (SD) | 5.08 (2.17) | 5.28 (2.10) | 4.91 (2.22) |
| DETAIL (parent-report) | Â Â | Â Â | Â Â |
| Â Â Min | 0 | 1 | 0 |
| Â Â Max | 9 | 9 | 9 |
| Â Â Mean (SD) | 4.01 (2.19) | 3.94 (1.89) | 4.06 (2.39) |
Create long format df
ftl = reshape(smallft, varying = list(c("AQ_chrep", "AQ_parrep"),
c("chAQ_fullscale", "parAQ_fullscale"),
c("chAQ_H1", "parAQ_H1"),
c("chAQ_H2", "parAQ_H2"),
c("chAQ_SBC1", "parAQ_SBC1"),
c("chAQ_SBC2", "parAQ_SBC2"),
c("chAQ_SBC3", "parAQ_SBC3"),
c("chAQ_SBC4", "parAQ_SBC4"),
c("chAQ_SBC5", "parAQ_SBC5")),
v.names =c("AQ","AQ_fullvar", "AQ_H1", "AQ_H2", "SOCIAL","ATT_SWITCH","COMMUNIC","IMAGIN","DETAIL"),
idvar = "ID",
timevar = "rater",
direction="long")
ftl <- ftl %>%
mutate(
rater = factor(rater,
levels = c (1,2),
labels = c("self", "parent")),
ID = factor(ID))Plot AQ distributions
AQ total analysis
# name dataset "d"
d = ftl
# set contrast for binary vars (had to repeat within lmer function)
contrasts(d$GENDER) <- c(-.5, .5)
contrasts(d$rater) <- c(-.5, .5)Analysis of AQ scores consisted of 3 consecutive models.
The first model (a) includes pT, sex, interaction pT-x-sex, pubertal factor, interaction pubertal factor-x-pT and the within-subject variable of rater (parent/teen)
Model (b) adjusts for maternal age at time of birth
Model (c) adjusts for the time of amniocentesis (gestational age in weeks)
These models were run separately due to the drop off in N due to missing data for maternal age and time of amnio.
# pubertal stage model
M1a.AQ.gauss <- lmer(AQ ~ ftz + GENDER + ftz:GENDER + rater
+ pub.stage.z + ftz:pub.stage.z
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
M1b.AQ.gauss <- lmer(AQ ~ ftz + GENDER + ftz:GENDER + rater
+ pub.stage.z + ftz:pub.stage.z
+ mat.age
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
M1c.AQ.gauss <- lmer(AQ ~ ftz + GENDER + ftz:GENDER + rater
+ pub.stage.z + ftz:pub.stage.z
+ mat.age + GA
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
# as above but with uncentred pubertal stage (for plots)
M1b.AQ.gauss_ <- lmer(AQ ~ ftz + GENDER + rater + pub.stage
+ ftz:GENDER
+ ftz:pub.stage
+ mat.age
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
stargazer(M1a.AQ.gauss, M1b.AQ.gauss, M1c.AQ.gauss, report=('vcstp'), digits = 2, digit.separator = "", order = c(), align = TRUE, type = "text", title="PUBERTAL STAGE MODEL")##
## PUBERTAL STAGE MODEL
## ====================================================
## Dependent variable:
## --------------------------------
## AQ
## (1) (2) (3)
## ----------------------------------------------------
## ftz -0.92 -0.32 -1.77
## (1.67) (1.86) (4.09)
## t = -0.55 t = -0.17 t = -0.43
## p = 0.59 p = 0.87 p = 0.67
##
## GENDER1 7.44 5.17 5.41
## (3.12) (3.39) (8.29)
## t = 2.38 t = 1.52 t = 0.65
## p = 0.02 p = 0.13 p = 0.52
##
## rater1 -3.61 -3.53 -3.09
## (0.83) (0.91) (1.25)
## t = -4.36 t = -3.86 t = -2.46
## p = 0.0001 p = 0.0002 p = 0.02
##
## pub.stage.z 2.50 1.42 -1.92
## (1.79) (2.01) (2.65)
## t = 1.40 t = 0.71 t = -0.72
## p = 0.17 p = 0.48 p = 0.47
##
## mat.age -0.01 -0.12
## (0.16) (0.20)
## t = -0.05 t = -0.58
## p = 0.96 p = 0.56
##
## GA -2.07
## (0.61)
## t = -3.37
## p = 0.001
##
## ftz:GENDER1 1.63 2.04 9.02
## (3.49) (3.74) (8.05)
## t = 0.47 t = 0.55 t = 1.12
## p = 0.64 p = 0.59 p = 0.27
##
## ftz:pub.stage.z 3.74 5.36 9.67
## (2.30) (2.48) (3.14)
## t = 1.63 t = 2.16 t = 3.08
## p = 0.11 p = 0.04 p = 0.003
##
## Constant 14.28 14.84 51.46
## (1.50) (5.56) (14.68)
## t = 9.52 t = 2.67 t = 3.51
## p = 0.00 p = 0.01 p = 0.0005
##
## ----------------------------------------------------
## Observations 170 139 94
## Log Likelihood -551.18 -447.11 -297.53
## Akaike Inf. Crit. 1120.35 914.22 617.05
## Bayesian Inf. Crit. 1148.57 943.57 645.03
## ====================================================
## Note: *p<0.1; **p<0.05; ***p<0.01# pubertal timing model
M2a.AQ.gauss <- lmer(AQ ~ ftz + GENDER + ftz:GENDER + rater
+ pub.timing + ftz:pub.timing
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
M2b.AQ.gauss <- lmer(AQ ~ ftz + GENDER + ftz:GENDER + rater
+ pub.timing + ftz:pub.timing
+ mat.age
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
M2c.AQ.gauss <- lmer(AQ ~ ftz + GENDER + ftz:GENDER + rater
+ pub.timing + ftz:pub.timing
+ mat.age + GA
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
stargazer(M2a.AQ.gauss, M2b.AQ.gauss, M2c.AQ.gauss, report=('vcstp'), digits = 2, digit.separator = "", order = c(), align = TRUE, type = "text", title="PUBERTAL TIMING MODEL")##
## PUBERTAL TIMING MODEL
## ===================================================
## Dependent variable:
## -------------------------------
## AQ
## (1) (2) (3)
## ---------------------------------------------------
## ftz -0.98 -0.41 -1.87
## (1.65) (1.82) (3.99)
## t = -0.60 t = -0.22 t = -0.47
## p = 0.56 p = 0.83 p = 0.64
##
## GENDER1 6.13 4.35 6.30
## (2.94) (3.17) (7.64)
## t = 2.08 t = 1.37 t = 0.82
## p = 0.04 p = 0.17 p = 0.41
##
## rater1 -3.67 -3.54 -3.13
## (0.83) (0.92) (1.26)
## t = -4.43 t = -3.86 t = -2.49
## p = 0.0000 p = 0.0002 p = 0.02
##
## pub.timing 2.13 1.88 -1.32
## (2.14) (2.21) (2.65)
## t = 1.00 t = 0.85 t = -0.50
## p = 0.32 p = 0.40 p = 0.62
##
## mat.age 0.03 -0.05
## (0.16) (0.20)
## t = 0.22 t = -0.27
## p = 0.83 p = 0.79
##
## GA -2.02
## (0.59)
## t = -3.44
## p = 0.001
##
## ftz:GENDER1 -0.57 -0.95 3.30
## (3.31) (3.52) (8.00)
## t = -0.17 t = -0.27 t = 0.41
## p = 0.87 p = 0.79 p = 0.68
##
## ftz:pub.timing 5.68 6.91 10.59
## (2.45) (2.60) (3.10)
## t = 2.32 t = 2.65 t = 3.41
## p = 0.03 p = 0.01 p = 0.001
##
## Constant 14.30 13.29 48.75
## (1.47) (5.51) (14.23)
## t = 9.72 t = 2.41 t = 3.43
## p = 0.00 p = 0.02 p = 0.001
##
## ---------------------------------------------------
## Observations 170 139 94
## Log Likelihood -550.07 -445.87 -296.63
## Akaike Inf. Crit. 1118.14 911.74 615.26
## Bayesian Inf. Crit. 1146.36 941.08 643.24
## ===================================================
## Note: *p<0.1; **p<0.05; ***p<0.01# Use model (b) to plot results to balance sample size with control of confounds
# EMMs for main effect of ft
main.predicted <- effect(c("ftz"), M1b.AQ.gauss,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25)),
se=TRUE, confidence.level = .95,
typical=mean)
main.pred <- as.data.frame(main.predicted)
# EMMs for main effect of ft for each sex
predicted.sex <- effect(c("ftz", "GENDER"), M1b.AQ.gauss,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25)),
se=TRUE, confidence.level = .95,
typical=mean)
pred.sex <- as.data.frame(predicted.sex)
pred.sex$GENDER = relevel(pred.sex$GENDER, "girl")
ggplot()+
geom_point(data = d,
aes(x = ftz, y = AQ, colour = GENDER, shape=rater))+
geom_ribbon(data = pred.sex, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=GENDER, fill=GENDER), alpha = 0.2) +
geom_line(data = pred.sex, aes(x = ftz, y = fit,
group = GENDER, colour=GENDER),size=1) +
geom_line(data = main.pred, aes(x =ftz, y = fit),size=1, linetype="dashed") +
labs(x = "Foetal Testosterone", y = "Latent AQ total (rater-averaged)") +
theme_minimal() 
predicted <- effect(c("ftz", "pub.stage"), M1b.AQ.gauss_,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.stage = seq(from=2, to=4, by=1)),
se=TRUE, confidence.level = .95,
typical=mean)
pred <- as.data.frame(predicted)
predicted.sex <- effect(c("ftz", "GENDER", "pub.stage"), M1b.AQ.gauss_,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.stage = seq(from=2, to=4, by=1)),
se=TRUE, confidence.level = .95,
typical=mean)
pred.sex <- as.data.frame(predicted.sex)
pred.sex$GENDER = relevel(pred.sex$GENDER, "girl")
gsex <- ggplot()+
geom_ribbon(data = pred.sex, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.stage, fill=pub.stage), alpha = 0.2) +
geom_line(data = pred.sex, aes(x = ftz, y = fit,
group = pub.stage, colour=pub.stage),size=1, linetype="dashed") +
labs(x = "", y = "", fill = "Pubertal Stage\n", colour="Pubertal Stage\n") +
theme_minimal()+
coord_cartesian(ylim = c(0,35)) +
scale_colour_gradient(low="orange", high="springgreen4")+
scale_fill_gradient(low="orange", high="springgreen4")+
facet_grid(~GENDER)
gfull<-ggplot()+
geom_ribbon(data = pred, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.stage, fill=pub.stage), alpha = 0.2) +
geom_line(data = pred, aes(x = ftz, y = fit,
group = pub.stage, colour=pub.stage),size=1, linetype="dashed") +
labs(x = "", y = "", fill = "Pubertal Stage\n", colour="Pubertal Stage\n") +
theme_minimal()+
coord_cartesian(ylim = c(0,35)) +
scale_colour_gradient(low="orange", high="springgreen4")+
scale_fill_gradient(low="orange", high="springgreen4")+
theme(legend.position = "none" )
gboth<-ggarrange(gfull, gsex,nrow=2, ncol=1)
annotate_figure(gboth,
bottom=text_grob("Prenatal Testosterone (PT)", size=16),
left = text_grob("\nAQ Total", rot=90, size=16))
predicted.sex <- effect(c("ftz", "GENDER", "pub.timing"), M2b.AQ.gauss,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.timing = seq(from=-0.5, to=0.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean)
pred.sex <- as.data.frame(predicted.sex)
pred.sex$GENDER = relevel(pred.sex$GENDER, "girl")
gsex<-ggplot()+
geom_ribbon(data = pred.sex, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.2) +
geom_line(data = pred.sex, aes(x = ftz, y = fit,
group = pub.timing, colour=pub.timing),size=1) +
labs(x = "", y = "", fill = "Pubertal Timing\n", colour="Pubertal Timing\n") +
theme_minimal()+
coord_cartesian(ylim = c(0,35)) +
scale_colour_gradient(low="orange", high="springgreen4")+
scale_fill_gradient(low="orange", high="springgreen4")+
facet_grid(~GENDER)
predicted<- effect(c("ftz", "pub.timing"), M2b.AQ.gauss,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.timing = seq(from=-0.5, to=0.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean)
pred <- as.data.frame(predicted)
gfull<-ggplot()+
geom_ribbon(data = pred, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.2) +
geom_line(data = pred, aes(x = ftz, y = fit,
group = pub.timing, colour=pub.timing),size=1) +
labs(x = "", y = "", fill = "Pubertal Timing\n", colour="Pubertal Timing\n") +
theme_minimal()+
coord_cartesian(ylim = c(0,35)) +
scale_colour_gradient(low="orange", high="springgreen4")+
scale_fill_gradient(low="orange", high="springgreen4")+
theme(legend.position = "none" )
gboth<-ggarrange(gfull, gsex,nrow=2, ncol=1)
annotate_figure(gboth,
bottom=text_grob("Prenatal Testosterone (PT)", size=16),
left = text_grob("\nAQ Total", rot=90, size=16))
PT may have a non-linear relationship with AQ. Explore if a quadratic pT term is correlated with AQ and whether plots show a strong non-linearity.
M1a.quad.lmm <- lmer(AQ ~ ftz + I(ftz^2) + GENDER + rater + pub.stage.z
+ ftz:GENDER + ftz:pub.stage.z
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
M1b.quad.lmm <- lmer(AQ ~ ftz + I(ftz^2) + GENDER + rater + pub.stage.z
+ ftz:GENDER + ftz:pub.stage.z + mat.age
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
M1c.quad.lmm <- lmer(AQ ~ ftz + I(ftz^2) + GENDER + rater + pub.stage.z
+ ftz:GENDER + ftz:pub.stage.z + mat.age + GA
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
stargazer(M1a.quad.lmm, M1b.quad.lmm, M1c.quad.lmm, digits = 2, report=('vcsp'), digit.separator = "", order = c(), align = TRUE, type = "text")##
## ====================================================
## Dependent variable:
## --------------------------------
## AQ
## (1) (2) (3)
## ----------------------------------------------------
## ftz -0.83 -0.03 -1.97
## (1.70) (1.88) (4.10)
## p = 0.63 p = 0.99 p = 0.64
##
## I(ftz2) -0.34 -0.70 -0.79
## (0.70) (0.71) (0.89)
## p = 0.63 p = 0.33 p = 0.38
##
## GENDER1 7.63 5.48 6.73
## (3.16) (3.41) (8.43)
## p = 0.02 p = 0.11 p = 0.43
##
## rater1 -3.61 -3.52 -3.09
## (0.83) (0.91) (1.26)
## p = 0.0001 p = 0.0002 p = 0.02
##
## pub.stage.z 2.63 1.69 -1.33
## (1.82) (2.03) (2.73)
## p = 0.15 p = 0.41 p = 0.63
##
## mat.age -0.005 -0.15
## (0.16) (0.20)
## p = 0.98 p = 0.47
##
## GA -2.18
## (0.63)
## p = 0.001
##
## ftz:GENDER1 2.59 4.10 11.92
## (4.02) (4.29) (8.71)
## p = 0.52 p = 0.34 p = 0.18
##
## ftz:pub.stage.z 3.67 5.35 9.42
## (2.31) (2.48) (3.16)
## p = 0.12 p = 0.04 p = 0.003
##
## Constant 14.26 14.72 53.91
## (1.51) (5.56) (14.96)
## p = 0.00 p = 0.01 p = 0.0004
##
## ----------------------------------------------------
## Observations 170 139 94
## Log Likelihood -550.50 -446.05 -296.33
## Akaike Inf. Crit. 1121.00 914.10 616.66
## Bayesian Inf. Crit. 1152.36 946.38 647.18
## ====================================================
## Note: *p<0.1; **p<0.05; ***p<0.01M2a.quad.lmm <- lmer(AQ ~ ftz + I(ftz^2) + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
M2b.quad.lmm <- lmer(AQ ~ ftz + I(ftz^2) + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing + mat.age
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
M2c.quad.lmm <- lmer(AQ ~ ftz + I(ftz^2) + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing + mat.age + GA
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
stargazer(M2a.quad.lmm, M2b.quad.lmm, M2c.quad.lmm, digits = 2, report=('vcsp'), digit.separator = "", order = c(), align = TRUE, type = "text")##
## ===================================================
## Dependent variable:
## -------------------------------
## AQ
## (1) (2) (3)
## ---------------------------------------------------
## ftz -0.94 -0.23 -1.96
## (1.67) (1.85) (4.04)
## p = 0.58 p = 0.91 p = 0.63
##
## I(ftz2) -0.17 -0.46 -0.17
## (0.69) (0.71) (0.91)
## p = 0.81 p = 0.52 p = 0.85
##
## GENDER1 6.19 4.47 6.61
## (2.97) (3.19) (7.85)
## p = 0.04 p = 0.17 p = 0.40
##
## rater1 -3.67 -3.53 -3.12
## (0.83) (0.92) (1.26)
## p = 0.0000 p = 0.0002 p = 0.02
##
## pub.timing 2.19 2.08 -1.17
## (2.16) (2.24) (2.81)
## p = 0.32 p = 0.36 p = 0.68
##
## mat.age 0.03 -0.06
## (0.16) (0.20)
## p = 0.83 p = 0.77
##
## GA -2.04
## (0.61)
## p = 0.001
##
## ftz:GENDER1 -0.07 0.45 4.08
## (3.89) (4.13) (8.94)
## p = 0.99 p = 0.92 p = 0.65
##
## ftz:pub.timing 5.62 6.69 10.39
## (2.48) (2.64) (3.29)
## p = 0.03 p = 0.02 p = 0.002
##
## Constant 14.29 13.27 49.23
## (1.48) (5.54) (14.65)
## p = 0.00 p = 0.02 p = 0.001
##
## ---------------------------------------------------
## Observations 170 139 94
## Log Likelihood -549.50 -445.08 -295.79
## Akaike Inf. Crit. 1118.99 912.17 615.58
## Bayesian Inf. Crit. 1150.35 944.44 646.10
## ===================================================
## Note: *p<0.1; **p<0.05; ***p<0.01# plot quadratic pT effects
predicted <- as.data.frame(effect(c("ftz", "pub.stage.z"), M1b.quad.lmm,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.stage.z = seq(from=-0.5, to=0.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean))
ggplot()+
geom_ribbon(data = predicted, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.stage.z, fill=pub.stage.z), alpha = 0.2) +
geom_line(data = predicted, aes(x = ftz, y = fit,
group = pub.stage.z, colour=pub.stage.z),size=1) +
labs(x = "Foetal Testosterone", y = "Latent AQ total (rater-averaged)", title = "Quadratic FT effects, across 3 pubertal stages") +
theme_minimal() 
predicted <- as.data.frame(effect(c("ftz", "pub.timing"), M2b.quad.lmm,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.timing = seq(from=-0.5, to=0.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean))
ggplot()+
geom_ribbon(data = predicted, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.2) +
geom_line(data = predicted, aes(x = ftz, y = fit,
group = pub.timing, colour=pub.timing),size=1) +
labs(x = "Foetal Testosterone", y = "Latent AQ total (rater-averaged)", title = "Quadratic FT effects, across 3 values of pubertal timing") +
theme_minimal() 
Sensitivity Analysis 1: simpler models
simplerAQ.a <- lmer(AQ ~ ftz + GENDER
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5)))
simplerAQ.b <- lmer(AQ ~ ftz*GENDER
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5)))
simplerAQ.c <- lmer(AQ ~ ftz*GENDER + rater
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
stargazer(simplerAQ.a, simplerAQ.b,simplerAQ.c,report=('vcsp'),
digits = 2, digit.separator = "",
order = c(), align = TRUE, type = "text")##
## =================================================
## Dependent variable:
## -----------------------------
## AQ
## (1) (2) (3)
## -------------------------------------------------
## ftz -0.96 -0.93 -0.92
## (0.90) (1.74) (1.76)
## p = 0.29 p = 0.60 p = 0.61
##
## GENDER1 6.20 6.14 6.42
## (1.82) (3.10) (3.14)
## p = 0.001 p = 0.05 p = 0.05
##
## rater1 -3.45
## (0.83)
## p = 0.0001
##
## ftz:GENDER1 -0.08 0.08
## (3.48) (3.53)
## p = 0.99 p = 0.99
##
## Constant 14.27 14.30 14.23
## (0.65) (1.55) (1.57)
## p = 0.00 p = 0.00 p = 0.00
##
## -------------------------------------------------
## Observations 171 171 171
## Log Likelihood -574.74 -572.58 -564.08
## Akaike Inf. Crit. 1159.48 1157.15 1142.17
## Bayesian Inf. Crit. 1175.19 1176.00 1164.16
## =================================================
## Note: *p<0.1; **p<0.05; ***p<0.01rows <- which(d$rater=="self")
d_self <- d[rows,]
rows <- which(d$rater=="parent")
d_parent <- d[rows,]
selfAQ1a <- lm(AQ ~ ftz + GENDER + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z ,
data = d_self,
contrasts=list(GENDER=c(-0.5,0.5)))
selfAQ1b <- lm(AQ ~ ftz + GENDER + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z + mat.age.z,
data = d_self,
contrasts=list(GENDER=c(-0.5,0.5)))
selfAQ1c<- lm(AQ ~ ftz + GENDER + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z + mat.age.z + GA.z,
data = d_self,
contrasts=list(GENDER=c(-0.5,0.5)))
selfAQ2a <- lm(AQ ~ ftz + GENDER + pub.timing
+ ftz:GENDER
+ ftz:pub.timing ,
data = d_self,
contrasts=list(GENDER=c(-0.5,0.5)))
selfAQ2b <- lm(AQ ~ ftz + GENDER + pub.timing
+ ftz:GENDER
+ ftz:pub.timing + mat.age.z,
data = d_self,
contrasts=list(GENDER=c(-0.5,0.5)))
selfAQ2c <- lm(AQ ~ ftz + GENDER + pub.timing
+ ftz:GENDER
+ ftz:pub.timing + mat.age.z + GA.z,
data = d_self,
contrasts=list(GENDER=c(-0.5,0.5)))
stargazer(selfAQ1a, selfAQ1b, selfAQ1c, report=('vcsp'), digits = 2, digit.separator = "",
order = c(), align = TRUE, type = "text")##
## ============================================================================
## Dependent variable:
## --------------------------------------------------------
## AQ
## (1) (2) (3)
## ----------------------------------------------------------------------------
## ftz -0.85 0.84 4.68
## (1.90) (1.95) (5.09)
## p = 0.66 p = 0.67 p = 0.37
##
## GENDER1 4.83 0.45 -8.83
## (3.53) (3.56) (10.22)
## p = 0.18 p = 0.90 p = 0.40
##
## pub.stage.z 1.49 0.14 -2.21
## (2.05) (2.16) (2.45)
## p = 0.47 p = 0.95 p = 0.38
##
## mat.age.z 0.04 0.05
## (0.17) (0.19)
## p = 0.80 p = 0.81
##
## GA.z -2.37
## (0.58)
## p = 0.0003
##
## ftz:GENDER1 0.89 0.81 -0.87
## (3.96) (3.91) (10.07)
## p = 0.83 p = 0.84 p = 0.94
##
## ftz:pub.stage.z 3.34 5.80 13.42
## (2.60) (2.52) (2.91)
## p = 0.21 p = 0.03 p = 0.0001
##
## Constant 16.41 17.14 21.12
## (1.71) (1.72) (4.97)
## p = 0.00 p = 0.00 p = 0.0002
##
## ----------------------------------------------------------------------------
## Observations 86 66 45
## R2 0.08 0.10 0.49
## Adjusted R2 0.03 0.01 0.40
## Residual Std. Error 6.62 (df = 80) 5.81 (df = 59) 4.90 (df = 37)
## F Statistic 1.46 (df = 5; 80) 1.11 (df = 6; 59) 5.15*** (df = 7; 37)
## ============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01stargazer(selfAQ2a, selfAQ2b, selfAQ2c, report=('vcsp'), digits = 2, digit.separator = "",
order = c(), align = TRUE, type = "text")##
## ============================================================================
## Dependent variable:
## --------------------------------------------------------
## AQ
## (1) (2) (3)
## ----------------------------------------------------------------------------
## ftz -0.90 0.75 5.06
## (1.87) (1.90) (4.94)
## p = 0.64 p = 0.70 p = 0.32
##
## GENDER1 4.16 0.39 -8.92
## (3.35) (3.34) (9.65)
## p = 0.22 p = 0.91 p = 0.37
##
## pub.timing 0.28 0.64 -2.01
## (2.49) (2.31) (2.43)
## p = 0.91 p = 0.79 p = 0.42
##
## mat.age.z 0.09 0.13
## (0.17) (0.18)
## p = 0.59 p = 0.49
##
## GA.z -2.29
## (0.55)
## p = 0.0002
##
## ftz:GENDER1 -1.03 -2.52 -9.84
## (3.75) (3.67) (9.92)
## p = 0.79 p = 0.50 p = 0.33
##
## ftz:pub.timing 5.73 7.60 14.41
## (2.79) (2.65) (2.86)
## p = 0.05 p = 0.01 p = 0.0001
##
## Constant 16.40 17.18 21.95
## (1.68) (1.67) (4.86)
## p = 0.00 p = 0.00 p = 0.0001
##
## ----------------------------------------------------------------------------
## Observations 86 66 45
## R2 0.11 0.14 0.53
## Adjusted R2 0.05 0.05 0.44
## Residual Std. Error 6.54 (df = 80) 5.68 (df = 59) 4.75 (df = 37)
## F Statistic 1.92 (df = 5; 80) 1.62 (df = 6; 59) 5.85*** (df = 7; 37)
## ============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01parentAQ1a <- lm(AQ ~ ftz + GENDER + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z ,
data = d_parent,
contrasts=list(GENDER=c(-0.5,0.5)))
parentAQ1b <- lm(AQ ~ ftz + GENDER + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z + mat.age.z,
data = d_parent,
contrasts=list(GENDER=c(-0.5,0.5)))
parentAQ1c <- lm(AQ ~ ftz + GENDER + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z + mat.age.z + GA.z,
data = d_parent,
contrasts=list(GENDER=c(-0.5,0.5)))
parentAQ2a <- lm(AQ ~ ftz + GENDER + pub.timing
+ ftz:GENDER
+ ftz:pub.timing ,
data = d_parent,
contrasts=list(GENDER=c(-0.5,0.5)))
parentAQ2b <- lm(AQ ~ ftz + GENDER + pub.timing
+ ftz:GENDER
+ ftz:pub.timing + mat.age.z,
data = d_parent,
contrasts=list(GENDER=c(-0.5,0.5)))
parentAQ2c <- lm(AQ ~ ftz + GENDER + pub.timing
+ ftz:GENDER
+ ftz:pub.timing + mat.age.z + GA.z,
data = d_parent,
contrasts=list(GENDER=c(-0.5,0.5)))
stargazer(parentAQ1a, parentAQ1b, parentAQ1c, report=('vcsp'), digits = 2, digit.separator = "", align = TRUE, type = "text")##
## ==============================================================================
## Dependent variable:
## ----------------------------------------------------------
## AQ
## (1) (2) (3)
## ------------------------------------------------------------------------------
## ftz -0.65 -1.09 -3.58
## (2.11) (2.42) (5.63)
## p = 0.77 p = 0.66 p = 0.53
##
## GENDER1 8.67 8.70 10.97
## (4.00) (4.43) (11.47)
## p = 0.04 p = 0.06 p = 0.35
##
## pub.stage.z 2.71 2.09 -1.27
## (2.35) (2.60) (4.07)
## p = 0.26 p = 0.43 p = 0.76
##
## mat.age.z -0.04 -0.25
## (0.21) (0.30)
## p = 0.84 p = 0.42
##
## GA.z -1.93
## (0.93)
## p = 0.05
##
## ftz:GENDER1 1.46 2.48 9.98
## (4.40) (4.89) (11.06)
## p = 0.74 p = 0.62 p = 0.38
##
## ftz:pub.stage.z 4.14 4.63 6.31
## (2.90) (3.27) (4.84)
## p = 0.16 p = 0.17 p = 0.20
##
## Constant 12.45 12.23 9.58
## (1.90) (2.10) (5.36)
## p = 0.00 p = 0.0000 p = 0.09
##
## ------------------------------------------------------------------------------
## Observations 84 73 49
## R2 0.20 0.18 0.22
## Adjusted R2 0.15 0.10 0.08
## Residual Std. Error 7.13 (df = 78) 7.59 (df = 66) 8.21 (df = 41)
## F Statistic 4.01*** (df = 5; 78) 2.39** (df = 6; 66) 1.62 (df = 7; 41)
## ==============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01stargazer(parentAQ2a, parentAQ2b, parentAQ2c, report=('vcsp'), digits = 2, digit.separator = "",
align = TRUE, type = "text")##
## ==============================================================================
## Dependent variable:
## ----------------------------------------------------------
## AQ
## (1) (2) (3)
## ------------------------------------------------------------------------------
## ftz -0.81 -1.19 -3.75
## (2.10) (2.40) (5.54)
## p = 0.71 p = 0.63 p = 0.51
##
## GENDER1 7.28 7.52 11.74
## (3.76) (4.17) (10.50)
## p = 0.06 p = 0.08 p = 0.28
##
## pub.timing 3.56 2.68 -0.65
## (2.60) (2.92) (4.17)
## p = 0.18 p = 0.37 p = 0.88
##
## mat.age.z -0.01 -0.20
## (0.21) (0.30)
## p = 0.98 p = 0.51
##
## GA.z -1.91
## (0.91)
## p = 0.05
##
## ftz:GENDER1 -0.72 -0.05 6.42
## (4.21) (4.66) (11.11)
## p = 0.87 p = 1.00 p = 0.57
##
## ftz:pub.timing 5.39 5.98 7.38
## (3.05) (3.48) (4.86)
## p = 0.09 p = 0.10 p = 0.14
##
## Constant 12.32 12.13 9.63
## (1.88) (2.08) (5.33)
## p = 0.00 p = 0.0000 p = 0.08
##
## ------------------------------------------------------------------------------
## Observations 84 73 49
## R2 0.22 0.19 0.23
## Adjusted R2 0.17 0.12 0.10
## Residual Std. Error 7.08 (df = 78) 7.54 (df = 66) 8.14 (df = 41)
## F Statistic 4.32*** (df = 5; 78) 2.58** (df = 6; 66) 1.74 (df = 7; 41)
## ==============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Sensitivity Analysis 2: sex-stratified models
dboys <- d[d$GENDER=="boy",]
# center variable within sex
dboys$pub.stage <- dboys$gen_puberty_score
dboys$pub.stage.z <- scale(dboys$gen_puberty_score, center = T, scale=F)
# run tests in males
boy.AQ1a <- lmer(AQ ~ ftz + rater + pub.stage.z
+ ftz:pub.stage.z
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))
boy.AQ1b <- lmer(AQ ~ ftz + rater + pub.stage.z
+ ftz:pub.stage.z + mat.age.z
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))
boy.AQ1c <- lmer(AQ ~ ftz + rater + pub.stage.z
+ ftz:pub.stage.z + mat.age.z + GA.z
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))
stargazer(boy.AQ1a, boy.AQ1b, boy.AQ1c, report=('vcsp'), digits = 2, digit.separator = "",
align = TRUE, type = "text")##
## =================================================
## Dependent variable:
## -----------------------------
## AQ
## (1) (2) (3)
## -------------------------------------------------
## ftz -0.91 -0.41 0.33
## (1.14) (1.21) (1.25)
## p = 0.43 p = 0.74 p = 0.80
##
## rater1 -1.78 -1.67 -1.76
## (1.24) (1.37) (1.55)
## p = 0.16 p = 0.23 p = 0.26
##
## pub.stage.z 3.88 -1.01 -6.40
## (3.22) (4.02) (3.80)
## p = 0.23 p = 0.81 p = 0.10
##
## mat.age.z 0.04 -0.24
## (0.24) (0.25)
## p = 0.87 p = 0.33
##
## GA.z -2.65
## (0.69)
## p = 0.0002
##
## ftz:pub.stage.z 2.27 9.07 16.11
## (4.45) (5.29) (4.79)
## p = 0.62 p = 0.09 p = 0.001
##
## Constant 17.32 16.69 16.15
## (1.25) (1.28) (1.21)
## p = 0.00 p = 0.00 p = 0.00
##
## -------------------------------------------------
## Observations 96 80 67
## Log Likelihood -323.68 -265.67 -214.48
## Akaike Inf. Crit. 661.36 547.35 446.95
## Bayesian Inf. Crit. 679.31 566.40 466.80
## =================================================
## Note: *p<0.1; **p<0.05; ***p<0.01boy.AQ2a <- lmer(AQ ~ ftz + rater + pub.timing
+ ftz:pub.timing
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))
boy.AQ2b <- lmer(AQ ~ ftz + rater + pub.timing
+ ftz:pub.timing + mat.age.z
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))
boy.AQ2c <- lmer(AQ ~ ftz + rater + pub.timing
+ ftz:pub.timing + mat.age.z + GA.z
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))
stargazer(boy.AQ2a, boy.AQ2b, boy.AQ2c, report=('vcsp'), digits = 2, digit.separator = "",
align = TRUE, type = "text")##
## =================================================
## Dependent variable:
## -----------------------------
## AQ
## (1) (2) (3)
## -------------------------------------------------
## ftz -1.46 -0.90 -0.84
## (1.11) (1.22) (1.32)
## p = 0.19 p = 0.47 p = 0.53
##
## rater1 -1.96 -1.71 -1.80
## (1.25) (1.37) (1.55)
## p = 0.12 p = 0.22 p = 0.25
##
## pub.timing -0.48 0.11 -4.64
## (4.12) (4.03) (3.79)
## p = 0.91 p = 0.98 p = 0.23
##
## mat.age.z 0.09 -0.22
## (0.24) (0.25)
## p = 0.71 p = 0.40
##
## GA.z -2.51
## (0.68)
## p = 0.0003
##
## ftz:pub.timing 8.61 9.53 14.66
## (4.75) (4.83) (4.45)
## p = 0.08 p = 0.05 p = 0.001
##
## Constant 17.47 16.65 16.69
## (1.19) (1.26) (1.23)
## p = 0.00 p = 0.00 p = 0.00
##
## -------------------------------------------------
## Observations 96 80 67
## Log Likelihood -322.81 -265.02 -214.60
## Akaike Inf. Crit. 659.61 546.04 447.21
## Bayesian Inf. Crit. 677.56 565.10 467.05
## =================================================
## Note: *p<0.1; **p<0.05; ***p<0.01## run pub stage models again without standardized version of pub.stage (for plots)
boy.AQ1a_ <- lmer(AQ ~ ftz + rater + pub.stage
+ ftz:pub.stage
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))
boy.AQ1b_ <- lmer(AQ ~ ftz + rater + pub.stage
+ ftz:pub.stage + mat.age.z
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))
boy.AQ1c_ <- lmer(AQ ~ ftz + rater + pub.stage
+ ftz:pub.stage + mat.age.z + GA.z
+ (1 | ID),
data = dboys,
contrasts=list(rater=c(-0.5,0.5)))dgirls <- d[d$GENDER=="girl",]
dgirls$pub.stage <- dgirls$gen_puberty_score
dgirls$pub.stage.z <- scale(dgirls$gen_puberty_score, center = T, scale=F)
# remove the girl with the highest FT just to see if shes leveraging results.
# rows <- which(dgirls$ID==153)
# dgirls <- dgirls[-c(rows),]
girl.AQ1a <- lmer(AQ ~ ftz + rater + pub.stage.z
+ ftz:pub.stage.z
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))
girl.AQ1b <- lmer(AQ ~ ftz + rater + pub.stage.z
+ ftz:pub.stage.z + mat.age.z
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))
girl.AQ1c <- lmer(AQ ~ ftz + rater + pub.stage.z
+ ftz:pub.stage.z + mat.age.z + GA.z
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))
stargazer(girl.AQ1a, girl.AQ1b, girl.AQ1c, report=('vcsp'), digits = 2, digit.separator = "",
align = TRUE, type = "text")##
## ===================================================
## Dependent variable:
## -------------------------------
## AQ
## (1) (2) (3)
## ---------------------------------------------------
## ftz -1.77 -2.29 6.71
## (2.26) (3.00) (6.70)
## p = 0.44 p = 0.45 p = 0.32
##
## rater1 -5.79 -5.89 -5.84
## (0.85) (0.96) (1.81)
## p = 0.00 p = 0.00 p = 0.002
##
## pub.stage.z -8.74 -10.51 -43.17
## (6.97) (9.70) (21.20)
## p = 0.21 p = 0.28 p = 0.05
##
## mat.age.z -0.13 0.36
## (0.18) (0.26)
## p = 0.47 p = 0.17
##
## GA.z 0.88
## (1.18)
## p = 0.46
##
## ftz:pub.stage.z -8.05 -9.50 -32.43
## (7.89) (10.60) (20.52)
## p = 0.31 p = 0.38 p = 0.12
##
## Constant 10.25 10.32 19.30
## (1.96) (2.58) (6.54)
## p = 0.0000 p = 0.0001 p = 0.004
##
## ---------------------------------------------------
## Observations 74 59 27
## Log Likelihood -204.87 -162.10 -64.79
## Akaike Inf. Crit. 423.74 340.20 147.58
## Bayesian Inf. Crit. 439.87 356.82 159.24
## ===================================================
## Note: *p<0.1; **p<0.05; ***p<0.01girl.AQ2a <- lmer(AQ ~ ftz + rater + pub.timing
+ ftz:pub.timing
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))
girl.AQ2b <- lmer(AQ ~ ftz + rater + pub.timing
+ ftz:pub.timing + mat.age.z
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))
girl.AQ2c <- lmer(AQ ~ ftz + rater + pub.timing
+ ftz:pub.timing + mat.age.z + GA.z
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))
stargazer(girl.AQ2a, girl.AQ2b, girl.AQ2c, report=('vcsp'), digits = 2, digit.separator = "",
align = TRUE, type = "text")##
## =================================================
## Dependent variable:
## -----------------------------
## AQ
## (1) (2) (3)
## -------------------------------------------------
## ftz -1.00 -0.96 9.88
## (2.18) (2.53) (6.85)
## p = 0.65 p = 0.71 p = 0.15
##
## rater1 -5.73 -5.91 -5.86
## (0.86) (0.96) (1.71)
## p = 0.00 p = 0.00 p = 0.001
##
## pub.timing -2.92 -10.08 -51.34
## (10.13) (11.51) (21.72)
## p = 0.78 p = 0.39 p = 0.02
##
## mat.age.z -0.12 0.60
## (0.18) (0.27)
## p = 0.51 p = 0.03
##
## GA.z 1.25
## (1.13)
## p = 0.27
##
## ftz:pub.timing -2.19 -8.57 -37.92
## (11.52) (12.87) (20.95)
## p = 0.85 p = 0.51 p = 0.08
##
## Constant 10.92 11.61 23.16
## (1.90) (2.15) (6.81)
## p = 0.00 p = 0.0000 p = 0.001
##
## -------------------------------------------------
## Observations 74 59 27
## Log Likelihood -205.22 -161.88 -63.49
## Akaike Inf. Crit. 424.44 339.77 144.98
## Bayesian Inf. Crit. 440.57 356.39 156.65
## =================================================
## Note: *p<0.1; **p<0.05; ***p<0.01## run another version of pub.stage models without standardized pub.stage (for plots)
girl.AQ1a_ <- lmer(AQ ~ ftz + rater + pub.stage
+ ftz:pub.stage
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))
girl.AQ1b_ <- lmer(AQ ~ ftz + rater + pub.stage
+ ftz:pub.stage + mat.age.z
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))
girl.AQ1c_ <- lmer(AQ ~ ftz + rater + pub.stage
+ ftz:pub.stage + mat.age.z + GA.z
+ (1 | ID),
data = dgirls,
contrasts=list(rater=c(-0.5,0.5)))# set colours
lcol <- "deepskyblue"
hcol <- "blue4"
lcol2 <- "turquoise"
hcol2 <- "darkcyan"
lcol3 <- "red"
hcol3 <- "darkred"
lcol4 <- "orange"
hcol4 <- "darkorange4"
lpub <-"orange"
hpub <- "springgreen4"
# M1. pubertal stage
boypred <- as.data.frame(effect(c("ftz", "pub.stage"), boy.AQ1b_,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.stage = seq(from=2, to=4, by=1)),
se=TRUE, confidence.level = .95,
typical=mean))
girlpred <- as.data.frame(effect(c("ftz", "pub.stage"), girl.AQ1b_,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.stage = seq(from=2, to=4, by=1)),
se=TRUE, confidence.level = .95,
typical=mean))
boy1 <- ggplot()+
#geom_point(data = dboys, aes(x=ftz, y=AQ, colour=pub.stage))+
geom_ribbon(data = boypred, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.stage, fill=pub.stage), alpha = 0.3) +
geom_line(data = boypred, aes(x = ftz, y = fit,
group=pub.stage, colour=pub.stage), size=1, linetype="dashed") +
labs(x = "", y = "", fill = "Male\nPubertal Stage", colour = "Male\nPubertal Stage") +
scale_colour_gradient(low=lpub, high=hpub)+
scale_fill_gradient(low=lpub, high=hpub)+
coord_cartesian(ylim = c(0,40)) +
theme_minimal() +
theme(legend.position = "none")+
annotate("text", label = "Males", x = 0.7, y = 37.5, size = 6, colour = "black")
girl1 <- ggplot()+
#geom_point(data = dgirls, aes(x=ftz, y=AQ, colour=pub.stage))+
geom_ribbon(data = girlpred, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.stage, fill=pub.stage), alpha = 0.3) +
geom_line(data = girlpred, aes(x = ftz, y = fit,
group=pub.stage, colour=pub.stage), size=1,linetype="dashed") +
labs(x = "", y = "", fill = "Female\nPubertal Stage", colour = "Female\nPubertal Stage") +
scale_colour_gradient(low=lpub, high=hpub)+
scale_fill_gradient(low=lpub, high=hpub)+
coord_cartesian(ylim = c(0,40)) +
theme_minimal() +
theme(legend.position = "none")+
annotate("text", label = "Females", x = 0.7, y = 37.5, size = 6, colour = "black")
p1<-ggarrange(girl1, boy1, nrow=1, ncol=2)
annotate_figure(p1,
left = text_grob("\nAQ Total", rot=90, size=14),
bottom = text_grob("Prenatal Testosterone (standardised)", size=14))
# M2. pubertal timing
boypred <- as.data.frame(effect(c("ftz", "pub.timing"), boy.AQ2b,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.timing = seq(from=-0.5, to=0.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean))
girlpred <- as.data.frame(effect(c("ftz", "pub.timing"), girl.AQ2b,
xlevels= list(ftz = seq(from=-1.25, to=3, by=0.25),
pub.timing = seq(from=-0.5, to=0.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean))
boy2 <- ggplot()+
geom_ribbon(data = boypred, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.3) +
geom_line(data = boypred, aes(x = ftz, y = fit,
group=pub.timing, colour=pub.timing), size=1) +
labs(x = "", y = "", fill = "Male\nPubertal Timing", colour = "Male\nPubertal Timing") +
scale_colour_gradient(low=lpub, high=hpub)+
scale_fill_gradient(low=lpub, high=hpub)+
coord_cartesian(ylim = c(0,40)) +
theme_minimal() +
theme(legend.position = "none")+
annotate("text", label = "Males", x = 0.7, y = 37.5, size = 6, colour = "black")
girl2 <- ggplot()+
geom_ribbon(data = girlpred, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.3) +
geom_line(data = girlpred, aes(x = ftz, y = fit,
group=pub.timing, colour=pub.timing), size=1) +
labs(x = "", y = "", fill = "Female\nPubertal Timing", colour = "Female\nPubertal Timing") +
scale_colour_gradient(low=lpub, high=hpub)+
scale_fill_gradient(low=lpub, high=hpub)+
coord_cartesian(ylim = c(0,40)) +
theme_minimal() +
theme(legend.position = "none")+
annotate("text", label = "Females", x = 0.7, y = 37.5, size = 6, colour = "black")
p3<-ggarrange(girl2, boy2, nrow=1, ncol=2)
annotate_figure(p3,
left = text_grob("\nAQ Total", rot=90, size=14),
bottom = text_grob("Prenatal Testosterone (standardised)", size=14))
Sensitivity Analysis 3: generalized linear mixed models
M1a.AQ.gamma <- glmer(AQ ~ ftz + GENDER + rater + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)),
family = Gamma(link=identity),
control = glmerControl(optimizer = "bobyqa" ,
optCtrl = list(maxfun= 100000)),
nAGQ = 20)
M1b.AQ.gamma <- glmer(AQ ~ ftz + GENDER + rater + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z + mat.age.z
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)),
family = Gamma(link=identity),
control = glmerControl(optimizer = "bobyqa" ,
optCtrl = list(maxfun= 100000)),
nAGQ = 20)
M1c.AQ.gamma <- glmer(AQ ~ ftz + GENDER + rater + pub.stage.z
+ ftz:GENDER
+ ftz:pub.stage.z + mat.age.z + GA.z
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)),
family = Gamma(link=identity),
control = glmerControl(optimizer = "bobyqa" ,
optCtrl = list(maxfun= 100000)),
nAGQ = 20)
stargazer(M1a.AQ.gamma, M1b.AQ.gamma, M1c.AQ.gamma, digits = 2, report = ('vcstp'), digit.separator = "", align = TRUE, type = "text")##
## ==============================================
## Dependent variable:
## ------------------------------
## AQ
## (1) (2) (3)
## ----------------------------------------------
## ftz -0.87 -0.79 -2.95
## (1.44) (1.60) (3.32)
## t = -0.60 t = -0.50 t = -0.89
## p = 0.55 p = 0.63 p = 0.38
##
## GENDER1 5.99 4.79 7.09
## (2.73) (2.96) (6.82)
## t = 2.19 t = 1.62 t = 1.04
## p = 0.03 p = 0.11 p = 0.30
##
## rater1 -4.67 -4.71 -4.01
## (0.69) (0.77) (0.98)
## t = -6.81 t = -6.11 t = -4.10
## p = 0.00 p = 0.00 p = 0.0001
##
## pub.stage.z 1.13 0.69 -1.66
## (1.64) (1.80) (2.70)
## t = 0.69 t = 0.38 t = -0.61
## p = 0.50 p = 0.71 p = 0.54
##
## mat.age.z -0.04 -0.11
## (0.14) (0.19)
## t = -0.29 t = -0.56
## p = 0.78 p = 0.58
##
## GA.z -1.52
## (0.59)
## t = -2.59
## p = 0.01
##
## ftz:GENDER1 1.75 1.96 9.69
## (3.03) (3.25) (6.56)
## t = 0.58 t = 0.60 t = 1.48
## p = 0.57 p = 0.55 p = 0.14
##
## ftz:pub.stage.z 2.86 3.44 5.98
## (2.09) (2.24) (3.25)
## t = 1.37 t = 1.54 t = 1.84
## p = 0.18 p = 0.13 p = 0.07
##
## Constant 13.10 13.24 10.83
## (1.30) (1.40) (3.15)
## t = 10.11 t = 9.46 t = 3.44
## p = 0.00 p = 0.00 p = 0.001
##
## ----------------------------------------------
## Observations 170 139 94
## ==============================================
## Note: *p<0.1; **p<0.05; ***p<0.01M2a.AQ.gamma <- glmer(AQ ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER
+ ftz:pub.timing
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)),
family = Gamma(link=identity),
control = glmerControl(optimizer = "bobyqa" ,
optCtrl = list(maxfun= 100000)),
nAGQ = 20)
M2b.AQ.gamma <- glmer(AQ ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER
+ ftz:pub.timing + mat.age.z
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)),
family = Gamma(link=identity),
control = glmerControl(optimizer = "bobyqa" ,
optCtrl = list(maxfun= 100000)),
nAGQ = 20)
M2c.AQ.gamma <- glmer(AQ ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER
+ ftz:pub.timing + mat.age.z + GA.z
+ (1 | ID),
data = d,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)),
family = Gamma(link=identity),
control = glmerControl(optimizer = "bobyqa" ,
optCtrl = list(maxfun= 100000)),
nAGQ = 20)
stargazer(M2a.AQ.gamma, M2b.AQ.gamma, M2c.AQ.gamma, digits = 2, report = ('vcstp'), digit.separator = "", align = TRUE, type = "text")##
## =============================================
## Dependent variable:
## ------------------------------
## AQ
## (1) (2) (3)
## ---------------------------------------------
## ftz -0.89 -0.85 -3.17
## (1.44) (1.59) (3.30)
## t = -0.62 t = -0.54 t = -0.96
## p = 0.54 p = 0.60 p = 0.34
##
## GENDER1 5.40 4.40 8.24
## (2.56) (2.76) (6.19)
## t = 2.11 t = 1.60 t = 1.33
## p = 0.04 p = 0.12 p = 0.19
##
## rater1 -4.66 -4.70 -4.05
## (0.68) (0.77) (0.97)
## t = -6.84 t = -6.12 t = -4.16
## p = 0.00 p = 0.00 p = 0.0001
##
## pub.timing 1.63 0.95 -1.18
## (1.98) (2.05) (2.85)
## t = 0.82 t = 0.46 t = -0.41
## p = 0.42 p = 0.65 p = 0.68
##
## mat.age.z -0.01 -0.08
## (0.14) (0.20)
## t = -0.09 t = -0.39
## p = 0.93 p = 0.71
##
## GA.z -1.51
## (0.58)
## t = -2.62
## p = 0.01
##
## ftz:GENDER1 0.17 0.04 6.42
## (2.88) (3.08) (6.60)
## t = 0.06 t = 0.01 t = 0.97
## p = 0.96 p = 0.99 p = 0.34
##
## ftz:pub.timing 3.71 4.54 6.61
## (2.32) (2.45) (3.36)
## t = 1.60 t = 1.85 t = 1.97
## p = 0.11 p = 0.07 p = 0.05
##
## Constant 13.06 13.21 10.85
## (1.29) (1.39) (3.17)
## t = 10.16 t = 9.51 t = 3.43
## p = 0.00 p = 0.00 p = 0.001
##
## ---------------------------------------------
## Observations 170 139 94
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01AQ subscale analysis
# Put all (standardized) outcome sub-scales together in DVlist
DVlist <- c("SOCIAL", "ATT_SWITCH", "COMMUNIC", "IMAGIN", "DETAIL")
X = length(DVlist)
contrasts(d$GENDER) <- c(-.5, .5)
contrasts(d$rater) <- c(-.5, .5)# pubertal stage model
M1a.outcomes.lme <- lapply(DVlist, function(x) {
lmer(eval(substitute(i ~ ftz + GENDER + rater + pub.stage.z
+ ftz:GENDER + ftz:pub.stage.z
+ (1|ID),
list(i = as.name(x)))),
data = d)})
M1b.outcomes.lme <- lapply(DVlist, function(x) {
lmer(eval(substitute(i ~ ftz + GENDER + rater + pub.stage.z
+ ftz:GENDER + ftz:pub.stage.z + mat.age.z
+ (1|ID),
list(i = as.name(x)))),
data = d)})
M1c.outcomes.lme <- lapply(DVlist, function(x) {
lmer(eval(substitute(i ~ ftz + GENDER + rater + pub.stage.z
+ ftz:GENDER + ftz:pub.stage.z + mat.age.z + GA.z
+ (1|ID),
list(i = as.name(x)))),
data = d)})
# pubertal timing models
M2a.outcomes.lme <- lapply(DVlist, function(x) {
lmer(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing
+ (1|ID),
list(i = as.name(x)))),
data = d)})
M2b.outcomes.lme <- lapply(DVlist, function(x) {
lmer(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing + mat.age.z
+ (1|ID),
list(i = as.name(x)))),
data = d)})
M2c.outcomes.lme <- lapply(DVlist, function(x) {
lmer(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing + mat.age.z + GA.z
+ (1|ID),
list(i = as.name(x)))),
data = d)})lpub <-"orange"
hpub <- "springgreen4"
DVlist.nice <- c("Social Skills", "Attention-Switching", "Communication", "Imagination", "Attention to Detail")
plot_list <- vector("list", length = length(DVlist.nice))
for (X in 1:length(DVlist.nice)) {
predicted <- effect(c("ftz", "pub.timing", "GENDER"), M2b.outcomes.lme[[X]],
xlevels= list(ftz = seq(from=-1, to=3.5, by=0.5),
pub.timing=seq(from=-.5, to=.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean)
pred.df <- as.data.frame(predicted)
pred.df$GENDER = relevel(pred.df$GENDER, "girl")
levels(pred.df$GENDER) <- c("Females", "Males")
pred.df$GENDER <- droplevels(pred.df$GENDER)
dv.col <- which(names(d)==DVlist[X])
p <- ggplot() +
geom_ribbon(data= pred.df, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.2) +
geom_line(data= pred.df,aes(x = ftz, y = fit, group = pub.timing, colour=pub.timing),size=1) +
labs(x = "", y = paste(DVlist.nice[X]),
title = paste(DVlist.nice[X])) +
scale_y_discrete(name = "", limits = c(0,2,4,6,8,10),
labels = c("0", "2", "4", "6", "8", "10"))+
coord_cartesian(ylim = c(0, 10)) +
scale_colour_gradient(low=lpub, high=hpub)+
scale_fill_gradient(low=lpub, high=hpub)+
theme_minimal()+
theme(legend.position = "none")+
facet_grid(~GENDER)
plot_list[[X]] <- p
}
gg<-ggarrange(plot_list[[1]], plot_list[[3]], nrow = 1, ncol=2)
annotate_figure(gg,
bottom = text_grob("Prenatal Testosterone (standardized)",
hjust=0.5, size=14))
DVlist.nice <- c("Social Skills", "Attention-Switching", "Communication", "Imagination", "Attention to Detail")
plot_list <- vector("list", length = length(DVlist.nice))
for (X in 1:length(DVlist.nice)) {
predicted <- effect(c("ftz", "pub.timing"), M2b.outcomes.lme[[X]],
xlevels= list(ftz = seq(from=-1, to=3.5, by=0.5),
pub.timing=seq(from=-.5, to=.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean)
pred.df <- as.data.frame(predicted)
p <- ggplot() +
geom_ribbon(data= pred.df, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.3) +
geom_line(data= pred.df,aes(x = ftz, y = fit, group = pub.timing, colour=pub.timing),size=1) +
labs(x = "", y = paste(DVlist.nice[X]),
title = paste(DVlist.nice[X]), fill = "Pubertal Timing", colour= "Pubertal Timing") +
scale_y_discrete(name = "", limits = c(0,2,4,6,8,10),
labels = c("0", "2", "4", "6", "8", "10"))+
scale_colour_gradient(low="orange", high="springgreen4")+
scale_fill_gradient(low="orange", high="springgreen4")+
coord_cartesian(ylim = c(1, 11)) +
theme_minimal()+
theme(legend.position = "none", plot.title = element_text(face = "bold"))
plot_list[[X]] <- p
}
gg_top<-ggarrange(plot_list[[1]], plot_list[[3]], nrow = 1, ncol=3)
gg_bottom<-ggarrange(plot_list[[2]], plot_list[[4]],plot_list[[5]], nrow = 1, ncol=3)
gg_all <- ggarrange(gg_top, gg_bottom, nrow=2, ncol=1)
annotate_figure(gg_all,
bottom = text_grob("Prenatal Testosterone (standardized)",
hjust=0.5, size=14))
Subscore Sensitivity Analysis 1: sex-stratify
contrasts(dboys$rater) <- c(-.5, .5)
contrasts(dgirls$rater) <- c(-.5, .5)
# NOTE: ft:rater is excluded
boyM2.lme <- lapply(DVlist, function(x) {
lmer(eval(substitute(i ~ ftz + rater + pub.timing
+ ftz:pub.timing + (1 | ID),
list(i = as.name(x)))),
data = dboys)})
girlM2.lme <- lapply(DVlist, function(x) {
lmer(eval(substitute(i ~ ftz + rater + pub.timing
+ ftz:pub.timing + (1 | ID),
list(i = as.name(x)))),
data = dgirls)})
stargazer(boyM2.lme[[1]], boyM2.lme[[2]], boyM2.lme[[3]], boyM2.lme[[4]], boyM2.lme[[5]], column.labels =
c("boy","boy", "boy", "boy", "boy"),
digits = 2, report=('vcsp'), digit.separator = "", order = c(), align = TRUE, type = "text")##
## ==================================================================
## Dependent variable:
## ----------------------------------------------
## SOCIAL ATT_SWITCH COMMUNIC IMAGIN DETAIL
## boy boy boy boy boy
## (1) (2) (3) (4) (5)
## ------------------------------------------------------------------
## ftz -0.17 -0.31 -0.37 -0.11 -0.50
## (0.35) (0.28) (0.29) (0.25) (0.29)
## p = 0.63 p = 0.27 p = 0.21 p = 0.67 p = 0.09
##
## rater1 -0.01 -1.08 -0.12 -0.02 -0.85
## (0.38) (0.39) (0.37) (0.37) (0.39)
## p = 0.99 p = 0.01 p = 0.75 p = 0.96 p = 0.04
##
## pub.timing -0.15 0.33 0.46 0.30 -1.59
## (1.31) (1.02) (1.07) (0.93) (1.09)
## p = 0.92 p = 0.75 p = 0.67 p = 0.76 p = 0.15
##
## ftz:pub.timing 3.62 0.81 1.85 1.01 1.46
## (1.51) (1.17) (1.23) (1.07) (1.25)
## p = 0.02 p = 0.49 p = 0.14 p = 0.35 p = 0.25
##
## Constant 2.82 4.08 2.80 2.94 4.79
## (0.38) (0.30) (0.31) (0.27) (0.32)
## p = 0.00 p = 0.00 p = 0.00 p = 0.00 p = 0.00
##
## ------------------------------------------------------------------
## Observations 96 96 96 96 96
## Log Likelihood -216.98 -205.85 -206.07 -199.64 -209.01
## Akaike Inf. Crit. 447.95 425.70 426.14 413.28 432.02
## Bayesian Inf. Crit. 465.90 443.65 444.09 431.23 449.97
## ==================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01stargazer(girlM2.lme[[1]], girlM2.lme[[2]], girlM2.lme[[3]], girlM2.lme[[4]], girlM2.lme[[5]], column.labels =
c("girl","girl", "girl", "girl", "girl"),
digits = 2, report=('vcsp'), digit.separator = "", order = c(), align = TRUE, type = "text")##
## ========================================================================
## Dependent variable:
## ----------------------------------------------------
## SOCIAL ATT_SWITCH COMMUNIC IMAGIN DETAIL
## girl girl girl girl girl
## (1) (2) (3) (4) (5)
## ------------------------------------------------------------------------
## ftz 0.68 -0.69 0.74 -0.15 -1.45
## (0.68) (0.67) (0.71) (0.54) (0.91)
## p = 0.32 p = 0.31 p = 0.30 p = 0.78 p = 0.12
##
## rater1 -1.00 -1.57 -0.67 -1.09 -1.35
## (0.32) (0.37) (0.35) (0.26) (0.39)
## p = 0.002 p = 0.0001 p = 0.06 p = 0.0001 p = 0.001
##
## pub.timing 0.75 -1.17 5.23 -2.84 -4.25
## (3.17) (3.17) (3.32) (2.52) (4.25)
## p = 0.82 p = 0.72 p = 0.12 p = 0.27 p = 0.32
##
## ftz:pub.timing 0.45 0.14 6.42 -2.43 -6.12
## (3.60) (3.59) (3.77) (2.86) (4.83)
## p = 0.91 p = 0.97 p = 0.09 p = 0.40 p = 0.21
##
## Constant 1.93 2.30 2.09 1.36 3.39
## (0.59) (0.59) (0.62) (0.47) (0.79)
## p = 0.002 p = 0.0001 p = 0.001 p = 0.004 p = 0.0001
##
## ------------------------------------------------------------------------
## Observations 74 74 74 74 74
## Log Likelihood -130.57 -135.47 -135.19 -116.03 -148.03
## Akaike Inf. Crit. 275.14 284.95 284.38 246.06 310.05
## Bayesian Inf. Crit. 291.27 301.07 300.50 262.19 326.18
## ========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01DVlist.nice <- c("Social Skills", "Attention-Switching", "Communication", "Imagination", "Attention to Detail")
boy_list <- vector("list", length = length(DVlist.nice))
girl_list <- vector("list", length = length(DVlist.nice))
for (X in 1:5) {
## BOYS ##
boypred <- as.data.frame(effect(c("ftz", "pub.timing"), boyM2.lme[[X]],
xlevels= list(ftz = seq(from=-1, to=3, by=0.25),
pub.timing=seq(from=-.5, to=.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean))
boyp <- ggplot() +
geom_ribbon(data= boypred, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.2) +
geom_line(data= boypred,aes(x = ftz, y = fit,
group = pub.timing, colour=pub.timing),size=1) +
labs(x = "", title=paste0(DVlist.nice[X])) +
scale_y_discrete(name = "", limits = c(1,3,5,7,9,11),
labels = c("0", "2", "4", "6", "8", "10"))+
coord_cartesian(ylim = c(1, 11)) +
scale_colour_gradient(low=lpub, high=hpub)+
scale_fill_gradient(low=lpub, high=hpub)+
theme_minimal()+
theme(legend.position = "none", plot.title = element_text(face = "bold"))
boy_list[[X]] <- boyp
## GIRLS ##
girlpred <- as.data.frame(effect(c("ftz", "pub.timing"), girlM2.lme[[X]],
xlevels= list(ftz = seq(from=-1, to=3.5, by=0.25),
pub.timing=seq(from=-.5, to=.5, by=0.5)),
se=TRUE, confidence.level = .95,
typical=mean))
girlp <- ggplot() +
geom_ribbon(data= girlpred, aes(x = ftz, ymin = fit - se, ymax = fit + se,
group=pub.timing, fill=pub.timing), alpha = 0.2) +
geom_line(data= girlpred,aes(x = ftz, y = fit,
group = pub.timing, colour=pub.timing),size=1) +
scale_y_discrete(name = "", limits = c(1,3,5,7,9,11),
labels = c("0", "2", "4", "6", "8", "10"))+
labs(x="")+
scale_colour_gradient(low=lpub, high=hpub)+
scale_fill_gradient(low=lpub, high=hpub)+
coord_cartesian(ylim = c(1, 11)) +
theme_minimal()+
theme(legend.position = "none",plot.title = element_text(face = "bold"))
girl_list[[X]] <- girlp
}
gboy<-ggarrange(boy_list[[5]], boy_list[[2]], boy_list[[4]], boy_list[[1]],boy_list[[3]],
nrow = 1, ncol=5, common.legend = T)
gboy<-annotate_figure(gboy,
left = text_grob("Males",
size=14, rot = 90))
ggirl<-ggarrange(girl_list[[5]],girl_list[[2]],girl_list[[4]],girl_list[[1]],girl_list[[3]],
nrow = 1, ncol=5, common.legend = T)
ggirl<-annotate_figure(ggirl,
left = text_grob("Females",
size=14, rot = 90))
gg <- ggarrange (gboy, ggirl, nrow=2, ncol=1, common.legend = T)
annotate_figure(gg,
bottom = text_grob("Foetal Testosterone (standardized)",
hjust=0.5, size=14))
Subscore Sensitivity Analysis 2: generalized linear mixed models
d$SOCIALplus = d$SOCIAL+1
d$ATT_SWITCHplus = d$ATT_SWITCH+1
d$COMMUNICplus = d$COMMUNIC+1
d$IMAGINplus = d$IMAGIN+1
d$DETAILplus = d$DETAIL+1
DVlist <- c("SOCIALplus") # GLMM only run for social skills as this was the only subscale significant in LMMs
# M2a.outcomes.glmm <- lapply(DVlist, function(x) {
# glmer(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
# + ftz:GENDER + ftz:pub.timing
# + (1|ID),
# list(i = as.name(x)))),
# data = d,
# family = Gamma(link=identity),
# control = glmerControl(optimizer = "bobyqa" ,
# optCtrl = list(maxfun= 100000)),
# nAGQ = 20)})
# did not converge
M2b.outcomes.glmm <- lapply(DVlist, function(x) {
glmer(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing + mat.age.z
+ (1|ID),
list(i = as.name(x)))),
data = d,
family = Gamma(link=identity),
control = glmerControl(optimizer = "bobyqa" ,
optCtrl = list(maxfun= 100000)),
nAGQ = 20)})
M2c.outcomes.glmm <- lapply(DVlist, function(x) {
glmer(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing + mat.age.z + GA.z
+ (1|ID),
list(i = as.name(x)))),
data = d,
family = Gamma(link=identity),
control = glmerControl(optimizer = "bobyqa" ,
optCtrl = list(maxfun= 100000)),
nAGQ = 20)})
# run M2a-c as glmmPQLs just to validate glmer findings
M2a.outcomes.pql <- lapply(DVlist, function(x) {
glmmPQL(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing,
list(i = as.name(x)))),
random = ~1|ID,
data = d,
family = Gamma(link="identity"),
verbose=0)})
M2b.outcomes.pql <- lapply(DVlist, function(x) {
glmmPQL(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing + mat.age.z,
list(i = as.name(x)))),
random = ~1|ID,
data = d,
family = Gamma(link="identity"),
verbose=0)})
M2c.outcomes.pql <- lapply(DVlist, function(x) {
glmmPQL(eval(substitute(i ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER + ftz:pub.timing + mat.age.z + GA.z,
list(i = as.name(x)))),
random = ~1|ID,
data = d,
family = Gamma(link="identity"),
verbose=0)})
# social skills
stargazer(M2b.outcomes.glmm[[1]], M2c.outcomes.glmm[[1]], report = ('vctsp'),
digits = 2, digit.separator = "", align = TRUE, type = "text")##
## ===========================================
## Dependent variable:
## ----------------------------
## SOCIALplus
## (1) (2)
## -------------------------------------------
## ftz 0.37 -0.20
## t = 0.74 t = -0.22
## (0.50) (0.93)
## p = 0.46 p = 0.83
##
## GENDER1 0.26 0.96
## t = 0.29 t = 0.55
## (0.88) (1.74)
## p = 0.78 p = 0.59
##
## rater1 -0.48 -0.21
## t = -2.39 t = -0.86
## (0.20) (0.24)
## p = 0.02 p = 0.40
##
## pub.timing 0.65 -0.65
## t = 1.09 t = -0.84
## (0.60) (0.77)
## p = 0.28 p = 0.40
##
## mat.age.z -0.001 0.04
## t = -0.03 t = 0.64
## (0.04) (0.06)
## p = 0.98 p = 0.53
##
## GA.z -0.12
## t = -0.71
## (0.18)
## p = 0.49
##
## ftz:GENDER1 -0.48 0.71
## t = -0.49 t = 0.38
## (0.98) (1.86)
## p = 0.63 p = 0.71
##
## ftz:pub.timing 1.52 3.01
## t = 2.15 t = 3.27
## (0.71) (0.92)
## p = 0.04 p = 0.002
##
## Constant 2.86 2.72
## t = 6.47 t = 3.05
## (0.44) (0.89)
## p = 0.00 p = 0.003
##
## -------------------------------------------
## Observations 139 94
## ===========================================
## Note: *p<0.1; **p<0.05; ***p<0.01Other plots
Again use model (b) to generate plots, compromise between control for confounds and sample size.
## Warning in stri_detect_regex(string, pattern, negate = negate, opts_regex =
## opts(pattern)): longer object length is not a multiple of shorter object length
## Warning in stri_detect_regex(string, pattern, negate = negate, opts_regex =
## opts(pattern)): longer object length is not a multiple of shorter object length
## Warning in stri_detect_regex(string, pattern, negate = negate, opts_regex =
## opts(pattern)): longer object length is not a multiple of shorter object length
## Warning in stri_detect_regex(string, pattern, negate = negate, opts_regex =
## opts(pattern)): longer object length is not a multiple of shorter object length
## Warning in stri_detect_regex(string, pattern, negate = negate, opts_regex =
## opts(pattern)): longer object length is not a multiple of shorter object length
Power Calculation
# note the rsim package will produce an error for every simulation (nsim) if there is missing data (https://github.com/pitakakariki/simr/issues/136)
tokeep <- c("ID", "GENDER", "rater", "ftz", "AQ", "pub.timing")
d_trim <- d[tokeep]
d_trim <- d_trim[!is.na(d_trim$AQ),]
d_trim <- d_trim[!is.na(d_trim$pub.timing),]
library(simr)
set.seed(123)
M3pow <- lmer(AQ ~ ftz + GENDER + rater + pub.timing
+ ftz:GENDER
+ ftz:pub.timing
+ (1 | ID),
data = d_trim,
contrasts=list(GENDER=c(-0.5,0.5), rater=c(-0.5,0.5)))
## power for main effect of pT
# Auyeung et al. (2009) showed that for every 1 nmol/L of FT, AQ scores went up 12pts (11.61-11.92). Thus, I would expect a 6 point (5.81-5.96) increase for every SD of FT (1 SD ~ 0.50 nmol/L). She also saw an ~11 pt increase in CAST scores (though they were transformed)
fixef(M3pow)["ftz"]<- 5
main <- powerSim(M3pow, fixed("ftz", "z"), nsim=100)## Simulating: | |Simulating: |= |Simulating: |== |Simulating: |=== |Simulating: |==== |Simulating: |===== |Simulating: |====== |Simulating: |======= |Simulating: |======== |Simulating: |========= |Simulating: |========== |Simulating: |=========== |Simulating: |============ |Simulating: |============= |Simulating: |============== |Simulating: |=============== |Simulating: |================ |Simulating: |================= |Simulating: |================== |Simulating: |=================== |Simulating: |==================== |Simulating: |===================== |Simulating: |====================== |Simulating: |======================= |Simulating: |======================== |Simulating: |========================= |Simulating: |========================== |Simulating: |=========================== |Simulating: |============================ |Simulating: |============================= |Simulating: |============================== |Simulating: |=============================== |Simulating: |================================ |Simulating: |================================= |Simulating: |================================== |Simulating: |=================================== |Simulating: |==================================== |Simulating: |===================================== |Simulating: |====================================== |Simulating: |======================================= |Simulating: |======================================== |Simulating: |========================================= |Simulating: |========================================== |Simulating: |=========================================== |Simulating: |============================================ |Simulating: |============================================= |Simulating: |============================================== |Simulating: |=============================================== |Simulating: |================================================ |Simulating: |================================================= |Simulating: |================================================== |Simulating: |=================================================== |Simulating: |==================================================== |Simulating: |===================================================== |Simulating: |====================================================== |Simulating: |======================================================= |Simulating: |======================================================== |Simulating: |========================================================= |Simulating: |========================================================== |Simulating: |=========================================================== |Simulating: |============================================================ |Simulating: |============================================================= |Simulating: |============================================================== |Simulating: |=============================================================== |Simulating: |================================================================ |Simulating: |================================================================= |Simulating: |==================================================================|main## Power for predictor 'ftz', (95% confidence interval):
## 88.00% (79.98, 93.64)
##
## Test: z-test
## Effect size for ftz is 5.0
##
## Based on 100 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 170
##
## Time elapsed: 0 h 0 m 6 s# To observe an effect of 6, we'd have a power of 98% (95% CI = 89-100)
# To observe an conservative effect of 5, we'd have a power of 82% (95% CI = 69-91)
## power for pT:sex interaction? try for various values
fixef(M3pow)["ftz:GENDER1"]<- 6
int1 <- powerSim(M3pow, fixed("ftz:GENDER1", "z"), nsim=500)## Simulating: | |Simulating: |= |Simulating: |== |Simulating: |=== |Simulating: |==== |Simulating: |===== |Simulating: |====== |Simulating: |======= |Simulating: |======== |Simulating: |========= |Simulating: |========== |Simulating: |=========== |Simulating: |============ |Simulating: |============= |Simulating: |============== |Simulating: |=============== |Simulating: |================ |Simulating: |================= |Simulating: |================== |Simulating: |=================== |Simulating: |==================== |Simulating: |===================== |Simulating: |====================== |Simulating: |======================= |Simulating: |======================== |Simulating: |========================= |Simulating: |========================== |Simulating: |=========================== |Simulating: |============================ |Simulating: |============================= |Simulating: |============================== |Simulating: |=============================== |Simulating: |================================ |Simulating: |================================= |Simulating: |================================== |Simulating: |=================================== |Simulating: |==================================== |Simulating: |===================================== |Simulating: |====================================== |Simulating: |======================================= |Simulating: |======================================== |Simulating: |========================================= |Simulating: |========================================== |Simulating: |=========================================== |Simulating: |============================================ |Simulating: |============================================= |Simulating: |============================================== |Simulating: |=============================================== |Simulating: |================================================ |Simulating: |================================================= |Simulating: |================================================== |Simulating: |=================================================== |Simulating: |==================================================== |Simulating: |===================================================== |Simulating: |====================================================== |Simulating: |======================================================= |Simulating: |======================================================== |Simulating: |========================================================= |Simulating: |========================================================== |Simulating: |=========================================================== |Simulating: |============================================================ |Simulating: |============================================================= |Simulating: |============================================================== |Simulating: |=============================================================== |Simulating: |================================================================ |Simulating: |================================================================= |Simulating: |==================================================================|int1## Power for predictor 'ftz:GENDER1', (95% confidence interval):
## 44.20% (39.79, 48.68)
##
## Test: z-test
## Effect size for ftz:GENDER1 is 6.0
##
## Based on 500 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 170
##
## Time elapsed: 0 h 0 m 27 s# power for effect=3 is 15% (95% CI: 12-18%)
# power for effect=5 is 31%
# power for effect=6 is 41%
# observed power is 5%
## power for pT:pub.timing interaction? try for various values
fixef(M3pow)["ftz:pub.timing"]<- 6
int2 <- powerSim(M3pow, fixed("ftz:pub.timing", "z"), nsim=200)## Simulating: | |Simulating: |= |Simulating: |== |Simulating: |=== |Simulating: |==== |Simulating: |===== |Simulating: |====== |Simulating: |======= |Simulating: |======== |Simulating: |========= |Simulating: |========== |Simulating: |=========== |Simulating: |============ |Simulating: |============= |Simulating: |============== |Simulating: |=============== |Simulating: |================ |Simulating: |================= |Simulating: |================== |Simulating: |=================== |Simulating: |==================== |Simulating: |===================== |Simulating: |====================== |Simulating: |======================= |Simulating: |======================== |Simulating: |========================= |Simulating: |========================== |Simulating: |=========================== |Simulating: |============================ |Simulating: |============================= |Simulating: |============================== |Simulating: |=============================== |Simulating: |================================ |Simulating: |================================= |Simulating: |================================== |Simulating: |=================================== |Simulating: |==================================== |Simulating: |===================================== |Simulating: |====================================== |Simulating: |======================================= |Simulating: |======================================== |Simulating: |========================================= |Simulating: |========================================== |Simulating: |=========================================== |Simulating: |============================================ |Simulating: |============================================= |Simulating: |============================================== |Simulating: |=============================================== |Simulating: |================================================ |Simulating: |================================================= |Simulating: |================================================== |Simulating: |=================================================== |Simulating: |==================================================== |Simulating: |===================================================== |Simulating: |====================================================== |Simulating: |======================================================= |Simulating: |======================================================== |Simulating: |========================================================= |Simulating: |========================================================== |Simulating: |=========================================================== |Simulating: |============================================================ |Simulating: |============================================================= |Simulating: |============================================================== |Simulating: |=============================================================== |Simulating: |================================================================ |Simulating: |================================================================= |Simulating: |==================================================================|int2## Power for predictor 'ftz:pub.timing', (95% confidence interval):
## 68.50% (61.57, 74.87)
##
## Test: z-test
## Effect size for ftz:pub.timing is 6.0
##
## Based on 200 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 170
##
## Time elapsed: 0 h 0 m 10 s# power for effect=3 is 20% (16-23%)
# power for effect=6 is 66%