---
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")
```