#callculate sample size required library(pwr) pwr.f2.test(u=2, f2=0.15, sig.level=0.01, power=0.90) #denominator df so round up and add no. of variables 89+3 #define variables as objects AnxietyT1 <- Final$"AnxietyT1" DepressionT1 <- Final$"DepressionT1" AnxietyT2 <- Final$"AnxietyT2" DepressionT2 <- Final$"DepressionT2" Resilience <- Final$"Resilience" Mindfulness <- Final$"Mindfulness" Fit1 <- Final$"Fit1" Fit2 <- Final$"Fit2" Age <- Final$`Age Group` CGAD <- Final$CGAD CPHQ <- Final$CPHQ AIMS <- Final$AIMS #set objects as required category Sex <- as.factor(Sex) Country <- as.factor(Country) Sport <- as.factor(Sport) sd(Fit1 <- as.factor(Fit1) Fit2 <- as.factor(Fit2) #Test for normality shapiro.test(Gad1) shapiro.test(Gad2) shapiro.test(PHQ1) shapiro.test(PHQ2) shapiro.test(Mind) shapiro.test(Risc) #Transform skewed dependent variables Log_Gad1 <- log(Gad1+1) Log_Gad2 <- log(Gad2+1) Log_PHQ1 <- log(PHQ1+1) Log_PHQ2 <- log(PHQ2+1) #calc correlations cor.test(DepressionT1, Mindfulness, method = "spearman") cor.test(DepressionT2, Mindfulness, method = "spearman") cor.test(AnxietyT1, Mindfulness, method = "spearman") cor.test(AnxietyT2, Mindfulness, method = "spearman") cor.test(AnxietyT1, Resilience, method = "spearman") cor.test(DepressionT1, Resilience, method = "spearman") cor.test(AnxietyT2, Resilience, method = "spearman") cor.test(DepressionT2, Resilience, method = "spearman") cor.test(DepressionT1, AnxietyT1, method = "spearman") cor.test(DepressionT2, AnxietyT2, method = "spearman") cor.test(DepressionT1, AnxietyT2, method = "spearman") cor.test(AnxietyT1, AnxietyT2, method = "spearman") cor.test(DepressionT2, AnxietyT1, method = "spearman") cor.test(DepressionT2, DepressionT1, method = "spearman") cor.test(Mindfulness, Resilience, method = "spearman") cor.test(AIMS, Mindfulness, method = "spearman") cor.test(AIMS, Resilience, method = "spearman") cor.test(AIMS, AnxietyT1, method = "spearman") cor.test(AIMS, AnxietyT2, method = "spearman") cor.test(AIMS, DepressionT1, method = "spearman") cor.test(AIMS, DepressionT2, method = "spearman") #prop tests prop.test(x = c(6, 2), c(23, 23)) prop.test(x = c(11, 13), c(160, 160)) #ttests #Between group means comparisons gadfit1 <- t.test(AnxietyT1~Fit1) t.test(AnxietyT1~Fit1) gadfit2 <- t.test(AnxietyT2~Fit2) t.test(AnxietyT2~Fit2) phqfit1 <- t.test(DepressionT1~Fit1) t.test(DepressionT1~Fit1) phqfit2 <- t.test(Depression2~Fit1) t.test(DepressionT2~Fit2) t.test(DepressionT1,DepressionT2,var.equal = T,paired = F) t_gf2 <-gadfit2$statistic[[1]] df_gf2 <- gadfit2$parameter[[1]] r_gf2 <- sqrt(t_gf2^2/(t_gf2^2+df_gf2)) r_gf2 t.test(AnxietyT1, AnxietyT2, alternative = "two.sided", var.equal = FALSE) t.test(DepressionT1, DepressionT2, alternative = "two.sided", var.equal = FALSE) #regression lm1 <- lm(AnxietyT1~ Mindfulness, data = Final) summary(lm1) pSD17 <- sd(Mindfulness) /sd(AnxietyT1) -0.03943*pSD17 lm2 <- lm(DepressionT1 ~ Mindfulness, data = Final) summary(lm2) lm3 <- lm(AnxietyT1 ~ Mindfulness+Resilience, data = Final) summary(lm3) lm4 <- lm(DepressionT1 ~ Mindfulness+Resilience, data = Final) summary(lm4) lm5 <- lm(Gad1~Risc, data = Final) summary(lm5) lm6 <- lm(PHQ1 ~ Risc, data = Final) summary(lm6) lm7 <- lm(Gad2 ~ Risc, data = Final) summary(lm7) lm8 <- lm(PHQ2 ~ Risc, data = Final) summary(lm8) #change regression lm_gadC <- lm(AnxietyT2 ~ AnxietyT1+Mindfulness, data = Final) summary(lm_gadC) lm_phqC <- lm(DepressionT2 ~ DepressionT1+Mindfulness, data = Final) summary(lm_phqC) plot(lm_phqC) lm_gadCR <- lm(AnxietyT2 ~ AnxietyT1+Resilience, data = Final) summary(lm_gadCR) lm_phqCR <- lm(DepressionT2 ~ DepressionT1+Resilience , data = Final) summary(lm_phqCR) #further builds lm_gadC1 <- lm(AnxietyT2 ~ AnxietyT1+Mindfulness+Fit1, data = Final) summary(lm_gadC1) lm_gadC2 <- lm(AnxietyT2 ~ AnxietyT1+Mindfulness+Fit1+Final$Transferability, data = Final) summary(lm_gadC2) Edu <- Final$Education lm_gadC3 <- lm(AnxietyT2 ~ AnxietyT1+Mindfulness+Fit1+Final$Transferability , data = Final) summary(lm_gadC3) lm_gadC4 <- lm(AnxietyT2 ~ AnxietyT1+Resilience+Fit1+Final$Transferability, data = Final) summary(lm_gadC4) lm_gadC5 <- lm(AnxietyT2 ~ AnxietyT1+Resilience+Fit1, data = Final) summary(lm_gadC5) plot(lm_gadC5) #plots for Dep #mod3 lm_phqC1 <- lm(DepressionT2 ~ DepressionT1+Mindfulness+Final$Transferability, data = Final) summary(lm_phqC1) lm_phqC2 <- lm(DepressionT2 ~ DepressionT1+Mindfulness+Final$Plan, data = Final) summary(lm_phqC2) lm_phqC3 <- lm(DepressionT2 ~ DepressionT1+Mindfulness+Sex, data = Final) summary(lm_phqC1) #mod3 lm_phqC4 <- lm(DepressionT2 ~ DepressionT1+Mindfulness+Final$Transferability, data = Final) summary(lm_phqC4) lm_phqC5 <- lm(DepressionT2 ~ DepressionT1+Mindfulness+AIMS, data = Final) summary(lm_phqC5) #mod4 Age <- Final$`Age Group` lm_phqC6 <- lm(DepressionT2 ~ DepressionT1+Resilience+Final$Transferability+Age+Sex, data = Final) summary(lm_phqC6) #basic theo models lm_gadCT1 <- lm(AnxietyT2 ~ AnxietyT1+Mindfulness+Age+Sex, data = Final) summary(lm_gadCT1) lm_gadCT2 <- lm(AnxietyT2 ~ AnxietyT1+Resilience+Age+Sex, data = Final) summary(lm_gadCT2) lm_phqCT1 <- lm(DepressionT2 ~ DepressionT1+Mindfulness+Age+Sex, data = Final) summary(lm_phqCT1) lm_phqCT2 <- lm(DepressionT2 ~ DepressionT1+Resilience+Age+Sex, data = Final) summary(lm_phqCT2) #advanced theo models Age <- Final$`Age Group` Trans <- Final$Transferability lm_gadCT1a <- lm(AnxietyT2 ~ AnxietyT1+Mindfulness+Fit1+Age+Sex, data = Final) summary(lm_gadCT1a) lm_gadCT1a <- lm(AnxietyT2 ~ AnxietyT1+Mindfulness+Fit1+Age+Sex+Trans, data = Final) summary(lm_gadCT1a) lm_gadCT2a <- lm(AnxietyT2 ~ AnxietyT1+Resilience+Age+Sex+Fit1, data = Final) summary(lm_gadCT2a) lm_gadCT2b <- lm(AnxietyT2 ~ AnxietyT1+Resilience+Age+Sex+Fit1+Trans, data = Final) summary(lm_gadCT2b) lm_phqCT1a <- lm(DepressionT2 ~ DepressionT1+Mindfulness+Age+Sex+Fit1, data = Final) summary(lm_phqCT1a) #stepwise models Mod1 <- lm(DepressionT2~1, data = Final) summary(Mod1) FitAll <- lm(DepressionT2~DepressionT1+Mindfulness+Resilience+Age+Final$Sex+AIMS++Fit2+Final$Fit1) DepStep1 <- step(Mod1, direction = "both", scope = formula(FitAll) ) summary(DepStep1) DepStep1 Mod2 <- lm(AnxietyT2~1, data = Final) FitAll1 <- lm(AnxietyT2~AnxietyT1+Mindfulness+Resilience+Age+Final$Sex+AIMS+Trans+Edu+Fit2+Final$Fit1) AnxStep1 <- step(Mod2, direction = "both", scope = formula(FitAll1) ) summary(AnxStep1) #again without transferability USED! Mod1a <- lm(DepressionT2~1, data = Final) summary(Mod1a) FitAll <- lm(DepressionT2~DepressionT1+Mindfulness+Resilience+Final$Sex+AIMS+Fit2+Final$Fit1+Final$`Age Group`+Final$Country) DepStep1 <- step(Mod1a, direction = "forward", scope = formula(FitAll) ) summary(DepStep1) Mod2a <- lm(AnxietyT2~1, data = Final) FitAll1 <- lm(AnxietyT2~AnxietyT1+Mindfulness+Resilience+Final$Sex+AIMS+Fit2+Final$Fit1+Final$Country+Final$`Age Group`) AnxStep1 <- step(Mod2a, direction = "forward", scope = formula(FitAll1) ) summary(AnxStep1) #standardised betas pooledSD <- sd(AnxietyT1) / sd(AnxietyT2) StandB <- pooledSD*0.45211 (sd(Mindfulness)/sd(AnxietyT2))*-0.11882 sd(Final$Fit1) # lm(scale(AnxietyT2) ~ scale(1), data=Final) install.packages("QuantPsyc") library(QuantPsyc) lm.beta(AnxStep1) lm.beta(DepStep1) #third try install.packages("MASS") library(MASS) # Fit the full model full.model <- lm(DepressionT2 ~1, data = Final) # Stepwise regression model step.model <- stepAIC(full.model, direction = "both", trace = FALSE) summary(step.model) step.model install.packages("leaps") library(leaps) Depmodel <- regsubsets(DepressionT2~1, data = Final, nvmax = 5, method = "seqrep") summary(Depmodel) #attempt 4 modelanx1 <- lm(AnxietyT2~AnxietyT1+Mindfulness+Resilience+Age+Final$Sex+AIMS+Edu+Fit2+Final$Fit1) modanx <- ols_step_forward_p(modelanx1, penter = 0.5 ) modanx modeldep1 <- lm(DepressionT2~DepressionT1+Mindfulness+Resilience+Age+Final$Sex+AIMS+Edu+Fit2+Final$Fit1) moddep <- ols_step_forward_p(modeldep1, penter = .07) moddep plot(moddep) #bivariate regression BVDep <- lm(DepressionT2~DepressionT1+Mindfulness) summary(BVDep) BVAnx <- lm(AnxietyT2~AnxietyT1+Mindfulness) summary(BVAnx) BVAnxa <- lm(AnxietyT2~AnxietyT1+Resilience) summary(BVAnxa) BVAnxb <- lm(AnxietyT2~AnxietyT1+Final$Fit1) summary(BVAnxb) BVAnxc <- lm(AnxietyT2~AnxietyT1+Mindfulness+Resilience+Final$Fit1) summary(BVAnxc) #Logistic regression RiskD <- Final1$RiskD as.factor(RiskD) modelD <- glm(RiskD ~ DepressionT1+Mindfulness,family=binomial,data=Final1) summary(modelD) modelA <- glm(RiskA~ AnxietyT1+Age, family = binomial, data = Final1) summary(modelA) ####### plot(residuals) termplot(lm_phqC,se=TRUE,partial.resid=TRUE, main="Auto-Regression Model", xlab="PHQ Score at T1", ylab="PHQ Score at T2") ggplot(data = Final, aes(x = PHQ2, y = PHQ1)) + geom_point() + geom_smooth(method="lm", formula= y ~ 1 + x, se=TRUE, fullrange=TRUE, color="red", size=2) + xlab("PHQ Score at T1") + ylab("PHQ Score at T2") + ggtitle("Auto regression model") library(ggplot2) lm_gadD <- lm(CGAD ~ Gad1, data = Cleaned_Final) summary(lm_gadD) #sem library(lavaan) library(semPlot) model <- ' GAD1~CAMSR+CDRisc CDRisc~CAMSR PHQ1~CAMSR+CDRisc GAD2~CAMSR+CDRisc PHQ2~CAMSR+CDRisc ' fit <- sem(model, data = Final) summary(fit, fit.measures = T, rsquare = T) semPaths(fit, whatLabels = "std", layout = "circle2") parameterEstimates(fit) #path analysis library(tidyverse) library(tidygraph) library(ggraph) library(lavaan) # Plot a fitted lavaan object ggsem <- function(fit, layout = "kk") { # Extract standardized parameters params <- lavaan::standardizedSolution(fit) # Edge properties param_edges <- params %>% filter(op %in% c("=~", "~", "~~"), lhs != rhs, pvalue < .10) %>% transmute(to = lhs, from = rhs, val = est.std, type = dplyr::case_when( op == "=~" ~ "loading", op == "~" ~ "regression", op == "~~" ~ "correlation", TRUE ~ NA_character_)) # Identify latent variables for nodes latent_nodes <- param_edges %>% filter(type == "loading") %>% distinct(to) %>% transmute(metric = to, latent = TRUE) # Node properties param_nodes <- params %>% filter(lhs == rhs) %>% transmute(metric = lhs, e = est.std) %>% left_join(latent_nodes) %>% mutate(latent = if_else(is.na(latent), FALSE, latent)) # Complete Graph Object param_graph <- tidygraph::tbl_graph(param_nodes, param_edges) # Plot ggraph(param_graph, layout = layout) + # Latent factor Nodes geom_node_point(aes(alpha = as.numeric(latent)), shape = 16, size = 5) + geom_node_point(aes(alpha = as.numeric(latent)), shape = 16, size = 4, color = "white") + # Observed Nodes geom_node_point(aes(alpha = as.numeric(!latent)), shape = 15, size = 5) + geom_node_point(aes(alpha = as.numeric(!latent)), shape = 15, size = 4, color = "white") + # Regression Paths (and text) geom_edge_link(aes(color = val, label = round(val, 2), alpha = as.numeric(type == "regression")), linetype = 1, angle_calc = "along", vjust = -.5, arrow = arrow(20, unit(.3, "cm"), type = "closed")) + # Factor Loadings (no text) geom_edge_link(aes(color = val , alpha = as.numeric(type == "loading")), linetype = 3, angle_calc = "along", arrow = arrow(20, unit(.3, "cm"), ends = "first", type = "closed")) + # Correlation Paths (no text) geom_edge_link(aes(color = val, alpha = as.numeric(type == "correlation")), linetype = 2, angle_calc = "along", arrow = arrow(20, unit(.3, "cm"), type = "closed", ends = "both")) + # Node names geom_node_text(aes(label = metric), nudge_y = .2, hjust = "along") + # Node residual error geom_node_text(aes(label = sprintf("%.2f", e)), nudge_y = -.1, size = 3) + # Scales and themes scale_alpha(guide = FALSE, range = c(0, 1)) + scale_edge_alpha(guide = FALSE, range = c(0, 1)) + scale_edge_colour_gradient2(guide = FALSE, low = "red", mid = "darkgray", high = "green") + scale_edge_linetype(guide = FALSE) + scale_size(guide = FALSE) + theme_graph() } model <- ' AnxietyT1~Mindfulness+Resilience Resilience~Mindfulness DepressionT1~Mindfulness+Resilience AnxietyT2~Mindfulness+Resilience DepressionT2~Mindfulness+Resilience ' fit <- sem(model, data = Final) ggsem(fit) ################################################# # Plot a fitted lavaan object ggsem <- function(fit, layout = "kk") { # Extract standardized parameters params <- lavaan::standardizedSolution(fit) # Edge properties param_edges <- params %>% filter(op %in% c("=~", "~", "~~"), lhs != rhs, pvalue < .10) %>% transmute(to = lhs, from = rhs, val = est.std, type = dplyr::case_when( op == "=~" ~ "loading", op == "~" ~ "regression", op == "~~" ~ "correlation", TRUE ~ NA_character_)) # Identify latent variables for nodes latent_nodes <- param_edges %>% filter(type == "loading") %>% distinct(to) %>% transmute(metric = to, latent = TRUE) # Node properties param_nodes <- params %>% filter(lhs == rhs) %>% transmute(metric = lhs, e = est.std) %>% left_join(latent_nodes) %>% mutate(latent = if_else(is.na(latent), FALSE, latent)) # Complete Graph Object param_graph <- tidygraph::tbl_graph(param_nodes, param_edges) # Plot ggraph(param_graph, layout = layout) + # Latent factor Nodes geom_node_point(aes(alpha = as.numeric(latent)), shape = 16, size = 5) + geom_node_point(aes(alpha = as.numeric(latent)), shape = 16, size = 4, color = "white") + # Observed Nodes geom_node_point(aes(alpha = as.numeric(!latent)), shape = 15, size = 5) + geom_node_point(aes(alpha = as.numeric(!latent)), shape = 15, size = 4, color = "white") + # Regression Paths (and text) geom_edge_link(aes(color = val, label = round(val, 2), alpha = as.numeric(type == "regression")), linetype = 1, angle_calc = "along", vjust = -.5, arrow = arrow(20, unit(.3, "cm"), type = "closed")) + # Factor Loadings (no text) geom_edge_link(aes(color = val , alpha = as.numeric(type == "loading")), linetype = 3, angle_calc = "along", arrow = arrow(20, unit(.3, "cm"), ends = "first", type = "closed")) + # Node names geom_node_text(aes(label = metric), nudge_y = .1, hjust = "along") + # Node residual error geom_node_text(aes(label = sprintf("%.2f", e)), nudge_y = -.1, size = 3) + # Scales and themes scale_alpha(guide = FALSE, range = c(0, 1)) + scale_edge_alpha(guide = FALSE, range = c(0, 1)) + scale_edge_colour_gradient2(guide = FALSE, low = "red", mid = "darkgray", high = "green") + scale_edge_linetype(guide = FALSE) + scale_size(guide = FALSE) + theme_graph() } model <- ' AnxietyT1~Mindfulness+Resilience Resilience~Mindfulness DepressionT1~Mindfulness+Resilience AnxietyT2~Mindfulness+Resilience DepressionT2~Mindfulness+Resilience ' #using ggsave ggplot(aes(model), data = Final) + geom_point() + geom_smooth(method=lm, fill = NA, fullrange=TRUE, color = "black") ggsave("test.tiff", units="in", width=5, height=4, dpi=300, compression = 'lzw') #using tiff() and dev.off tiff('test.tiff', units="in", width=5, height=4, res=300, compression = 'lzw') ggplot(aes(x, y), data = ddata) + geom_point() + geom_smooth(method=lm, fill = NA, fullrange=TRUE, color = "black") tiff("test.tiff", units="in", width=5, height=5, res=300) fit <- sem(model, data = Final) ggsem(fit) #make a bar chart for figure 1 library(ggplot2) Anx <- Figure1$Anx Dep <- Figure1$Dep as.factor(Anx) as.factor(Dep) Anxiety <- Figure1$Level Time <- Figure1$Condition ggplot(Figure1, aes(fill=Time, y=Anx, x = reorder(Anxiety, - Anx))) + geom_bar(position="dodge", stat = "identity") ggplot(Figure1, aes(x=reorder(Anxiety, - Anx), y=Anx, fill=Time)) + geom_bar(position=position_dodge(), stat="identity") + geom_errorbar(aes(ymin=Anx-ci, ymax=Anx+ci), width=.2, # Width of the error bars position=position_dodge(.9)) #scatterplot # Libraries library(GGally) # Plot To Use! ggparcoord(data = Figure1_scat, columns = c(1:2), groupColumn = 2, scale = "globalminmax", boxplot = FALSE, title = "Change in Anxiety (Gad7 score) from T1 to T2") #Scatterplot coloured by change ggparcoord(data = Figure1_scat, columns = c(1:2), groupColumn = "Change", scale = "globalminmax", boxplot = FALSE, title = "Change in Anxiety (Gad7 score) from T1 to T2", alphaLines = 0.5)) #Scatterplot highlighted change library(ggplot2) library(GGally) library(dplyr) as.list(Figure1_scat$Change) Reduced <-within(Figure1_scat, Change<-if_else(cut=="Reduce", "other")) ggparcoord(data = Figure1_scat, columns = c(1:2), groupColumn = "Change", scale = "globalminmax", boxplot = FALSE, title = "Change in Anxiety (Gad7 score) from T1 to T2", scale_color_manual(values=c("maroon","gray"))) #New analysis shapiro.test(CGad) hist(CGad) CGad <- Final$CGAD cor.test(Mindfulness, CGad) cor.test(Resilience, CGad) shapiro.test(CDep) hist(CDep) CDep <- Final$CPHQ cor.test(Mindfulness, CDep, method = "pearson") cor.test(Resilience, CDep, method = "pearson") plot(Mindfulness, CGad) #ANOVA analysis install.packages("compute.es") library(compute.es) install.packages("car") library(car) install.packages("ggplot") library(ggplot2) install.packages("multcomp") library(multcomp) install.packages("pastecs") library(pastecs) by(Final$AnxietyT1, Final$`Age Group`, stat.desc) by(Final$AnxietyT2, Final$`Age Group`, stat.desc) by(Final$DepressionT1, Final$`Age Group`, stat.desc) by(Final$DepressionT2, Final$`Age Group`, stat.desc) agem <- aov(DepressionT2~Final$`Age Group`, data = Final) summary(agem) agem1 <- aov(AnxietyT2~Final$`Age Group`, data = Final) summary(agem1) par(1,1) plot(agem)