## ----setup, include=FALSE------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----init----------------------------------------------------------------
library(readr)
dat <- read_csv("C:/Users/dkxdk/OneDrive - Leo Pharma A S/stat/R course/dat_mmrm2.csv")

## ----dat structure-------------------------------------------------------
summary(dat)

## ----redness data--------------------------------------------------------
dat_red <- dat[dat$param == "Redness",]
head(dat_red)

## ----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()

## ----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 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)

## ----lme4----------------------------------------------------------------
library(lme4)
model1 <- lmer(aval ~ avisit*trt01p + (1 | subjid), data = dat_red)
summary(model1)
anova(model1)
AIC(model1)

## ----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)

summary(model2)
intervals(model2, which = "var-cov")
anova(model2)
AIC(model2)


#sighat <- extract.lme.cov(model2, data = sim)
#sighat2 <- extract.lme.cov2(model2, data = sim)

# model with unstructured covariance matrix does not converge
# model3 <- lme(aval ~ avisit*trt01p,
#               random= ~ 1|subjid,
#               correlation=corSymm(),
#               data=dat_red)




## ----valid---------------------------------------------------------------
plot(model2)
plot(model2, resid(., type="p") ~ fitted(.) | trt01p)
qqnorm(model2)

## ----diff lsmeans--------------------------------------------------------
library(emmeans)
emm <- emmeans(model2, pairwise ~ trt01p| avisit)
emm$contrasts
emm_options()