Less Volume, More Creativity

The new support for documentation creation in RStudio is great.

• These slides are HTML, but I created them in RMarkdown (+ a little bit of HTML fiddling)

• Other slides will be PDF created using beamer – but again, I wrote them in RMarkdown

• A single RMarkdown file can generate PDF, HTML, or Word

  • no need to know HTML, LateX or Word
  • but if you do, you can take advantage

One key to successfully introducing R is finding a set of commands that is

• small: fewer is better
• coherent: commands should be as similar as possible
• powerful: can do what needs doing

It is not enough to use R, it must be used elegantly.

Two examples of this principle:

• the mosaic package
• the dplyr package (Hadley Wickham)

Goal: a minimal set of R commands for Intro Stats

Result: Minimal R Vignette

Much of the work on the mosaic package has been motivated by

• The Less Volume, More Creativity approach
• The Minimal R goal

R is case sensitive

• many students are not case sensitive

Arrows and Tab

• up/down arrows scroll through history
• TAB completion can simplify typing

If all else fails, try ESC

• If you see a + prompt, it means R is waiting for more input
• If this is unintentional, you probably have a typo
• ESC will get you pack to the command prompt

Other versions:

# simpler version
goal( ~ x, data = mydata )          
# fancier version
goal( y ~ x | z , data = mydata )   
# unified version
goal( formula , data = mydata )     

What do you want R to do? (goal)

What must R know to do that?

What do you want R to do? (goal)

• This determines the function to use

What must R know to do that?

• This determines the inputs to the function
• Must identify the variables and data frame

What is the Goal?

What does R need to know?

What is the Goal?

• a scatter plot

What does R need to know?

• which variable goes where
• which data set

What is the Goal?

• a scatter plot (xyplot())

What does R need to know?

• which variable goes where (births ~ dayofyear)
• which data set (data=Births78)
  • use ?Births78 for documentation

xyplot( births ~ dayofyear, data=Births78) 

Two Questions?

1. Command: bwplot()

2. The data: HELPrct

   • Variables: age, substance
   • use ?HELPrct for info about data

bwplot( age ~ substance, data=HELPrct)

1. Command: bwplot()

1. The data: HELPrct
   • Variables: age, substance
   • use ?HELPrct for info about data

bwplot( substance ~ age, data=HELPrct )

histogram( ~ age, data=HELPrct) 

Note: When there is one variable it is on the right side of the formula.

One Variable

  histogram( ~age, data=HELPrct ) 
densityplot( ~age, data=HELPrct ) 
     bwplot( ~age, data=HELPrct ) 
     qqmath( ~age, data=HELPrct ) 
freqpolygon( ~age, data=HELPrct ) 
   bargraph( ~sex, data=HELPrct )

Two Variables

xyplot(  i1 ~ age,       data=HELPrct ) 
bwplot( age ~ substance, data=HELPrct ) 
bwplot( substance ~ age, data=HELPrct ) 

• i1 average number of drinks (standard units) consumed per day in the past 30 days (measured at baseline)

One variable

• histogram(), qqmath(), densityplot(), freqpolygon(), bargraph()

Two Variables

• xyplot(), bwplot()

Create a plot of your own choosing with one of these data sets

names(KidsFeet)    # 4th graders' feet
?KidsFeet

names(Utilities)   # utility bill data
?Utilities

names(NHANES)      # body shape, etc.
?NHANES

• Add groups =group to overlay.
• Use y ~ x | z to create multipanel plots.

densityplot( ~ age | sex, data=HELPrct,  
               groups=substance,  
               auto.key=TRUE)   

• titles
• axis labels
• colors
• sizes
• transparency
• etc, etc.

My approach:

• Let the students ask or
• Let the data analysis drive

xyplot( births ~ dayofyear, data=Births78,  
  groups=dayofyear %% 7, type='l',
  auto.key=list(columns=4, lines=T, points=F),
  par.settings=list(
    superpose.line=list( lty=1 ) ))

Big idea: Replace plot name with summary name

• Nothing else changes

histogram( ~ age, data=HELPrct )
     mean( ~ age, data=HELPrct )

[1] 35.65

The mosaic package includes formula aware versions of mean(), sd(), var(), min(), max(), sum(), IQR(), …

