Chapter 19, Section 19.3

 

Study of Log C-Peptide Concentration in Children with Diabetes

Mixed Model Representation of Penalized Spline for Log C-Peptide Concentration

 

 

 

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")