rms2-fev.Rmd

---
title: 'RMS Chapter 2 Example: Analysis of FEV Dataset'
author: "Frank Harrell"
date: "`r Sys.Date()`"
output:
  html_document:
    df_print: paged
    toc: yes
  html_notebook:
    highlight: textmate
    toc: yes
    toc_float:
      collapsed: yes
minor: curve fitting
major: regression modeling; rms
---
```{r setup,results='hide'}
require(Hmisc)
knitrSet(lang='markdown')
```
# Descriptive Statistics
```{r importdesc,w=7,h=6,results='asis'}
require(rms)
options(prType='html')
getHdata(FEV)
html(contents(FEV))
options(grType='plotly')  # use plotly interactive graphics
plot(describe(FEV))
require(plotly, quietly=TRUE)
ggplotly(ggplot(FEV, aes(x=age, y=fev, color=height)) + geom_point() + facet_grid(smoke ~ sex))
dd <- datadist(FEV); options(datadist='dd')
htmlVerbatim(dd)    # in Hmisc package
```

# Regression Models of Increasing Complexity
FEV is log-transformed to satisfy normality and equal variance assumptions.

## Single Predictor, Linear Effect
```{r lin,w=6.5,mfrow=c(1,2),results='asis'}
f <- ols(log(fev) ~ age, data=FEV)
f
r <- as.numeric(resid(f))
with(FEV, plot(age, r))
qqnorm(r)
qqline(r)
```
```{r linplot,results='asis'}
# Show algebraic form of fitted model
g <- Function(f)
htmlVerbatim(g, exp(g()), exp(g(age=10)))   # Evaluate the fitted model, antilog to get original scale
latex(f, md=TRUE)
ggplot(Predict(f), addlayer=geom_point(aes(x=age, y=log(fev), color=sex), FEV))
```

## Restricted Cubic Spline With 3 Default Knots
The knots are at these quantiles of age: 0.1, 0.5, 0.9.
```{r rcs3,results='asis'}
f <- ols(log(fev) ~ rcs(age, 3), data=FEV, x=TRUE)
print(f)
g <- Function(f)
htmlVerbatim(g(10:20), g(), g)
ggplot(Predict(f))
ha <- function(f) print(anova(f), dec.ms=2, dec.ss=2)
ha(f)
```

## RCS With 5 Default Knots
```{r rcs5,results='asis'}
f <- ols(log(fev) ~ rcs(age, 5), data=FEV)
f
Function(f)
ha(f)
ggplot(Predict(f))
```

## RCS with 5 Knots, Additive (Non-Interacting) Sex Effect
```{r rcs5sex,w=5,results='asis'}
f <- ols(log(fev) ~ rcs(age, 5) + sex, data=FEV)
Function(f)
ha(f)
ggplot(Predict(f, age, sex),
       addlayer=geom_point(aes(x=age, y=log(fev), color=sex), FEV))
```

## RCS in Age Interacting With Sex
The following model allows for different shapes of effects for males and females.
```{r rcsia,mfrow=c(1,2),w=6.5,results='asis'}
f <- ols(log(fev) ~ rcs(age, 5) * sex, data=FEV)
r <- as.numeric(resid(f))
plot(fitted(f), r)
qqnorm(r); qqline(r)
Function(f)
ha(f)
```
```{r rcsiap}
ggplot(Predict(f, age, sex))
```
Plot predicted median FEV by anti-logging predicted values.
```{r rcsiapm}
ggplot(Predict(f, age, sex, fun=exp), ylab='Estimated Median FEV')
```

## RCS with Age * Sex and Additive Nonlinear Effect of Height
We show the joint age and height effects using a color image plot
```{r rcsiah,results='asis'}
f <- ols(log(fev) ~ rcs(age, 5) * sex + rcs(height, 5), data=FEV)
ha(f)
ggplot(Predict(f, age, sex), adj.subtitle=FALSE)
```
```{r rcsbplot}
p <- Predict(f, age, height)
bplot(p)
```
```{r rcspl1}
p <- with(p, list(age=unique(age), height=unique(height),
                  Yhat=matrix(yhat, ncol=200)))
with(p, plot_ly(x=age, y=height, z=Yhat, type='heatmap'))
```
```{r rcspl2}
with(p, plot_ly(x=age, y=height, z=Yhat, type='surface'))
```

## Model Ignoring Sex and Height But Including Smoking History
By plotting spike histograms showing the smoking-specific age distribution we see that there are almost no very young children who have smoked.  This limits the power to detect an age by smoking interaction.
```{r amoke,w=5}
f <- ols(log(fev) ~ rcs(age, 5) + smoke, data=FEV)
ggplot(Predict(f, age, smoke), rdata=FEV)
```

# Computing Environment
```{r rsession,echo=FALSE}
si <- sessionInfo(); si$loadedOnly <- NULL
print(si, locale=FALSE)
```
```{r cite,results='asis',echo=FALSE}
print(citation(), style='text')
```