BIO 201 Homework 4

## Settings for RMarkdown http://yihui.name/knitr/options#chunk_options
opts_chunk$set(comment = "", warning = FALSE, message = FALSE, tidy = FALSE, 
    echo = TRUE, fig.width = 7, fig.height = 7)
options(width = 116, scipen = 10)

setwd("~/statistics/bio201/")

library(ggplot2)

References

1. 6.23-6.24

## 6.23
dat.6.23 <- read.table(header = TRUE, text = "
variable        base.mean       base.sd         fu.mean         fu.sd
Cr              0.97            0.22            1.00            0.19
K               4.43            0.64            4.49            0.71
Pi              1.68            0.47            1.57            0.40
PAIS            36.50           16.08           23.27           13.79
")

## 95% CI by t-distribution (If n > 200, normal approximation is possible. Rosner page 172)
dat.6.23ans <- within(dat.6.23, {
    base.95.ci.upper <- base.mean + qt(p = 0.975, df = 102 - 1) * base.sd / sqrt(102)
    base.95.ci.lower <- base.mean + qt(p = 0.025, df = 102 - 1) * base.sd / sqrt(102)

    fu.95.ci.upper   <- fu.mean   + qt(p = 0.975, df =  69 - 1) * fu.sd   / sqrt(69)
    fu.95.ci.lower   <- fu.mean   + qt(p = 0.025, df =  69 - 1) * fu.sd   / sqrt(69)
})

## Baseline
dat.6.23ans[, grep("base", names(dat.6.23ans))]
  base.mean base.sd base.95.ci.lower base.95.ci.upper
1      0.97    0.22           0.9268            1.013
2      4.43    0.64           4.3043            4.556
3      1.68    0.47           1.5877            1.772
4     36.50   16.08          33.3416           39.658
## Follow-up
dat.6.23ans[, grep("fu", names(dat.6.23ans))]
  fu.mean fu.sd fu.95.ci.lower fu.95.ci.upper
1    1.00  0.19         0.9544          1.046
2    4.49  0.71         4.3194          4.661
3    1.57  0.40         1.4739          1.666
4   23.27 13.79        19.9573         26.583

## 6.24
## Not done

2. 6.77-6.80

## 6.77
## Not done

## 6.78
binom.test(x = 156, n = 1996)$conf.int
[1] 0.06676 0.09081
attr(,"conf.level")
[1] 0.95

## 6.79
FPR <- 156 / 1996
1 - (1 - FPR)^5
[1] 0.3343

## 6.80
## Not sure. Here doing same calculation with 95% confidence limit values.
1 - (1 - binom.test(156, 1966)$conf.int)^5
[1] 0.2960 0.3834
attr(,"conf.level")
[1] 0.95

3. 6.86-6.87

library(foreign)
bone.den <- read.dta("./BONEDEN.DAT.dta")

bone.den <- within(bone.den, {
    diff.ls <- ls2 - ls1
    diff.fn <- fn2 - fn1
    diff.fs <- fs2 - fs1
})

## page 176 example
t.test(bone.den$diff.ls)$conf
[1] -0.063721 -0.007986
attr(,"conf.level")
[1] 0.95
## 6.86 Femoral neck
t.test(bone.den$diff.fn)$conf
[1] -0.03013  0.02866
attr(,"conf.level")
[1] 0.95
## 6.87 Femoral shaft
t.test(bone.den$diff.fs)$conf
[1] -0.066704  0.005728
attr(,"conf.level")
[1] 0.95

4. 6.97-6.101

## 6.97
121.85 + qt(p = c(0.025, 0.975), df = 9) * 5.75 / sqrt(10)
[1] 117.7 126.0
sbp <- c(121,109,117.5,125,125,129,123,118.5,123.5,127)
t.test(sbp)$conf
[1] 117.7 126.0
attr(,"conf.level")
[1] 0.95

## 6.98: Not hypertensive

## 6.99
0.4
[1] 0.4

## 6.100
binom.test(4, 10)$conf
[1] 0.1216 0.7376
attr(,"conf.level")
[1] 0.95

## 6.101: Inconclusive.

5. 6.102-6.103

## 6.102
res.pois <- poisson.test(x = 45, T = 10112, r = 2.73/1000)
res.pois

    Exact Poisson test

data:  45 time base: 10112 
number of events = 45, time base = 10112, p-value = 0.002155
alternative hypothesis: true event rate is not equal to 0.00273 
95 percent confidence interval:
 0.003246 0.005955 
sample estimates:
event rate 
   0.00445 


## 6.103
attributes(res.pois$conf) <- NULL
res.pois$conf * 1000
[1] 3.246 5.955