##############################################
# Analysis of Longitudinal Data (SuSe 2017)  #
#                                            #
# Solution for lab session 2                 #  
# Date: 18.05.2017                           #
##############################################


##############
# Exercise 1 #
##############

# Read data
setwd("../Vorher_online/")

library(nlme)
load("rats.long.RData")
str(rats.long)

# Transformation of TIME:
rats.long$logT     <- log(1+(rats.long$TIME-45)/10) 
# take out NAs
rats.long          <- rats.long[complete.cases(rats.long),]



## 1a)
# Fit a linear mixed model with a linear Trend in logT and 
# a subjekt-specific random intercept
lmm_RI   <- lme(RESPONSE ~ logT, random = ~1|SUBJECT, data=rats.long)
summary(lmm_RI)

# (i)
# --> theoretic considerations

# (ii)
# variance-covariance matrix of the random effects
getVarCov(lmm_RI)  

# --> This special case of an LMM (i.e. only with a random intercept) can
#     using the marginal point of view be interpreted as 'uniform
#     correlation model' (Proof was shown in lab). 
#     Note: The difference to the general 'uniform correlation model'
#     is that the correlations here are always positive whereas
#     in general they also can be negative.
#     We get:

(d.squared       <- getVarCov(lmm_RI)[1])
(sigma.squared   <- lmm_RI$sigma^2)
(cor.withinGroup <- d.squared/(d.squared + sigma.squared))

# --> The marginal correlation between two measurements of the same rat
#     is thus 0.73.

# (iii) 
# --> Theoretic considerations
#     Thus, the conditional correlation is zero (conditional independence!).

## 1c)

# subset of rats for rats with at least 3 measurements (with less measurements no (useful) model could be fitted)
obs           <- table(rats.long$SUBJECT)
toomanyNAs    <- levels(rats.long$SUBJECT)[obs < 3]
rats_enough          <- rats.long[!(rats.long$SUBJECT %in% toomanyNAs),]

# transform rats_enough to groupedData object
rats_enough<-groupedData(RESPONSE ~ logT|SUBJECT,data=rats_enough)

# fit seperate linear models with linear trend in logT
lm_separate<-lmList(RESPONSE~logT,data=rats_enough)
plot(intervals(lm_separate))

# --> confidence intervals for intercept overlap less, thus a random intercept is recommended,
#     but also a random slope might be neccesary as there are considerable variations in logT, too.

# (i)
# --> Correlation between measurements of the same subject is not accounted for. 
#     No possibility to obtain estimates for group effects.
#     Often there are not enough measurements.

# (ii)
lmm_RIRS    <- lme(RESPONSE ~ logT, random = ~ 1 + logT|SUBJECT, data=rats_enough)
summary(lmm_RIRS)

# Note: only 43 rats (R-output: Number of groups) left as those with less than 3 measurements are taken out.

# (iv)
# Get covariance matrix of random effects
(D  <-  getVarCov(lmm_RIRS))

# Correlation of random intercepts and random slopes
cov2cor(D)
# --> The correlation between random intercept and random slope 
#     is approximately 0.0002.



# (v)
plot(compareFits(coef(lm_separate), coef(lmm_RIRS)))

# --> The subject-specific difference in the trend are
#     (almost totally) removed in the linear mixed model 
#     ('shrinkage-effect').
#     Shrikage-effect: The subject-specific response
#     profiles in the linear mixed model are shrunk by
#     the normality assumption of the random effects towards
#     the global mean.
#     This is a stabilization of the parameter estimations
#     and can be seen as a compromise between a model
#     with a fixed effect for each rat and a model for the
#     pooled data where the longitudinal structure is
#     not taken into account.
#     Theoretical background: see lecture.

plot(comparePred(lm_separate, lmm_RIRS))
# --> The differences in the subject-specific regression
#     lines are often very clear. However, the estimated subject-
#     specific trend for the animals in the linear mixed model
#     is almost identical for all rats.


# also have a look at the residual variance
(lmm_RIRS$sigma^2)
D
# --> The residual variance is much larger than the variance of the 
#     random slope. Thus, the slope is shrunk to the population mean.
#     The variance of the random intercept is larger which explains
#     why it is not shrunk as much.