Fract.Poly: Fit fractional polynomials

Description Usage Arguments Value Author(s) References Examples

View source: R/Fract.Poly.R

Description

Fits regression models with m terms of the form X^{p}, where the exponents p are selected from a small predefined set S of both integer and non-integer values.

Usage

1
Fract.Poly(Covariate, Outcome, S=c(-2,-1,-0.5,0,0.5,1,2,3), Max.M=5, Dataset)

Arguments

Covariate

The covariate to be considered in the models.

Outcome

The outcome to be considered in the models.

S

The set S from which each power p^{m} is selected. Default S={-2, -1, -0.5, 0, 0.5, 1, 2, 3}.

Max.M

The maximum order M to be considered for the fractional polynomial. This value can be 5 at most. When M=5, then fractional polynomials of order 1 to 5 are considered. Default Max.M=5.

Dataset

A data.frame that should consist of multiple lines per subject ('long' format).

Value

Results.M1

The results (powers and AIC values) of the fractional polynomials of order 1.

Results.M2

The results (powers and AIC values) of the fractional polynomials of order 2.

Results.M3

The results (powers and AIC values) of the fractional polynomials of order 3.

Results.M4

The results (powers and AIC values) of the fractional polynomials of order 4.

Results.M5

The results (powers and AIC values) of the fractional polynomials of order 5.

Author(s)

Wim Van der Elst, Geert Molenberghs, Ralf-Dieter Hilgers, & Nicole Heussen

References

Van der Elst, W., Molenberghs, G., Hilgers, R., & Heussen, N. (2015). Estimating the reliability of repeatedly measured endpoints based on linear mixed-effects models. A tutorial. Submitted.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
# Open data
data(Example.Data)

# Fit fractional polynomials, mox. order = 3
FP <- Fract.Poly(Covariate = Time, Outcome = Outcome, 
Dataset = Example.Data, Max.M=3)

# Explore results
summary(FP)
# best fitting model (based on AIC) for m=3, 
# powers:  p_{1}=3,  p_{2}=3, and  p_{3}=2

# Fit model and compare with observed means

    # plot of mean 
Spaghetti.Plot(Dataset = Example.Data, Outcome = Outcome,
Time = Time, Id=Id, Add.Profiles = FALSE, Lwd.Me=1, 
ylab="Mean Outcome")

    # Coding of predictors (note that when p_{1}=p_{2},
    # beta_{1}*X ** {p_{1}} + beta_{2}*X ** {p_{1}} * log(X)
    #  and when p=0, X ** {0}= log(X)    )
term1 <- Example.Data$Time**3
term2 <- (Example.Data$Time**3) * log(Example.Data$Time) 
term3 <- Example.Data$Time**2 

    # fit model
Model <- lm(Outcome~term1+term2+term3, data=Example.Data)
Model$coef # regression weights (beta's)  

    # make prediction for time 1 to 47
term1 <- (1:47)**3 
term2 <- ((1:47)**3) * log(1:47)   
term3 <- (1:47)**2 
    # compute predicted values
pred <- Model$coef[1] + (Model$coef[2] * term1) + 
  (Model$coef[3] * term2) + (Model$coef[4] * term3)
   # Add predicted values to plot
lines(x = 1:47, y=pred,  lty=2)
legend("topright", c("Observed", "Predicted"), lty=c(1, 2))

Example output

Loading required package: nlme
Best fitting model for m=1: 
  power1      AIC 
  -0.500 2715.026 

Best fitting model for m=2: 
  power1   power2      AIC 
  -2.000   -2.000 2714.357 

Best fitting model for m=3: 
  power1   power2   power3      AIC 
   2.000    3.000    3.000 2702.398 
 (Intercept)        term1        term2        term3 
147.80836947   0.08290322  -0.01802568  -0.67134790 

CorrMixed documentation built on May 2, 2019, 3:26 p.m.