lme()
in the nlme
package can do this for you. We need to fit a basic linear mixed effects model, with a fixed effect for the intercept (i.e. mean), as well as a subject-level random effect for this intercept. Once the model is fit, the confidence intervals can be obtained with the lme method for the generic function intervals()
. Here is an example:
require(ggplot2)
require(nlme)
n = 12
# Generate data table
set.seed(12345)
ID = factor(rep(1:n, len=2*n))
value = 5 + rnorm(n)[as.numeric(ID)] + rnorm(2*n)
data = data.frame(ID, value)
# Plot data
dev.new(width=4, height=4)
qplot(ID, value, data=data, geom='boxplot')
# Fit linear mixed-effects model with an intercept term, as well as a random
# subject-level effect on intercept
fit = lme(value ~ 1, random = ~ 1 | ID, data=data)
summary(fit)
intervals(fit)
Raw data:
> summary(fit)
Linear mixed-effects model fit by REML
Data: data
AIC BIC logLik
80.90274 84.30923 -37.45137
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 0.9126617 0.871738
Fixed effects: value ~ 1
Value Std.Error DF t-value p-value
(Intercept) 5.423818 0.3179249 12 17.06006 0
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.29134118 -0.78398894 -0.09942166 0.66779227 1.48250943
Number of Observations: 24
Number of Groups: 12
> intervals(fit)
Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 4.731119 5.423818 6.116517
attr(,"label")
[1] "Fixed effects:"
Random Effects:
Level: ID
lower est. upper
sd((Intercept)) 0.4827668 0.9126617 1.72537
Within-group standard error:
lower est. upper
0.5841779 0.8717380 1.3008490