Harold Nelson

11/25/2018

First we’ll get some data on birthweights from the MASS package.

This dataframe was created to study the effect of various risk factors on the birthweight of a baby.

Here is the basic data dictionary and a source citation copied from the documentation in the MASS package.

low: indicator of birth weight less than 2.5 kg.

age: mother’s age in years.

lwt: mother’s weight in pounds at last menstrual period.

race: mother’s race (1 = white, 2 = black, 3 = other).

smoke: smoking status during pregnancy.

ptl: number of previous premature labours.

ht: history of hypertension.

ui: presence of uterine irritability.

ftv: number of physician visits during the first trimester.

bwt: birth weight in grams.

Source

Hosmer, D.W. and Lemeshow, S. (1989) Applied Logistic Regression. New York: Wiley

Once MASS has been downloaded to the computer with the install.packages() command, you should comment out the first line in the chunk below, as I have done.

```
#install.packages("MASS")
library(MASS)
str(birthwt)
```

```
## 'data.frame': 189 obs. of 10 variables:
## $ low : int 0 0 0 0 0 0 0 0 0 0 ...
## $ age : int 19 33 20 21 18 21 22 17 29 26 ...
## $ lwt : int 182 155 105 108 107 124 118 103 123 113 ...
## $ race : int 2 3 1 1 1 3 1 3 1 1 ...
## $ smoke: int 0 0 1 1 1 0 0 0 1 1 ...
## $ ptl : int 0 0 0 0 0 0 0 0 0 0 ...
## $ ht : int 0 0 0 0 0 0 0 0 0 0 ...
## $ ui : int 1 0 0 1 1 0 0 0 0 0 ...
## $ ftv : int 0 3 1 2 0 0 1 1 1 0 ...
## $ bwt : int 2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
```

To see if race might be related to birthweight we can use the tapply() function. There are three arguments we need to specify for tapply() to give us useful information.

The first argument is the variable of interest, which is uaually quantitative. In this case it is birthwt$bwt.

The second argument is either a categorical variable or a numerical variable with a small number of values. In this case, it is birthwt$race.

The third argument is the name of a function. This function will be applied to subsets of the first argument defined by the values of the second argument. Let’s do this and examine the results.

`tapply(birthwt$bwt,birthwt$race,summary)`

```
## $`1`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1021 2585 3062 3103 3651 4990
##
## $`2`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1135 2370 2849 2720 3057 3860
##
## $`3`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 709 2313 2835 2805 3274 4054
```

We can see the usual summary statistics of birthweigt for the three subsets of our dataframe. We can make the obvious comparisons based on the values of the mean or the median. Race 1 has the highest mean birthweight, etc. Of course, what we’re doing falls far short of formal statistical inference. We’re just “kicking the tires” of our data.

We also might want to get other summary statistics, such as the standard deviation, range, or interquartile range not provided by the function summary(). It’s easy to do this by just re-running the code above with the corresponding function names. Let’s do that.

`print("sd")`

`## [1] "sd"`

`tapply(birthwt$bwt,birthwt$race,sd)`

```
## 1 2 3
## 727.8861 638.6839 722.1944
```

`print("IQR")`

`## [1] "IQR"`

`tapply(birthwt$bwt,birthwt$race,IQR)`

```
## 1 2 3
## 1066.25 686.50 961.00
```

Can we get graphical results? Let’s try to get histograms.

`tapply(birthwt$bwt,birthwt$race,hist)`

```
## $`1`
## $breaks
## [1] 1000 1500 2000 2500 3000 3500 4000 4500 5000
##
## $counts
## [1] 1 6 16 20 20 25 6 2
##
## $density
## [1] 2.083333e-05 1.250000e-04 3.333333e-04 4.166667e-04 4.166667e-04
## [6] 5.208333e-04 1.250000e-04 4.166667e-05
##
## $mids
## [1] 1250 1750 2250 2750 3250 3750 4250 4750
##
## $xname
## [1] "X[[i]]"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## $`2`
## $breaks
## [1] 1000 1500 2000 2500 3000 3500 4000
##
## $counts
## [1] 1 2 8 7 6 2
##
## $density
## [1] 7.692308e-05 1.538462e-04 6.153846e-04 5.384615e-04 4.615385e-04
## [6] 1.538462e-04
##
## $mids
## [1] 1250 1750 2250 2750 3250 3750
##
## $xname
## [1] "X[[i]]"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## $`3`
## $breaks
## [1] 500 1000 1500 2000 2500 3000 3500 4000 4500
##
## $counts
## [1] 1 2 6 16 11 19 11 1
##
## $density
## [1] 2.985075e-05 5.970149e-05 1.791045e-04 4.776119e-04 3.283582e-04
## [6] 5.671642e-04 3.283582e-04 2.985075e-05
##
## $mids
## [1] 750 1250 1750 2250 2750 3250 3750 4250
##
## $xname
## [1] "X[[i]]"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
```

