R Notebook

Solutions

suppressMessages(require(rms, quietly = TRUE, warn.conflicts = FALSE))
require(splines, quietly = TRUE)
require(plotly, quietly = TRUE, warn.conflicts = FALSE)

Challenge 1

The first time you install a package within a new R session, you must select the CRAN mirror. Usually the best choice is the repository that is geographically closest to you. Once selected subsequent package installations will use this same server.

What command allows you to select a new CRAN mirror?

chooseCRANmirror()

How do you manually specify the default repository?

options(repos = "http://debian.mc.vanderbilt.edu/R/CRAN/")

Challenge 2

gender <- c('M','M','F','M','F','F','M','F','M')
age <- c(34, 64, 38, 63, 40, 73, 27, 51, 47)
smoker <- c('no','yes','no','no','yes','no','no','no','yes')
exercise <- factor(c('moderate','frequent','some','some','moderate','none',
                     'none','moderate','moderate'),
                    levels=c('none','some','moderate','frequent'), ordered=TRUE
)
los <- c(4,8,1,10,6,3,9,4,8)
x <- data.frame(gender, age, smoker, exercise, los)

1. Create a linear model using x, estimating the association between los and all remaining variables

lm(los ~ gender + age + smoker + exercise, data=x)

Call:
lm(formula = los ~ gender + age + smoker + exercise, data = x)

Coefficients:
(Intercept)      genderM          age    smokeryes   exercise.L   exercise.Q   exercise.C  
   0.588144     4.508675     0.033377     2.966623    -2.749852    -0.710942     0.002393  

2. Create a new model, this time predicting los by gender; show the model summary

f <- lm(los ~ gender)
summary(f)

Call:
lm(formula = los ~ gender)

Residuals:
   Min     1Q Median     3Q    Max 
  -3.8   -0.5    0.2    1.2    2.5 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)    3.500      1.099   3.186   0.0154 *
genderM        4.300      1.474   2.917   0.0224 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.197 on 7 degrees of freedom
Multiple R-squared:  0.5487,    Adjusted R-squared:  0.4842 
F-statistic:  8.51 on 1 and 7 DF,  p-value: 0.02243

3. What is the estimate for the intercept? What is the estimate for gender?

coef(f)

(Intercept)     genderM 
        3.5         4.3 

4. Re-calculate the standard errors, by taking the square root of the diagonal of the variance-covariance matrix of the summary of the linear model

sqrt(diag(vcov(summary(f))))

(Intercept)     genderM 
   1.098701    1.474061 

5. Predict los with the following new data set

newdat <- data.frame(gender=c('F','M','F'))
predict(f, newdat)

  1   2   3 
3.5 7.8 3.5 

6. Sum the square of the residuals of the model. Compare this to passing the model to the deviance function.

sum(resid(f)^2)

[1] 33.8

deviance(f)

[1] 33.8

7. Create a subset of x by taking all records where gender is ‘M’ and assigning it to the variable men. Do the same for the variable women.

men <- x[x[,'gender'] == 'M',]
women <- x[x[,'gender'] == 'F',]

8. Call the t.test function, where the first argument is los for women and the second argument is los for men. Add the argument var.equal and set it to TRUE. Does this match the p-value computed in the model summary?

t.test(women[,'los'], men[,'los'], var.equal=TRUE)

    Two Sample t-test

data:  women[, "los"] and men[, "los"]
t = -2.9171, df = 7, p-value = 0.02243
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -7.7856014 -0.8143986
sample estimates:
mean of x mean of y 
      3.5       7.8 

Challenge 3

x <- seq(0, 3.5*pi, length=50)
dat <- data.frame(x = x, y = 4 * sin(x) + rnorm(length(x)))

1. Create a scatterplot of dat.

plot(dat)

2. f <- lm(y ~ ns(x,1), data = dat) will create a linear model where y is predicted with a natural cubic spline of x; in this case it uses a spline with one degree of freedom. Create five models, increasing the degrees of freedom from 1 to 5.

f1 <- lm(y ~ ns(x,1), data = dat)
f2 <- lm(y ~ ns(x,2), data = dat)
f3 <- lm(y ~ ns(x,3), data = dat)
f4 <- lm(y ~ ns(x,4), data = dat)
f5 <- lm(y ~ ns(x,5), data = dat)

3. Generate predicted values for the five models with the given data set.

newdat <- data.frame(x=seq(x[1]*0.9, x[50]*1.1, length=50))

p1 <- predict(f1, newdat)
p2 <- predict(f2, newdat)
p3 <- predict(f3, newdat)
p4 <- predict(f4, newdat)
p5 <- predict(f5, newdat)

4. Recreate the scatterplot of dat, but add lines for the predicted values (try it with base, ggplot, and plotly).

