R Notebook

Descriptive Statistics

Histograms of volumes and clusters

Boxplots sied by side of lobes

Boxplots side by side of vascular territories

library(ggplot2)
library(tidyverse)
library(viridis)
library(hrbrthemes)#for the violin plot
library(lme4)
library(reshape2)
library(splines)
library(modelr)
library(corrplot)
library(ggeffects)#ggpredict
#options(na.action = na.warn)
library(dplyr)

data <- read.csv("D:\\Dropbox/PhD/ukb/cross_merged_lobes.csv")

data$eid <- as.factor(data$eid)
data$sex <- as.factor(data$sex)


# descriptives
data$age_at_scan <-as.numeric(data$age_at_scan)
data$pvs_vol_mm3 <-as.numeric(data$pvs_vol_mm3)
data$overall_health_rating <-as.numeric(data$overall_health_rating)
data$overall_health_rating <- replace(data$overall_health_rating, which(data$overall_health_rating < 1), NA)
#data$overall_health_rating <- factor(data$overall_health_rating,levels=c(1,2,3,4),ordered=TRUE)
data$sex <- as.factor(data$sex)
data$diagnosed_diabetes <- factor(data$diagnosed_diabetes,c(0,1))
data$sex <- factor(data$sex,c(0,1))

summary(data$sex)
summary(data$age_at_scan)

data$typical_daily_alcohol <- replace(data$typical_daily_alcohol, which(data$typical_daily_alcohol < 0), NA)
data$typical_daily_alcohol <- factor(data$typical_daily_alcohol,c(1,2,3,4,5))
data$alcohol_frequency <- replace(data$alcohol_frequency, which(data$alcohol_frequency < 0), NA)
data$alcohol_frequency <- factor(data$alcohol_frequency)

summary(data$alcohol_frequency)
#hist(data$alcohol_frequency)


data$sleep_duration_sr <- replace(data$sleep_duration_sr, which(data$sleep_duration_sr < 0), NA)
summary(data$sleep_duration_sr)


data$freq_depressed_twoweeks <- replace(data$freq_depressed_twoweeks, which(data$freq_depressed_twoweeks < 0), NA)

summary(data$freq_depressed_twoweeks)
summary(data$height)
summary(data$weight)
summary(data$BMI)
summary(data$townsend_deprivation_index)

data$snoring_sr <- replace(data$snoring_sr, which(data$snoring_sr < 0), NA)
data$insomnia_sr <- replace(data$insomnia_sr, which(data$insomnia_sr < 0), NA)
data$difficulty_getting_up_morning <- replace(data$difficulty_getting_up_morning, which(data$difficulty_getting_up_morning < 0), NA)
data$avg_total_household_income <- replace(data$avg_total_household_income, which(data$avg_total_household_income < 0), NA)

summary(data$snoring_sr)
summary(data$insomnia_sr)
summary(data$difficulty_getting_up_morning)
summary(data$avg_total_household_income)
summary(data$ipaq_activity_group)

data$snoring_sr <- factor(data$snoring_sr)
data$insomnia_sr <- factor(data$insomnia_sr)
data$ipaq_activity_group <- factor(data$ipaq_activity_group)
data$avg_total_household_income <- factor(data$avg_total_household_income)

summary(data$age_at_scan)
summary(data$sex)
summary(data$pvs_vol_mm3)
summary(data$pvs_clusters)
summary(data$age_at_scan)

# create new vars

data$front_pvs_vol <- data$frontal_r + data$frontal_l
data$par_pvs_vol <- data$parietal_r + data$parietal_l
data$temp_pvs_vol <- data$temporal_r + data$temporal_l
data$occ_pvs_vol <- data$occipital_r + data$occipital_l
data$cerebel_pvs_vol <- data$cerebellum_r + data$cerebellum_l
#data$brainstem

data$anterior_pvs_vol <- data$ACAL + data$ACAR + data$MLSL + data$MLSR
data$middle_pvs_vol <- data$LLSL+data$LLSR+data$MCAFL+data$MCAFR+
                        data$MCAPL+data$MCAPR+data$MCATL+data$MCAOR+
                        data$MCAOL+data$MCAOR.1+data$MCAIL+data$MCAIR
data$post_pvs_vol <- data$PCATL+data$PCATR+data$PCAOL+data$PCAOR+
                      data$PCTPL+data$PCTPR+data$ACTPL+data$ACTPR
data$vb_pvs_vol <- data$BL+data$BR+data$SCL+data$SCR+data$ICL+data$ICR
data$vents_pvs_vol <- data$LVL+data$LVR

data$ACA_pvs_vol <- data$ACAL+data$ACAR
data$LLS_pvs_vol <- data$LLSL + data$LLSR
data$MCAF_pvs_vol <- data$MCAFL + data$MCAFR
data$MCAP_pvs_vol <- data$MCAPL + data$MCAPR

data$MCAT_pvs_vol <- data$MCATL+data$MCAOR
data$MCA0_pvs_vol <- data$MCAOL+data$MCAOR.1

# boxplots
# format and add in descriptive variables: median, mean, sd, min , max
boxplot(data$pvs_vol_mm3,
         main = "Boxplots",
         ylab = bquote("PVS Volume"~(mm^3)),
         xlab =c("PVS Volume"))

boxplot(data$pvs_clusters,
         main = "Boxplots",
         ylab = bquote("PVS Clusters"),
         xlab =c("PVS Clusters"))

# histograms
vol_density <- hist(data$pvs_vol_mm3,breaks=50,col="#87CEEB",xlab=bquote("PVS Volume"~(mm^3)),main="",cex.lab=1.3)
clu_density <- hist(data$pvs_clusters, breaks=50, col="#87CEEB",xlab="PVS Clusters",main="",cex.lab=1.3)

summary(data$pvs_vol_mm3)
summary(data$pvs_clusters)

boxplot(data$front_pvs_vol,data$par_pvs_vol,data$temp_pvs_vol,
        data$occ_pvs_vol,data$cerebel_pvs_vol,
         main = "Lobar Distributions",
         ylab = bquote("PVS Volume"~(mm^3)),
         names =c("Frontal", "Parietal", "Temporal", "Occipital Lobe","Cerebellar"))

# violin plot for wm lobes
subset_data_lobes <- data[c("front_pvs_vol","par_pvs_vol","temp_pvs_vol")]#,"occ_pvs_vol","cerebel_pvs_vol")]
meltData_lobes <- melt(subset_data_lobes)
boxplot(data=meltData_lobes, value~variable)

