In Brief

The mosaic package provides facilities for putting model fitting into a “mathematical function” framework.

Use with lm( ), glm( )

Example:

Mortality at 20 years of follow-up for smokers and non-smokers.

sample( Whickham, size=5 )
##      outcome smoker age orig.ids
## 665    Alive     No  22      665
## 781    Alive     No  32      781
## 1062    Dead     No  67     1062
## 956     Dead     No  72      956
## 200    Alive     No  51      200
smod <- glm( outcome=='Alive' ~ age*smoker, family='binomial',
             data=Whickham)
f <- makeFun(smod) # EXTRACT THE FUNCTION

Evaluation … For a 60-year old …

f( age=60, smoker=c('Yes'))
##      1 
## 0.5129

Plot data and the functional form

xyplot( outcome=='Alive' ~ age | smoker, data=Whickham,
        pch=20, alpha=.2 )

plotFun( f(age=age, smoker="No") ~ age, add=TRUE)
plotFun( f(age=age, smoker="Yes") ~ age, add=TRUE, 
         col='red', lwd=2)

plot of chunk unnamed-chunk-6

lm()

hmod <- lm( height ~ father, data=subset(Galton,sex=="M"))

Extract the function … just like resid( ) or fitted( )

hfun <- makeFun(hmod)

Plot it … first data, then the function.

xyplot(height ~ father | sex, data=Galton)
plotFun( hfun(x) ~ x, add=TRUE)

plot of chunk unnamed-chunk-9

Multiple variables:

mod2 <- lm( height ~ father + sex, data=Galton)
hfun2 <- makeFun(mod2)
xyplot(height ~ father | sex, data=Galton)
plotFun( hfun2(x, sex="F") ~ x, add=TRUE)
plotFun( hfun2(x, sex="M") ~ x, add=TRUE, col='red')

plot of chunk unnamed-chunk-10

Or …

mod3 <- lm( height ~ father*mother + sex, data=Galton)
hfun3 <- makeFun(mod3)
plotFun( hfun3(father=x, mother=y, sex="F") ~ x&y,
         x.lim=c(60,80),y.lim=c(57,73)) 
plotPoints( mother ~ father, data=Galton, add=TRUE,
            pch=20,col='red')

plot of chunk unnamed-chunk-11

Confidence Intervals

hfun2( father=66, sex=c('M','F'), interval='confidence')
##     fit   lwr   upr
## 1 67.87 67.59 68.16
## 2 62.70 62.40 62.99
hfun2( father=66, sex=c('M','F'), interval='prediction')
##     fit   lwr   upr
## 1 67.87 63.40 72.35
## 2 62.70 58.22 67.18

NOTE: ggplot( ) does this well.

ggplot(data=Galton, aes(x=father, y=height)) + geom_point()  + 
  aes(colour=sex)  + stat_smooth(method=lm) + 
  theme(legend.position="none") + labs(title="") 

plot of chunk unnamed-chunk-13

Convenient interface via mosaic:mPlot().