Also provides favstats() to compute our favorites.

favstats( ~ age, data=HELPrct )

 min Q1 median Q3 max  mean   sd   n missing
  19 30     35 40  60 35.65 7.71 453       0

tally( ~ sex, data=HELPrct)

female   male 
   107    346 

tally( ~ substance, data=HELPrct)

alcohol cocaine  heroin 
    177     152     124 

Three ways to think about this. All do the same thing.

sd(   age ~ substance, data=HELPrct )
sd( ~ age | substance, data=HELPrct )
sd( ~ age, groups=substance, data=HELPrct )

alcohol cocaine  heroin 
  7.652   6.693   7.986 

tally( sex ~ substance, data=HELPrct )

        substance
sex      alcohol cocaine heroin
  female  0.2034  0.2697 0.2419
  male    0.7966  0.7303 0.7581

tally( ~ sex + substance, data=HELPrct )

        substance
sex      alcohol cocaine heroin
  female      36      41     30
  male       141     111     94

mean( age ~ substance | sex, data=HELPrct )

  A.F   C.F   H.F   A.M   C.M   H.M     F     M 
39.17 34.85 34.67 37.95 34.36 33.05 36.25 35.47 

mean( age ~ substance | sex, data=HELPrct, 
      .format="table" )

  substance sex  mean
1         A   F 39.17
2         A   M 37.95
3         C   F 34.85
4         C   M 34.36
5         H   F 34.67
6         H   M 33.05

• I've abbreviated the names to make things fit on slide
• Also works for median(), min(), max(), sd(), var(), favstats(), etc.

• single and multiple variable graphical summaries
• single and multiple variabble numerical summaries
• linear models

  mean( age ~ sex, data=HELPrct )
bwplot( age ~ sex, data=HELPrct ) 
    lm( age ~ sex, data=HELPrct )

female   male 
 36.25  35.47 

(Intercept)     sexmale 
    36.2523     -0.7841 

The mosaic package includes some other things, too

