> library(foreign)
> ds <- read.dta("c:/diabetes.dta")
> attach(ds)
> bf1 <- (age-5)*I(age > 5)
> bf2 <- (age-6)*I(age > 6)
> bf3 <- (age-7)*I(age > 7)
> bf4 <- (age-8)*I(age > 8)
> bf5 <- (age-9)*I(age > 9)
> bf6 <- (age-10)*I(age > 10)
> bf7 <- (age-11)*I(age > 11)
> bf8 <- (age-12)*I(age > 12)
> bf9 <- (age-13)*I(age > 13)
> bf10 <- (age-14)*I(age > 14)
> library(nlme)
> Const <- factor(rep(1,length(logc)))
> model1 <- lme(logc ~ age, random=list(Const=pdIdent(~-1 + bf1 + bf2 +
+ bf3 + bf4 + bf5 + bf6 + bf7 + bf8 + bf9 + bf10)))
> summary(model1)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
-99.15086 -92.29657 53.57543
Random effects:
Formula: ~-1 + bf1 + bf2 + bf3 + bf4 + bf5 + bf6 + bf7 + bf8 + bf9 + bf10 | Const
Structure: Multiple of an Identity
bf1 bf2 bf3 bf4 bf5 bf6
StdDev: 0.01020370 0.01020370 0.01020370 0.01020370 0.01020370 0.01020370
bf7 bf8 bf9 bf10 Residual
StdDev: 0.01020370 0.01020370 0.01020370 0.01020370 0.05540213
Fixed effects: logc ~ age
Value Std.Error DF t-value p-value
(Intercept) 0.5246126 0.03173165 41 16.532786 0e+00
age 0.0268907 0.00692908 41 3.880842 4e-04
Correlation:
(Intr)
age -0.883
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.83142861 -0.40195966 0.03754746 0.45272629 2.37011118
Number of Observations: 43
Number of Groups: 1
> ranef(model1)
bf1 bf2 bf3 bf4 bf5 bf6
1 -0.007727653 -0.00820119 -0.006290257 -0.003866588 -0.001573929 -0.001148187
bf7 bf8 bf9 bf10
1 -0.0004626088 0.0006163229 0.001050489 0.001375959
> predicted <- predict(model1)
> plot(age[order(age)], predicted[order(age)], type="l", xlab="Age (in years)",
+ ylab="Predicted Log C-Peptide Concentation",
+ main="Predicted Log C-Peptide Concentration versus Age")