library(nlme)


################################################################################
# Exercise 2
# Test for fixed and random effects
################################################################################

### Simulate data ###
set.seed(123)
Ji <- rep(6,30)
n <- sum(Ji)
I <- length(Ji)
ID <- rep(1:I, Ji)
gender <- 1*(ID>15)
dosis <- rep(c(10, 20, 40, 80, 120, 160), I)
RI <- rnorm(I, 10)
RS <- rnorm(I, sd = 0.09) 
eps <- rnorm(n, sd = 2)

fun <- function(x){
  -2*sqrt(x) - x/20
}

SBD <- 150 + fun(dosis) - 2 * gender - dosis * gender/100 + 
  RI[ID] + RS[ID] * dosis + eps


####### as an alternative, simulate other data to see greater p-values 
## gender, dosis do not interact
## very small effect of dosis
## very samll random slopes
if(FALSE){
  SBD <- 150 + 0.01 * fun(dosis) - gender - 0 * dosis * gender/100 + 
    RI[ID] + 0.1 * RS[ID] * dosis + eps
}


blutdruck <- data.frame(ID=ID, gender = gender, dosis = dosis, SBD = SBD)


###### (2a)
# create dummy-variables 
blutdruck$genderFemaleDosis <- 0
blutdruck$genderFemaleDosis[blutdruck$gender==0] <- (blutdruck$dosis)[blutdruck$gender==0]

blutdruck$genderMaleDosis <- 0
blutdruck$genderMaleDosis[blutdruck$gender==1] <- (blutdruck$dosis)[blutdruck$gender==1]

# treat gender as factor
blutdruck$gender <- factor(blutdruck$gender)

# Model with random intercept and different slope for male and female  
fm1 <- lme(SBD ~ -1 + gender + genderFemaleDosis + genderMaleDosis, random=~1|ID, data=blutdruck)
summary(fm1) 


###### (2b)
# H0: dosis has no influence for males or females
# use Wald-test, as LQ-Test not possible as estimation was done using REML
# Calculate the p-value of the Wald-test for the effect of genderFemaleDosis and genderMaleDosis
hatbeta <- fixef(fm1)[3:4]
tw <- hatbeta %*% solve(vcov(fm1)[3:4,3:4]) %*% hatbeta
1 - pchisq(tw, df=2)
# reject H0 on 5% level

# EXTRA: using anova() on one model gives F-test of effects for single variables
# see ?anova.lme
anova(fm1)

### do a likelihood-ratio test (LRT), models have to be fitted using ML
fm1ML1 <- lme(SBD ~ -1 + gender + genderFemaleDosis + genderMaleDosis, random=~1|ID, data=blutdruck, method="ML")
fm1ML0 <- lme(SBD ~ -1 + gender, random=~1|ID, data=blutdruck, method="ML")
# LRT for models fitted by ML
anova(fm1ML0, fm1ML1)

## or LRT by hand
tLQT <- 2*fm1ML1$logLik - 2*fm1ML0$logLik
1 - pchisq(tLQT, df=2)

### try LQT for models fitted with REML
fm1REML0 <- lme(SBD ~ -1 + gender, random=~1|ID, data=blutdruck, method="REML")
anova(fm1, fm1REML0) ## gives a warning, as it ist not meaningful!

#### use package pbkrtest, for F-test and for parametric bootstrap 
library(pbkrtest)

### (1) use F-test with approximated df
# An approximate F-test based on the Kenward-Roger approach.
?KRmodcomp

## you have to fit the models using lmer to use KRmodcomp()
library(lme4)
fm1lmer <- lmer(SBD ~ -1 + gender + genderFemaleDosis + genderMaleDosis + (1|ID), data=blutdruck)
fm0lmer <- lmer(SBD ~ -1 + gender + (1|ID), data=blutdruck)

# F-test with approximated df, based on Kenward-Roger approximation
KRmodcomp(fm1lmer, fm0lmer)

### (2) use parametric bootstrap methods
# (takes about 1 minute)
PBmodcomp(fm1lmer, fm0lmer)



###### (2c)
# H0: The slope in dose is the same for females and males
# Calculate the p-value of the Wald-test for the hypothesis 
L <- t(c(1,-1)) # matrix with one row
hatbeta <- fixef(fm1)[3:4]
tw2 <- t(hatbeta)%*%t(L) %*% solve(L%*%vcov(fm1)[3:4,3:4]%*%t(L)) %*% L%*%hatbeta
1 - pchisq(tw2, df=2)
# do not reject H0 on 5% level

###### (2d)
# Test on variance component of random intercepts:
# H0: var(b_{0i})=\tau^2=0

### exact test of random intercepts 
library(RLRsim)
?exactRLRT
# it is enough to specify m, as only one variance component is tested
exactRLRT(m = fm1) 
# reject H0 on 5% level


### Test statistic to approximate p-values
# model under H0
fm1Fix <- lm(SBD ~ -1 + gender + genderFemaleDosis + genderMaleDosis, data=blutdruck)
# approximate likelihood-ratio test 
T_RLRTc <- 2 * logLik(fm1, REML = TRUE)[1] - 2 * logLik(fm1Fix, REML = TRUE)[1]
# approximate p-value with mixture of chi^2 distributions, q=1
0.5 *(1-pchisq(T_RLRTc, 1)) + 0.5 *(1-pchisq(T_RLRTc, 0)) 
# reject H0 on 5% level

###### (2e)
# compute conditional Akaike information (cAIC)
# do the model fit using lme4
library(lme4)
# model from (a)
fm1lmer <- lmer(SBD ~ -1 + gender + genderFemaleDosis + genderMaleDosis + (1|ID), data=blutdruck)
# model from (e) with random slopes 
fm2lmer <- lmer(SBD ~ -1 + gender + genderFemaleDosis + genderMaleDosis + (1 + dosis|ID), data=blutdruck)
# ## try another optimizer to check results, see
# # ?lmerControl
# fm2lmerB <- update(fm2lmer, control=lmerControl(optimizer="Nelder_Mead"))
# require(optimx)
# fm2lmerC <- update(fm2lmer, control=lmerControl(optimizer="optimx", optCtrl=list(method="nlminb")))


# compute the cAIC
library(cAIC4)
?cAIC
(c1 <- cAIC(fm1lmer, method="steinian"))
(c2 <- cAIC(fm2lmer, method="steinian"))

## preferred model is the one with the minimum AIC value
c1$caic
c2$caic

###### (2f)
# Test on variance component of random slope

# Model with random intercept (m0, model under H0)
fm1 <- lme(SBD ~ -1 + gender + genderFemaleDosis + genderMaleDosis, random=~1|ID, data=blutdruck)
summary(fm1)

# Model with random intercept and random slope (mA, model under alternative)
fm2 <- lme(SBD ~ -1 + gender + genderFemaleDosis + genderMaleDosis, random=~1+dosis|ID, data=blutdruck)
summary(fm2)

# Test statistic to approximate p-values
T_RLRTd <- 2 * logLik(fm2, REML = TRUE)[1] - 2 * logLik(fm1, REML = TRUE)[1]
# approximate p-value with mixture of chi^2 distributions, 
# q=2 as the variance component of random slopes 
# and the covariance between random intercept and random slope are tested to be 0
0.5 *(1-pchisq(T_RLRTd, 2)) + 0.5 *(1-pchisq(T_RLRTd, 1)) 
# reject H0 on 5% level