Homework 4 PUBlIC HEALTH 251D
4.12.1
Stat <- rep(c("Dead","Survived"), 4)
trt <- rep(c("Tolbutamide", "Placebo"), 4)
ag <- c(rep("<55",4),rep("55+",4))
freq <- c(8, 98, 5, 115, 22, 76, 16, 69)
(df <- data.frame(Status = Stat,
Treatment = trt,
Agegrp = ag,
Freq = freq))
## Status Treatment Agegrp Freq
## 1 Dead Tolbutamide <55 8
## 2 Survived Placebo <55 98
## 3 Dead Tolbutamide <55 5
## 4 Survived Placebo <55 115
## 5 Dead Tolbutamide 55+ 22
## 6 Survived Placebo 55+ 76
## 7 Dead Tolbutamide 55+ 16
## 8 Survived Placebo 55+ 69
4.12.2
fn <- c("Mark", "Michael", "Nick")
ln <- c("Mccoy", "Koehler", "Glennister")
af <- c("friend", "brother", "friend")
em <- c("mmccoy@gmail.com", "na", "nglennister@sdsu.edu")
dte <- c(rep("10/17/18",3))
(df2 <- data.frame(Firstname = fn,
Lastname = ln,
Affiliation = af,
Email = em,
Date = dte))
## Firstname Lastname Affiliation Email Date
## 1 Mark Mccoy friend mmccoy@gmail.com 10/17/18
## 2 Michael Koehler brother na 10/17/18
## 3 Nick Glennister friend nglennister@sdsu.edu 10/17/18
4.12.3
aids <- read.table(url("http://www.medepi.net/data/aids.txt"), header=TRUE)
aids$cases <- as.numeric(aids$cases)
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following object is masked _by_ '.GlobalEnv':
##
## Stat
aids$cases[is.na(aids$cases)] <- 0
ggplot(aids, aes(x=year, y = cases)) + geom_line()
4.12.4
measles <- read.table(url("http://www.medepi.net/data/measles.txt"), header=T)
plot(measles$year, measles$cases, type="l", lwd=2, col="skyblue1", main="Normal scale")
plot(measles$year, log(measles$cases), type="l", lwd=2, col="skyblue1", main="log scale")
4.12.5
hepb <- read.table(url("http://www.medepi.net/data/hepb.txt"), header=T, na.strings = ".")
matplot(hepb$year, cbind(hepb$cases,aids$cases),
type = "l", lwd = 2, xlab = "Year", ylab = "Cases",
main = "Reported cases of Hepatitis B and AIDS, United States, 1980--2003")
legend(1980, 100000, legend = c("Hepatitis B", "AIDS"),
lwd = 2, lty = 1:2, col = 1:2)
4.12.6
edat <- read.table("http://www.medepi.net/data/evans.txt", header = TRUE, sep = "")
str(edat)
## 'data.frame': 609 obs. of 12 variables:
## $ id : int 21 31 51 71 74 91 111 131 141 191 ...
## $ chd: int 0 0 1 0 0 0 1 0 0 0 ...
## $ cat: int 0 0 1 1 0 0 0 0 0 0 ...
## $ age: int 56 43 56 64 49 46 52 63 42 55 ...
## $ chl: int 270 159 201 179 243 252 179 217 176 250 ...
## $ smk: int 0 1 1 1 1 1 1 0 1 0 ...
## $ ecg: int 0 0 1 0 0 0 1 0 0 1 ...
## $ dbp: int 80 74 112 100 82 88 80 92 76 114 ...
## $ sbp: int 138 128 164 200 145 142 128 135 114 182 ...
## $ hpt: int 0 0 1 1 0 0 0 0 0 1 ...
## $ ch : int 0 0 1 1 0 0 0 0 0 0 ...
## $ cc : int 0 0 201 179 0 0 0 0 0 0 ...
xtabs(~chd, data = edat)
## chd
## 0 1
## 538 71
edat$chd2 <- factor(edat$chd, levels = 0:1, labels = c("No", "Yes"))
xtabs(~chd2, data = edat)
## chd2
## No Yes
## 538 71
xtabs(~cat, data = edat)
## cat
## 0 1
## 487 122
edat$cat2 <- factor(edat$cat, levels = 0:1, labels = c("Normal", "High"))
xtabs(~cat2, data = edat)
## cat2
## Normal High
## 487 122
xtabs(~smk, data = edat)
## smk
## 0 1
## 222 387
edat$smk2 <- factor(edat$smk, levels = 0:1, labels = c("Never", "Ever"))
xtabs(~smk2, data = edat)
## smk2
## Never Ever
## 222 387
xtabs(~ecg, data = edat)
## ecg
## 0 1
## 443 166
edat$ecg2 <- factor(edat$ecg, levels = 0:1, labels = c("Normal", "Abnormal"))
xtabs(~ecg2, data = edat)
## ecg2
## Normal Abnormal
## 443 166
xtabs(~hpt, data = edat)
## hpt
## 0 1
## 354 255
edat$hpt2 <- factor(edat$hpt, levels = 0:1, labels = c("No", "Yes"))
xtabs(~hpt2, data = edat)
## hpt2
## No Yes
## 354 255
quantile(edat$age)
## 0% 25% 50% 75% 100%
## 40 46 52 60 76
edat$age4 <- cut(edat$age, quantile(edat$age), right = FALSE, include.lowest = TRUE)
xtabs(~age4, data = edat)
## age4
## [40,46) [46,52) [52,60) [60,76]
## 134 158 158 159
hptnew <- rep(NA, nrow(edat))
normal <- edat$sbp < 120 & edat$dbp < 80
hptnew[normal] <- 1
prehyp <- (edat$sbp >= 120 & edat$sbp < 140) | (edat$dbp >= 80 & edat$dbp < 90)
hptnew[prehyp] <- 2
stage1 <- (edat$sbp >= 140 & edat$sbp < 160) | (edat$dbp >= 90 & edat$dbp < 100)
hptnew[stage1] <- 3
stage2 <- edat$sbp >= 160 | edat$dbp >= 100
hptnew[stage2] <- 4
edat$hpt4 <- factor(hptnew, levels = 1:4,
labels=c("Normal", "PreHTN", "HTN.Stage1", "HTN.Stage2"))
xtabs(~hpt4, data = edat) # test
## hpt4
## Normal PreHTN HTN.Stage1 HTN.Stage2
## 56 165 177 211
#### (4)
xtabs(~hpt2 + hpt4, data = edat)
## hpt4
## hpt2 Normal PreHTN HTN.Stage1 HTN.Stage2
## No 56 165 133 0
## Yes 0 0 44 211
4.12.7
wdat <- read.table("http://www.medepi.net/data/wnv/wnv2004raw.txt",
header = TRUE, sep = ",", as.is = TRUE, na.strings = c(".","Unknown"))
str(wdat)
## 'data.frame': 779 obs. of 8 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ county : chr "San Bernardino" "San Bernardino" "San Bernardino" "San Bernardino" ...
## $ age : int 40 64 19 12 12 17 61 74 71 26 ...
## $ sex : chr "F" "F" "M" "M" ...
## $ syndrome : chr "WNF" "WNF" "WNF" "WNF" ...
## $ date.onset : chr "05/19/2004" "05/22/2004" "05/22/2004" "05/16/2004" ...
## $ date.tested: chr "06/02/2004" "06/16/2004" "06/16/2004" "06/16/2004" ...
## $ death : chr "No" "No" "No" "No" ...
wdat$date.onset2 <- as.Date(wdat$date.onset, format = "%m/%d/%Y")
wdat$date.tested2 <- as.Date(wdat$date.tested, format = "%m/%d/%Y")
#write.table(wdat, "~/Desktop/wnvdat.txt",
#sep = ",", row.names = FALSE)
4.12.8
odat <- read.table("http://www.medepi.net/data/oswego/oswego.txt", sep = "",
header = TRUE, na.strings = ".")
str(odat)
## 'data.frame': 75 obs. of 21 variables:
## $ id : int 2 3 4 6 7 8 9 10 14 16 ...
## $ age : int 52 65 59 63 70 40 15 33 10 32 ...
## $ sex : Factor w/ 2 levels "F","M": 1 2 1 1 2 1 1 1 2 1 ...
## $ meal.time : Factor w/ 7 levels "10:00 PM","11:00 AM",..: 6 3 3 5 5 5 1 4 5 7 ...
## $ ill : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 2 2 2 ...
## $ onset.date : Factor w/ 3 levels "4/18","4/19",..: 2 2 2 1 1 2 2 1 2 2 ...
## $ onset.time : Factor w/ 18 levels "1:00 AM","10:00 PM",..: 9 9 9 5 5 10 1 6 10 4 ...
## $ baked.ham : Factor w/ 2 levels "N","Y": 2 2 2 2 2 1 1 2 1 2 ...
## $ spinach : Factor w/ 2 levels "N","Y": 2 2 2 2 2 1 1 2 1 2 ...
## $ mashed.potato : Factor w/ 3 levels "N","NA","Y": 3 3 1 1 3 1 1 3 1 1 ...
## $ cabbage.salad : Factor w/ 2 levels "N","Y": 1 2 1 2 1 1 1 1 1 1 ...
## $ jello : Factor w/ 2 levels "N","Y": 1 1 1 2 2 1 1 1 1 1 ...
## $ rolls : Factor w/ 2 levels "N","Y": 2 1 1 1 2 1 1 2 1 2 ...
## $ brown.bread : Factor w/ 2 levels "N","Y": 1 1 1 1 2 1 1 2 1 1 ...
## $ milk : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
## $ coffee : Factor w/ 2 levels "N","Y": 2 2 2 1 2 1 1 1 1 2 ...
## $ water : Factor w/ 2 levels "N","Y": 1 1 1 2 2 1 1 2 1 1 ...
## $ cakes : Factor w/ 2 levels "N","Y": 1 1 2 1 1 1 2 1 1 2 ...
## $ vanilla.ice.cream : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 1 2 2 2 ...
## $ chocolate.ice.cream: Factor w/ 3 levels "N","NA","Y": 1 3 3 1 1 3 3 3 3 3 ...
## $ fruit.salad : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
head(odat)
## id age sex meal.time ill onset.date onset.time baked.ham spinach
## 1 2 52 F 8:00 PM Y 4/19 12:30 AM Y Y
## 2 3 65 M 6:30 PM Y 4/19 12:30 AM Y Y
## 3 4 59 F 6:30 PM Y 4/19 12:30 AM Y Y
## 4 6 63 F 7:30 PM Y 4/18 10:30 PM Y Y
## 5 7 70 M 7:30 PM Y 4/18 10:30 PM Y Y
## 6 8 40 F 7:30 PM Y 4/19 2:00 AM N N
## mashed.potato cabbage.salad jello rolls brown.bread milk coffee water
## 1 Y N N Y N N Y N
## 2 Y Y N N N N Y N
## 3 N N N N N N Y N
## 4 N Y Y N N N N Y
## 5 Y N Y Y Y N Y Y
## 6 N N N N N N N N
## cakes vanilla.ice.cream chocolate.ice.cream fruit.salad
## 1 N Y N N
## 2 N Y Y N
## 3 Y Y Y N
## 4 N Y N N
## 5 N Y N N
## 6 N Y Y N
mdt <- paste("4/18/1940", odat$meal.time)
meal.dt <- strptime(mdt, "%m/%d/%Y %I:%M %p")
odt <- paste(paste(odat$onset.date,"/1940",sep = ""), odat$onset.time)
onset.dt <- strptime(odt, "%m/%d/%Y %I:%M %p")
hist(onset.dt, breaks = 30, freq = TRUE)
min.obs.pos <- which(onset.dt == min(onset.dt, na.rm = TRUE))
min.obs.pos
## [1] 33
max.obs.pos <- which(onset.dt == max(onset.dt, na.rm = TRUE))
max.obs.pos
## [1] 10
odat[min.obs.pos, ]
## id age sex meal.time ill onset.date onset.time baked.ham spinach
## 33 52 8 M 11:00 AM Y 4/18 3:00 PM N N
## mashed.potato cabbage.salad jello rolls brown.bread milk coffee water
## 33 N N N N N N N N
## cakes vanilla.ice.cream chocolate.ice.cream fruit.salad
## 33 N Y Y N
odat[max.obs.pos, ]
## id age sex meal.time ill onset.date onset.time baked.ham spinach
## 10 16 32 F NA Y 4/19 10:30 AM Y Y
## mashed.potato cabbage.salad jello rolls brown.bread milk coffee water
## 10 N N N Y N N Y N
## cakes vanilla.ice.cream chocolate.ice.cream fruit.salad
## 10 Y Y Y N
onset.ct <- as.POSIXct(onset.dt)
odat2 <- odat[order(odat$ill, onset.ct), ]; odat2
## id age sex meal.time ill onset.date onset.time baked.ham spinach
## 47 1 11 M NA N NA NA N N
## 48 5 13 F NA N NA NA N N
## 49 11 65 M NA N NA NA Y Y
## 50 12 38 F NA N NA NA Y Y
## 51 13 62 F NA N NA NA Y Y
## 52 15 25 M NA N NA NA Y Y
## 53 19 11 M NA N NA NA Y Y
## 54 23 64 M NA N NA NA N N
## 55 25 65 F NA N NA NA Y Y
## 56 28 62 M NA N NA NA Y Y
## 57 30 17 M 10:00 PM N NA NA N N
## 58 34 40 M NA N NA NA Y Y
## 59 35 35 F NA N NA NA Y Y
## 60 37 36 M NA N NA NA Y N
## 61 41 54 F NA N NA NA Y Y
## 62 45 20 M 10:00 PM N NA NA N N
## 63 46 17 M NA N NA NA Y Y
## 64 50 9 F NA N NA NA N N
## 65 51 50 M NA N NA NA Y Y
## 66 53 35 F NA N NA NA N N
## 67 56 11 F NA N NA NA N N
## 68 61 37 M NA N NA NA N N
## 69 62 24 F NA N NA NA Y Y
## 70 63 69 F NA N NA NA N Y
## 71 64 7 M NA N NA NA Y Y
## 72 67 11 F 7:30 PM N NA NA Y Y
## 73 68 17 M 7:30 PM N NA NA Y Y
## 74 69 36 F NA N NA NA N N
## 75 73 14 F 10:00 PM N NA NA N N
## 33 52 8 M 11:00 AM Y 4/18 3:00 PM N N
## 20 31 35 M NA Y 4/18 9:00 PM Y Y
## 23 36 35 F NA Y 4/18 9:15 PM Y Y
## 26 40 68 M NA Y 4/18 9:30 PM Y N
## 29 44 58 M NA Y 4/18 9:30 PM Y Y
## 16 24 3 M NA Y 4/18 9:45 PM N Y
## 17 26 59 F NA Y 4/18 9:45 PM N Y
## 13 20 33 F NA Y 4/18 10:00 PM Y Y
## 12 18 36 M NA Y 4/18 10:15 PM Y Y
## 4 6 63 F 7:30 PM Y 4/18 10:30 PM Y Y
## 5 7 70 M 7:30 PM Y 4/18 10:30 PM Y Y
## 32 49 52 F NA Y 4/18 10:30 PM Y Y
## 36 57 74 M NA Y 4/18 10:30 PM Y Y
## 8 10 33 F 7:00 PM Y 4/18 11:00 PM Y Y
## 15 22 7 M NA Y 4/18 11:00 PM Y Y
## 19 29 37 F NA Y 4/18 11:00 PM Y Y
## 35 55 25 M NA Y 4/18 11:00 PM Y N
## 46 75 45 F NA Y 4/18 11:00 PM Y Y
## 24 38 57 F NA Y 4/18 11:30 PM Y Y
## 39 60 53 F 7:30 PM Y 4/18 11:30 PM Y Y
## 34 54 48 F NA Y 4/19 12:00 AM Y Y
## 44 72 18 F 7:30 PM Y 4/19 12:00 AM Y Y
## 1 2 52 F 8:00 PM Y 4/19 12:30 AM Y Y
## 2 3 65 M 6:30 PM Y 4/19 12:30 AM Y Y
## 3 4 59 F 6:30 PM Y 4/19 12:30 AM Y Y
## 11 17 62 F NA Y 4/19 12:30 AM N N
## 30 47 62 F NA Y 4/19 12:30 AM Y Y
## 41 66 8 F NA Y 4/19 12:30 AM Y N
## 42 70 21 F NA Y 4/19 12:30 AM Y N
## 7 9 15 F 10:00 PM Y 4/19 1:00 AM N N
## 14 21 13 F 10:00 PM Y 4/19 1:00 AM N N
## 18 27 15 F 10:00 PM Y 4/19 1:00 AM N N
## 21 32 15 M 10:00 PM Y 4/19 1:00 AM N N
## 22 33 50 F 10:00 PM Y 4/19 1:00 AM N N
## 25 39 16 F 10:00 PM Y 4/19 1:00 AM N N
## 31 48 20 F 7:00 PM Y 4/19 1:00 AM N N
## 37 58 12 F 10:00 PM Y 4/19 1:00 AM N N
## 40 65 17 F 10:00 PM Y 4/19 1:00 AM N N
## 43 71 60 M 7:30 PM Y 4/19 1:00 AM N N
## 6 8 40 F 7:30 PM Y 4/19 2:00 AM N N
## 9 14 10 M 7:30 PM Y 4/19 2:00 AM N N
## 28 43 72 F NA Y 4/19 2:00 AM Y Y
## 45 74 52 M NA Y 4/19 2:15 AM Y N
## 27 42 77 M NA Y 4/19 2:30 AM N N
## 38 59 44 F 7:30 PM Y 4/19 2:30 AM Y Y
## 10 16 32 F NA Y 4/19 10:30 AM Y Y
## mashed.potato cabbage.salad jello rolls brown.bread milk coffee water
## 47 N N N N N N N N
## 48 N N N N N N N N
## 49 Y N Y Y N N N N
## 50 Y N N Y N N Y N
## 51 N Y Y Y Y N N Y
## 52 Y Y Y Y Y Y Y Y
## 53 NA Y N Y N N N Y
## 54 N N N N N N N N
## 55 Y Y Y N Y N Y N
## 56 N Y N Y Y N Y Y
## 57 N N N N N N N N
## 58 N N N Y Y N Y Y
## 59 Y N N Y Y N Y Y
## 60 Y Y N Y Y N Y N
## 61 Y N N Y N N Y N
## 62 N N N N N N N N
## 63 Y N N Y N N N Y
## 64 N N N N N N N N
## 65 Y Y Y Y Y Y Y Y
## 66 N N N N N N N N
## 67 N N N N N N N N
## 68 N N N N N N N N
## 69 Y N N Y N N Y N
## 70 Y N Y N Y N N Y
## 71 Y Y Y Y N N N Y
## 72 Y Y N Y N N Y Y
## 73 Y Y N Y N N Y N
## 74 N N N N N N N N
## 75 N N N N N N N N
## 33 N N N N N N N N
## 20 Y N Y Y Y N Y N
## 23 Y Y N Y Y N Y N
## 26 Y Y N N Y N Y N
## 29 Y N N N Y Y Y N
## 16 Y N N Y N N N Y
## 17 Y Y N Y Y N N Y
## 13 Y Y Y Y N N Y Y
## 12 N Y N Y Y N N N
## 4 N Y Y N N N N Y
## 5 Y N Y Y Y N Y Y
## 32 Y Y N Y N N Y N
## 36 Y Y Y Y Y N Y N
## 8 Y N N Y Y N N Y
## 15 Y Y Y Y Y N N Y
## 19 Y N Y Y Y N Y N
## 35 Y N N Y Y N N Y
## 46 Y Y Y Y Y N Y N
## 24 N Y Y Y Y N Y N
## 39 Y Y Y N Y N Y Y
## 34 Y Y Y Y Y Y Y N
## 44 Y Y Y N N N N Y
## 1 Y N N Y N N Y N
## 2 Y Y N N N N Y N
## 3 N N N N N N Y N
## 11 N N N N N N N N
## 30 N N N Y N N N Y
## 41 Y Y Y N N N N N
## 42 N Y Y N N N N N
## 7 N N N N N N N N
## 14 N N N N N N N N
## 18 N N N N N N N N
## 21 N N N N N N N N
## 22 N N N N N N N N
## 25 N N N N N N N N
## 31 N N N N N N N N
## 37 N N N N N N N N
## 40 N N N N N N N N
## 43 N N N N N N N N
## 6 N N N N N N N N
## 9 N N N N N N N N
## 28 N Y Y N Y N Y N
## 45 Y N Y Y Y N Y Y
## 27 N N N N N N N N
## 38 Y N N Y N N N Y
## 10 N N N Y N N Y N
## cakes vanilla.ice.cream chocolate.ice.cream fruit.salad
## 47 N N Y N
## 48 N N Y N
## 49 N Y N N
## 50 N Y Y Y
## 51 N N N N
## 52 Y Y N N
## 53 N N Y N
## 54 N Y N N
## 55 Y Y Y N
## 56 Y N Y N
## 57 Y Y Y N
## 58 Y N Y Y
## 59 N N Y N
## 60 N N Y N
## 61 Y N Y N
## 62 Y Y Y N
## 63 N Y Y N
## 64 Y N Y N
## 65 Y N Y N
## 66 N Y Y N
## 67 N N Y N
## 68 N N Y N
## 69 N N N N
## 70 Y N Y N
## 71 Y N Y N
## 72 N N Y N
## 73 Y Y N N
## 74 N N Y N
## 75 Y Y N N
## 33 N Y Y N
## 20 Y Y N Y
## 23 N Y N N
## 26 N Y N N
## 29 N Y NA Y
## 16 Y Y N N
## 17 Y Y N N
## 13 Y Y Y N
## 12 N Y N N
## 4 N Y N N
## 5 N Y N N
## 32 N Y Y N
## 36 Y Y N N
## 8 N Y Y N
## 15 Y Y Y N
## 19 Y Y N N
## 35 Y Y Y N
## 46 Y Y N Y
## 24 Y Y Y N
## 39 Y Y Y N
## 34 Y Y Y N
## 44 Y Y Y N
## 1 N Y N N
## 2 N Y Y N
## 3 Y Y Y N
## 11 N Y N N
## 30 N Y N N
## 41 Y Y Y N
## 42 N Y Y N
## 7 Y N Y N
## 14 Y Y N N
## 18 Y Y Y N
## 21 Y Y N N
## 22 N Y N N
## 25 Y N Y N
## 31 N Y N N
## 37 Y Y Y N
## 40 Y Y Y N
## 43 Y Y N N
## 6 N Y Y N
## 9 N Y Y N
## 28 Y Y Y N
## 45 Y Y Y N
## 27 N Y N Y
## 38 Y N Y N
## 10 Y Y Y N
incub.dt <- onset.dt - meal.dt
library(MASS)
truehist(as.numeric(incub.dt), nbins = 7, prob = FALSE, ylab = 'Frequency',
col = "skyblue", xlab = "Incubation Period (hours)")
mean(incub.dt, na.rm = TRUE)
## Time difference of 4.295455 hours
median(incub.dt, na.rm = TRUE)
## Time difference of 4 hours
range(incub.dt, na.rm = TRUE)
## Time differences in hours
## [1] 3 7