# Violin basic
meltData_lobes %>%
  ggplot( aes(factor(variable), value, fill=variable),) +
    geom_violin() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    #theme_ipsum() +
    theme(
      legend.position="none",
      plot.title = element_text(size=15,hjust=0.5),
      text = element_text(size = 15)
      ) +
  ggtitle("Violin chart of Lobar Distributions") + 
  scale_x_discrete(breaks=c("front_pvs_vol","par_pvs_vol","temp_pvs_vol"),
                   labels=c("Frontal", "Parietal", "Temporal")) +
  ylab(bquote("PVS Volume"~(mm^3))) + xlab("")
  #geom_jitter(color="black", size=0.01, alpha=0.1) +
    


boxplot(data$anterior_pvs_vol,data$middle_pvs_vol,data$post_pvs_vol,
        data$vb_pvs_vol,data$vents_pvs_vol,
         main = "PVS Distribution by Vascular Territory",
         ylab = bquote("PVS Volume"~(mm^3)),
         names =c("ACA", "MCA", "PCA", "VB","LV"))
#VB=Vertebro-Basilar
#LV=Lateral Ventricular

# breakdown of ACA  territory
boxplot(data$ACAL+data$ACAR,data$MLSL+data$MLSR,
         main = "Subdivisions of the Anterior Territory",
         ylab = bquote("PVS Volume"~(mm^3)),
         names =c("ACA","Medial Lenticulostriate"))

# breakdown of MCA territory
boxplot(data$LLSL+data$LLSR,data$MCAFL+data$MCAFR,
        data$MCAPL+data$MCAPR,data$MCATL+data$MCAOR,
        data$MCAOL+data$MCAOR.1,data$MCAIL+data$MCAIR,
         main = "Subdivisions of the Middle Territory",
         ylab = bquote("PVS Volume"~(mm^3)),
         names =c("LLS","Frontal Pars","Parietal Pars","Temporal Pars","Occipital Pars","Insular Pars"),
        cex.axis=0.7)
       

#mca_df <- subset(data, select=c('LLS_pvs_vol','MCAF_pvs_vol','MCAP_pvs_vol'))
#mca_df <- melt(data)
# creating the modified dataframe
#data_mod < - melt(data_frame, id.vars='eid', 
#                  measure.vars=c('col2', 'col3', 'col4'))


subset_data <- data[c("ACA_pvs_vol","LLS_pvs_vol","MCAF_pvs_vol","MCAP_pvs_vol",'MCAT_pvs_vol','MCA0_pvs_vol')]
meltData <- melt(subset_data)
boxplot(data=meltData, value~variable)
#coudl i also use gather?
#https://www.statology.org/gather-function-in-r/

#p <- ggplot(meltData, aes(factor(variable), value)) 
#p + geom_boxplot() + facet_wrap(~variable, scale="free")

# Violin basic
meltData %>%
  ggplot( aes(factor(variable), value, fill=variable)) +
    geom_violin() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    #theme_ipsum() +
    theme(
      legend.position="none",
      plot.title = element_text(size=14,hjust=0.5),
      text = element_text(size = 11)
    ) +
    ggtitle("PVS Volumes within Vascular Sub-Divisions") +
  scale_x_discrete(breaks=c("ACA_pvs_vol","LLS_pvs_vol","MCAF_pvs_vol",'MCAP_pvs_vol','MCAT_pvs_vol','MCA0_pvs_vol'),
                   labels=c("ACA", "LLS", "MCAF",'MCAP','MCAT','MCAO'))+
  #geom_jitter(color="black", size=0.01, alpha=0.1) ++
  ylab(bquote("PVS Volume"~(mm^3))) + xlab("")



# Focus on the ACA in anteiror territory
#data$ACAL+data$ACAR
# Focus on LL, Frontal pars and parietal pars in Middle territory
#data$LLSL+data$LLSR, data$MCAFL+data$MCAFR, data$MCAPL+data$MCAPR

Assocaitions with confounds

gamlj::gamljGlm( formula = LLSL ~ sex + age_at_scan + height + weight + I(age_at_scan^2) + height:age_at_scan + weight:age_at_scan, data = data, plotHAxis = age_at_scan, plotSepLines = weight, scaling=c(age_at_scan = “none”, height = “centered”, weight = “centered”))

modelling

# same with lm as glm
#model_break1 <- lm(pvs_vol_mm3 ~ age_at_scan + I(age_at_scan^2) + sex + sleep_duration_sr , data=data)

model <- lm(pvs_vol_mm3 ~ ns(age_at_scan,3) + sex + ns(pollution_score,3) +
              overall_health_rating + 
              overall_health_rating:age_at_scan +
              alcohol_frequency+ns(townsend_deprivation_index,1)+
              freq_depressed_twoweeks+ns(height,3)+ns(weight,1)+# ns(BMI,1)+
              diagnosed_diabetes+diagnosed_diabetes:age_at_scan+
              insomnia_sr+insomnia_sr:age_at_scan+snoring_sr+snoring_sr:age_at_scan+
              avg_total_household_income+ipaq_activity_group,data=data)