plot(dat, xlim=range(newdat[,'x']))
lines(newdat[,'x'], p1, col=1)

lines(newdat[,'x'], p2, col=2)
lines(newdat[,'x'], p3, col=3)

lines(newdat[,'x'], p4, col=4)
lines(newdat[,'x'], p5, col=5)

predat <- data.frame(x=newdat[,'x'], y=c(p1, p2, p3, p4, p5), model=factor(rep(1:5, each=50)))
p <- ggplot(predat) + aes(x = x, y = y, color = model) + geom_line() + geom_point(data = dat, color='black')
p

ggplotly(p)

Challenge 4

1. Retrieve the titanic3 data set.

getHdata(titanic3)

2. Describe the data.

html(describe(titanic3))

titanic3

14 Variables   1309 Observations

---

pclass

n	missing	distinct
1309	0	3

 Value        1st   2nd   3rd
 Frequency    323   277   709
 Proportion 0.247 0.212 0.542
 

---

survived: Survived

n	missing	distinct	Info	Sum	Mean	Gmd
1309	0	2	0.708	500	0.382	0.4725

---

name: Name

n	missing	distinct
1309	0	1307

lowest :	Abbing, Mr. Anthony	Abbott, Master. Eugene Joseph	Abbott, Mr. Rossmore Edward	Abbott, Mrs. Stanton (Rosa Hunt	Abelseth, Miss. Karen Marie
highest:	Zabour, Miss. Hileni	Zabour, Miss. Thamine	Zakarian, Mr. Mapriededer	Zakarian, Mr. Ortin	Zimmerman, Mr. Leo

---

sex

n	missing	distinct
1309	0	2

 Value      female   male
 Frequency     466    843
 Proportion  0.356  0.644
 

---

age: Age Year

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1046	263	98	0.999	29.88	16.06	5	14	21	28	39	50	57

lowest : 0.1667 0.3333 0.4167 0.6667 0.7500 , highest: 70.5000 71.0000 74.0000 76.0000 80.0000

---

sibsp: Number of Siblings/Spouses Aboard

n	missing	distinct	Info	Mean	Gmd
1309	0	7	0.67	0.4989	0.777

 Value          0     1     2     3     4     5     8
 Frequency    891   319    42    20    22     6     9
 Proportion 0.681 0.244 0.032 0.015 0.017 0.005 0.007
 

---

parch: Number of Parents/Children Aboard

n	missing	distinct	Info	Mean	Gmd
1309	0	8	0.549	0.385	0.6375

 Value          0     1     2     3     4     5     6     9
 Frequency   1002   170   113     8     6     6     2     2
 Proportion 0.765 0.130 0.086 0.006 0.005 0.005 0.002 0.002
 

---

ticket: Ticket Number

n	missing	distinct
1309	0	929

lowest :	110152	110413	110465	110469	110489
highest:	W./C. 6607	W./C. 6608	W./C. 6609	W.E.P. 5734	WE/P 5735

---

fare: Passenger Fare British Pound (¤)

        n  missing distinct     Info     Mean      Gmd      .05      .10      .25 
     1308        1      281        1     33.3    38.61    7.225    7.567    7.896 
      .50      .75      .90      .95 
   14.454   31.275   78.051  133.650 
 

lowest : 0.0000 3.1708 4.0125 5.0000 6.2375 , highest: 227.5250 247.5208 262.3750 263.0000 512.3292

---

cabin

n	missing	distinct
1309	0	187

lowest : A10 A11 A14 A16 , highest: F33 F38 F4 G6 T

---

embarked

n	missing	distinct
1307	2	3

 Value        Cherbourg  Queenstown Southampton
 Frequency          270         123         914
 Proportion       0.207       0.094       0.699
 

---

boat

n	missing	distinct
1309	0	28

lowest : 1 10 11 12 , highest: A B C C D D

---

body: Body Identification Number

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
121	1188	121	1	160.8	113	16	35	72	155	256	297	307

lowest : 1 4 7 9 14 , highest: 312 314 322 327 328

---

home.dest: Home/Destination

n	missing	distinct
745	564	368

lowest :	Aberdeen / Portland, OR	Albany, NY	Altdorf, Switzerland	Amenia, ND	Antwerp, Belgium / Stanton, OH
highest:	Worcester, England	Worcester, MA	Yoevil, England / Cottage Grove	Youngstown, OH	Zurich, Switzerland

---

3. naclus can be used to examine missing data. Use it on the dataset and plot the results.

plot(naclus(titanic3))

4. Describe age, sex, pclass, and embarked with the summaryM function while stratifying on survived. Show the output as a table and plot.

s <- summaryM(age + sex + pclass + embarked ~ survived, data=titanic3)
plot(s)

html(s)

Descriptive Statistics (N=1309).
	N	0 N=809	1 N=500
Age Year	1046	21 28 39	20 28 38
sex : male	1309	0.84 (682)	0.32 (161)
pclass : 1st	1309	0.15 (123)	0.40 (200)
2nd		0.20 (158)	0.24 (119)
3rd		0.65 (528)	0.36 (181)
embarked : Cherbourg	1307	0.15 (120)	0.30 (150)
Queenstown		0.10 ( 79)	0.09 ( 44)
Southampton		0.75 (610)	0.61 (304)
a b c represent the lower quartile a, the median b, and the upper quartile c for continuous variables.   N is the number of non-missing values. Numbers after proportions are frequencies.

5. Use a logistic regression model (lrm) with the formula survived ~ rcs(age,5)*sex+pclass. rcs(age,5) indicates that age will be represented by a restricted cubic spline with five knots. Save the fit object as f and print it.

dd <- datadist(titanic3)
options(datadist = "dd")
f <- lrm(survived ~ rcs(age,5)*sex+pclass, data=titanic3)
f

Frequencies of Missing Values Due to Each Variable
survived      age      sex   pclass 
       0      263        0        0 

Logistic Regression Model
 
 lrm(formula = survived ~ rcs(age, 5) * sex + pclass, data = titanic3)
 
 
                       Model Likelihood     Discrimination    Rank Discrim.    
                          Ratio Test           Indexes           Indexes       
 Obs          1046    LR chi2     465.62    R2       0.485    C       0.854    
  0            619    d.f.            11    g        1.964    Dxy     0.708    
  1            427    Pr(> chi2) <0.0001    gr       7.125    gamma   0.714    
 max |deriv| 8e-08                          gp       0.347    tau-a   0.343    
                                            Brier    0.144                     
 
                   Coef    S.E.   Wald Z Pr(>|Z|)
 Intercept          2.6792 0.5103  5.25  <0.0001 
 age                0.0015 0.0404  0.04  0.9699  
 age'              -0.0341 0.2444 -0.14  0.8891  
 age''              0.2303 1.7373  0.13  0.8945  
 age'''            -0.3131 2.5798 -0.12  0.9034  
 sex=male           0.2338 0.6343  0.37  0.7124  
 pclass=2nd        -1.3745 0.2434 -5.65  <0.0001 
 pclass=3rd        -2.3129 0.2382 -9.71  <0.0001 
 age * sex=male    -0.1892 0.0561 -3.37  0.0007  
 age' * sex=male    0.7244 0.3320  2.18  0.0291  
 age'' * sex=male  -4.1784 2.2962 -1.82  0.0688  
 age''' * sex=male  5.1827 3.3417  1.55  0.1209  
 

6. Obtain predictions for a data.frame that contains variables for age, sex, and pclass. Provide your own values for three or four observations. Pass the predicted values to plogis to perform a logistic transformation.

dat <- expand.grid(age=30, sex=c('female','male'), pclass=c('1st','2nd'))
plogis(predict(f, dat))

        1         2         3         4 
0.9301634 0.5549355 0.7711294 0.2397830 

7. Plot the nomogram of the fit object. Add fun=plogis as an argument to nomogram.

plot(nomogram(f, fun=plogis))

8. Use ggplot and Predict to plot the predictors age, sex, and pclass. Include fun=plogis as an argument to Predict.

ggplot(Predict(f, age, sex, pclass, fun=plogis))

Challenge 5

The CDC has many available data sets. Take a look at the demographics data file for the NHANES National Youth Fitness Survey. This is an XPT file, or SAS export file. The code book can be viewed here.

1. Download, then import this dataset.

dat <- foreign::read.xport("~/Downloads/Y_DEMO.XPT")

2. Keep only the following variables: RIAGENDR, RIDAGEYR, INDHHIND2, DMDHHSIZ. Describe the data.

dat <- dat[,c('RIAGENDR','RIDAGEYR','INDHHIN2','DMDHHSIZ')]
html(describe(dat))

dat

4 Variables   1640 Observations

---

RIAGENDR

n	missing	distinct	Info	Mean	Gmd
1640	0	2	0.75	1.498	0.5003

 Value          1     2
 Frequency    823   817
 Proportion 0.502 0.498
 

---

RIDAGEYR

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1640	0	13	0.994	8.978	4.255	3	4	6	9	12	14	15

 Value          3     4     5     6     7     8     9    10    11    12    13    14
 Frequency    118   123   127   140   134   129   110   132   117   140   129   136
 Proportion 0.072 0.075 0.077 0.085 0.082 0.079 0.067 0.080 0.071 0.085 0.079 0.083
                 
 Value         15
 Frequency    105
 Proportion 0.064
 

---

INDHHIN2

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1634	6	16	0.988	10.98	8.935	2	3	5	8	14	15	15

 Value          1     2     3     4     5     6     7     8     9    10    12    13
 Frequency     44    57   100   101   138   171   117   101   110    82    63    23
 Proportion 0.027 0.035 0.061 0.062 0.084 0.105 0.072 0.062 0.067 0.050 0.039 0.014
                                   
 Value         14    15    77    99
 Frequency    158   325    26    18
 Proportion 0.097 0.199 0.016 0.011
 

---

DMDHHSIZ

n	missing	distinct	Info	Mean	Gmd
1640	0	6	0.937	4.612	1.38

 Value          2     3     4     5     6     7
 Frequency     50   233   572   384   251   150
 Proportion 0.030 0.142 0.349 0.234 0.153 0.091
 

---

3. Turn RIAGENDR into a factor variable with the proper value labels.

dat[,'RIAGENDR'] <- factor(dat[,'RIAGENDR'], labels=c('Male','Female'))

4. Modify INDHHIN2 by collapsing:

• 1,2,3,4,13 into 1
• 5,6,7,8 into 2
• 9,10,14 into 3
• 15 into 4
• 12,77,99 into NA

Then turn it into an ordered factor variable with the labels $0-$19999, $20000-$54999, $55000-$99999, $100000+.

tmpinc <- dat[,'INDHHIN2']
tmpinc[tmpinc %in% c(1,2,3,4,13)] <- 1
tmpinc[tmpinc %in% c(5,6,7,8)] <- 2
tmpinc[tmpinc %in% c(9,10,14)] <- 3
tmpinc[tmpinc == 15] <- 4
tmpinc[!(tmpinc %in% c(1,2,3,4))] <- NA
dat[,'INDHHIN2'] <- factor(tmpinc,
                           labels=c('$0-$19999', '$20000-$54999', '$55000-$99999', '$100000+'),
                           ordered=TRUE)

5. Show the table of gender against age, and household income against size.

with(dat, table(RIAGENDR, RIDAGEYR))

        RIDAGEYR
RIAGENDR  3  4  5  6  7  8  9 10 11 12 13 14 15
  Male   60 65 62 75 64 62 56 66 54 74 62 63 60
  Female 58 58 65 65 70 67 54 66 63 66 67 73 45

with(dat, table(INDHHIN2, DMDHHSIZ))

               DMDHHSIZ
INDHHIN2          2   3   4   5   6   7
  $0-$19999      18  52  82  92  47  34
  $20000-$54999  21  81 184  75 105  61
  $55000-$99999   8  52 147  91  34  18
  $100000+        2  35 132  97  38  21

Challenge 6

Suppose you have a sample of N=100 from the Gamma distribution with shape=2 and scale=1. Calculate the 75th percentile with a basic bootstrapped 95% confidence interval for the upper quartile (75th percentile).

N <- 100
G <- rgamma(N, 2, 1)
# point estimate
PE <- quantile(G, 0.75)
# Bootstrap the 75th percentile
loops <- 1000
x <- numeric(loops)
for(i in seq(loops)) {
  x[i] <- quantile(sample(G, replace=TRUE), 0.75)
}
lb <- quantile(x, 0.025)
ub <- quantile(x, 0.975)
myCI <- unname(2 * PE - c(ub, lb))
myCI

[1] 2.015070 3.145078

Did the CI work?

truth <- quantile(rgamma(1e6, 2, 1), 0.75)
truth > myCI[1] & truth < myCI[2]

 75% 
TRUE 

boot package solution

require(boot, quietly = TRUE, warn.conflicts = FALSE)
boot.ci(boot(G, function(i, j) quantile(i[j], 0.75), R = 1000))

bootstrap variances needed for studentized intervals

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot(G, function(i, j) quantile(i[j], 0.75), 
    R = 1000))

Intervals : 
Level      Normal              Basic         
95%   ( 2.195,  3.317 )   ( 2.015,  3.161 )  

Level     Percentile            BCa          
95%   ( 2.375,  3.521 )   ( 2.375,  3.521 )  
Calculations and Intervals on Original Scale

Challenge 7

https://fivethirtyeight.com/features/can-you-solve-these-colorful-puzzles/

You play a game with four balls: One ball is red, one is blue, one is green and one is yellow. They are placed in a box. You draw a ball out of the box at random and note its color. Without replacing the first ball, you draw a second ball and then paint it to match the color of the first. Replace both balls, and repeat the process. The game ends when all four balls have become the same color. What is the expected number of turns to finish the game?

Simulate the answer.

draw <- function(x) sample(x)[c(1,1,3,4)]
test <- function(ii=c('r','b','g','y')) {
  i <- 1
  while(length(unique(ii<-draw(ii))) > 1) i <- i+1
  i
}
N <- 1000
ans <- numeric(N)
for(i in seq(N)) {
  ans[i] <- test()
}
mean(ans)

[1] 9.001