--- title: "Mindset-OA_boswellm" output: pdf_document: default html_document: default --- ```{r setup, include=FALSE, warning=FALSE} knitr::opts_chunk$set(echo = FALSE) ``` ## The code is the analysis associated with the following paper: Boswell M., et al. 2021. Mindset is Associated with Future Physical Activity and Management Strategies in Individuals with Knee Osteoarthritis: A Repeated Cross-Sectional Survey. Annals of Physical and Rehab Med. ```{r} # Load libraries library(psych) require(stddiff) library(plyr) library(sjPlot) library(sjmisc) library(lmSupport) library(stats) library(jtools) library(ggstance) library(Rmisc) library(tidyr) library(dplyr) library(lme4) library(psy) library(nFactors) library(ltm) library(irr) ``` ```{r} # Set working directory setwd("/Users/melissa/Google Drive/NMBL/OA/Mindset/Osteoarthritis and Mindset/Study 2 - Survey/SimTK") # Load Data OA_T1_0 <- read.csv("OA_T1.csv") # OA group at T1 (baseline) noOA_T1_0 <- read.csv("noOA_T1.csv") # No OA group at T1 (baseline) OA_T2_0 <- read.csv("OA_T2.csv") # OA group at T2 (three weeks later) noOA_T2_0 <- read.csv("noOA_T2.csv") # No OA group at T2 (three weeks later) # Remove extra headings from survey OA_T1_0 <- OA_T1_0[-c(1), ] noOA_T1_0 <- noOA_T1_0[-c(1), ] OA_T2_0 <- OA_T2_0[-c(1), ] noOA_T2_0 <- noOA_T2_0[-c(1), ] # Make new data frame OA_T1 <- OA_T1_0 noOA_T1 <- noOA_T1_0 OA_T2 <- OA_T2_0 noOA_T2 <- noOA_T2_0 # Make Numeric OA_T1[] <- sapply(OA_T1, as.numeric) noOA_T1[] <- sapply(noOA_T1, as.numeric) OA_T2[] <- sapply(OA_T2, as.numeric) noOA_T2[] <- sapply(noOA_T2, as.numeric) # Bring back aids OA_T1$aid <- OA_T1_0$aid noOA_T1$aid <- noOA_T1_0$aid OA_T2$aid <- OA_T2_0$aid noOA_T2$aid <- noOA_T2_0$aid ``` ```{r} # Survey measure summary scores # Re-code and Summarize PASE scale # OA_T1 OA_T1$PASE_2a[OA_T1$PASE_2a == 0] <- 0.00 OA_T1$PASE_3a[OA_T1$PASE_3a == 0] <- 0.00 OA_T1$PASE_4a[OA_T1$PASE_4a == 0] <- 0.00 OA_T1$PASE_5a[OA_T1$PASE_5a == 0] <- 0.00 OA_T1$PASE_6a[OA_T1$PASE_6a == 0] <- 0.00 OA_T1$PASE_2a[OA_T1$PASE_2a == 1] <- 1.50 OA_T1$PASE_3a[OA_T1$PASE_3a == 1] <- 1.50 OA_T1$PASE_4a[OA_T1$PASE_4a == 1] <- 1.50 OA_T1$PASE_5a[OA_T1$PASE_5a == 1] <- 1.50 OA_T1$PASE_6a[OA_T1$PASE_6a == 1] <- 1.50 OA_T1$PASE_2a[OA_T1$PASE_2a == 2] <- 3.50 OA_T1$PASE_3a[OA_T1$PASE_3a == 2] <- 3.50 OA_T1$PASE_4a[OA_T1$PASE_4a == 2] <- 3.50 OA_T1$PASE_5a[OA_T1$PASE_5a == 2] <- 3.50 OA_T1$PASE_6a[OA_T1$PASE_6a == 2] <- 3.50 OA_T1$PASE_2a[OA_T1$PASE_2a == 3] <- 6.00 OA_T1$PASE_3a[OA_T1$PASE_3a == 3] <- 6.00 OA_T1$PASE_4a[OA_T1$PASE_4a == 3] <- 6.00 OA_T1$PASE_5a[OA_T1$PASE_5a == 3] <- 6.00 OA_T1$PASE_6a[OA_T1$PASE_6a == 3] <- 6.00 OA_T1$PASE_2b[OA_T1$PASE_2b == 1] <- 0.50 OA_T1$PASE_3b[OA_T1$PASE_3b == 1] <- 0.50 OA_T1$PASE_4b[OA_T1$PASE_4b == 1] <- 0.50 OA_T1$PASE_5b[OA_T1$PASE_5b == 1] <- 0.50 OA_T1$PASE_6b[OA_T1$PASE_6b == 1] <- 0.50 OA_T1$PASE_2b[OA_T1$PASE_2b == 2] <- 1.50 OA_T1$PASE_3b[OA_T1$PASE_3b == 2] <- 1.50 OA_T1$PASE_4b[OA_T1$PASE_4b == 2] <- 1.50 OA_T1$PASE_5b[OA_T1$PASE_5b == 2] <- 1.50 OA_T1$PASE_6b[OA_T1$PASE_6b == 2] <- 1.50 OA_T1$PASE_2b[OA_T1$PASE_2b == 3] <- 3.00 OA_T1$PASE_3b[OA_T1$PASE_3b == 3] <- 3.00 OA_T1$PASE_4b[OA_T1$PASE_4b == 3] <- 3.00 OA_T1$PASE_5b[OA_T1$PASE_5b == 3] <- 3.00 OA_T1$PASE_6b[OA_T1$PASE_6b == 3] <- 3.00 OA_T1$PASE_2b[OA_T1$PASE_2b == 4] <- 5.00 OA_T1$PASE_3b[OA_T1$PASE_3b == 4] <- 5.00 OA_T1$PASE_4b[OA_T1$PASE_4b == 4] <- 5.00 OA_T1$PASE_5b[OA_T1$PASE_5b == 4] <- 5.00 OA_T1$PASE_6b[OA_T1$PASE_6b == 4] <- 5.00 # Deal with NAs OA_T1$PASE_2b[is.na(OA_T1$PASE_2b)] <- 0 OA_T1$PASE_3b[is.na(OA_T1$PASE_3b)] <- 0 OA_T1$PASE_4b[is.na(OA_T1$PASE_4b)] <- 0 OA_T1$PASE_5b[is.na(OA_T1$PASE_5b)] <- 0 OA_T1$PASE_6b[is.na(OA_T1$PASE_6b)] <- 0 OA_T1$PASE_2tot <- OA_T1$PASE_2a * OA_T1$PASE_2b/7 OA_T1$PASE_3tot <- OA_T1$PASE_3a * OA_T1$PASE_3b/7 OA_T1$PASE_4tot <- OA_T1$PASE_4a * OA_T1$PASE_4b/7 OA_T1$PASE_5tot <- OA_T1$PASE_5a * OA_T1$PASE_5b/7 OA_T1$PASE_6tot <- OA_T1$PASE_6a * OA_T1$PASE_6b/7 # OA_T2 OA_T2$PASE_2a[OA_T2$PASE_2a == 0] <- 0.00 OA_T2$PASE_3a[OA_T2$PASE_3a == 0] <- 0.00 OA_T2$PASE_4a[OA_T2$PASE_4a == 0] <- 0.00 OA_T2$PASE_5a[OA_T2$PASE_5a == 0] <- 0.00 OA_T2$PASE_6a[OA_T2$PASE_6a == 0] <- 0.00 OA_T2$PASE_2a[OA_T2$PASE_2a == 1] <- 1.50 OA_T2$PASE_3a[OA_T2$PASE_3a == 1] <- 1.50 OA_T2$PASE_4a[OA_T2$PASE_4a == 1] <- 1.50 OA_T2$PASE_5a[OA_T2$PASE_5a == 1] <- 1.50 OA_T2$PASE_6a[OA_T2$PASE_6a == 1] <- 1.50 OA_T2$PASE_2a[OA_T2$PASE_2a == 2] <- 3.50 OA_T2$PASE_3a[OA_T2$PASE_3a == 2] <- 3.50 OA_T2$PASE_4a[OA_T2$PASE_4a == 2] <- 3.50 OA_T2$PASE_5a[OA_T2$PASE_5a == 2] <- 3.50 OA_T2$PASE_6a[OA_T2$PASE_6a == 2] <- 3.50 OA_T2$PASE_2a[OA_T2$PASE_2a == 3] <- 6.00 OA_T2$PASE_3a[OA_T2$PASE_3a == 3] <- 6.00 OA_T2$PASE_4a[OA_T2$PASE_4a == 3] <- 6.00 OA_T2$PASE_5a[OA_T2$PASE_5a == 3] <- 6.00 OA_T2$PASE_6a[OA_T2$PASE_6a == 3] <- 6.00 OA_T2$PASE_2b[OA_T2$PASE_2b == 1] <- 0.50 OA_T2$PASE_3b[OA_T2$PASE_3b == 1] <- 0.50 OA_T2$PASE_4b[OA_T2$PASE_4b == 1] <- 0.50 OA_T2$PASE_5b[OA_T2$PASE_5b == 1] <- 0.50 OA_T2$PASE_6b[OA_T2$PASE_6b == 1] <- 0.50 OA_T2$PASE_2b[OA_T2$PASE_2b == 2] <- 1.50 OA_T2$PASE_3b[OA_T2$PASE_3b == 2] <- 1.50 OA_T2$PASE_4b[OA_T2$PASE_4b == 2] <- 1.50 OA_T2$PASE_5b[OA_T2$PASE_5b == 2] <- 1.50 OA_T2$PASE_6b[OA_T2$PASE_6b == 2] <- 1.50 OA_T2$PASE_2b[OA_T2$PASE_2b == 3] <- 3.00 OA_T2$PASE_3b[OA_T2$PASE_3b == 3] <- 3.00 OA_T2$PASE_4b[OA_T2$PASE_4b == 3] <- 3.00 OA_T2$PASE_5b[OA_T2$PASE_5b == 3] <- 3.00 OA_T2$PASE_6b[OA_T2$PASE_6b == 3] <- 3.00 OA_T2$PASE_2b[OA_T2$PASE_2b == 4] <- 5.00 OA_T2$PASE_3b[OA_T2$PASE_3b == 4] <- 5.00 OA_T2$PASE_4b[OA_T2$PASE_4b == 4] <- 5.00 OA_T2$PASE_5b[OA_T2$PASE_5b == 4] <- 5.00 OA_T2$PASE_6b[OA_T2$PASE_6b == 4] <- 5.00 # Deal with NAs OA_T2$PASE_2b[is.na(OA_T2$PASE_2b)] <- 0 OA_T2$PASE_3b[is.na(OA_T2$PASE_3b)] <- 0 OA_T2$PASE_4b[is.na(OA_T2$PASE_4b)] <- 0 OA_T2$PASE_5b[is.na(OA_T2$PASE_5b)] <- 0 OA_T2$PASE_6b[is.na(OA_T2$PASE_6b)] <- 0 OA_T2$PASE_2tot <- OA_T2$PASE_2a * OA_T2$PASE_2b/7 OA_T2$PASE_3tot <- OA_T2$PASE_3a * OA_T2$PASE_3b/7 OA_T2$PASE_4tot <- OA_T2$PASE_4a * OA_T2$PASE_4b/7 OA_T2$PASE_5tot <- OA_T2$PASE_5a * OA_T2$PASE_5b/7 OA_T2$PASE_6tot <- OA_T2$PASE_6a * OA_T2$PASE_6b/7 # noOA_T1 noOA_T1$PASE_2a[noOA_T1$PASE_2a == 0] <- 0.00 noOA_T1$PASE_3a[noOA_T1$PASE_3a == 0] <- 0.00 noOA_T1$PASE_4a[noOA_T1$PASE_4a == 0] <- 0.00 noOA_T1$PASE_5a[noOA_T1$PASE_5a == 0] <- 0.00 noOA_T1$PASE_6a[noOA_T1$PASE_6a == 0] <- 0.00 noOA_T1$PASE_2a[noOA_T1$PASE_2a == 1] <- 1.50 noOA_T1$PASE_3a[noOA_T1$PASE_3a == 1] <- 1.50 noOA_T1$PASE_4a[noOA_T1$PASE_4a == 1] <- 1.50 noOA_T1$PASE_5a[noOA_T1$PASE_5a == 1] <- 1.50 noOA_T1$PASE_6a[noOA_T1$PASE_6a == 1] <- 1.50 noOA_T1$PASE_2a[noOA_T1$PASE_2a == 2] <- 3.50 noOA_T1$PASE_3a[noOA_T1$PASE_3a == 2] <- 3.50 noOA_T1$PASE_4a[noOA_T1$PASE_4a == 2] <- 3.50 noOA_T1$PASE_5a[noOA_T1$PASE_5a == 2] <- 3.50 noOA_T1$PASE_6a[noOA_T1$PASE_6a == 2] <- 3.50 noOA_T1$PASE_2a[noOA_T1$PASE_2a == 3] <- 6.00 noOA_T1$PASE_3a[noOA_T1$PASE_3a == 3] <- 6.00 noOA_T1$PASE_4a[noOA_T1$PASE_4a == 3] <- 6.00 noOA_T1$PASE_5a[noOA_T1$PASE_5a == 3] <- 6.00 noOA_T1$PASE_6a[noOA_T1$PASE_6a == 3] <- 6.00 noOA_T1$PASE_2b[noOA_T1$PASE_2b == 1] <- 0.50 noOA_T1$PASE_3b[noOA_T1$PASE_3b == 1] <- 0.50 noOA_T1$PASE_4b[noOA_T1$PASE_4b == 1] <- 0.50 noOA_T1$PASE_5b[noOA_T1$PASE_5b == 1] <- 0.50 noOA_T1$PASE_6b[noOA_T1$PASE_6b == 1] <- 0.50 noOA_T1$PASE_2b[noOA_T1$PASE_2b == 2] <- 1.50 noOA_T1$PASE_3b[noOA_T1$PASE_3b == 2] <- 1.50 noOA_T1$PASE_4b[noOA_T1$PASE_4b == 2] <- 1.50 noOA_T1$PASE_5b[noOA_T1$PASE_5b == 2] <- 1.50 noOA_T1$PASE_6b[noOA_T1$PASE_6b == 2] <- 1.50 noOA_T1$PASE_2b[noOA_T1$PASE_2b == 3] <- 3.00 noOA_T1$PASE_3b[noOA_T1$PASE_3b == 3] <- 3.00 noOA_T1$PASE_4b[noOA_T1$PASE_4b == 3] <- 3.00 noOA_T1$PASE_5b[noOA_T1$PASE_5b == 3] <- 3.00 noOA_T1$PASE_6b[noOA_T1$PASE_6b == 3] <- 3.00 noOA_T1$PASE_2b[noOA_T1$PASE_2b == 4] <- 5.00 noOA_T1$PASE_3b[noOA_T1$PASE_3b == 4] <- 5.00 noOA_T1$PASE_4b[noOA_T1$PASE_4b == 4] <- 5.00 noOA_T1$PASE_5b[noOA_T1$PASE_5b == 4] <- 5.00 noOA_T1$PASE_6b[noOA_T1$PASE_6b == 4] <- 5.00 # Deal with Nas noOA_T1$PASE_2b[is.na(noOA_T1$PASE_2b)] <- 0 noOA_T1$PASE_3b[is.na(noOA_T1$PASE_3b)] <- 0 noOA_T1$PASE_4b[is.na(noOA_T1$PASE_4b)] <- 0 noOA_T1$PASE_5b[is.na(noOA_T1$PASE_5b)] <- 0 noOA_T1$PASE_6b[is.na(noOA_T1$PASE_6b)] <- 0 # Compute activity frequency noOA_T1$PASE_2tot <- noOA_T1$PASE_2a * noOA_T1$PASE_2b/7 noOA_T1$PASE_3tot <- noOA_T1$PASE_3a * noOA_T1$PASE_3b/7 noOA_T1$PASE_4tot <- noOA_T1$PASE_4a * noOA_T1$PASE_4b/7 noOA_T1$PASE_5tot <- noOA_T1$PASE_5a * noOA_T1$PASE_5b/7 noOA_T1$PASE_6tot <- noOA_T1$PASE_6a * noOA_T1$PASE_6b/7 ## Question 10 OA_T1$PASE_10tot <- OA_T1$PASE_10a/7 OA_T2$PASE_10tot <- OA_T2$PASE_10a/7 noOA_T1$PASE_10tot <- noOA_T1$PASE_10a/7 ## If 10b = 0 then 10 total = 0 OA_T1$PASE_10tot[which(OA_T1$PASE_10b == 0)] <- 0 OA_T2$PASE_10tot[which(OA_T2$PASE_10b == 0)] <- 0 noOA_T1$PASE_10tot[which(noOA_T1$PASE_10b == 0)] <- 0 ## PASE Total Score OA_T1$PASE_T1 <- OA_T1$PASE_2tot*20 + OA_T1$PASE_3tot*21 + OA_T1$PASE_4tot*23 + OA_T1$PASE_5tot*23 + OA_T1$PASE_6tot*30 + OA_T1$PASE_7*25 + OA_T1$PASE_8*25 + OA_T1$PASE_9_a*30 + OA_T1$PASE_9_b*36 + OA_T1$PASE_9_c*20 + OA_T1$PASE_9_d*35 + OA_T1$PASE_10tot*21 OA_T2$PASE_T2 <- OA_T2$PASE_2tot*20 + OA_T2$PASE_3tot*21 + OA_T2$PASE_4tot*23 + OA_T2$PASE_5tot*23 + OA_T2$PASE_6tot*30 + OA_T2$PASE_7*25 + OA_T2$PASE_8*25 + OA_T2$PASE_9_a*30 + OA_T2$PASE_9_b*36 + OA_T2$PASE_9_c*20 + OA_T2$PASE_9_d*35 + OA_T2$PASE_10tot*21 noOA_T1$PASE_T1 <- noOA_T1$PASE_2tot*20 + noOA_T1$PASE_3tot*21 + noOA_T1$PASE_4tot*23 + noOA_T1$PASE_5tot*23 + noOA_T1$PASE_6tot*30 + noOA_T1$PASE_7*25 + noOA_T1$PASE_8*25 + noOA_T1$PASE_9_a*30 + noOA_T1$PASE_9_b*36 + noOA_T1$PASE_9_c*20 + noOA_T1$PASE_9_d*35 + noOA_T1$PASE_10tot*21 # Mindset about the MPH of Health - Exercise OA_T1$MPH_T1 <- rowMeans(OA_T1[match(c("PES_1", "PES_2", "PES_3", "PES_4", "PES_5", "PES_6", "PES_7"), names(OA_T1))]) OA_T2$MPH_T2 <- rowMeans(OA_T2[match(c("PES_1", "PES_2", "PES_3", "PES_4", "PES_5", "PES_6", "PES_7"), names(OA_T2))]) noOA_T1$MPH_T1 <- rowMeans(noOA_T1[match(c("PES_1", "PES_2", "PES_3", "PES_4", "PES_5", "PES_6", "PES_7"), names(noOA_T1))]) # Global Health (higher better) ## Re-code Reverse Questions OA_T1$Global10R[OA_T1$Global10 == 1] <- 5 OA_T1$Global10R[OA_T1$Global10 == 2] <- 4 OA_T1$Global10R[OA_T1$Global10 == 3] <- 3 OA_T1$Global10R[OA_T1$Global10 == 4] <- 2 OA_T1$Global10R[OA_T1$Global10 == 5] <- 1 OA_T1$Global08R[OA_T1$Global08 == 1] <- 5 OA_T1$Global08R[OA_T1$Global08 == 2] <- 4 OA_T1$Global08R[OA_T1$Global08 == 3] <- 3 OA_T1$Global08R[OA_T1$Global08 == 4] <- 2 OA_T1$Global08R[OA_T1$Global08 == 5] <- 1 OA_T2$Global10R[OA_T2$Global10 == 1] <- 5 OA_T2$Global10R[OA_T2$Global10 == 2] <- 4 OA_T2$Global10R[OA_T2$Global10 == 3] <- 3 OA_T2$Global10R[OA_T2$Global10 == 4] <- 2 OA_T2$Global10R[OA_T2$Global10 == 5] <- 1 OA_T2$Global08R[OA_T2$Global08 == 1] <- 5 OA_T2$Global08R[OA_T2$Global08 == 2] <- 4 OA_T2$Global08R[OA_T2$Global08 == 3] <- 3 OA_T2$Global08R[OA_T2$Global08 == 4] <- 2 OA_T2$Global08R[OA_T2$Global08 == 5] <- 1 noOA_T1$Global10R[noOA_T1$Global10 == 1] <- 5 noOA_T1$Global10R[noOA_T1$Global10 == 2] <- 4 noOA_T1$Global10R[noOA_T1$Global10 == 3] <- 3 noOA_T1$Global10R[noOA_T1$Global10 == 4] <- 2 noOA_T1$Global10R[noOA_T1$Global10 == 5] <- 1 noOA_T1$Global08R[noOA_T1$Global08 == 1] <- 5 noOA_T1$Global08R[noOA_T1$Global08 == 2] <- 4 noOA_T1$Global08R[noOA_T1$Global08 == 3] <- 3 noOA_T1$Global08R[noOA_T1$Global08 == 4] <- 2 noOA_T1$Global08R[noOA_T1$Global08 == 5] <- 1 ## Re-code Pain Question OA_T1$Globalpain[OA_T1$Global07 == 1] <- 5 OA_T1$Globalpain[OA_T1$Global07 == 2] <- 4 OA_T1$Globalpain[OA_T1$Global07 == 3] <- 4 OA_T1$Globalpain[OA_T1$Global07 == 4] <- 4 OA_T1$Globalpain[OA_T1$Global07 == 5] <- 3 OA_T1$Globalpain[OA_T1$Global07 == 6] <- 3 OA_T1$Globalpain[OA_T1$Global07 == 7] <- 3 OA_T1$Globalpain[OA_T1$Global07 == 8] <- 2 OA_T1$Globalpain[OA_T1$Global07 == 9] <- 2 OA_T1$Globalpain[OA_T1$Global07 == 10] <- 2 OA_T1$Globalpain[OA_T1$Global07 == 11] <- 1 OA_T2$Globalpain[OA_T2$Global07 == 1] <- 5 OA_T2$Globalpain[OA_T2$Global07 == 2] <- 4 OA_T2$Globalpain[OA_T2$Global07 == 3] <- 4 OA_T2$Globalpain[OA_T2$Global07 == 4] <- 4 OA_T2$Globalpain[OA_T2$Global07 == 5] <- 3 OA_T2$Globalpain[OA_T2$Global07 == 6] <- 3 OA_T2$Globalpain[OA_T2$Global07 == 7] <- 3 OA_T2$Globalpain[OA_T2$Global07 == 8] <- 2 OA_T2$Globalpain[OA_T2$Global07 == 9] <- 2 OA_T2$Globalpain[OA_T2$Global07 == 10] <- 2 OA_T2$Globalpain[OA_T2$Global07 == 11] <- 1 noOA_T1$Globalpain[noOA_T1$Global07 == 1] <- 5 noOA_T1$Globalpain[noOA_T1$Global07 == 2] <- 4 noOA_T1$Globalpain[noOA_T1$Global07 == 3] <- 4 noOA_T1$Globalpain[noOA_T1$Global07 == 4] <- 4 noOA_T1$Globalpain[noOA_T1$Global07 == 5] <- 3 noOA_T1$Globalpain[noOA_T1$Global07 == 6] <- 3 noOA_T1$Globalpain[noOA_T1$Global07 == 7] <- 3 noOA_T1$Globalpain[noOA_T1$Global07 == 8] <- 2 noOA_T1$Globalpain[noOA_T1$Global07 == 9] <- 2 noOA_T1$Globalpain[noOA_T1$Global07 == 10] <- 2 noOA_T1$Globalpain[noOA_T1$Global07 == 11] <- 1 OA_T1$Global_T1 <- rowMeans(OA_T1[match(c("Global01", "Global02", "Global03", "Global04", "Global05", "Global06", "Global08R", "Global09", "Global10R", "Globalpain"), names(OA_T1))]) OA_T2$Global_T2 <- rowMeans(OA_T2[match(c("Global01", "Global02", "Global03", "Global04", "Global05", "Global06", "Global08R", "Global09", "Global10R", "Globalpain"), names(OA_T2))]) noOA_T1$Global_T1 <- rowMeans(noOA_T1[match(c("Global01", "Global02", "Global03", "Global04", "Global05", "Global06", "Global08R", "Global09", "Global10R", "Globalpain"), names(noOA_T1))]) # sex assigned at birth OA_T1$sex <- OA_T1$Dem_4a noOA_T1$sex <- noOA_T1$Dem_4a # Age OA_T1$age <- OA_T1$Dem_1 noOA_T1$age <- noOA_T1$Dem_1 # BMI OA_T1$BMI <- OA_T1$Dem_3/(((OA_T1$Dem_2_4)*12+OA_T1$Dem_2_5)^2)*703 noOA_T1$BMI <- noOA_T1$Dem_3/(((noOA_T1$Dem_2_4)*12+noOA_T1$Dem_2_5)^2)*703 # WOMAC # WOMAC is a total score that is calculated automatically in qualtrics OA_T1$WOMAC_T1 <- OA_T1$SC0 OA_T2$WOMAC_T2 <- OA_T2$SC0 # T2 label to MPH items OA_T2$PES_1_T2 <- OA_T2$PES_1 OA_T2$PES_2_T2 <- OA_T2$PES_2 OA_T2$PES_3_T2 <- OA_T2$PES_3 OA_T2$PES_4_T2 <- OA_T2$PES_4 OA_T2$PES_5_T2 <- OA_T2$PES_5 OA_T2$PES_6_T2 <- OA_T2$PES_6 OA_T2$PES_7_T2 <- OA_T2$PES_7 ``` # Simplify data frames ```{r} # Match aid case OA_T2$aid = tolower(OA_T2$aid) # Data frame of only data used in analyses dfOA_T1 <- OA_T1[, c("sex", "age", "BMI", "Global_T1", "PASE_T1", "MPH_T1", "PES_1", "PES_2" , "PES_3", "PES_4", "PES_5", "PES_6", "PES_7", "aid")] dfOA_T2 <- OA_T2[, c("Global_T2", "PASE_T2", "WOMAC_T2", "MPH_T2", "PES_1_T2", "PES_2_T2" , "PES_3_T2", "PES_4_T2", "PES_5_T2", "PES_6_T2", "PES_7_T2", "aid")] dfNoOA_T1 <- noOA_T1[, c("sex", "age", "BMI", "Global_T1", "PASE_T1", "MPH_T1", "aid")] ``` # Visualize the data ```{r} nOA_T1 = 1:nrow(dfOA_T1) nOA_T2 = 1:nrow(dfOA_T2) nNoOA_T1 = 1:nrow(dfNoOA_T1) par(mfrow=c(1,2)) plot(nOA_T1, dfOA_T1$age, xlab="n", ylab="Age", main="OA T1", pch=19) plot(nNoOA_T1, dfNoOA_T1$age, xlab="n", ylab="Age", main="No OA T1", pch=19) par(mfrow=c(1,2)) plot(nOA_T1, dfOA_T1$BMI, xlab="n", ylab="BMI", main="OA T1", pch=19) plot(nNoOA_T1, dfNoOA_T1$BMI, xlab="n", ylab="BMI", main="No OA T1", pch=19) par(mfrow=c(1,3)) plot(nOA_T1, dfOA_T1$Global_T1, xlab="n", ylab="Global_T1", main="OA T1", pch=19) plot(nNoOA_T1, dfNoOA_T1$Global_T1, xlab="n", ylab="Global_T1", main="No OA T1", pch=19) plot(nOA_T2, dfOA_T2$Global_T2, xlab="n", ylab="Global_T2", main="OA T2", pch=19) par(mfrow=c(1,3)) plot(nOA_T1, dfOA_T1$PASE_T1, xlab="n", ylab="PASE_T1", main="OA T1", pch=19) plot(nNoOA_T1, dfNoOA_T1$PASE_T1, xlab="n", ylab="PASE_T1", main="No OA T1", pch=19) plot(nOA_T2, dfOA_T2$PASE_T2, xlab="n", ylab="PASE_T2", main="OA T2", pch=19) par(mfrow=c(1,3)) plot(nOA_T1, dfOA_T1$MPH_T1, xlab="n", ylab="MPH_T1", main="OA T1", pch=19) plot(nNoOA_T1, dfNoOA_T1$MPH_T1, xlab="n", ylab="MPH_T1", main="No OA T1", pch=19) plot(nOA_T2, dfOA_T2$MPH_T2, xlab="n", ylab="MPH_T2", main="OA T2", pch=19) par(mfrow=c(1,2)) plot(nOA_T1, OA_T1$WOMAC_T1, xlab="n", ylab="WOMAC_T1", main="OA T1", pch=19) plot(nOA_T2, dfOA_T2$WOMAC_T2, xlab="n", ylab="WOMAC_T2", main="OA T2", pch=19) ``` # Descriptive Tables ```{r} describe(dfOA_T1) describe(dfOA_T2) describe(dfNoOA_T1) ``` # Standardized Mean Difference (between OA and non-OA) ```{r} dfOA_T1_SMD <- OA_T1[, c("sex", "age", "BMI", "Global_T1", "PASE_T1", "MPH_T1", "aid")] dfNoOA_T1_SMD <- noOA_T1[, c("sex", "age", "BMI", "Global_T1", "PASE_T1", "MPH_T1", "aid")] dfOA_T1_SMD$ifOA = 1 dfNoOA_T1_SMD$ifOA = 0 dfSMD = rbind(dfOA_T1_SMD, dfNoOA_T1_SMD) stddiff.numeric(data=dfSMD, gcol=8, vcol=(1:6)) ``` # Merge OA dataframes ```{r} # Add WOMAC dfOA_T1$WOMAC_T1 <- OA_T1$WOMAC_T1 # Add time point dfOA_T1$TP = 1 dfOA_T2$TP = 2 # Combine data frames based on aid df <- join(dfOA_T1, dfOA_T2, by="aid", type="inner") ``` # Standardized Mean Difference (between OA who answered at T2 and who did not) ```{r} # Participants not at T2 OA_NoT2_SMD <- OA_T1[!OA_T1$aid %in% OA_T2$aid, , drop = FALSE] OA_NoT2_SMD <- OA_NoT2_SMD[, c("sex", "age", "BMI", "Global_T1", "WOMAC_T1", "PASE_T1", "MPH_T1", "aid")] # Particiants at T2 OA_T2_SMD <- df[, c("sex", "age", "BMI", "Global_T1", "WOMAC_T1", "PASE_T1", "MPH_T1", "aid")] # Add time point OA_NoT2_SMD$TP = 1 OA_T2_SMD$TP = 2 dfSMD_TP = rbind(OA_NoT2_SMD, OA_T2_SMD) stddiff.numeric(data=dfSMD_TP, gcol=9, vcol=(1:7)) ``` # Correlation matrix with significance fuction ```{r, echo = FALSE} # x is a matrix containing the data # method : correlation method. "pearson"" or "spearman"" is supported # removeTriangle : remove upper or lower triangle # results : if "html" or "latex" # the results will be displayed in html or latex format corstars <-function(x, method=c("pearson", "spearman"), removeTriangle=c("upper", "lower"), result=c("none", "html", "latex")){ #Compute correlation matrix require(Hmisc) x <- as.matrix(x) correlation_matrix<-rcorr(x, type=method[1]) R <- correlation_matrix$r # Matrix of correlation coeficients p <- correlation_matrix$P # Matrix of p-value ## Define notions for significance levels; spacing is important. mystars <- ifelse(p < .0001, "****", ifelse(p < .001, "*** ", ifelse(p < .01, "** ", ifelse(p < .05, "* ", " ")))) ## trunctuate the correlation matrix to two decimal R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1] ## build a new matrix that includes the correlations with their apropriate stars Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x)) diag(Rnew) <- paste(diag(R), " ", sep="") rownames(Rnew) <- colnames(x) colnames(Rnew) <- paste(colnames(x), "", sep="") ## remove upper triangle of correlation matrix if(removeTriangle[1]=="upper"){ Rnew <- as.matrix(Rnew) Rnew[upper.tri(Rnew, diag = TRUE)] <- "" Rnew <- as.data.frame(Rnew) } ## remove lower triangle of correlation matrix else if(removeTriangle[1]=="lower"){ Rnew <- as.matrix(Rnew) Rnew[lower.tri(Rnew, diag = TRUE)] <- "" Rnew <- as.data.frame(Rnew) } ## remove last column and return the correlation matrix Rnew <- cbind(Rnew[1:length(Rnew)-1]) if (result[1]=="none") return(Rnew) else{ if(result[1]=="html") print(xtable(Rnew), type="html") else print(xtable(Rnew), type="latex") } } ``` # Correlations at T1 ```{r, echo = FALSE} dfOA_T1_corr <- data.frame(dfOA_T1$MPH_T1, dfOA_T1$Global_T1, dfOA_T1$WOMAC_T1, dfOA_T1$sex, dfOA_T1$age, dfOA_T1$BMI) colnames(dfOA_T1_corr) <- c( "MPH", "Global-10", "WOMAC", "sex", "age", "BMI") corstars(dfOA_T1_corr, method = "pearson") ``` # Correlations at T2 ```{r, echo = FALSE} dfOA_T2_corr <- data.frame(df$PASE_T2, df$PASE_T1, df$MPH_T1, df$Global_T1, df$WOMAC_T1, df$MPH_T2, df$Global_T2, df$WOMAC_T2, df$sex, df$age, df$BMI) colnames(dfOA_T2_corr) <- c( "PASE T2", "PASE T1", "MPH T1" , "Global T1", "WOMAC T1", "MPH T2", "Global T2", "WOMAC T2", "sex", "age", "BMI") corstars(dfOA_T2_corr, method = "pearson") ``` # Regression Analyses: MPH mindset ```{r, echo = FALSE} mod1 <- lm(dfSMD$MPH_T1 ~ scale(dfSMD$age) + dfSMD$sex + scale(dfSMD$BMI) + scale(dfSMD$Global_T1) + scale(dfSMD$PASE_T1) + dfSMD$ifOA) summ(mod1,confint = TRUE, digits = 3) ``` # Regression Analyses: Future physical activity ```{r, echo = FALSE} mod3 <- lm(df$PASE_T2 ~ scale(df$age) + df$sex + scale(df$BMI) + scale(df$Global_T1) + scale(df$WOMAC_T1) + scale(df$MPH_T1)) summ(mod3,confint = TRUE, digits = 3) res <- summary(mod3) ``` # Regression Analyses: Future physical activity, without outlier ```{r, echo = FALSE} df_outlier <- df[df$PASE_T2 < (mean(df$PASE_T2) + 3*sd(df$PASE_T2)),] mod4 <- lm(df_outlier$PASE_T2 ~ scale(df_outlier$age) + df_outlier$sex + scale(df_outlier$BMI) + scale(df_outlier$Global_T1) + scale(df_outlier$WOMAC_T1) + scale(df_outlier$MPH_T1)) summ(mod4,confint = TRUE, digits = 3) res <- summary(mod4) ``` ```{r, echo = FALSE} # Load Data management1 <- read.csv("Coding_Management_1.csv") # Load rater 1 coding results management2 <- read.csv("Coding_Management_2.csv") # Load rater 2 coding results management_final <- read.csv("Coding_Management_final.csv") # Load final decision ``` # Inter-rater Reliability: Cohen's Kappa for 2 Raters (Weights: unweighted) ```{r, echo = FALSE} library(irr) meds <- data.frame("rater1" = management1$Medication.or.Shots, "rater2" = management2$Medication.or.Shots) exercise <- data.frame("rater1" = management1$Exercise, "rater2" = management2$Exercise) K_meds = kappa2(meds[, c("rater1", "rater2")], weight = "unweighted") K_meds K_exercise = kappa2(exercise[, c("rater1", "rater2")], weight = "unweighted") K_exercise ``` # Column Sums ```{r, echo = FALSE} # Column Sums managementSums <- colSums(Filter(is.numeric, management_final)) managementSums ``` # Box Plots ```{r, echo = FALSE} # Code vs. Mindset library(ggplot2) library(ggpubr) library(gplots) med_only <- dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 1 & management_final[,6] == 0)]] med_ex <- dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 1 & management_final[,6] == 1)]] ex_only <- dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 0 & management_final[,6] == 1)]] diffs <- management_final[,2] - management_final[,6] a <- data.frame(group = "Med", value = med_only) b <- data.frame(group = "Med & ex", value = med_ex) c <- data.frame(group = "Ex", value = ex_only) plot.data <- rbind(a, b, c) plot.data$group <- factor(plot.data$group , levels=c("Med", "Med & ex", "Ex")) tiff("boxplot.tiff", units="in", width=7, height=5, res=300) ggplot(plot.data, aes(x=group, y=value, fill=group)) + geom_boxplot(fill=col2hex("cadetblue3")) + geom_jitter(color="black", size=1.0) + labs(x="Management Strategy", y = "MPH") + theme_classic(base_size = 15) dev.off() ``` # T-test for differences ```{r, echo = FALSE} # Medication vs. meds and exercise x <- t.test(dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 1 & management_final[,6] == 0)]], dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 1 & management_final[,6] == 1)]]) x x$conf.int # Meds and exercise vs. Exercise x <- t.test(dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 0 & management_final[,6] == 1)]], dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 1 & management_final[,6] == 1)]]) x x$conf.int # Meds vs. Exercise x <- t.test(dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 1 & management_final[,6] == 0)]], dfOA_T1$MPH_T1[row(management_final)[which(management_final[,2] == 0 & management_final[,6] == 1)]]) x x$conf.int ``` # MPH Factor Analysis, ICC, and and Pearson Correlation ```{r, echo = FALSE} cronbach.alpha(data.frame(OA_T1$PES_1, OA_T1$PES_2, OA_T1$PES_3, OA_T1$PES_4, OA_T1$PES_5, OA_T1$PES_6, OA_T1$PES_7), CI = TRUE,) cronbach.alpha(data.frame(OA_T2$PES_1, OA_T2$PES_2, OA_T2$PES_3, OA_T2$PES_4, OA_T2$PES_5, OA_T2$PES_6, OA_T2$PES_7), CI = TRUE,) df_icc <- data.frame(df$MPH_T1, df$MPH_T2) icc( df_icc, model = "twoway", type = "agreement", unit = "single" ) MPH <- data.frame(OA_T1$PES_1, OA_T1$PES_2, OA_T1$PES_3, OA_T1$PES_4, OA_T1$PES_5, OA_T1$PES_6, OA_T1$PES_7) n.factors <- 1 fit <- factanal(MPH, n.factors, # number of factors to extract rotation="varimax") # 'varimax' is an ortho rotation load <- fit$loadings[,1] # [,1:3] load scree.plot(fit$correlation) ev <- eigen(cor(MPH)) # get eigenvalues ev ap <- parallel(subject=nrow(MPH),var=ncol(MPH), rep=100, cent=.05) nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) plotnScree(nS) MPH_FA <- factanal(data.frame(OA_T1$PES_1, OA_T1$PES_2, OA_T1$PES_3, OA_T1$PES_4, OA_T1$PES_5, OA_T1$PES_6, OA_T1$PES_7), factors = 1) MPH_FA corPMI <- cor.test(df$MPH_T1, df$MPH_T2, method = "pearson") corPMI ``` # MPH Factor Correlations with PASE at T1 ```{r, echo = FALSE} df_MPH_corr <- data.frame(OA_T1$PASE_T1, OA_T1$PES_1, OA_T1$PES_2, OA_T1$PES_3, OA_T1$PES_4, OA_T1$PES_5, OA_T1$PES_6, OA_T1$PES_7) colnames(df_MPH_corr) <- c( "PASE T1", "Q1", "Q2" , "Q3", "Q4", "Q5", "Q6", "Q7") corstars(df_MPH_corr, method = "pearson") ``` # MPH Factor Correlations with PASE at T2 ```{r, echo = FALSE} df_MPH_corr <- data.frame(df$PASE_T2, df$PES_1, df$PES_2, df$PES_3, df$PES_4, df$PES_5, df$PES_6, df$PES_7) colnames(df_MPH_corr) <- c( "PASE T2", "Q1", "Q2" , "Q3", "Q4", "Q5", "Q6", "Q7") corstars(df_MPH_corr, method = "pearson") ```