#INCLUDE GRAPH
age_vol <- ggpredict(model, terms=c('age_at_scan [all]','sex'))
age_vol$group <- str_replace(age_vol$group, "0", "Female")
age_vol$group <- str_replace(age_vol$group, "1", "Male")
age_plot <- plot(age_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

#INCLUDE GRAPH
height_vol <- ggpredict(model, terms=c('height [all]'))
height_plot <- plot(height_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "Height (cm)")+ theme(text = element_text(size = 15))
#INCLUDE GRAPH
weight_vol <- ggpredict(model, terms=c('weight [all]'))
weight_plot <- plot(weight_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "Weight (kg)")+ theme(text = element_text(size = 15))

#bmi_vol <- ggpredict(model, terms=c('BMI [all]'))
#bmi_plot <- plot(bmi_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
#  labs(y= bquote("PVS Volume"~(mm^3)), x =bquote("BMI"~(kg^2/m)))+ theme(text = element_text(size = 15))

#INCLUDE GRAPH
alco_vol <- ggpredict(model, terms=c('alcohol_frequency'))
#alco_vol$group <- str_replace_all(alco_vol$x, c(`0`="`Never`",`1`='`"0-1 per Month"`',`3`='`"2-3 per Week"`',`4`='`"4+ per Week"`',`2`='`"2-4 per Month"`'))

alco_vol$x <- recode_factor(alco_vol$x, `0` = "Never", 
                                `1` = "0-1 per Month", `2`='2-4 per Month',`3`='2-3 per Week',`4`='4+ per Week')

alco_plot <- plot(alco_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=FALSE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x ='Frequency of Alcohol Consumption')+ theme(text = element_text(size = 15))

#INCLUDE GRAPH
ipaq_vol <- ggpredict(model, terms=c('ipaq_activity_group'))
ipaq_plot <- plot(ipaq_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "IPAQ Activity Group")+ theme(text = element_text(size = 15))

insom_vol <- ggpredict(model, terms=c("age_at_scan [all]",'insomnia_sr'))
insom_plot <- plot(insom_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "age")+ theme(text = element_text(size = 15))

depres_vol <- ggpredict(model, terms=c("age_at_scan [all]",'freq_depressed_twoweeks'))
depres_plot <- plot(depres_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "age")+ theme(text = element_text(size = 15))

heath_vol <- ggpredict(model, terms=c("age_at_scan [all]",'overall_health_rating'))
health_plot <- plot(heath_vol, ci=FALSE,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "age")+ theme(text = element_text(size = 15))

diab_vol <- ggpredict(model, terms=c("age_at_scan [all]",'diagnosed_diabetes'))
diab_plot <- plot(diab_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "age")+ theme(text = element_text(size = 15))

income_vol <- ggpredict(model, terms=c('avg_total_household_income'))
income_plot <- plot(income_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "income group")+ theme(text = element_text(size = 15))

ipaq_vol <- ggpredict(model, terms=c("age_at_scan [all]",'ipaq_activity_group'))
ipaq_plot <- plot(ipaq_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))
income_vol <- ggpredict(model, terms=c("age_at_scan [all]",'avg_total_household_income'))
income_plot <- plot(income_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

#results:
#age2,sex,poll^2,daily alcohol^1,freq_depressed_twoweeks,height^3,weight^, diabetea,insomina^, age:diabetes
# not sig: health rating,twonsend,sleep duration, snoring, income, geeting up in mroonig, ipaq, 
# not sig interctions with age: mornnig get up, income, ipaq, insomnia, snoring, health rating

#####clu
##model_clu <- lm(pvs_clusters ~ ns(age_at_scan,3) + sex + ns(pollution_score,3) + 
##              overall_health_rating + 
##              overall_health_rating:age_at_scan +
##              ns(typical_daily_alcohol,3)+ns(townsend_deprivation_index,1)+
##              freq_depressed_twoweeks+ns(height,3)+ns(weight,1)+#ns(BMI,3)+
##              diagnosed_diabetes+diagnosed_diabetes:age_at_scan+
##              insomnia_sr+insomnia_sr:age_at_scan+snoring_sr+
##              avg_total_household_income+
##              ipaq_activity_group,data=data)
# not sig:sleep dur, health rating, alcohol (contradict), townsend, snoring, income
#sig: diabetes, age, insomnia (contradict), getting up^, 
# sig interaction with age: getting up^(cont vol), insomnia (contradict vol), diabetes



model_aca <- lm(ACA_pvs_vol ~ ns(age_at_scan,3) + sex + sex:age_at_scan + ns(pollution_score,3) + 
              ns(sleep_duration_sr,1) + overall_health_rating + 
              overall_health_rating:age_at_scan + ns(hrs, 1) + 
              ns(typical_daily_alcohol,1)+ns(townsend_deprivation_index,1)+
              freq_depressed_twoweeks+ns(height,3)+ns(weight,1)+#ns(BMI,3)+
              diagnosed_diabetes+diagnosed_diabetes:age_at_scan+
              insomnia_sr+insomnia_sr:age_at_scan+snoring_sr+snoring_sr:age_at_scan+
              avg_total_household_income+avg_total_household_income:age_at_scan+
              difficulty_getting_up_morning+difficulty_getting_up_morning:age_at_scan,data=data)

model_lls <- lm(LLS_pvs_vol ~ ns(age_at_scan,3) + sex + sex:age_at_scan + ns(pollution_score,3) + 
              ns(sleep_duration_sr,1) + overall_health_rating + 
              overall_health_rating:age_at_scan + ns(hrs, 1) + 
              ns(typical_daily_alcohol,1)+ns(townsend_deprivation_index,1)+
              freq_depressed_twoweeks+ns(height,3)+ns(weight,1)+#ns(BMI,3)+
              diagnosed_diabetes+diagnosed_diabetes:age_at_scan+
              insomnia_sr+insomnia_sr:age_at_scan+snoring_sr+snoring_sr:age_at_scan+
              avg_total_household_income+avg_total_household_income:age_at_scan+
              difficulty_getting_up_morning+difficulty_getting_up_morning:age_at_scan,data=data)

#without colineraities between time of day, age, insominia, difficulty getting up, typica daily alcohol, weight
model_hrs_vol<- lm(pvs_vol_mm3 ~ sex +
              ns(sleep_duration_sr,2) + ns(hrs, 2) +sex:sleep_duration_sr,data=data)
summary(model_hrs_vol)

Call:
lm(formula = pvs_vol_mm3 ~ sex + ns(sleep_duration_sr, 2) + ns(hrs, 
    2) + sex:sleep_duration_sr, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-2412.1  -893.4  -281.7   600.4  8223.4 

Coefficients: (1 not defined because of singularities)
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                2085.72     123.93  16.829  < 2e-16 ***
sex1                        201.30     109.00   1.847  0.06478 .  
ns(sleep_duration_sr, 2)1   609.78     189.06   3.225  0.00126 ** 
ns(sleep_duration_sr, 2)2  -172.89     169.99  -1.017  0.30915    
ns(hrs, 2)1                 218.13      47.06   4.635 3.58e-06 ***
ns(hrs, 2)2                 -67.22      30.62  -2.195  0.02816 *  
sex0:sleep_duration_sr      -36.69      15.08  -2.434  0.01496 *  
sex1:sleep_duration_sr          NA         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1267 on 31173 degrees of freedom
  (70 observations deleted due to missingness)
Multiple R-squared:  0.03305,   Adjusted R-squared:  0.03286 
F-statistic: 177.6 on 6 and 31173 DF,  p-value: < 2.2e-16

sleep_vol <- ggpredict(model_hrs_vol, terms=c("sleep_duration_sr [all]",'sex'))

Warning: prediction from a rank-deficient fit may be misleading

sleep_plot <- plot(sleep_vol, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "Sleep Duration (Hours)")+ theme(text = element_text(size = 15))

model_alco_vol<- lm(pvs_vol_mm3 ~ sex + ns(age_at_scan,3)+ns(typical_daily_alcohol,1)+typical_daily_alcohol:age_at_scan,data=data)
summary(model_alco_vol)

Call:
lm(formula = pvs_vol_mm3 ~ sex + ns(age_at_scan, 3) + ns(typical_daily_alcohol, 
    1) + typical_daily_alcohol:age_at_scan, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-2515.1  -859.6  -270.8   574.1  8280.2 

Coefficients: (1 not defined because of singularities)
                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                         922.377    247.689   3.724 0.000197 ***
sex1                                413.066     18.712  22.075  < 2e-16 ***
ns(age_at_scan, 3)1                 832.775     83.234  10.005  < 2e-16 ***
ns(age_at_scan, 3)2                1495.096    212.267   7.043 1.94e-12 ***
ns(age_at_scan, 3)3                 481.330    131.136   3.670 0.000243 ***
ns(typical_daily_alcohol, 1)        806.890    366.054   2.204 0.027516 *  
typical_daily_alcohol1:age_at_scan    9.555      4.674   2.044 0.040944 *  
typical_daily_alcohol2:age_at_scan    7.609      3.602   2.112 0.034664 *  
typical_daily_alcohol3:age_at_scan    4.958      2.533   1.957 0.050309 .  
typical_daily_alcohol4:age_at_scan    3.143      1.559   2.016 0.043843 *  
typical_daily_alcohol5:age_at_scan       NA         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1249 on 19340 degrees of freedom
  (11900 observations deleted due to missingness)
Multiple R-squared:  0.07539,   Adjusted R-squared:  0.07496 
F-statistic: 175.2 on 9 and 19340 DF,  p-value: < 2.2e-16

alco <- ggpredict(model_alco_vol, terms=c("age_at_scan [all]",'typical_daily_alcohol'))

Warning: prediction from a rank-deficient fit may be misleading

alco_plot <- plot(alco, ci=95,add.data=FALSE, show.title=FALSE,show.legend=TRUE)+
  labs(y= bquote("PVS Volume"~(mm^3)), x = "Age")+ theme(text = element_text(size = 15))


### plots
###vol <- ggpredict(model, terms=c("age_at_scan [all]"))
###vol_plot <- plot(vol,ci=95,add.data=FALSE,show.title=FALSE)+
###  labs(y= bquote("PVS Volume"~(mm^3)), x = "Age")+ theme(text = element_text(size = 15))

###time <- ggpredict(model, terms=c("hrs [all]"))
###time_plot <- plot(time,ci=95,add.data=FALSE,show.title=FALSE)+
###  labs(y= bquote("PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

# mean ci plots (dot plots with bar for each of the 4 health groups)
###dx_vol <- ggpredict(model, terms=c("overall_health_rating [all]"))
###dx_vol_plot <- plot(dx_vol,ci=95,add.data=FALSE,show.title=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453'))+
###  labs(y= bquote("PVS Volume"~(mm^3)), x = "Health Rating")+ theme(text = element_text(size = 15))+
###  geom_line(aes(group = "overall_health_rating"), linetype = "dashed")

# mean ci plots (dot plots with bar for each of the pollution)
###pollution_pred <- ggpredict(model, terms=c("pollution_score"))
###pollution_plot <- plot(pollution_pred,ci=95,add.data=FALSE,show.title=FALSE)+
###  labs(y= bquote("PVS Volume"~(mm^3)), x = "Pollution Level")+ theme(text = element_text(size = 15))+
###  geom_line(aes(group = "pollution_score"), linetype = "dashed")


model1 <- glm(ACA_pvs_vol ~ age_at_scan + I(age_at_scan^2) + sex, data=data)
model2 <- glm(LLS_pvs_vol ~ age_at_scan + I(age_at_scan^2) + sex, data=data)
model3 <- glm(MCAF_pvs_vol ~ age_at_scan + I(age_at_scan^2) + sex, data=data)
model4 <- glm(MCAP_pvs_vol ~ age_at_scan + I(age_at_scan^2) + sex, data=data)

summary(model1)

Call:
glm(formula = ACA_pvs_vol ~ age_at_scan + I(age_at_scan^2) + 
    sex, data = data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-989.8  -362.4  -114.7   245.0  3297.9  

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -2.193e+03  1.972e+02 -11.125   <2e-16 ***
age_at_scan       7.626e+01  6.144e+00  12.413   <2e-16 ***
I(age_at_scan^2) -4.717e-01  4.742e-02  -9.948   <2e-16 ***
sex1              1.246e+02  5.841e+00  21.329   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 263121)

    Null deviance: 8830994910  on 31249  degrees of freedom
Residual deviance: 8221477571  on 31246  degrees of freedom
AIC: 478701

Number of Fisher Scoring iterations: 2

summary(model2)

Call:
glm(formula = LLS_pvs_vol ~ age_at_scan + I(age_at_scan^2) + 
    sex, data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-489.72   -92.36   -10.46    80.39   919.65  

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      669.84587   52.34623  12.796  < 2e-16 ***
age_at_scan       -4.62955    1.63122  -2.838  0.00454 ** 
I(age_at_scan^2)   0.03304    0.01259   2.625  0.00868 ** 
sex1              45.19143    1.55088  29.139  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 18547.86)

    Null deviance: 595572578  on 31249  degrees of freedom
Residual deviance: 579546427  on 31246  degrees of freedom
AIC: 395818

Number of Fisher Scoring iterations: 2

summary(model3)

Call:
glm(formula = MCAF_pvs_vol ~ age_at_scan + I(age_at_scan^2) + 
    sex, data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-459.19  -181.68   -66.59   111.82  2049.99  

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -1148.3412   102.7209 -11.179   <2e-16 ***
age_at_scan         38.3972     3.2010  11.995   <2e-16 ***
I(age_at_scan^2)    -0.2405     0.0247  -9.737   <2e-16 ***
sex1                76.2932     3.0434  25.069   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 71423.46)

    Null deviance: 2391628056  on 31249  degrees of freedom
Residual deviance: 2231697452  on 31246  degrees of freedom
AIC: 437952

Number of Fisher Scoring iterations: 2

summary(model4)

Call:
glm(formula = MCAP_pvs_vol ~ age_at_scan + I(age_at_scan^2) + 
    sex, data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-596.38  -226.84   -67.02   157.12  2363.68  

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -1.762e+03  1.223e+02  -14.41   <2e-16 ***
age_at_scan       6.526e+01  3.812e+00   17.12   <2e-16 ***
I(age_at_scan^2) -4.696e-01  2.942e-02  -15.96   <2e-16 ***
sex1              9.652e+01  3.624e+00   26.63   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 101265.3)

    Null deviance: 3306020500  on 31249  degrees of freedom
Residual deviance: 3164136448  on 31246  degrees of freedom
AIC: 448862

Number of Fisher Scoring iterations: 2

mod_vol <- lm(pvs_vol_mm3 ~ ns(age_at_scan, 3), data = data)
mod_aca <- lm(ACA_pvs_vol ~ ns(age_at_scan, 3), data = data)
mod_lls <- lm(LLS_pvs_vol ~ ns(age_at_scan, 3), data = data)
mod_mcaf <- lm(MCAF_pvs_vol ~ ns(age_at_scan, 3), data = data)
mod_mcap <- lm(MCAP_pvs_vol ~ ns(age_at_scan, 3), data = data)


grid <- data %>% 
  data_grid(age_at_scan = seq_range(age_at_scan, n = 50, expand = 0.1)) %>% 
  gather_predictions(mod_vol,mod_aca,mod_lls,mod_mcaf,mod_mcap)#, .pred = "pvs_vol_mm3")

Error in data_grid(., age_at_scan = seq_range(age_at_scan, n = 50, expand = 0.1)) : 
  unused argument (age_at_scan = seq_range(age_at_scan, n = 50, expand = 0.1))

# i hate r
# the errors dont fking make sense
model1 <- glm(ACA_pvs_vol ~ age_at_scan + I(age_at_scan^2) + sex, data=data)
model2 <- glm(LLS_pvs_vol ~ age_at_scan + I(age_at_scan^2) + sex, data=data)
model3 <- glm(MCAF_pvs_vol ~ age_at_scan + I(age_at_scan^2) + sex, data=data)
model4 <- glm(MCAP_pvs_vol ~ age_at_scan + I(age_at_scan^2) + sex, data=data)

summary(model1)
summary(model2)
summary(model3)
summary(model4)

mod_vol <- lm(pvs_vol_mm3 ~ ns(age_at_scan, 3), data = data)
mod_aca <- lm(ACA_pvs_vol ~ ns(age_at_scan, 3), data = data)
mod_lls <- lm(LLS_pvs_vol ~ ns(age_at_scan, 3), data = data)
mod_mcaf <- lm(MCAF_pvs_vol ~ ns(age_at_scan, 3), data = data)
mod_mcap <- lm(MCAP_pvs_vol ~ ns(age_at_scan, 3), data = data)

grid <- data %>% 
  data_grid(age_at_scan = seq_range(age_at_scan, n = 50, expand = 0.1)) %>% 
  gather_predictions(mod_vol,mod_aca,mod_lls,mod_mcaf,mod_mcap, .pred = "pvs_vol_mm3")

ggplot(data, aes(age_at_scan, pred)) + 
  #geom_point() +
  #geom_line(data = subset(grid,model=="mod_vol"),colour = "red") +
  geom_line(data = subset(grid,model=="mod_aca"),colour = "blue")+
  #geom_line(data = subset(grid,model=="mod_lls"),colour = "purple")+
  geom_line(data = subset(grid,model=="mod_mcaf"),colour = "green")+
  geom_line(data = subset(grid,model=="mod_mcap"),colour = "orange")# +
  #facet_wrap(~ model)
ggplot(data, aes(age_at_scan, pred)) + 
  geom_line(data = subset(grid,model=="mod_lls"),colour = "purple")

LLS Model and Plots Below

lls_model <- lm(LLS_pvs_vol ~ ns(age_at_scan,2) + sex + ns(pollution_score,3) + 
              ns(sleep_duration_sr,1) + overall_health_rating + 
              overall_health_rating:age_at_scan + ns(hrs, 1) + overall_health_rating:hrs + 
                sleep_duration_sr:hrs + diagnosed_diabetes + diagnosed_diabetes:age_at_scan, data=data)

lls_vol <- ggpredict(lls_model, terms=c("age_at_scan [all]"))
lls_vol_plot <- plot(lls_vol,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("LLS PVS Volume"~(mm^3)), x = "Age")+ theme(text = element_text(size = 15))

lls_time <- ggpredict(lls_model, terms=c("hrs [all]"))
lls_time_plot <- plot(lls_time,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("LLS PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

# mean ci plots (dot plots with bar for each of the 4 health groups)
lls_dx_vol <- ggpredict(lls_model, terms=c("overall_health_rating [all]"))
lls_dx_vol_plot <- plot(lls_dx_vol,ci=95,add.data=FALSE,show.title=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453'))+
  labs(y= bquote("LLS PVS Volume"~(mm^3)), x = "Health Rating")+ theme(text = element_text(size = 15))+
  geom_line(aes(group = "overall_health_rating"), linetype = "dashed")

# mean ci plots (dot plots with bar for each of the pollution)
lls_pollution_pred <- ggpredict(lls_model, terms=c("pollution_score"))
lls_pollution_plot <- plot(lls_pollution_pred,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("LLS PVS Volume"~(mm^3)), x = "Pollution Level")+ theme(text = element_text(size = 15))+
  geom_line(aes(group = "pollution_score"), linetype = "dashed")


# over time plots
lls_prediction_a_vol <- ggpredict(lls_model, terms=c("age_at_scan [all]","overall_health_rating"))
lls_prediction_a_vol_plot<-plot(lls_prediction_a_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453', 'pink'),line.size =1.2,show.legend=TRUE,show.title=FALSE)+
  labs(y= bquote("LLS PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))


lls_prediction_time_vol <- ggpredict(lls_model, terms=c("hrs [all]","overall_health_rating"))
lls_prediction_time_vol_plot<-plot(lls_prediction_time_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453', 'pink'),line.size =1.2,show.legend=FALSE,show.title=FALSE)+
  labs(y= bquote("LLS PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

lls_sleep <- ggpredict(lls_model, terms=c("hrs [all]","sleep_duration_sr"))
lls_prediction_sleep_vol_plot<-plot(lls_sleep,ci=TRUE,add.data=FALSE,line.size =1.2,show.legend=TRUE,show.title=FALSE)+
  labs(y= bquote("LLS PVS Volume"~(mm^3)), x = "Time")+ theme(text = element_text(size = 15))


# same as before.... lls plot against health rating
# no inflection point as i saw in jamovi
simple_lls_model <- lm(LLS_pvs_vol ~ ns(age_at_scan,2) + overall_health_rating + overall_health_rating:age_at_scan, data=data)

simp_lls_prediction_a_vol <- ggpredict(simple_lls_model, terms=c("age_at_scan [all]","overall_health_rating"))
simp_lls_prediction_a_vol_plot <-plot(simp_lls_prediction_a_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453','pink'),
                                       line.size=1.2,show.legend=FALSE,show.title=FALSE)+
                                    labs(y= bquote("LLS PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

# look at left only
left_lls_model <- lm(LLSL ~ ns(age_at_scan,2) + overall_health_rating + 
                       overall_health_rating:age_at_scan +sex+ 
                       diagnosed_diabetes + diagnosed_diabetes:age_at_scan, data=data)

left_lls_prediction_a_vol <- ggpredict(left_lls_model, terms=c("age_at_scan [all]","overall_health_rating"))
left_lls_prediction_a_vol_plot <-plot(left_lls_prediction_a_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453','pink'),
                                       line.size=1.2,show.legend=FALSE,show.title=FALSE)+
                                    labs(y= bquote("Left LLS PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

# look at right only
right_lls_model <- lm(LLSR ~ ns(age_at_scan,2) + overall_health_rating + 
                        overall_health_rating:age_at_scan +sex+ 
                        diagnosed_diabetes + diagnosed_diabetes:age_at_scan, data=data)

right_lls_prediction_a_vol <- ggpredict(right_lls_model, terms=c("age_at_scan [all]","overall_health_rating"))
right_lls_prediction_a_vol_plot <-plot(right_lls_prediction_a_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453','pink'),
                                       line.size=1.2,show.legend=FALSE,show.title=FALSE)+
                                    labs(y= bquote("Right LLS PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))


# LLS and diabetes

llls_pred_agediab <- ggpredict(left_lls_model, terms=c("age_at_scan [all]","diagnosed_diabetes"))
left_lls_pred_a_diab_vol_plot <-plot(llls_pred_agediab,ci=TRUE,add.data=FALSE,
                                       line.size=1.2,show.legend=TRUE,show.title=FALSE,colors=c('blue','red'))+
                                    labs(y= bquote("Left LLS PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

rlls_pred_agediab <- ggpredict(right_lls_model, terms=c("age_at_scan [all]","diagnosed_diabetes"))
right_lls_pred_a_diab_vol_plot <-plot(rlls_pred_agediab,ci=TRUE,add.data=FALSE,
                                       line.size=1.2,show.legend=TRUE,show.title=FALSE,colors=c('blue','red'))+
                                    labs(y= bquote("Right LLS PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

Do the same shit but for the ACA

#data$ACA_pvs_vol <- data$ACAL+data$ACAR
aca_model <- lm(ACA_pvs_vol ~ ns(age_at_scan,2) + sex + ns(pollution_score,3) + 
              ns(sleep_duration_sr,3) + overall_health_rating + 
              overall_health_rating:age_at_scan + ns(hrs, 3) + overall_health_rating:hrs + 
                sleep_duration_sr:hrs + diagnosed_diabetes + diagnosed_diabetes:age_at_scan, data=data)

aca_vol <- ggpredict(aca_model, terms=c("age_at_scan [all]"))
aca_vol_plot <- plot(aca_vol,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Age")+ theme(text = element_text(size = 15))

aca_time <- ggpredict(aca_model, terms=c("hrs [all]"))
aca_time_plot <- plot(aca_time,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

# mean ci plots (dot plots with bar for each of the 4 health groups)
# only when it is converted to factors, when as continous numerilca, it's a line plot
aca_dx_vol <- ggpredict(aca_model, terms=c("overall_health_rating [all]"))
aca_dx_vol_plot <- plot(aca_dx_vol,ci=95,add.data=FALSE,show.title=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453'))+
  labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Health Rating")+ theme(text = element_text(size = 15))+
  geom_line(aes(group = "overall_health_rating"), linetype = "dashed")

# mean ci plots (dot plots with bar for each of the pollution)
aca_pollution_pred <- ggpredict(aca_model, terms=c("pollution_score"))
aca_pollution_plot <- plot(aca_pollution_pred,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Pollution Level")+ theme(text = element_text(size = 15))+
  geom_line(aes(group = "pollution_score"), linetype = "dashed")


# over time plots
aca_prediction_a_vol <- ggpredict(aca_model, terms=c("age_at_scan [all]","overall_health_rating"))
aca_prediction_a_vol_plot<-plot(aca_prediction_a_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453', 'pink'),line.size =1.2,show.legend=TRUE,show.title=FALSE)+
  labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))


aca_prediction_time_vol <- ggpredict(aca_model, terms=c("hrs [all]","overall_health_rating"))
aca_prediction_time_vol_plot<-plot(aca_prediction_time_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453', 'pink'),line.size =1.2,show.legend=FALSE,show.title=FALSE)+
  labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

aca_sleep <- ggpredict(aca_model, terms=c("hrs [all]","sleep_duration_sr"))
aca_prediction_sleep_vol_plot<-plot(aca_sleep,ci=TRUE,add.data=FALSE,line.size =1.2,show.legend=TRUE,show.title=FALSE)+
  labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Time")+ theme(text = element_text(size = 15))

aca_pred_agediab <- ggpredict(aca_model, terms=c("age_at_scan [all]","diagnosed_diabetes"))
aca_pred_a_diab_vol_plot <-plot(aca_pred_agediab,ci=TRUE,add.data=FALSE,
                                       line.size=1.2,show.legend=TRUE,show.title=FALSE,colors=c('blue','red'))+
                                    labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

# now left and right

laca_model <- lm(ACAL ~ ns(age_at_scan,2) + sex + ns(pollution_score,3) + 
              ns(sleep_duration_sr,3) + overall_health_rating + 
              overall_health_rating:age_at_scan + ns(hrs, 3) + overall_health_rating:hrs + 
                sleep_duration_sr:hrs + diagnosed_diabetes + diagnosed_diabetes:age_at_scan, data=data)

laca_vol <- ggpredict(laca_model, terms=c("age_at_scan [all]"))
laca_vol_plot <- plot(laca_vol,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("LACA PVS Volume"~(mm^3)), x = "Age")+ theme(text = element_text(size = 15))

laca_time <- ggpredict(laca_model, terms=c("hrs [all]"))
laca_time_plot <- plot(aca_time,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("LACA PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

# mean ci plots (dot plots with bar for each of the 4 health groups)
# only when it is converted to factors, when as continous numerilca, it's a line plot
laca_dx_vol <- ggpredict(laca_model, terms=c("overall_health_rating [all]"))
laca_dx_vol_plot <- plot(laca_dx_vol,ci=95,add.data=FALSE,show.title=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453'))+
  labs(y= bquote("LACA PVS Volume"~(mm^3)), x = "Health Rating")+ theme(text = element_text(size = 15))+
  geom_line(aes(group = "overall_health_rating"), linetype = "dashed")

# mean ci plots (dot plots with bar for each of the pollution)
laca_pollution_pred <- ggpredict(laca_model, terms=c("pollution_score"))
laca_pollution_plot <- plot(laca_pollution_pred,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("LACA PVS Volume"~(mm^3)), x = "Pollution Level")+ theme(text = element_text(size = 15))+
  geom_line(aes(group = "pollution_score"), linetype = "dashed")


# over time plots
laca_prediction_a_vol <- ggpredict(laca_model, terms=c("age_at_scan [all]","overall_health_rating"))
laca_prediction_a_vol_plot<-plot(laca_prediction_a_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453', 'pink'),line.size =1.2,show.legend=TRUE,show.title=FALSE)+
  labs(y= bquote("LACA PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))


laca_prediction_time_vol <- ggpredict(laca_model, terms=c("hrs [all]","overall_health_rating"))
laca_prediction_time_vol_plot<-plot(laca_prediction_time_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453', 'pink'),line.size =1.2,show.legend=FALSE,show.title=FALSE)+
  labs(y= bquote("LACA PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

laca_sleep <- ggpredict(laca_model, terms=c("hrs [all]","sleep_duration_sr"))
laca_prediction_sleep_vol_plot<-plot(laca_sleep,ci=TRUE,add.data=FALSE,line.size =1.2,show.legend=TRUE,show.title=FALSE)+
  labs(y= bquote("LACA PVS Volume"~(mm^3)), x = "Time")+ theme(text = element_text(size = 15))

laca_pred_agediab <- ggpredict(laca_model, terms=c("age_at_scan [all]","diagnosed_diabetes"))
laca_pred_a_diab_vol_plot <-plot(laca_pred_agediab,ci=TRUE,add.data=FALSE,
                                       line.size=1.2,show.legend=TRUE,show.title=FALSE,colors=c('blue','red'))+
                                    labs(y= bquote("ACA PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))


raca_model <- lm(ACAR ~ ns(age_at_scan,2) + sex + ns(pollution_score,3) + 
              ns(sleep_duration_sr,3) + overall_health_rating + 
              overall_health_rating:age_at_scan + ns(hrs, 3) + overall_health_rating:hrs + 
                sleep_duration_sr:hrs + diagnosed_diabetes + diagnosed_diabetes:age_at_scan, data=data)

raca_vol <- ggpredict(raca_model, terms=c("age_at_scan [all]"))
raca_vol_plot <- plot(raca_vol,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("RACA PVS Volume"~(mm^3)), x = "Age")+ theme(text = element_text(size = 15))

raca_time <- ggpredict(raca_model, terms=c("hrs [all]"))
raca_time_plot <- plot(raca_time,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("RACA PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

# mean ci plots (dot plots with bar for each of the 4 health groups)
# only when it is converted to factors, when as continous numerilca, it's a line plot
raca_dx_vol <- ggpredict(raca_model, terms=c("overall_health_rating [all]"))
raca_dx_vol_plot <- plot(raca_dx_vol,ci=95,add.data=FALSE,show.title=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453'))+
  labs(y= bquote("RACA PVS Volume"~(mm^3)), x = "Health Rating")+ theme(text = element_text(size = 15))+
  geom_line(aes(group = "overall_health_rating"), linetype = "dashed")

# mean ci plots (dot plots with bar for each of the pollution)
raca_pollution_pred <- ggpredict(raca_model, terms=c("pollution_score"))
raca_pollution_plot <- plot(raca_pollution_pred,ci=95,add.data=FALSE,show.title=FALSE)+
  labs(y= bquote("RACA PVS Volume"~(mm^3)), x = "Pollution Level")+ theme(text = element_text(size = 15))+
  geom_line(aes(group = "pollution_score"), linetype = "dashed")


# over time plots
raca_prediction_a_vol <- ggpredict(raca_model, terms=c("age_at_scan [all]","overall_health_rating"))
raca_prediction_a_vol_plot<-plot(raca_prediction_a_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453', 'pink'),line.size =1.2,show.legend=TRUE,show.title=FALSE)+
  labs(y= bquote("RACA PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))


raca_prediction_time_vol <- ggpredict(raca_model, terms=c("hrs [all]","overall_health_rating"))
raca_prediction_time_vol_plot<-plot(raca_prediction_time_vol,ci=TRUE,add.data=FALSE,colors=c("#79A6EA",'#A8A8A8','#E8B453', 'pink'),line.size =1.2,show.legend=FALSE,show.title=FALSE)+
  labs(y= bquote("RACA PVS Volume"~(mm^3)), x = "Time of Scan")+ theme(text = element_text(size = 15))

raca_sleep <- ggpredict(raca_model, terms=c("hrs [all]","sleep_duration_sr"))
raca_prediction_sleep_vol_plot<-plot(raca_sleep,ci=TRUE,add.data=FALSE,line.size =1.2,show.legend=TRUE,show.title=FALSE)+
  labs(y= bquote("RACA PVS Volume"~(mm^3)), x = "Time")+ theme(text = element_text(size = 15))

raca_pred_agediab <- ggpredict(raca_model, terms=c("age_at_scan [all]","diagnosed_diabetes"))
raca_pred_a_diab_vol_plot <-plot(raca_pred_agediab,ci=TRUE,add.data=FALSE,
                                       line.size=1.2,show.legend=TRUE,show.title=FALSE,colors=c('blue','red'))+
                                    labs(y= bquote("RACA PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

# plots the old way (honours project)
#pred_age_vol <- ggpredict(vol_lmm15, terms=c("Age [all]"))
#pred_age_vol_plot<-plot(pred_age_vol,ci=TRUE,add.data=TRUE, color ='blue',line.size =1,show.legend=FALSE,show.title=FALSE)+
#  labs(y= bquote("PVS Volume"~(mm^3)), x = "Age (years)")+ theme(text = element_text(size = 15))

# correlation plots

subset_data <- data[c("pvs_vol_mm3","pvs_clusters",'age_at_scan','sex',
                      "hip_circ","waist_circ","height","weight","BMI",
                      "pollution_score","townsend_deprivation_index","reason_ltfu",
                      "ipaq_activity_group","num_treatments_medications",
                      "overall_health_rating",
                      "avg_total_household_income",
                      "sleep_duration_sr","difficulty_getting_up_morning","chronotype","nap_during_day",
                      "insomnia_sr","snoring_sr", "narcolepsy_sr",
                      "seen_gp_anx_dep","freq_tiredness_twoweeks",
                      "freq_tension_twoweeks","freq_unenthusiam_twoweeks",
                      "freq_depressed_twoweeks","mood_swings",
                      
                      "miserable","irritable","fed_up",
                      "nervous","worrier_anxious","tension_highlystrung",
                      
                      "neuroticism_score",
                      "typical_daily_alcohol","alcohol_frequency",'hrs')]

subset_data <- sapply(subset_data, as.numeric)

my_cor <- cor(subset_data, use="complete.obs")
#subset_data <- lapply(subset_data,as.numeric)
corr_plot <- corrplot(my_cor, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45,
         tl.cex = 0.8)#number.cex = 0.5
#is.corr = FALSE, method = "square"

Didn’t work in the correlation plot “long_term_illness”,“weight_change_last_year”,“job_night_shift”,“job_shift_work”, “age_at_death”,###### ,“trouble_falling_asleep”,“oversleeping”, “diagnosed_diabetes”,“diagnosed_cancer”, “other_serious_condition”,“number_live_births” “contraceptive_pill_use”,“age_stroke”, “number_depression_episodes”,“seen_psych_anx_dep”,“sensitive_hurtfeelings”, “worry_after_embarressment”,“nerves”,“able_walk_cycle_tenmins”, “age_first_depression”,“age_last_depression”

“pollution_score”, “waist_circ”,“hip_circ”,“height”,“weight”,“BMI”,

“suburb”,“ethnicity”,

“ipaq_activity_group”,“num_treatments_medications”,“townsend_deprivation_index”, “reason_ltfu”,“date_ltfu”,

“overall_health_rating”,“long_term_illness”,“weight_change_last_year”, “avg_total_household_income”,“job_night_shift”,“job_shift_work”,

“age_at_death”,“age_stroke”,“number_depression_episodes”,

“sleep_duration_sr”,“difficulty_getting_up_morning”,“chronotype”,“nap_during_day” “insomnia_sr”,“snoring_sr”, “narcolepsy_sr”, “diagnosed_diabetes”,“diagnosed_cancer”, “other_serious_condition”,“number_live_births”. “contraceptive_pill_use”, “seen_psych_anx_dep”, “seen_gp_anx_dep”,“freq_tiredness_twoweeks”, “freq_tension_twoweeks”,“freq_unenthusiam_twoweeks”, “freq_depressed_twoweeks”,“mood_swings”, “miserable”,“irritable”, “sensitive_hurtfeelings”,“fed_up”, “nervous”,“worrier_anxious”, “tension_highlystrung”,“worry_after_embarressment”, “nerves”,“able_walk_cycle_tenmins”, “neuroticism_score”,“typical_daily_alcohol”, “alcohol_frequency”,“age_first_depression”, “age_last_depression”,“trouble_falling_asleep”, “oversleeping”,

“time”,“year”,“month”, “day”,“hrs”, “mnths”,“pollution_score”,“season”

time_of_scan
age_at_death
age_stroke

sleep_duration_sr
insomnia_sr
nap_during_day
snoring_sr
narcolepsy_sr

overall_health_rating
long_term_illness
diagnosed_diabetes

BMI
height
weight
waist_circ
hip_circ

Associations with disease

gamlj::gamljGlm( formula = LLSL ~ sex + age_at_scan + overall_health_rating + I(age_at_scan^2) + overall_health_rating:age_at_scan, data = data, plotHAxis = age_at_scan, plotSepLines = overall_health_rating, scaling=c(age_at_scan = “none”, sleep_duration_sr = “none”, overall_health_rating = “centered”, typical_daily_alcohol = “centered”, diagnosed_diabetes = “centered”))

ukb startup code by Lachy

Create Core Dataset

# packages

library(dplyr)
library(readr)
library(purrr)
library(tidyr)

# read UKB data

d <- readRDS("/home/lcri0006/bc41_scratch/UKB_data/ukb51956.rds")

# read data dictionary
# dict has "Field_id", "New_name", 'coding_units', 'Detail'
datadict <- read_csv("/home/lcri0006/bc41/core_ukb/core_data_dictionary.csv")


# select core variables

core <- d %>% 
  select(all_of(datadict$Field_id))

# rename variables 

names(core) <- datadict$New_name

# count missing 

map_df(core, ~sum(is.na(.))) %>% 
  pivot_longer(everything()) %>% 
  arrange(desc(value)) %>% 
  print(n = 150)

# save core dataset
write_csv(core, "/fs02/bc41/core_ukb/core_ukb.csv")

ukb Lachy code

cut ukb data

library(tidyverse)

# r formatted data storage
# tab delimited file
d <- readRDS("/home/lcri0006/bc41/UKB_data/ukb51956.rds")

# get bulk data ids
# mri ID's 

mri_id <- d %>% 
  select(eid, `20252-2.0`) %>% 
  filter(!is.na(`20252-2.0`))

write.table(mri_id, "/home/lcri0006/bc41/UKB_data/mri_id.txt",
            quote = F, col.names = F, row.names = F)

# Actigraphy ID's

acc_id <- d %>% 
  select(eid, `90001-0.0`) %>% 
  filter(!is.na(`90001-0.0`))

write.table(acc_id, "/home/lcri0006/bc41/UKB_data/accel_id.txt",
            quote = F, col.names = F, row.names = F)

# first 100 rows of actigraphy data (test)

test_acc <- acc_id %>% 
  slice(1:100)

write.table(test_acc, "/home/lcri0006/bc41/UKB_data/accel_test.txt",
            quote = F, col.names = F, row.names = F)