Confidence intervals with Lme4

# Silly example revisited: How high do people jump after 1 and 2 cups of coffee?

### Overall settings
participants <- 10
trialspercondition <- 10

### Set parameters for each condition
condition1m <- 50 # Mean Condition 1
condition2m <- 60 # Mean Condition 2
sd <- 10 # sd for both conditions
indvar <- 0.1 # Individual differences
trialerror <- 5 # How much noise to we have per trial?

### Make dataframe with participants and condition
sim <- matrix(NA,(participants*trialspercondition),3)
sim2 <- matrix(NA,(participants*trialspercondition),3)

### Loop over different conditions and participants
k = 0 # This is for increasing each round, emulating trials
for(i in 1:participants){

  #### Individual differences
  indv <- rnorm(1, 1, indvar) #
  for(j in 1:trialspercondition){
    sim[j+k,1] <- i
    sim[j+k,2] <- rnorm(1, condition1m*indv, sd) + rnorm(1,0,trialerror) # add little error
    sim[j+k,3] <- 1
    sim2[j+k,1] <- i
    sim2[j+k,2] <- rnorm(1, condition2m*indv, sd) + rnorm(1,0,trialerror) # add little error
    sim2[j+k,3] <- 2}
  k = k+trialspercondition
}

# Combine in data frame
sims <- data.frame(rbind(sim,sim2))
names(sims) <- c("Participant","Jump","Coffee")
sims$Participant <- factor(sims$Participant)
sims$Coffee <- factor(sims$Coffee)
head(sims)
summary(sims)

### Load lme4 which is a great place to start for doing repeated analysis
library(lme4)

# Make Null model
m0 <- lmer(Jump~(1|Participant), data= sims)
summary(m0)

# Make model with coffee
m1 <- lmer(Jump~ - 1+ Coffee + (1|Participant), data= sims)
summary(m1)

# Compare if a model with coffee is "worth the extra price" or a simpler model is preffered
anova(m0,m1,refit=FALSE)

# Calculate confidence interval (SHOULD CONTAIN THE TRUE VALUES OF SIMULATION MOST OF THE TIME)
ci95 <- confint(m1)
round(c(summary(m1)$coefficients[1,1],ci95[3,1:2]),2)
round(c(summary(m1)$coefficients[2,1],ci95[4,1:2]),2)

# Plot true vs captured
plot(as.factor(1:2),c(condition1m,condition2m), pch=19,xlim=c(0,3), ylim=c(40,70), lwd=3,
     main ="Comparing model to Truth", ylab="Height Jumped", xlab="Cups of Coffee") # TRUTH
points(1:2,c(summary(m1)$coefficients[,1]), pch=19, col="red") # Model Estimate
points(1:2,ci95[3:4,1], pch=6, col="red") # Lower
points(1:2,ci95[3:4,2], pch=2, col="red") # Higher