Julia — Oct 20, 2013, 9:13 PM
#Julia Wei
#Homework Chapter 3
#PH 251D
#10.20.2013
#3.1.
Status <- c("Dead", "Survived", "Dead", "Survived", "Dead", "Survived", "Dead", "Survived")
Treatment <- c("Tolbutamide", "Tolbutamide", "Placebo", "Placebo", "Tolbutamide", "Tolbutamide", "Placebo", "Placebo")
Agegrp <- c("<55", "<55", "<55", "<55", "55+", "55+", "55+", "55+")
Freq <- c(8, 98, 5, 115, 22, 76, 16, 69)
df1 <- data.frame(Status, Treatment, Agegrp, Freq)
df1
Status Treatment Agegrp Freq
1 Dead Tolbutamide <55 8
2 Survived Tolbutamide <55 98
3 Dead Placebo <55 5
4 Survived Placebo <55 115
5 Dead Tolbutamide 55+ 22
6 Survived Tolbutamide 55+ 76
7 Dead Placebo 55+ 16
8 Survived Placebo 55+ 69
#3.2.
first.name <- c("Julia", "Michelle", "Kelly")
last.name <- c("Wei", "Vo", "Adamski")
affiliation <- c("UC Berkeley", "UC Berkeley", "UC Berkeley")
first.email.address <- c("julia.c.wei@gmail.com", "mnessvi@gmail.com", "kelly.adamski@gmail.com")
second.email.address <- c("julia-wei@berkeley.edu", "mtvo@berkeley.edu", "knadamski@berkeley.edu")
date.today <- c("10/19/2013", "10/19/2013", "10/19/2013")
df2 <- data.frame(first.name, last.name, affiliation, first.email.address, second.email.address, date.today)
df2
first.name last.name affiliation first.email.address
1 Julia Wei UC Berkeley julia.c.wei@gmail.com
2 Michelle Vo UC Berkeley mnessvi@gmail.com
3 Kelly Adamski UC Berkeley kelly.adamski@gmail.com
second.email.address date.today
1 julia-wei@berkeley.edu 10/19/2013
2 mtvo@berkeley.edu 10/19/2013
3 knadamski@berkeley.edu 10/19/2013
#3.3.
AIDSdata <- read.table("http://www.medepi.net/data/aids.txt", header=TRUE, sep="", na.strings = ".")
AIDSdata
cases year
1 NA 1980
2 NA 1981
3 NA 1982
4 NA 1983
5 4445 1984
6 8249 1985
7 12932 1986
8 21070 1987
9 31001 1988
10 33722 1989
11 41595 1990
12 43672 1991
13 45472 1992
14 103691 1993
15 78279 1994
16 71547 1995
17 66885 1996
18 58492 1997
19 46521 1998
20 45104 1999
21 40758 2000
22 41868 2001
23 42745 2002
24 44232 2003
plot(AIDSdata$year, AIDSdata$cases, type = "l", xlab = "year", lwd = 2, ylab = "AIDS cases", main = "AIDS cases in the US from 1984-2003")
#3.4.
measlesdata <- read.table("http://www.medepi.net/data/measles.txt", header=TRUE, sep="")
measlesdata
year cases
1 1950 319000
2 1951 530000
3 1952 683000
4 1953 449000
5 1954 683000
6 1955 555000
7 1956 612000
8 1957 487000
9 1958 763000
10 1959 406000
11 1960 442000
12 1961 424000
13 1962 482000
14 1963 385000
15 1964 458000
16 1965 262000
17 1966 204000
18 1967 63000
19 1968 22000
20 1969 26000
21 1970 47351
22 1971 75290
23 1972 32275
24 1973 26690
25 1974 22094
26 1975 24374
27 1976 41126
28 1977 57345
29 1978 26871
30 1979 13597
31 1980 13506
32 1981 3124
33 1982 1714
34 1983 1497
35 1984 2587
36 1985 2822
37 1986 6282
38 1987 3655
39 1988 3396
40 1989 18193
41 1990 27786
42 1991 9643
43 1992 2237
44 1993 312
45 1994 963
46 1995 309
47 1996 508
48 1997 138
49 1998 100
50 1999 100
51 2000 86
52 2001 116
plot(measlesdata$year, measlesdata$cases, type = "l", lwd = 2, xlab = "year", ylab = "measles cases", main = "Measles cases in the US from 1950-2001")
plot(measlesdata$year, measlesdata$cases, type = "l", xlab = "year", log = "y", ylab = "log measles cases", main = "Measles cases in the US from 1950-2001 (semi-logarithmic scale)")
#3.5.
hepatitisbdata <- read.table("http://www.medepi.net/data/hepb.txt", header=TRUE, sep="")
hepatitisbdata
cases year
1 19015 1980
2 21152 1981
3 22177 1982
4 24318 1983
5 26115 1984
6 26611 1985
7 26107 1986
8 25916 1987
9 23177 1988
10 23419 1989
11 21102 1990
12 18003 1991
13 16126 1992
14 13361 1993
15 12517 1994
16 10805 1995
17 10637 1996
18 10416 1997
19 10258 1998
20 7694 1999
21 8036 2000
22 7843 2001
23 7996 2002
24 7526 2003
matplot(hepatitisbdata$year, cbind(hepatitisbdata$cases,AIDSdata$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)
#3.6.
#a.
EvansStudyData <- read.table("http://www.medepi.net/data/evans.txt", header=TRUE, sep="")
EvansStudyData$recodechd <- factor(EvansStudyData$chd, levels = 0:1, labels = c("No", "Yes"))
table(EvansStudyData$recodechd)
No Yes
538 71
EvansStudyData$recodecat <- factor(EvansStudyData$cat, levels = 0:1, labels = c("Normal", "High"))
table(EvansStudyData$recodecat)
Normal High
487 122
EvansStudyData$recodesmk <- factor(EvansStudyData$smk, levels = 0:1, labels = c("Never smoked", "Ever smoked"))
table(EvansStudyData$recodesmk)
Never smoked Ever smoked
222 387
EvansStudyData$recodeecg <- factor(EvansStudyData$ecg, levels = 0:1, labels = c("No abnormality", "Abnormality"))
table(EvansStudyData$recodeecg)
No abnormality Abnormality
443 166
EvansStudyData$recodehpt <- factor(EvansStudyData$hpt, levels = 0:1, labels = c("No", "Yes"))
table(EvansStudyData$recodehpt)
No Yes
354 255
#b.
EvansStudyData$discretizeage <- cut(EvansStudyData$age, quantile(EvansStudyData$age), right = FALSE, include.lowest = TRUE)
table(EvansStudyData$discretizeage)
[40,46) [46,52) [52,60) [60,76]
134 158 158 159
#c.
reclasshpt <- rep(NA, nrow(EvansStudyData))
normal <- EvansStudyData$sbp<120 & EvansStudyData$dbp<80
reclasshpt[normal] <- 1
prehypertension <- EvansStudyData$sbp>=120 & EvansStudyData$sbp<140 | EvansStudyData$dbp>=80 & EvansStudyData$dbp<90
reclasshpt[prehypertension] <- 2
hypertension.stage1 <- EvansStudyData$sbp>=140 & EvansStudyData$sbp<160 | EvansStudyData$dbp>=90 & EvansStudyData$dbp<100
reclasshpt[hypertension.stage1] <- 3
hypertension.stage2 <- EvansStudyData$sbp>= 160 | EvansStudyData$dbp>=100
reclasshpt[hypertension.stage2] <- 4
EvansStudyData$newhptclass <- factor(reclasshpt, levels=1:4, labels = c("normal", "prehypertension", "hypertension stage 1", "hypertension stage 2"))
table(EvansStudyData$newhptclass)
normal prehypertension hypertension stage 1
56 165 177
hypertension stage 2
211
#d.
table("old hypertension"=EvansStudyData$recodehpt, "new hypertension"=EvansStudyData$newhptclass)
new hypertension
old hypertension normal prehypertension hypertension stage 1
No 56 165 133
Yes 0 0 44
new hypertension
old hypertension hypertension stage 2
No 0
Yes 211
#3.7.
wnvdata <- read.table("http://www.medepi.net/data/wnv/wnv2004raw.txt", header=TRUE, sep=",", as.is=TRUE, na.strings=c(".", "Unknown"))
wnvdata$date.onset.international <- as.Date(wnvdata$date.onset, format = "%m/%d/%Y")
wnvdata$date.tested.international <- as.Date(wnvdata$date.tested, format = "%m/%d/%Y")
write.table(wnvdata, "c:/Users/Julia/Documents/wnvdata.txt", sep=".", row.names=FALSE)
#3.8.
#a.
oswego <- read.table("http://www.medepi.net/data/oswego.txt", header=TRUE, sep="", na.strings=".")
str(oswego)
'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/ 6 levels "10:00 PM","11:00 AM",..: 6 3 3 5 5 5 1 4 5 NA ...
$ ill : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 2 2 2 ...
$ onset.date : Factor w/ 2 levels "4/18","4/19": 2 2 2 1 1 2 2 1 2 2 ...
$ onset.time : Factor w/ 17 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/ 2 levels "N","Y": 2 2 1 1 2 1 1 2 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/ 2 levels "N","Y": 1 2 2 1 1 2 2 2 2 2 ...
$ fruit.salad : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
head(oswego)
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
mealdatetime <- paste("4/18/1940", oswego$meal.time)
mealdatetime2 <- strptime(mealdatetime, "%m/%d/%Y %I:%M %p")
onsetdatetime <- paste(paste(oswego$onset.date, "/1940", sep=""), oswego$onset.time)
onsetdatetime2 <- strptime(onsetdatetime, "%m/%d/%Y %I:%M %p")
hist(onsetdatetime2, breaks = 30, freq = TRUE)
#b.
minimum <- which(onsetdatetime2==min(onsetdatetime2,na.rm=T))
minimum
[1] 33
maximum <- which(onsetdatetime2==max(onsetdatetime2,na.rm=T))
maximum
[1] 10
oswego[minimum,]
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
oswego[maximum,]
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
#c.
onset.ct <- as.POSIXct(onsetdatetime2)
oswego2 <- oswego[order(oswego$ill, onset.ct),]
oswego2
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
#d.
incubationdistribution <- onsetdatetime2 - mealdatetime2
library(MASS)
truehist(as.numeric(incubationdistribution), nbins = 7, prob = FALSE, col = "skyblue", xlab = "Incubation Period (hours)")
mean(incubationdistribution, na.rm = TRUE)
Time difference of 4.295 hours
median(incubationdistribution, na.rm = TRUE)
Time difference of 4 hours
range(incubationdistribution, na.rm = TRUE)
Time differences in hours
[1] 3 7