Add code

20210823_dFF-prediction-BMD_TV_v3.Rmd

+---
+title: "Predicting change in FF in Becker"
+author: "Thom TJ Veeger"
+date: "23/08/2021"
+output: 
+  html_document:
+    toc: true
+---
+
+# Preamble
+
+Set options, read the data.
+
+```{r set-options, echo=FALSE, cache=FALSE}
+options(width=200)
+```
+
+```{r loading data, warning=FALSE, message=FALSE}
+
+library(foreign)
+library(visreg)
+library(ggplot2)
+library(nlme)
+library(lme4)     # for the mixed model
+library(lmerTest)
+library(bootpredictlme4) # for confidence intervals visreg lmer
+library(lattice)
+
+dataset=read.spss("/Users/Thom/Documents/PhD/Drafts/BMD negative issue/R analysis/20210518_T2stdPHPDE-dFF-50thres-long_TV_v2.sav",to.data.frame=TRUE)
+dataset2=read.spss("/Users/Thom/Documents/PhD/Drafts/BMD negative issue/R analysis/dFF&FFbase_upper&lower_long.sav",to.data.frame=TRUE)
+```
+
+
+# Linear mixed model FFbase to predict ∆FF (upper and lower leg data)
+
+```{r simple regression with only FFbase and upperleg, warning=FALSE, message=FALSE}
+model = lmer(dFF ~ poly(FFbase,4) + (1|Subject) + (1|Muscle),data=dataset2, control=lmerControl(optCtrl=list(maxfun=1e8)),REML=FALSE)
+
+summary(model)
+CI=confint(model)
+CI
+
+tmp=visreg(model,plot=FALSE)
+fig3 <- plot(tmp, gg=TRUE,ylab="∆FF",xlab="FFbase",fill.par=list(fill="#ADD8E6"),points.par = list(size = 0.3))+theme(axis.text=element_text(size=10),axis.title=element_text(size=10, face = 'bold'))
+
+```
+
+# Linear mixed effects model FFbase, stdT2, PDE/ATP to predict ∆FF (lower leg data)
+
+```{r linear mixed effects model, warning=FALSE, message=FALSE}
+
+# Check for outliers and influential points
+plot(cooks.distance(model1))
+dataset$stdT2_Z = dataset %>% mutate(zscore = (stdT2 - mean(stdT2))/sd(stdT2))
+plot(dataset$stdT2_Z$stdT2)
+
+
+# Removal of datapoint 8 and running the model
+dataset_filt = dataset[-c(8),]
+model1 = lmer(dFF ~ poly(FFbase,4) + stdT2 + PDE_ATP + (1|Subject) + (1|Muscle),data=dataset_filt, control=lmerControl(optCtrl=list(maxfun=1e8)),REML=FALSE)
+summary(model1)
+r.squaredGLMM(model1)
+
+CI=confint(model1)
+CI
+
+options(bootnsim=1000)
+tmp2=visreg(model1,plot=FALSE)
+
+fig4_1 <- plot(tmp2[[1]],gg=TRUE,ylab="∆FF",xlab="FFbase",fill.par=list(fill="#ADD8E6"),points.par = list(size = 0.3))+theme(axis.text=element_text(size=10),axis.title=element_text(size=10,face="bold"))
+fig4_1
+
+fig4_2 <- plot(tmp2[[2]],gg=TRUE,ylab="∆FF",xlab=expression(bold(paste(stdT[2],phantom(l),'[ms]'))),fill.par=list(fill="#ADD8E6"),points.par = list(size = 0.3))+theme(axis.text=element_text(size=10),axis.title=element_text(size=10,face="bold"))
+fig4_2
+
+fig4_3 <- plot(tmp2[[3]],gg=TRUE,ylab="∆FF",xlab="PDE/ATP",fill.par=list(fill="#ADD8E6"),points.par = list(size = 0.3))+theme(axis.text=element_text(size=10),axis.title=element_text(size=10,face="bold"))
+fig4_3
+
+```