The labels are very poor, but in this situation, with only a few values of the second argument, we know that the first histogram is for race 1, etc.

Does the first argument have to be quantitative? Let’s try to get a table of the smoke variable for each value of the race variable.

`tapply(birthwt$smoke,birthwt$race,table)`

```
## $`1`
##
## 0 1
## 44 52
##
## $`2`
##
## 0 1
## 16 10
##
## $`3`
##
## 0 1
## 55 12
```

Yes, this does work, but there are easier ways to do this, which produce more readable output. The simple base R table() command with two arguments is one way.

`table(birthwt$race,birthwt$smoke)`

```
##
## 0 1
## 1 44 52
## 2 16 10
## 3 55 12
```

I started out to just clarify the use of tapply(), but I saw some related things, which I felt compelled to discuss, primarily to point out that tapply() is not the best tool for all situations.

What I want you to do is apply the usage rules of tapply() to a few examples.

Use tapply() to produce summary statistics, the standard deviation, and the IQR, of the variable bwt for the subsets of birthwt defined by the variable smoke. Note that this requires three calls to tapply(). Do this yourself before you look at the answer on the next slide.

`print("Basic Summary")`

`## [1] "Basic Summary"`

`tapply(birthwt$bwt,birthwt$smoke,summary)`

```
## $`0`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1021 2509 3100 3056 3622 4990
##
## $`1`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 709 2370 2776 2772 3246 4238
```

`print("Standard Deviation")`

`## [1] "Standard Deviation"`

`tapply(birthwt$bwt,birthwt$smoke,sd)`

```
## 0 1
## 752.6566 659.6349
```

`print("IQR")`

`## [1] "IQR"`

`tapply(birthwt$bwt,birthwt$smoke,IQR)`

```
## 0 1
## 1112.50 875.25
```

Use tapply() to produce summary statistics, the standard deviation, and the IQR, of the variable ftv (first trimester visits) for the subsets of birthwt defined by the variable low (low birthweight). Note that this requires three calls to tapply(). Do this yourself before you look at the answer on the next slide.

`print("Basic Summary")`

`## [1] "Basic Summary"`

`tapply(birthwt$ftv,birthwt$low,summary)`

```
## $`0`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.8385 1.0000 6.0000
##
## $`1`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.6949 1.0000 4.0000
```

`print("Standard Deviation")`

`## [1] "Standard Deviation"`

`tapply(birthwt$ftv,birthwt$low,sd)`

```
## 0 1
## 1.069694 1.038139
```

`print("IQR")`

`## [1] "IQR"`

`tapply(birthwt$ftv,birthwt$low,IQR)`

```
## 0 1
## 1 1
```

Use tapply() to produce summary statistics, the standard deviation, and the IQR, of the variable age for the subsets of birthwt defined by the variable low. Note that this requires three calls to tapply(). Do this yourself before you look at the answer on the next slide.

`print("Basic Summary")`

`## [1] "Basic Summary"`

`tapply(birthwt$age,birthwt$low,summary)`

```
## $`0`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.00 19.00 23.00 23.66 28.00 45.00
##
## $`1`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.00 19.50 22.00 22.31 25.00 34.00
```

`print("Standard Deviation")`

`## [1] "Standard Deviation"`

`tapply(birthwt$age,birthwt$low,sd)`

```
## 0 1
## 5.584522 4.511496
```

`print("IQR")`

`## [1] "IQR"`

`tapply(birthwt$age,birthwt$low,IQR)`

```
## 0 1
## 9.0 5.5
```