chapter6.Rmd

#Analysis of longitudinal data

Dependencies:
```{r message=FALSE}
library(dplyr)
library(tidyr)
library(GGally)
library(ggplot2)
library(lme4)
```

Load BPRS data
```{r}
BPRS <- read.csv("https://raw.githubusercontent.com/KimmoVehkalahti/MABS/master/Examples/data/BPRS.txt", sep  =" ", header = T)
```


Factor treatment and subject:
```{r}
BPRS$treatment <- factor(BPRS$treatment)
BPRS$subject <- factor(BPRS$subject)
```


Convert to long form and extract the week number:
```{r}
BPRSL <- BPRS %>% gather(key = weeks, value = bprs, -treatment, -subject)
BPRSL <- BPRSL %>% mutate(week = as.integer(substr(weeks,5,5)))
```

```{r}
ggplot(BPRSL, aes(x = week, y = bprs, linetype = subject)) +
  geom_line() +
  scale_linetype_manual(values = rep(1:10, times=4)) +
  facet_grid(. ~ treatment, labeller = label_both) +
  theme(legend.position = "none") + 
  scale_y_continuous(limits = c(min(BPRSL$bprs), max(BPRSL$bprs)))
```


Standardize bprs:
```{r}
BPRSL <- BPRSL %>%
  group_by(week) %>%
  mutate(stdbprs = (bprs - mean(bprs))/sd(bprs) ) %>%
  ungroup()
```

Glimpse of the data:
```{r}
glimpse(BPRSL)
```

Plot with the standardised bprs:
```{r}
ggplot(BPRSL, aes(x = week, y = stdbprs, linetype = subject)) +
  geom_line() +
  scale_linetype_manual(values = rep(1:10, times=4)) +
  facet_grid(. ~ treatment, labeller = label_both) +
  scale_y_continuous(name = "standardized bprs")
```

Number of weeks, baseline (week 0) included
```{r}
n <- BPRSL$week %>% unique() %>% length()
```

Summary of the data with mean and standard error of bprs by treatment and week:
```{r}
BPRSS <- BPRSL %>%
  group_by(treatment, week) %>%
  summarise( mean = mean(bprs), se = sd(bprs)/sqrt(n) ) %>%
  ungroup()
```

Glimpse of the data
```{r}
glimpse(BPRSS)
```

Mean profiles plot:
```{r}
ggplot(BPRSS, aes(x = week, y = mean, linetype = treatment, shape = treatment)) +
  geom_line() +
  scale_linetype_manual(values = c(1,2)) +
  geom_point(size=3) +
  scale_shape_manual(values = c(1,2)) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se, linetype="1"), width=0.3) +
  theme(legend.position = c(0.8,0.8)) +
  scale_y_continuous(name = "mean(bprs) +/- se(bprs)")
```

Create a summary data by treatment and subject with mean as the summary variable.
```{r}
BPRSL8S <- BPRSL %>%
  filter(week > 0) %>%
  group_by(treatment, subject) %>%
  summarise( mean=mean(bprs) ) %>%
  ungroup()
```

Glimpse the data
```{r}
glimpse(BPRSL8S)
```

Draw a boxplot of the mean versus treatment
```{r}
ggplot(BPRSL8S, aes(x = treatment, y = mean)) + geom_boxplot()
```

Add the baseline from the original data as a new variable to the summary data
```{r}
BPRSL8S2 <- BPRSL8S %>%
  mutate(baseline = BPRS$week0)
```

Perform a two-sample t-test
```{r}
t.test(mean ~ treatment, data = BPRSL8S2, var.equal = TRUE)
```

Fit the linear model with the mean as the response
```{r}
fit <- lm(mean ~ baseline + treatment, data = BPRSL8S2)
```

Compute the analysis of variance table for the fitted model
```{r}
anova(fit)
```

Load RATSL data
```{r}
RATS <- read.table("https://raw.githubusercontent.com/KimmoVehkalahti/MABS/master/Examples/data/rats.txt", header = TRUE, sep = '\t')
```

Convert data to long form:
```{r}
RATSL <- RATS %>%
  gather(key = WD, value = Weight, -ID, -Group) %>%
  mutate(Time = as.integer(substr(WD,3,4))) 
```

Regression model of RATSL data:
```{r}
RATS_reg <- lm(Weight ~ Time + Group, data = RATSL)
```

Summary:
```{r}
summary(RATS_reg)
```

Random intercept model:
```{r}
RATS_ref <- lmer(Weight ~ Time + Group + (1 | ID), data = RATSL, REML = FALSE)
```

Summary of the model:
```{r}
summary(RATS_ref)
```

Random intercept model and random slope model:
```{r}
RATS_ref1 <- lmer(Weight ~ Time + Group + (Time | ID), data = RATSL, REML = FALSE)
```

Anova test on the models:
```{r}
anova(RATS_ref1, RATS_ref)
```