---
title: "Second part R course"
author: "Alexandra Kuznetsova"
date: "19 nov 2018"
output:
html_document: default
pdf_document: default
word_document: default
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Data for the second part
The aim of the (hypothetical) trial was to investigate the effect of three different treatments in terms of thickness and redness of the skin in psoriasis patients.
There were 10 subjects
per treatment arm and the effectiveness of the drug was measured in terms of redness and thickness of the skin. Thickness of the skin in the lesion was measured with a special device, whereas redness was evaluated by an investigator: a score ranging from 0 to 30. The data has a form of ADaM data ADXD with paramcd having two values: "Redness", "Thickness" and aval contains the values of these parameters. ADY variable contains days and AVISIT contains visit days in a character format.
## Import and plotting the data
```{r init}
library(readr)
dat <- read_csv("C:/Users/dkxdk/OneDrive - Leo Pharma A S/stat/R course/dat_mmrm2.csv")
str(dat)
```
Now we check the data structure and the variable formats.
```{r dat structure}
str(dat)
```
We are going to look at the redness of the skin. First we subset the data.
```{r redness data}
dat_red <- dat[dat$param == "Redness",]
head(dat_red)
```
THen we plot redness scores versus days for each subject, colouring by treatment group :
```{r plot}
library(ggplot2)
dat_red$subjid<- as.factor(dat_red$subjid)
dat_red$trt01p <- as.factor(dat_red$trt01p)
ggplot(dat_red, aes(x=ADY, y=aval, group=subjid, colour=trt01p)) +
geom_line()
```
It is also of interest to plot treatment group average time profiles:
```{r plot average}
require(plyr)
dat_red$avisit <- as.factor(dat_red$avisit)
mns <- ddply(dat_red, ~ trt01p + avisit + ADY, summarize,
redness = mean(aval))
ggplot(mns, aes(x=ADY, y=redness, group=trt01p, colour=trt01p)) +
geom_point() + geom_line()
```
Plot data points with a fitted curve:
```{r plot smooth}
ggplot(dat_red, aes(x=ADY, y=aval, colour = trt01p, group = trt01p)) +
geom_point(shape=1) + # Use hollow circles
geom_smooth(method = "lm", formula = y ~ splines::ns(x, 2), se = FALSE)
```
It seems like there is a difference in treatment groups. It is also clear that there is a between subjects variation.
## Modelling of the data
First we fit data with a linear mixed model with one random subject effect by using **lme4** package.
The covariance matrix $V$ has the following form:
\[ V_{y_{i_1}, y_{i_2}} =
\begin{cases}
0 & \quad \text{if } subj_{i1} \neq subj_{i2} \text{ and } i_1 \neq i_2\\
\sigma_{subj}^2 & \quad \text{if } subj_{i1} = subj_{i2} \text{ and } i_1 \neq i_2 \\
\sigma_{subj}^2+ \sigma^2 & \quad \text{if } i_1 = i_2
\end{cases}
\]
In this model two measurements on the subject are correlated, but equally correlated (no matter how far observations are taken).
This type of model can be fit in R in different ways. Using the **nlme** or **lme4** packages. With **lme4**:
```{r lme4}
library(lme4)
model1 <- lmer(aval ~ avisit*trt01p + (1 | subjid), data = dat_red)
summary(model1)
anova(model1)
```
Now we will try to fit mixed model with repeated measurements:
$$redness \sim N(\mu, V)$$
with
$$\mu_i = \mu + \alpha(treat_i) + \beta(day_i) + \gamma(treat_i, day_i)$$
$$\text{Redness}_{ijk} = \mu + \alpha_i + \beta_j + \gamma_{ij} + \delta_{k} + \epsilon_{ijk}$$
where $\alpha$ stands for treatment effect, $$\beta$$ for day factor, $$\gamma$$ - treatment -by-day interaction, $$\delta$$ is a random subject effect
and covariance $V$ having the following form:
\[ V_{i_1, i_2} =
\begin{cases}
0 & \quad \text{if } subj_{i1} \neq subj_{i2} \text{ and } i_1 \neq i_2\\
\nu^2 + \tau^2*exp(\frac{-(day_{i_1} - day_{i_2})^2}{\rho^2}) & \quad \text{if } subj_{i1} = subj_{i2} \text{ and } i_1 \neq i_2 \\
\nu^2 + \tau^2 + \sigma^2 & \quad \text{if } i_1 = i_2
\end{cases}
\]
In this model two observations very close to each other have covariance $\nu^2 + \tau^2$, two observations very far have covariance $\nu^2$
This type of model can be fit with **nlme** package
```{r nlme ancova}
library(nlme)
?corClasses
str(dat_red)
model2 <- lme(aval ~ avisit*trt01p,
random= ~ 1|subjid,
correlation=corGaus(form= ~ ADY|subjid, nugget=TRUE),
data=dat_red)
model2
anova(model2)
#sighat <- extract.lme.cov(model2, data = sim)
#sighat2 <- extract.lme.cov2(model2, data = sim)
# model with unstructured covariance structure does not converge
# model3 <- lme(aval ~ avisit*trt01p,
# random= ~ 1|subjid,
# correlation=corSymm(),
# data=dat_red)
```
Validation plots
```{r valid}
plot(model2)
plot(model2, resid(., type="p") ~ fitted(.) | trt01p)
qqnorm(model2)
```
Differences of least squares means at Day 8
```{r diff lsmeans}
library(emmeans)
emm <- emmeans(model2, pairwise ~ trt01p| avisit)
emm$contrasts
emm_options()
```