• Data sets (you've already seen some of them)
• xtras: xchisq.test(), xpnorm(), xqqmath()
• mplot()
  • mplot(HELPrct) interactive plot creation
  • replacements for plot() in some situations
• simplified histogram() controls (e.g., width)
• simplified ways to add onto lattice plots

xpnorm( 700, mean=500, sd=100)

If X ~ N(500,100), then 

    P(X <= 700) = P(Z <= 2) = 0.9772
    P(X >  700) = P(Z >  2) = 0.0228

[1] 0.9772

xpnorm( c(300, 700), mean=500, sd=100)

If X ~ N(500,100), then 

    P(X <= 300) = P(Z <= -2) = 0.0228
    P(X <= 700) = P(Z <= 2) = 0.9772
    P(X >  300) = P(Z >  -2) = 0.9772
    P(X >  700) = P(Z >  2) = 0.0228

[1] 0.02275 0.97725

xchisq.test(phs)

    Pearson's Chi-squared test with Yates' continuity correction

data:  phs
X-squared = 24.43, df = 1, p-value = 7.71e-07

   104.00   10933.00 
(  146.52) (10890.48)
[12.34]  [ 0.17] 
<-3.51>  < 0.41> 

   189.00   10845.00 
(  146.48) (10887.52)
[12.34]  [ 0.17] 
< 3.51>  <-0.41> 

key:
    observed
    (expected)
    [contribution to X-squared]
    <residual>

Modeling is really the starting point for the mosaic design.

• linear models (lm() and glm()) defined the template
• lattice graphics use the template (so we chose lattice)
• we added functionality so numerical summaries can be done with the same template
• additional things added to make modeling easier for beginners.

model <- lm(width ~ length * sex, 
            data=KidsFeet)
Width <- makeFun(model)
Width( length=25, sex="B")

    1 
9.168 

Width( length=25, sex="G")

    1 
8.939 

xyplot( width ~ length, data=KidsFeet, 
        groups=sex, auto.key=TRUE )
plotFun( Width(length, sex="B") ~ length, 
         col=1, add=TRUE)
plotFun( Width(length, sex="G") ~ length, 
         col=2, add=TRUE)

Lady Tasting Tea

• Often used on first day of class

• Story

  • woman claims she can tell whether milk has been poured into tea or vice versa.
  • Question: How do we test this claim?

Use rflip() to simulate flipping coins

rflip()

Flipping 1 coin [ Prob(Heads) = 0.5 ] ...

H

Number of Heads: 1 [Proportion Heads: 1]

Faster if we flip multiple coins at once:

rflip(10)

Flipping 10 coins [ Prob(Heads) = 0.5 ] ...

H H H T T T H H H T

Number of Heads: 6 [Proportion Heads: 0.6]

• easier to consider heads = correct; tails = incorrect than to compare with a given pattern
  • this switch bothers me more than it bothers my students

rflip(10) simulates 1 lady tasting 10 cups 1 time.

We can do that many times to see how guessing ladies do:

do(2) * rflip(10)

   n heads tails prop
1 10     4     6  0.4
2 10     6     4  0.6

• do() is clever about what it remembers
• 2 isn't many – we'll do many next.

Ladies <- do(5000) * rflip(10)
head(Ladies, 1)

   n heads tails prop
1 10     5     5  0.5

histogram( ~ heads, data=Ladies, width=1 )

tally( ~(heads >= 9) , data=Ladies)

 TRUE FALSE 
   54  4946 

tally( ~(heads >= 9) , data=Ladies, 
                       format="prop")

  TRUE  FALSE 
0.0108 0.9892 

tally( ~(heads >= 9) , data=Ladies, 
                       format="prop")

  TRUE  FALSE 
0.0108 0.9892 

 prop( ~(heads >= 9), data=Ladies)

  TRUE 
0.0108 

1. Do it for your data
2. Do it for “random” data
3. Do it lots of times for “random” data

• definition of “random” is important, but can often be handled by shuffle() or resample()

diffmean(age ~ sex, data=HELPrct)

diffmean 
 -0.7841 

do(1) * 
  diffmean(age ~ shuffle(sex), data=HELPrct)

  diffmean
1    1.357

Null <- do(5000) * 
  diffmean(age ~ shuffle(sex), data=HELPrct)

prop( ~(abs(diffmean) > 0.7841), data=Null )

  TRUE 
0.3566 

histogram(~ diffmean, data=Null, v=-.7841) 

Bootstrap <- do(5000) * 
  diffmean(age~sex, data= resample(HELPrct))

histogram( ~diffmean, data=Bootstrap, 
                      v=-.7841, glwd=4 )

cdata(~diffmean, data=Bootstrap, p=.95)

      low        hi central.p 
  -2.4639    0.8782    0.9500 

confint(Bootstrap, method="quantile")

      name  lower  upper level   method
1 diffmean -2.464 0.8782  0.95 quantile

confint(Bootstrap)  # default uses st. err

      name  lower  upper level method estimate margin.of.error
1 diffmean -2.455 0.8731  0.95 stderr  -0.7908           1.664

do(1) * lm(width ~ length, data=KidsFeet)

  Intercept length  sigma r.squared
1     2.862 0.2479 0.3963     0.411

do(3) * lm( width ~ shuffle(length), data=KidsFeet)

  Intercept   length  sigma r.squared
1     8.259  0.02966 0.5149  0.005883
2    10.388 -0.05644 0.5109  0.021293
3     8.476  0.02087 0.5157  0.002912

do(1) * 
  lm(width ~ length + sex, data=KidsFeet)

  Intercept length    sexG  sigma r.squared
1     3.641  0.221 -0.2325 0.3849    0.4595

do(3) * 
  lm( width ~ length + shuffle(sex), 
                       data=KidsFeet)

  Intercept length      sexG  sigma r.squared
1     2.933 0.2458 -0.034556 0.4014    0.4122
2     2.863 0.2479  0.003523 0.4018    0.4110
3     2.876 0.2526 -0.261422 0.3781    0.4783

Null <- do(5000) * 
  lm( width ~ length + shuffle(sex), 
                       data=KidsFeet)
histogram( ~ sexG, data=Null, 
           v=-0.2325, glwd=4)

histogram(~sexG, data=Null, 
           v=-0.2325, glwd=4)

prop(~ (sexG <= -0.2325), data=Null)

 TRUE 
0.036