##############################################
# Analysis of Longitudinal Data (SuSe 2017)  #
#                                            #
# Solution for lab session 3                 #  
# Date: 30.05.2017                           #
##############################################



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

## 1a)

# Load package nlme
library(nlme)

# First glance at the data 
?Orthodont
str(Orthodont)
getGroups(Orthodont)

# --> Every two years, the distance of two spots at the jaw was measured using x-ray.
# --> nfnGroupedData: numeric covariate, single level of grouping
#     (here: each child is a group)
# --> Other possibilities:
#     nffGroupedData: factor covariate, single level of nesting
#     nmGroupedData: multiple levels of nesting
#     see: ?groupedData

## 1b)
# Transformation of age
Orthodont$newage<-Orthodont$age-11
mean(Orthodont$newage)

# --> By this transformation, the data are centered and thus the mean of the transformed data is zero. 
#     This has the consequence that the intercept can be interpreted as the value for the mean age (here 11 years).

## 1c)
# Model with random intercept
m_RI <- lme(distance ~ newage + Sex,random = ~1|Subject,data=Orthodont)
summary(m_RI)

## 1d)
#i) Estimation with ML instead of REML
m_RI_ml <- update(m_RI,method="ML")
summary(m_RI_ml)

#ii)
# Compare estimated covariance matrix of the random effects (ML and REML)
getVarCov(m_RI)
getVarCov(m_RI_ml)

# --> The REML estimation variance is larger than that of resulting 
#     from the ML estimation which is biased downwards.

#iii)
# Compare the estimated fixed effects (ML and REML)
fixef(m_RI)
fixef(m_RI_ml)

# --> Estimated fixed effects are identical. This is not always the case!

## 1e)
# Add interaction of age and sex of the children
m_RI_int<-update(m_RI,distance~ newage*Sex)
  
#i) 
# --> Yes, as it can be assumed that the effect of sex plays a greater role when children grow older.

#ii)
# --> The likelihood ratio test cannot be used with REML estimation.
#     Reason: Likelihood is based on error contrast which depend on X.
#     As here, fixed effects differ: REML likelihoods are not comparable.
#     Per default, summary.lme() gives out the results of a approximative t-test which can be used for ML and REML.
#     Other alternatives: approximative Wald-test or approximative F-test (for general hypotheses).

## 1f)

#i)
# Marginal modell with same fixed effects with unstructured correlation structure  
m_marg_unstruc<-gls(distance ~ newage*Sex,correlation=corSymm(form = ~1|Subject),
                      data=Orthodont)
summary(m_marg_unstruc)

# Correlation of the marginal model
(cor_marg_unstruc<-cov2cor(getVarCov(m_marg_unstruc)))

# Marginal correlation of the lme 
# Compare lecture slides on compound symmetry

# --> the correlation is d^2/(d^2+\sigma^2)
# extract variance of random intercepts
var_b<-getVarCov(m_RI_int)[1]

(cor_marg_lme<-var_b/(var_b+m_RI_int$sigma^2))

# --> Marginal correlation of two measurements of the same child in m_RI_int is always the same. 
#     This is because we only have random intercepts and assume conditional independence -> compound symmetry.
#     The entries in correlation matrix (4 measurements per child!) of the marginal model
#     are similar to the marginal correlation of the lme.
#     -> it is thus reasonable to replace the unstructured correlation by a simpler form.

#ii)
# Marginal model with simpler correlation structure 
m_marg_const<-update(m_marg_unstruc,correlation=corCompSymm(form=~1|Subject))
summary(m_marg_const)

# Correlation
(cov2cor(getVarCov(m_marg_const)))

#iii)
# --> No, as this is only the marginal perspective of the linear mixed model.
#     The lmm implies a marginal model but the reverse is not true.
#     The fixed effects have the same interpretation (in the case of LMM NOT GLMM!).
#     Random intercept model with conditional independence can be seen as a special 
#     case of the marginal uniform correlation model (from the marginal perspective).
#     But: uniform correlation model is more general and negative correlations are allowed.