CorrMixed (version 0.1-1)

Fract.Poly: Fit fractional polynomials

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

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.M1The results (powers and AIC values) of the fractional polynomials of order $1$.
  • Results.M2The results (powers and AIC values) of the fractional polynomials of order $2$.
  • Results.M3The results (powers and AIC values) of the fractional polynomials of order $3$.
  • Results.M4The results (powers and AIC values) of the fractional polynomials of order $4$.
  • Results.M5The results (powers and AIC values) of the fractional polynomials of order $5$.

References

Van der Elst, W., Molenberghs, G., Hilgers, R., & Heussen, N. (2015). Correlation in continuous monitoring of vital parameters - estimating reliability using linear mixed-effects models. Submitted.

Examples

Run this code
# 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))

Run the code above in your browser using DataCamp Workspace