【独家发布】R In Action:Data Analysis and Graphics with R

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

[size=+1]R in Action
Data Analysis and Graphics with R
Second edition of this book is now available

Robert I. Kabacoff

August 2011 | 472 pages
ISBN: 9781935182399

http://www.manning.com/kabacoff/

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Graphics Analysis GRAPHIC Analysi alysis 2011

1. #------------------------------------------------------#
   
2. # R in Action: Chapter 1                               #
   
3. #------------------------------------------------------#
   
4. 
   
5. # In the following code q() is commented out so that
   
6. # you don't quit the session
   
7. 
   
8. # Listing 1.1 - A Sample R session
   
9. 
   
10. age <- c(1, 3, 5, 2, 11, 9, 3, 9, 12, 3)
    
11. weight <- c(4.4, 5.3, 7.2, 5.2, 8.5, 7.3, 6, 10.4,
    
12.     10.2, 6.1)
    
13. mean(weight)
    
14. sd(weight)
    
15. cor(age, weight)
    
16. plot(age, weight)
    
17. # q()
    
18. 
    
19. # Listing 1.2 - An example of commands used to manage
    
20. # the R Workspace. Change the next line to one of your
    
21. # directories
    
22. 
    
23. setwd("C:/myprojects/project1")
    
24. options()
    
25. options(digits=3)
    
26. x <- runif(20)
    
27. summary(x)
    
28. hist(x)
    
29. savehistory()
    
30. save.image()
    
31. # q()
    
32. 
    
33. # Listing 1.3 - Working with a new package
    
34. 
    
35. help.start()
    
36. install.packages("vcd")
    
37. help(package = "vcd")
    
38. library(vcd)
    
39. help(Arthritis)
    
40. Arthritis
    
41. example(Arthritis)
    
42. # q()

复制代码

1. #--------------------------------------------------------#
   
2. # R in Action: Chapter 2                                 #
   
3. #--------------------------------------------------------#
   
4. 
   
5. #  Creating vectors
   
6. 
   
7. a <- c(1, 2, 5, 3, 6, -2, 4)
   
8. b <- c("one", "two", "three")
   
9. c <- c(TRUE, TRUE, TRUE, FALSE, TRUE, FALSE)
   
10. 
    
11. # Using vector subscripts
    
12. 
    
13. a <- c(1, 2, 5, 3, 6, -2, 4)
    
14. a[3]
    
15. a[c(1, 3, 5)]
    
16. a[2:6]
    
17. 
    
18. # Listing 2.1 - Creating Matrices
    
19. 
    
20. y <- matrix(1:20, nrow = 5, ncol = 4)
    
21. y
    
22. cells <- c(1, 26, 24, 68)
    
23. rnames <- c("R1", "R2")
    
24. cnames <- c("C1", "C2")
    
25. mymatrix <- matrix(cells, nrow = 2, ncol = 2, byrow = TRUE,
    
26.     dimnames = list(rnames, cnames))
    
27. mymatrix
    
28. mymatrix <- matrix(cells, nrow = 2, ncol = 2, byrow = FALSE,
    
29.     dimnames = list(rnames, cnames))
    
30. mymatrix
    
31. 
    
32. # Listing 2.2 - Using matrix subscripts
    
33. 
    
34. x <- matrix(1:10, nrow = 2)
    
35. x
    
36. x[2, ]
    
37. x[, 2]
    
38. x[1, 4]
    
39. x[1, c(4, 5)]
    
40. 
    
41. # Listing 2.3 - Creating an array
    
42. 
    
43. dim1 <- c("A1", "A2")
    
44. dim2 <- c("B1", "B2", "B3")
    
45. dim3 <- c("C1", "C2", "C3", "C4")
    
46. z <- array(1:24, c(2, 3, 4), dimnames = list(dim1,
    
47.     dim2, dim3))
    
48. z
    
49. 
    
50. # Listing 2.4 - Creating a dataframe
    
51. 
    
52. patientID <- c(1, 2, 3, 4)
    
53. age <- c(25, 34, 28, 52)
    
54. diabetes <- c("Type1", "Type2", "Type1", "Type1")
    
55. status <- c("Poor", "Improved", "Excellent", "Poor")
    
56. patientdata <- data.frame(patientID, age, diabetes,
    
57.     status)
    
58. patientdata
    
59. 
    
60. # Listing 2.5 - Specifying elements of a dataframe
    
61. 
    
62. patientdata[1:2]
    
63. patientdata[c("diabetes", "status")]
    
64. patientdata$age
    
65. 
    
66. # Listing 2.6 - Using factors
    
67. 
    
68. patientID <- c(1, 2, 3, 4)
    
69. age <- c(25, 34, 28, 52)
    
70. diabetes <- c("Type1", "Type2", "Type1", "Type1")
    
71. status <- c("Poor", "Improved", "Excellent", "Poor")
    
72. diabetes <- factor(diabetes)
    
73. status <- factor(status, order = TRUE)
    
74. patientdata <- data.frame(patientID, age, diabetes,
    
75.     status)
    
76. str(patientdata)
    
77. summary(patientdata)
    
78. 
    
79. #  Listing 2.7 - Creating a list
    
80. 
    
81. g <- "My First List"
    
82. h <- c(25, 26, 18, 39)
    
83. j <- matrix(1:10, nrow = 5)
    
84. k <- c("one", "two", "three")
    
85. mylist <- list(title = g, ages = h, j, k)
    
86. mylist

复制代码

1. #-----------------------------------------------------#
   
2. # R in Action: Chapter 3                              #
   
3. # requires that the Hmisc package has been installed  #
   
4. # install.packages('Hmisc')                           #
   
5. #-----------------------------------------------------#
   
6. 
   
7. # pause after each graph
   
8. par(ask = TRUE)
   
9. 
   
10. # --Section 3.1--
    
11. 
    
12. attach(mtcars)
    
13. plot(wt, mpg)
    
14. abline(lm(mpg ~ wt))
    
15. title("Regression of MPG on Weight")
    
16. detach(mtcars)
    
17. 
    
18. # --Section 3.2--
    
19. 
    
20. dose <- c(20, 30, 40, 45, 60)
    
21. drugA <- c(16, 20, 27, 40, 60)
    
22. drugB <- c(15, 18, 25, 31, 40)
    
23. plot(dose, drugA, type = "b")
    
24. 
    
25. # --Section 3.3--
    
26. 
    
27. opar <- par(no.readonly = TRUE)
    
28. par(lty = 2, pch = 17)
    
29. plot(dose, drugA, type = "b")
    
30. par(opar)
    
31. 
    
32. plot(dose, drugA, type = "b", lty = 2, pch = 17)
    
33. plot(dose, drugA, type = "b", lyt = 3, lwd = 3, pch = 15,
    
34.     cex = 2)
    
35. 
    
36. n <- 10
    
37. mycolors <- rainbow(n)
    
38. pie(rep(1, n), labels = mycolors, col = mycolors)
    
39. mygrays <- gray(0:n/n)
    
40. pie(rep(1, n), labels = mygrays, col = mygrays)
    
41. 
    
42. # Listing 3.1 - Using graphical parameters to control
    
43. # graph appearance
    
44. 
    
45. dose <- c(20, 30, 40, 45, 60)
    
46. drugA <- c(16, 20, 27, 40, 60)
    
47. drugB <- c(15, 18, 25, 31, 40)
    
48. opar <- par(no.readonly = TRUE)
    
49. par(pin = c(2, 3))
    
50. par(lwd = 2, cex = 1.5)
    
51. par(cex.axis = 0.75, font.axis = 3)
    
52. plot(dose, drugA, type = "b", pch = 19, lty = 2, col = "red")
    
53. plot(dose, drugB, type = "b", pch = 23, lty = 6, col = "blue",
    
54.     bg = "green")
    
55. par(opar)
    
56. 
    
57. # --Section 3.4--
    
58. 
    
59. plot(dose, drugA, type = "b", col = "red", lty = 2,
    
60.     pch = 2, lwd = 2, main = "Clinical Trials for Drug A",
    
61.     sub = "This is hypothetical data",
    
62.     xlab = "Dosage", ylab = "Drug Response", xlim = c(0, 60),
    
63.     ylim = c(0, 70))
    
64. 
    
65. 
    
66. # Listing 3.2 - An Example of Custom Axes
    
67. 
    
68. x <- c(1:10)
    
69. y <- x
    
70. z <- 10/x
    
71. opar <- par(no.readonly = TRUE)
    
72. par(mar = c(5, 4, 4, 8) + 0.1)
    
73. 
    
74. plot(x, y, type = "b", pch = 21, col = "red", yaxt = "n",
    
75.     lty = 3, ann = FALSE)
    
76. lines(x, z, type = "b", pch = 22, col = "blue", lty = 2)
    
77. axis(2, at = x, labels = x, col.axis = "red", las = 2)
    
78. axis(4, at = z, labels = round(z, digits = 2), col.axis = "blue",
    
79.     las = 2, cex.axis = 0.7, tck = -0.01)
    
80. mtext("y=1/x", side = 4, line = 3, cex.lab = 1, las = 2,
    
81.     col = "blue")
    
82. title("An Example of Creative Axes", xlab = "X values",
    
83.     ylab = "Y=X")
    
84. par(opar)
    
85. 
    
86. # Listing 3.3 - Comparing Drug A and Drug B response by dose
    
87. 
    
88. dose <- c(20, 30, 40, 45, 60)
    
89. drugA <- c(16, 20, 27, 40, 60)
    
90. drugB <- c(15, 18, 25, 31, 40)
    
91. opar <- par(no.readonly = TRUE)
    
92. par(lwd = 2, cex = 1.5, font.lab = 2)
    
93. plot(dose, drugA, type = "b", pch = 15, lty = 1, col = "red",
    
94.     ylim = c(0, 60), main = "Drug A vs. Drug B", xlab = "Drug Dosage",
    
95.     ylab = "Drug Response")
    
96. lines(dose, drugB, type = "b", pch = 17, lty = 2,
    
97.     col = "blue")
    
98. abline(h = c(30), lwd = 1.5, lty = 2, col = "grey")
    
99. library(Hmisc)
    
100. minor.tick(nx = 3, ny = 3, tick.ratio = 0.5)
     
101. legend("topleft", inset = 0.05, title = "Drug Type",
     
102.     c("A", "B"), lty = c(1, 2), pch = c(15, 17), col = c("red",
     
103.         "blue"))
     
104. par(opar)
     
105. 
     
106. # --Section 3.4.5--
     
107. 
     
108. # Example of labeling points
     
109. 
     
110. attach(mtcars)
     
111. plot(wt, mpg, main = "Milage vs. Car Weight", xlab = "Weight",
     
112.     ylab = "Mileage", pch = 18, col = "blue")
     
113. text(wt, mpg, row.names(mtcars), cex = 0.6, pos = 4,
     
114.     col = "red")
     
115. detach(mtcars)
     
116. 
     
117. # View font families
     
118. opar <- par(no.readonly = TRUE)
     
119. par(cex = 1.5)
     
120. plot(1:7, 1:7, type = "n")
     
121. text(3, 3, "Example of default text")
     
122. text(4, 4, family = "mono", "Example of mono-spaced text")
     
123. text(5, 5, family = "serif", "Example of serif text")
     
124. par(opar)
     
125. 
     
126. # --Section 3.5--
     
127. 
     
128. # combining graphs
     
129. 
     
130. # Figure 3.14
     
131. attach(mtcars)
     
132. opar <- par(no.readonly = TRUE)
     
133. par(mfrow = c(2, 2))
     
134. plot(wt, mpg, main = "Scatterplot of wt vs. mpg")
     
135. plot(wt, disp, main = "Scatterplot of wt vs disp")
     
136. hist(wt, main = "Histogram of wt")
     
137. boxplot(wt, main = "Boxplot of wt")
     
138. par(opar)
     
139. detach(mtcars)
     
140. 
     
141. # Figure 3.15
     
142. attach(mtcars)
     
143. opar <- par(no.readonly = TRUE)
     
144. par(mfrow = c(3, 1))
     
145. hist(wt)
     
146. hist(mpg)
     
147. hist(disp)
     
148. par(opar)
     
149. detach(mtcars)
     
150. 
     
151. # Figure 3.16
     
152. attach(mtcars)
     
153. layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE))
     
154. hist(wt)
     
155. hist(mpg)
     
156. hist(disp)
     
157. detach(mtcars)
     
158. 
     
159. # Figure 3.17
     
160. attach(mtcars)
     
161. layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE),
     
162.     widths = c(3, 1), heights = c(1, 2))
     
163. hist(wt)
     
164. hist(mpg)
     
165. hist(disp)
     
166. detach(mtcars)
     
167. 
     
168. # Listing 3.4 - Fine placement of figures in a graph
     
169. 
     
170. opar <- par(no.readonly = TRUE)
     
171. par(fig = c(0, 0.8, 0, 0.8))
     
172. plot(mtcars$wt, mtcars$mpg, xlab = "Miles Per Gallon",
     
173.     ylab = "Car Weight")
     
174. par(fig = c(0, 0.8, 0.55, 1), new = TRUE)
     
175. boxplot(mtcars$wt, horizontal = TRUE, axes = FALSE)
     
176. par(fig = c(0.65, 1, 0, 0.8), new = TRUE)
     
177. boxplot(mtcars$mpg, axes = FALSE)
     
178. mtext("Enhanced Scatterplot", side = 3, outer = TRUE,
     
179.     line = -3)
     
180. par(opar)

复制代码

1. #---------------------------------------------------------#
   
2. # R in Action: Chapter 4                                  #
   
3. # requires that the reshape and sqldf packages have       #
   
4. # been installed                                          #
   
5. # install.packages(c('reshape', 'sqldf'))                 #
   
6. #---------------------------------------------------------#
   
7. 
   
8. # Listing 4.1 - Creating the leadership data frame
   
9. 
   
10. manager <- c(1, 2, 3, 4, 5)
    
11. date <- c("10/24/08", "10/28/08", "10/1/08", "10/12/08",
    
12.     "5/1/09")
    
13. gender <- c("M", "F", "F", "M", "F")
    
14. age <- c(32, 45, 25, 39, 99)
    
15. q1 <- c(5, 3, 3, 3, 2)
    
16. q2 <- c(4, 5, 5, 3, 2)
    
17. q3 <- c(5, 2, 5, 4, 1)
    
18. q4 <- c(5, 5, 5, NA, 2)
    
19. q5 <- c(5, 5, 2, NA, 1)
    
20. leadership <- data.frame(manager, date, gender, age,
    
21.     q1, q2, q3, q4, q5, stringsAsFactors = FALSE)
    
22.    
    
23. # the individual vectors are no longer needed
    
24. rm(manager, date, gender, age, q1, q2, q3, q4, q5)
    
25. 
    
26. # Listing 4.2 - Creating new variables
    
27. 
    
28. mydata <- data.frame(x1 = c(2, 2, 6, 4), x2 = c(3,
    
29.     4, 2, 8))
    
30. mydata$sumx <- mydata$x1 + mydata$x2
    
31. mydata$meanx <- (mydata$x1 + mydata$x2)/2
    
32. 
    
33. attach(mydata)
    
34. mydata$sumx <- x1 + x2
    
35. mydata$meanx <- (x1 + x2)/2
    
36. detach(mydata)
    
37. 
    
38. mydata <- transform(mydata, sumx = x1 + x2, meanx = (x1 +
    
39.     x2)/2)
    
40. 
    
41. # --Section 4.3--
    
42. 
    
43. # Recoding variables
    
44. 
    
45. leadership$agecat[leadership$age > 75] <- "Elder"
    
46. leadership$agecat[leadership$age > 45 &
    
47.     leadership$age <= 75] <- "Middle Aged"
    
48. leadership$agecat[leadership$age <= 45] <- "Young"
    
49. 
    
50. # or more compactly
    
51. 
    
52. leadership <- within(leadership, {
    
53.     agecat <- NA
    
54.     agecat[age > 75] <- "Elder"
    
55.     agecat[age >= 55 & age <= 75] <- "Middle Aged"
    
56.     agecat[age < 55] <- "Young"
    
57. })
    
58. 
    
59. # --Section 4.4--
    
60. 
    
61. # Renaming variables with the reshape package
    
62. library(reshape)
    
63. rename(leadership, c(manager = "managerID", date = "testDate"))
    
64. 
    
65. # --Section 4.5--
    
66. 
    
67. # Applying the is.na() function
    
68. is.na(leadership[, 6:10])
    
69. 
    
70. # recode 99 to missing for the variable age
    
71. leadership[leadership$age == 99, "age"] <- NA
    
72. leadership
    
73. 
    
74. # Using na.omit() to delete incomplete observations
    
75. newdata <- na.omit(leadership)
    
76. newdata
    
77. 
    
78. # --Section 4.6--
    
79. 
    
80. mydates <- as.Date(c("2007-06-22", "2004-02-13"))
    
81. 
    
82. # Converting character values to dates
    
83. 
    
84. strDates <- c("01/05/1965", "08/16/1975")
    
85. dates <- as.Date(strDates, "%m/%d/%Y")
    
86. 
    
87. myformat <- "%m/%d/%y"
    
88. leadership$date <- as.Date(leadership$date, myformat)
    
89. 
    
90. # Useful date functions
    
91. 
    
92. Sys.Date()
    
93. date()
    
94. 
    
95. today <- Sys.Date()
    
96. format(today, format = "%B %d %Y")
    
97. format(today, format = "%A")
    
98. 
    
99. # Calculations with with dates
    
100. 
     
101. startdate <- as.Date("2004-02-13")
     
102. enddate <- as.Date("2009-06-22")
     
103. days <- enddate - startdate
     
104. 
     
105. # Date functions and formatted printing
     
106. 
     
107. today <- Sys.Date()
     
108. format(today, format = "%B %d %Y")
     
109. dob <- as.Date("1956-10-10")
     
110. format(dob, format = "%A")
     
111. 
     
112. # --Section 4.7--
     
113. 
     
114. # Listing 4.5 - Converting from one data type to another
     
115. 
     
116. a <- c(1, 2, 3)
     
117. a
     
118. is.numeric(a)
     
119. is.vector(a)
     
120. a <- as.character(a)
     
121. a
     
122. is.numeric(a)
     
123. is.vector(a)
     
124. is.character(a)
     
125. 
     
126. # --Section 4.8--
     
127. 
     
128. # Sorting a dataset
     
129. 
     
130. attach(leadership)
     
131. newdata <- leadership[order(age), ]
     
132. newdata
     
133. detach(leadership)
     
134. 
     
135. attach(leadership)
     
136. newdata <- leadership[order(gender, -age), ]
     
137. newdata
     
138. detach(leadership)
     
139. 
     
140. # -- Section 4.10--
     
141. 
     
142. # Selecting variables
     
143. 
     
144. newdata <- leadership[, c(6:10)]
     
145. 
     
146. myvars <- c("q1", "q2", "q3", "q4", "q5")
     
147. newdata <- leadership[myvars]
     
148. 
     
149. myvars <- paste("q", 1:5, sep = "")
     
150. newdata <- leadership[myvars]
     
151. 
     
152. # Dropping variables
     
153. 
     
154. myvars <- names(leadership) %in% c("q3", "q4")
     
155. newdata <- leadership[!myvars]
     
156. 
     
157. newdata <- leadership[c(-7, -8)]
     
158. 
     
159. # You could use the following to delete q3 and q4
     
160. # from the leadership dataset (commented out so
     
161. # the rest of the code in this file will work)
     
162. #
     
163. # leadership$q3 <- leadership$q4 <- NULL
     
164. 
     
165. # Selecting observations
     
166. 
     
167. # Listing 4.6 - Selecting Observations
     
168. 
     
169. newdata <- leadership[1:3, ]
     
170. 
     
171. newdata <- leadership[which(leadership$gender == "M" &
     
172.     leadership$age > 30), ]
     
173. 
     
174. attach(leadership)
     
175. newdata <- leadership[which(leadership$gender == "M" &
     
176.     leadership$age > 30), ]
     
177. detach(leadership)
     
178. 
     
179. # Selecting observations based on dates
     
180. 
     
181. leadership$date <- as.Date(leadership$date, "%m/%d/%y")
     
182. startdate <- as.Date("2009-01-01")
     
183. enddate <- as.Date("2009-10-31")
     
184. newdata <- leadership[leadership$date >= startdate &
     
185.     leadership$date <= enddate, ]
     
186. 
     
187. # Using the subset() function
     
188. 
     
189. newdata <- subset(leadership, age >= 35 | age < 24,
     
190.     select = c(q1, q2, q3, q4))
     
191. newdata <- subset(leadership, gender == "M" & age >
     
192.     25, select = gender:q4)
     
193. 
     
194. # --Section 4.11--
     
195. 
     
196. # Listing 4.7 - Using SQL statements to manipulate data frames
     
197. 
     
198. library(sqldf)
     
199. newdf <- sqldf("select * from mtcars where carb=1 order by mpg",
     
200.     row.names = TRUE)
     
201. newdf <- sqldf("select avg(mpg) as avg_mpg, avg(disp) as avg_disp,
     
202.     gear from mtcars where cyl in (4, 6) group by gear")
     
203.    

复制代码

1. #----------------------------------------------------#
   
2. # R in Action: Chapter 5                             #
   
3. #----------------------------------------------------#
   
4. 
   
5. # Listing 5.1 -  Calculating the mean and
   
6. # standard deviation
   
7. 
   
8. x <- c(1, 2, 3, 4, 5, 6, 7, 8)
   
9. mean(x)
   
10. sd(x)
    
11. n <- length(x)
    
12. meanx <- sum(x)/n
    
13. css <- sum((x - meanx)**2)            
    
14. sdx <- sqrt(css / (n-1))
    
15. meanx
    
16. sdx
    
17. 
    
18. # Listing 5.2 - Generating pseudo-random numbers from
    
19. # a uniform distribution
    
20. 
    
21. runif(5)
    
22. runif(5)
    
23. set.seed(1234)                                                     
    
24. runif(5)
    
25. set.seed(1234)                                                      
    
26. runif(5)
    
27. 
    
28. # Listing 5.3 - Generating data from a multivariate
    
29. # normal distribution
    
30. 
    
31. library(MASS)
    
32. options(digits=3)
    
33. set.seed(1234)
    
34. 
    
35. mean <- c(230.7, 146.7, 3.6)                                          
    
36. sigma <- matrix( c(15360.8, 6721.2, -47.1,                              
    
37.                     6721.2, 4700.9, -16.5,
    
38.                      -47.1,  -16.5,   0.3), nrow=3, ncol=3)
    
39. 
    
40. mydata <- mvrnorm(500, mean, sigma)                                    
    
41. mydata <- as.data.frame(mydata)                                         
    
42. names(mydata) <- c("y", "x1", "x2")                                       
    
43. 
    
44. dim(mydata)                                                            
    
45. head(mydata, n=10)   
    
46. 
    
47. # Listing 5.4 - Applying functions to data objects
    
48. 
    
49. a <- 5
    
50. sqrt(a)
    
51. b <- c(1.243, 5.654, 2.99)
    
52. round(b)
    
53. c <- matrix(runif(12), nrow=3)
    
54. c
    
55. log(c)
    
56. mean(c)
    
57. 
    
58. # Listing 5.5 - Applying a function to the rows
    
59. # (columns) of a matrix
    
60. 
    
61. mydata <- matrix(rnorm(30), nrow=6)
    
62. mydata
    
63. apply(mydata, 1, mean)     
    
64. apply(mydata, 2, mean)
    
65. apply(mydata, 2, mean, trim=.4)   
    
66. 
    
67. # Listing 5.6 - A solution to the learning example
    
68. 
    
69. options(digits=2)
    
70. 
    
71. Student <- c("John Davis", "Angela Williams",
    
72.     "Bullwinkle Moose", "David Jones",
    
73.     "Janice Markhammer", "Cheryl Cushing",
    
74.     "Reuven Ytzrhak", "Greg Knox", "Joel England",
    
75.     "Mary Rayburn")
    
76. Math <- c(502, 600, 412, 358, 495, 512, 410, 625, 573, 522)
    
77. Science <- c(95, 99, 80, 82, 75, 85, 80, 95, 89, 86)
    
78. English <- c(25, 22, 18, 15, 20, 28, 15, 30, 27, 18)
    
79. roster <- data.frame(Student, Math, Science, English,
    
80.     stringsAsFactors=FALSE)
    
81.    
    
82. z <- scale(roster[,2:4])
    
83. score <- apply(z, 1, mean)                                            
    
84. roster <- cbind(roster, score)
    
85. 
    
86. y <- quantile(score, c(.8,.6,.4,.2))                                   
    
87. roster$grade[score >= y[1]] <- "A"                                    
    
88. roster$grade[score < y[1] & score >= y[2]] <- "B"
    
89. roster$grade[score < y[2] & score >= y[3]] <- "C"
    
90. roster$grade[score < y[3] & score >= y[4]] <- "D"
    
91. roster$grade[score < y[4]] <- "F"
    
92. 
    
93. name <- strsplit((roster$Student), " ")                                
    
94. lastname <- sapply(name, "[", 2)
    
95. firstname <- sapply(name, "[", 1)
    
96. 
    
97. roster <- cbind(firstname,lastname, roster[,-1])
    
98. roster <- roster[order(lastname,firstname),]
    
99. 
    
100. roster
     
101. 
     
102. # Listing 5.7 - A switch example
     
103. 
     
104. feelings <- c("sad", "afraid")
     
105. for (i in feelings)
     
106.     print(
     
107.       switch(i,
     
108.              happy  = "I am glad you are happy",
     
109.              afraid = "There is nothing to fear",
     
110.              sad    = "Cheer up",
     
111.              angry  = "Calm down now"
     
112.     )
     
113.   )
     
114.   
     
115. # Listing 5.8 - mystats(): a user-written function
     
116. # for summary statistics
     
117. 
     
118. mystats <- function(x, parametric=TRUE, print=FALSE) {
     
119.   if (parametric) {
     
120.     center <- mean(x); spread <- sd(x)
     
121.   } else {
     
122.     center <- median(x); spread <- mad(x)
     
123.   }
     
124.   if (print & parametric) {
     
125.     cat("Mean=", center, "\n", "SD=", spread, "\n")
     
126.   } else if (print & !parametric) {
     
127.     cat("Median=", center, "\n", "MAD=", spread, "\n")
     
128.   }
     
129.   result <- list(center=center, spread=spread)
     
130.   return(result)
     
131. }
     
132. 
     
133. # trying it out
     
134. set.seed(1234)
     
135. x <- rnorm(500)
     
136. y <- mystats(x)
     
137. y <- mystats(x, parametric=FALSE, print=TRUE)
     
138. 
     
139. # Another switch example
     
140. mydate <- function(type="long") {
     
141.     switch(type,
     
142.     long =  format(Sys.time(), "%A %B %d %Y"),
     
143.     short = format(Sys.time(), "%m-%d-%y"),
     
144.     cat(type, "is not a recognized type\n"))
     
145. }
     
146. mydate("long")
     
147. mydate("short")
     
148. mydate()
     
149. mydate("medium")
     
150. 
     
151. # Listing 5.9 - Transposing a dataset
     
152. 
     
153. cars <- mtcars[1:5, 1:4]      
     
154. cars
     
155. t(cars)
     
156. 
     
157. # Listing 5.10 - Aggregating data
     
158. 
     
159. options(digits=3)
     
160. attach(mtcars)
     
161. aggdata <-aggregate(mtcars, by=list(cyl,gear),
     
162.     FUN=mean, na.rm=TRUE)
     
163. aggdata

复制代码

1. #----------------------------------------------------------------#
   
2. # R in Action: Chapter 6                                         #
   
3. # requires that the vcd, plotrix, sm, vioplot packages have been #
   
4. # installed                                                      #
   
5. # install.packages(c('vcd', 'plotrix', 'sm', 'vioplot'))         #
   
6. #----------------------------------------------------------------#
   
7. 
   
8. # pause after each graph
   
9. par(ask = TRUE)
   
10. 
    
11. # save original graphic settings
    
12. opar <- par(no.readonly = TRUE)
    
13. 
    
14. # Load vcd package
    
15. library(vcd)
    
16. 
    
17. # Get cell counts for improved variable
    
18. counts <- table(Arthritis$Improved)
    
19. counts
    
20. 
    
21. # Listing 6.1 - Simple bar plot
    
22. 
    
23. # simple bar plot
    
24. barplot(counts, main = "Simple Bar Plot", xlab = "Improvement",
    
25.     ylab = "Frequency")
    
26. 
    
27. # horizontal bar plot
    
28. barplot(counts, main = "Horizontal Bar Plot", xlab = "Frequency",
    
29.     ylab = "Improvement", horiz = TRUE)
    
30. 
    
31. 
    
32. # get counts for Improved by Treatment table
    
33. counts <- table(Arthritis$Improved, Arthritis$Treatment)
    
34. counts
    
35. 
    
36. # Listing 6.2 - Stacked and groupde bar plots
    
37. 
    
38. # stacked barplot
    
39. barplot(counts, main = "Stacked Bar Plot", xlab = "Treatment",
    
40.     ylab = "Frequency", col = c("red", "yellow", "green"),
    
41.     legend = rownames(counts))
    
42. 
    
43. # grouped barplot
    
44. barplot(counts, main = "Grouped Bar Plot", xlab = "Treatment",
    
45.     ylab = "Frequency", col = c("red", "yellow", "green"),
    
46.     legend = rownames(counts),
    
47.     beside = TRUE)
    
48. 
    
49. 
    
50. # Listing 6.3 - Mean bar plots
    
51. 
    
52. states <- data.frame(state.region, state.x77)
    
53. means <- aggregate(states$Illiteracy,
    
54.     by = list(state.region),
    
55.     FUN = mean)
    
56. means
    
57. 
    
58. means <- means[order(means$x), ]
    
59. means
    
60. 
    
61. barplot(means$x, names.arg = means$Group.1)
    
62. title("Mean Illiteracy Rate")
    
63. 
    
64. # Listing 6.4 - Fitting labels in bar plots
    
65. 
    
66. par(mar = c(5, 8, 4, 2))
    
67. par(las = 2)
    
68. counts <- table(Arthritis$Improved)
    
69. 
    
70. barplot(counts, main = "Treatment Outcome", horiz = TRUE,
    
71.     cex.names = 0.8, names.arg = c("No Improvement",
    
72.     "Some Improvement", "Marked Improvement"))
    
73. 
    
74. # Section --6.1.5 Spinograms--
    
75. 
    
76. library(vcd)
    
77. attach(Arthritis)
    
78. counts <- table(Treatment, Improved)
    
79. spine(counts, main = "Spinogram Example")
    
80. detach(Arthritis)
    
81. 
    
82. 
    
83. # Listing 6.5 - Pie charts
    
84. 
    
85. par(mfrow = c(2, 2))
    
86. slices <- c(10, 12, 4, 16, 8)
    
87. lbls <- c("US", "UK", "Australia", "Germany", "France")
    
88. 
    
89. pie(slices, labels = lbls, main = "Simple Pie Chart")
    
90. 
    
91. pct <- round(slices/sum(slices) * 100)
    
92. lbls2 <- paste(lbls, " ", pct, "%", sep = "")
    
93. pie(slices, labels = lbls2, col = rainbow(length(lbls)),
    
94.     main = "Pie Chart with Percentages")
    
95. 
    
96. library(plotrix)
    
97. pie3D(slices, labels = lbls, explode = 0.1, main = "3D Pie Chart ")
    
98. 
    
99. mytable <- table(state.region)
    
100. lbls <- paste(names(mytable), "\n", mytable, sep = "")
     
101. pie(mytable, labels = lbls,
     
102.     main = "Pie Chart from a Table\n (with sample sizes)")
     
103. 
     
104. # restore original graphic parameters
     
105. par(opar)
     
106. 
     
107. # fan plots
     
108. library(plotrix)
     
109. slices <- c(10, 12, 4, 16, 8)
     
110. lbls <- c("US", "UK", "Australia", "Germany", "France")
     
111. fan.plot(slices, labels = lbls, main = "Fan Plot")
     
112. 
     
113. 
     
114. # Listing 6.6 - Histograms
     
115. 
     
116. par(mfrow = c(2, 2))
     
117. 
     
118. hist(mtcars$mpg)
     
119. 
     
120. hist(mtcars$mpg, breaks = 12, col = "red",
     
121.     xlab = "Miles Per Gallon",
     
122.     main = "Colored histogram with 12 bins")
     
123. 
     
124. hist(mtcars$mpg, freq = FALSE, breaks = 12, col = "red",
     
125.     xlab = "Miles Per Gallon",
     
126.     main = "Histogram, rug plot, density curve")
     
127. rug(jitter(mtcars$mpg))
     
128. lines(density(mtcars$mpg), col = "blue", lwd = 2)
     
129. 
     
130. # Histogram with Superimposed Normal Curve
     
131. # (Thanks to Peter Dalgaard)
     
132. x <- mtcars$mpg
     
133. h <- hist(x, breaks = 12, col = "red",
     
134.     xlab = "Miles Per Gallon",
     
135.     main = "Histogram with normal curve and box")
     
136. xfit <- seq(min(x), max(x), length = 40)
     
137. yfit <- dnorm(xfit, mean = mean(x), sd = sd(x))
     
138. yfit <- yfit * diff(h$mids[1:2]) * length(x)
     
139. lines(xfit, yfit, col = "blue", lwd = 2)
     
140. box()
     
141. 
     
142. # restore original graphic parameters
     
143. par(opar)
     
144. 
     
145. # Listing 6.7 - Kernel density plot
     
146. 
     
147. par(mfrow = c(2, 1))
     
148. d <- density(mtcars$mpg)
     
149. 
     
150. plot(d)
     
151. 
     
152. d <- density(mtcars$mpg)
     
153. plot(d, main = "Kernel Density of Miles Per Gallon")
     
154. polygon(d, col = "red", border = "blue")
     
155. rug(mtcars$mpg, col = "brown")
     
156. 
     
157. # restore original graphic parameters
     
158. par(opar)
     
159. 
     
160. # Listing 6.8 - Comparing kernel density plots
     
161. 
     
162. par(lwd = 2)
     
163. library(sm)
     
164. attach(mtcars)
     
165. 
     
166. cyl.f <- factor(cyl, levels = c(4, 6, 8),
     
167.     labels = c("4 cylinder", "6 cylinder", "8 cylinder"))
     
168. 
     
169. sm.density.compare(mpg, cyl, xlab = "Miles Per Gallon")
     
170. title(main = "MPG Distribution by Car Cylinders")
     
171. 
     
172. colfill <- c(2:(2 + length(levels(cyl.f))))
     
173. cat("Use mouse to place legend...", "\n\n")
     
174. legend(locator(1), levels(cyl.f), fill = colfill)
     
175. detach(mtcars)
     
176. par(lwd = 1)
     
177. 
     
178. # --Section 6.5--
     
179. 
     
180. boxplot(mpg ~ cyl, data = mtcars,
     
181.     main = "Car Milage Data",
     
182.     xlab = "Number of Cylinders",
     
183.     ylab = "Miles Per Gallon")
     
184. 
     
185. boxplot(mpg ~ cyl, data = mtcars, notch = TRUE,
     
186.     varwidth = TRUE, col = "red",
     
187.     main = "Car Mileage Data",
     
188.     xlab = "Number of Cylinders",
     
189.     ylab = "Miles Per Gallon")
     
190. 
     
191. # Listing 6.9 - Box plots for two crossed factors
     
192. 
     
193. mtcars$cyl.f <- factor(mtcars$cyl, levels = c(4, 6,
     
194.     8), labels = c("4", "6", "8"))
     
195. 
     
196. mtcars$am.f <- factor(mtcars$am, levels = c(0, 1),
     
197.     labels = c("auto", "standard"))
     
198. 
     
199. boxplot(mpg ~ am.f * cyl.f, data = mtcars,
     
200.     varwidth = TRUE, col = c("gold", "darkgreen"),
     
201.     main = "MPG Distribution by Auto Type",
     
202.     xlab = "Auto Type")
     
203. 
     
204. # Listing 6.10 - Violin plots
     
205. 
     
206. library(vioplot)
     
207. x1 <- mtcars$mpg[mtcars$cyl == 4]
     
208. x2 <- mtcars$mpg[mtcars$cyl == 6]
     
209. x3 <- mtcars$mpg[mtcars$cyl == 8]
     
210. vioplot(x1, x2, x3,
     
211.     names = c("4 cyl", "6 cyl", "8 cyl"),
     
212.     col = "gold")
     
213. title("Violin Plots of Miles Per Gallon")
     
214. 
     
215. # --Section 6.6--
     
216. 
     
217. dotchart(mtcars$mpg, labels = row.names(mtcars),
     
218.     cex = 0.7,
     
219.     main = "Gas Milage for Car Models",
     
220.     xlab = "Miles Per Gallon")
     
221. 
     
222. # Listing 6.11 - sorted colored grouped dot chart
     
223. 
     
224. x <- mtcars[order(mtcars$mpg), ]
     
225. x$cyl <- factor(x$cyl)
     
226. x$color[x$cyl == 4] <- "red"
     
227. x$color[x$cyl == 6] <- "blue"
     
228. x$color[x$cyl == 8] <- "darkgreen"
     
229. dotchart(x$mpg, labels = row.names(x), cex = 0.7,
     
230.     pch = 19, groups = x$cyl,
     
231.     gcolor = "black", color = x$color,
     
232.     main = "Gas Milage for Car Models\ngrouped by cylinder",
     
233.     xlab = "Miles Per Gallon")
     
234.    

复制代码

1. #--------------------------------------------------------------------#
   
2. # R in Action: Chapter 7                                             #
   
3. # requires that the npmc, ggm, gmodels, vcd, Hmisc,                  #
   
4. # pastecs, psych, doBy, and reshape packages have been installed     #
   
5. # install.packages(c('npmc', 'ggm', 'gmodels', 'vcd', 'Hmisc',       #
   
6. #     'pastecs', 'psych', 'doBy', 'reshape'))                        #
   
7. #---------------------------------------------------------------------
   
8. 
   
9. vars <- c("mpg", "hp", "wt")
   
10. head(mtcars[vars])
    
11. 
    
12. # Listing 7.1 - descriptive stats via summary
    
13. 
    
14. summary(mtcars[vars])
    
15. 
    
16. # Listing 7.2 - descriptive stats via sapply()
    
17. 
    
18. mystats <- function(x, na.omit = FALSE) {
    
19.     if (na.omit)
    
20.         x <- x[!is.na(x)]
    
21.     m <- mean(x)
    
22.     n <- length(x)
    
23.     s <- sd(x)
    
24.     skew <- sum((x - m)^3/s^3)/n
    
25.     kurt <- sum((x - m)^4/s^4)/n - 3
    
26.     return(c(n = n, mean = m, stdev = s, skew = skew, kurtosis = kurt))
    
27. }
    
28. 
    
29. sapply(mtcars[vars], mystats)
    
30. 
    
31. # Listing 7.3 - Descriptive statistics (Hmisc package)
    
32. 
    
33. library(Hmisc)
    
34. describe(mtcars[vars])
    
35. 
    
36. # Listing 7.4 - Descriptive statistics (pastecs package)
    
37. 
    
38. library(pastecs)
    
39. stat.desc(mtcars[vars])
    
40. 
    
41. # Listing 7.5 - Descriptive statistics (psych package)
    
42. 
    
43. library(psych)
    
44. describe(mtcars[vars])
    
45. 
    
46. # Listing 7.6 - Descriptive statistics by group with aggregate()
    
47. 
    
48. aggregate(mtcars[vars], by = list(am = mtcars$am), mean)
    
49. aggregate(mtcars[vars], by = list(am = mtcars$am), sd)
    
50. 
    
51. # Listing 7.7 - Descriptive statistics by group via by()
    
52. 
    
53. dstats <- function(x)(c(mean=mean(x), sd=sd(x)))
    
54. by(mtcars[vars], mtcars$am, dstats)
    
55. 
    
56. # Listing 7.8 Summary statists by group (doBy package)
    
57. 
    
58. library(doBy)
    
59. summaryBy(mpg + hp + wt ~ am, data = mtcars, FUN = mystats)
    
60. 
    
61. # Listing 7.9 - Summary statistics by group (psych package)
    
62. 
    
63. library(psych)
    
64. describe.by(mtcars[vars], mtcars$am)
    
65. 
    
66. # Listing 1.10 Summary statistics by group (reshape package)
    
67. 
    
68. library(reshape)
    
69. dstats <- function(x) (c(n = length(x), mean = mean(x),
    
70.     sd = sd(x)))
    
71. dfm <- melt(mtcars, measure.vars = c("mpg", "hp",
    
72.     "wt"), id.vars = c("am", "cyl"))
    
73. cast(dfm, am + cyl + variable ~ ., dstats)
    
74. 
    
75. # Section --7.2--
    
76. 
    
77. # get Arthritis data
    
78. library(vcd)
    
79. 
    
80. # one way table
    
81. 
    
82. mytable <- with(Arthritis, table(Improved))
    
83. mytable
    
84. prop.table(mytable)
    
85. prop.table(mytable)*100
    
86. 
    
87. 
    
88. # two way table
    
89. 
    
90. mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
    
91. mytable
    
92. margin.table(mytable, 1)
    
93. prop.table(mytable, 1)
    
94. margin.table(mytable, 2)
    
95. prop.table(mytable, 2)
    
96. prop.table(mytable)
    
97. addmargins(mytable)
    
98. admargins(prop.table(mytable))
    
99. addmargins(prop.table(mytable, 1), 2)
    
100. addmargins(prop.table(mytable, 2, 1)
     
101. 
     
102. # Listing 7.11 - Two-way table using CrossTable
     
103. 
     
104. library(gmodels)
     
105. CrossTable(Arthritis$Treatment, Arthritis$Improved)
     
106. 
     
107. # Listing 7.12 - Three-way contingency table
     
108. 
     
109. mytable <- xtabs(~ Treatment+Sex+Improved, data=Arthritis)
     
110. mytable
     
111. ftable(mytable)
     
112. margin.table(mytable, 1)
     
113. margin.table(mytable, 2)
     
114. margin.table(mytable, 3)
     
115. margin.table(mytable, c(1,3))
     
116. ftable(prop.table(mytable, c(1, 2)))
     
117. ftable(addmargins(prop.table(mytable, c(1, 2)), 3))
     
118. 
     
119. gtable(addmargins(prop.table(mytable, c(1, 2)), 3)) * 100
     
120. 
     
121. # Listing 7.13 - Chis-square test of independence
     
122. 
     
123. library(vcd)
     
124. mytable <- xtabs(~Treatment+Improved, data=Arthritis)
     
125. chisq.test(mytable)
     
126. mytable <- xtabs(~Improved+Sex, data=Arthritis)
     
127. chisq.test(mytable)
     
128. 
     
129. # Fisher's exact test
     
130. 
     
131. mytable <- xtabs(~Treatment+Improved, data=Arthritis)
     
132. fisher.test(mytable)
     
133. 
     
134. # Chochran-Mantel-Haenszel test
     
135. 
     
136. mytable <- xtabs(~Treatment+Improved+Sex, data=Arthritis)
     
137. mantelhaen.test(mytable)
     
138. 
     
139. # Listing 7.14 - Measures of association for a two-way table
     
140. 
     
141. library(vcd)
     
142. mytable <- xtabs(~Treatment+Improved, data=Arthritis)
     
143. assocstats(mytable)
     
144. 
     
145. 
     
146. # Listing 7.15 - converting a table into a flat file via table2flat
     
147. 
     
148. table2flat <- function(mytable) {
     
149.     df <- as.data.frame(mytable)
     
150.     rows <- dim(df)[1]
     
151.     cols <- dim(df)[2]
     
152.     x <- NULL
     
153.     for (i in 1:rows) {
     
154.         for (j in 1:df$Freq[i]) {
     
155.             row <- df[i, c(1:(cols - 1))]
     
156.             x <- rbind(x, row)
     
157.         }
     
158.     }
     
159.     row.names(x) <- c(1:dim(x)[1])
     
160.     return(x)
     
161. }
     
162. 
     
163. # Listing 7.16 - Using table2flat with published data
     
164. 
     
165. treatment <- rep(c("Placebo", "Treated"), 3)
     
166. improved <- rep(c("None", "Some", "Marked"), each = 2)
     
167. Freq <- c(29, 13, 7, 7, 7, 21)
     
168. mytable <- as.data.frame(cbind(treatment, improved, Freq))
     
169. mydata <- table2flat(mytable)
     
170. head(mydata)
     
171. 
     
172. # Listing 7.17 - Covariances and correlations
     
173. 
     
174. states <- state.x77[, 1:6]
     
175. cov(states)
     
176. cor(states)
     
177. cor(states, method="spearman")
     
178. 
     
179. x <- states[, c("Population", "Income", "Illiteracy", "HS Grad")]
     
180. y <- states[, c("Life Exp", "Murder")]
     
181. cor(x, y)
     
182. 
     
183. # partial correlation of population and murder rate, controlling
     
184. # for income, illiteracy rate, and HS graduation rate
     
185. 
     
186. library(ggm)
     
187. pcor(c(1, 5, 2, 3, 6), cov(states))
     
188. 
     
189. # Listing 7.18 - Testing correlations for significance
     
190. 
     
191. cor.test(states[, 3], states[, 5])
     
192. 
     
193. # Listing 7.19 - Correlation matrix and tests of significance via corr.test
     
194. 
     
195. library(psych)
     
196. corr.test(states, use = "complete")
     
197. 
     
198. # --Section 7.4--
     
199. 
     
200. # independent t-test
     
201. 
     
202. library(MASS)
     
203. t.test(Prob ~ So, data=UScrime)
     
204. 
     
205. # dependent t-test
     
206. 
     
207. library(MASS)
     
208. sapply(UScrime[c("U1", "U2")], function(x) (c(mean = mean(x),
     
209.     sd = sd(x))))
     
210. with(UScrime, t.test(U1, U2, paired = TRUE))
     
211. 
     
212. # --Section 7.5--
     
213. 
     
214. # Wilcoxon two group comparison
     
215. 
     
216. with(UScrime, by(Prob, So, median))
     
217. wilcox.test(Prob ~ So, data=UScrime)
     
218. sapply(UScrime[c("U1", "U2")], median)
     
219. with(UScrime, wilcox.test(U1, U2, paired = TRUE))
     
220. 
     
221. # Kruskal Wallis test
     
222. 
     
223. states <- as.data.frame(cbind(state.region, state.x77))
     
224. kruskal.test(Illiteracy ~ state.region, data=states)
     
225. 
     
226. # Listing 7.20 - Nonparametric multiple comparisons
     
227. 
     
228. class <- state.region
     
229. var <- state.x77[, c("Illiteracy")]
     
230. mydata <- as.data.frame(cbind(class, var))
     
231. rm(class,var)
     
232. library(npmc)
     
233. summary(npmc(mydata), type = "BF")
     
234. aggregate(mydata, by = list(mydata$class), median)

复制代码

1. #------------------------------------------------------------------------#
   
2. # R in Action: Chapter 8                                                 #
   
3. # requires that the car, gvlma, MASS, leaps packages have been installed #
   
4. # install.packages(c('car', 'gvlma', 'MASS', 'leaps'))                   #
   
5. #------------------------------------------------------------------------#
   
6. 
   
7. # pause on each graph
   
8. par(ask = TRUE)
   
9. 
   
10. # save current graphical parameters
    
11. opar <- par(no.readonly = TRUE)
    
12. 
    
13. # Listing 8.1 - simple linear regression
    
14. 
    
15. fit <- lm(weight ~ height, data = women)
    
16. summary(fit)
    
17. women$weight
    
18. fitted(fit)
    
19. residuals(fit)
    
20. 
    
21. # scatter plot of height by weight
    
22. 
    
23. plot(women$height, women$weight, main = "Women Age 30-39",
    
24.     xlab = "Height (in inches)", ylab = "Weight (in pounds)")
    
25. # add the line of best fit
    
26. abline(fit)
    
27. 
    
28. # Listing 8.2 - Polynomial regression
    
29. 
    
30. fit2 <- lm(weight ~ height + I(height^2), data = women)
    
31. summary(fit2)
    
32. 
    
33. plot(women$height, women$weight, main = "Women Age 30-39",
    
34.     xlab = "Height (in inches)", ylab = "Weight (in lbs)")
    
35. lines(women$height, fitted(fit2))
    
36. 
    
37. # scatterplot for women data
    
38. 
    
39. library(car)
    
40. scatterplot(weight ~ height, data = women, spread = FALSE,
    
41.     lty.smooth = 2, pch = 19, main = "Women Age 30-39", xlab = "Height (inches)",
    
42.     ylab = "Weight (lbs.)")
    
43. 
    
44. # Listing 8.3 - Examining bivariate relationship
    
45. 
    
46. states <- as.data.frame(state.x77[, c("Murder", "Population",
    
47.     "Illiteracy", "Income", "Frost")])
    
48.    
    
49. cor(states)
    
50. 
    
51. library(car)
    
52. scatterplotMatrix(states, spread = FALSE, lty.smooth = 2,
    
53.     main = "Scatterplot Matrix")
    
54.    
    
55. # Listing 8.4 - Multiple linear regression
    
56. 
    
57. fit <- lm(Murder ~ Population + Illiteracy + Income +
    
58.     Frost, data = states)
    
59.    
    
60. # Listing 8.5 Multiple linear regression with a significant
    
61. # interaction term
    
62. 
    
63. fit <- lm(mpg ~ hp + wt + hp:wt, data = mtcars)
    
64. summary(fit)
    
65. 
    
66. library(effects)
    
67. plot(effect("hp:wt", fit, list(wt = c(2.2, 3.2, 4.2))),
    
68.     multiline = TRUE)
    
69.    
    
70. # --Section 8.3--
    
71. 
    
72. fit <- lm(Murder ~ Population + Illiteracy + Income +
    
73.     Frost, data=states)
    
74. confint(fit)
    
75. 
    
76. # simple regression diagnostics
    
77. 
    
78. fit <- lm(weight ~ height, data = women)
    
79. par(mfrow = c(2, 2))
    
80. plot(fit)
    
81. par(opar)
    
82. 
    
83. # regression diagnostics for quadratic fit
    
84. 
    
85. newfit <- lm(weight ~ height + I(height^2), data = women)
    
86. par(mfrow = c(2, 2))
    
87. plot(newfit)
    
88. par(opar)
    
89. 
    
90. # regression diagnostics for quadratic fit
    
91. # with deleted observations
    
92. 
    
93. newfit <- lm(weight ~ height + I(height^2), data = women[-c(13, 15),])
    
94. par(mfrow = c(2, 2))
    
95. plot(newfit)
    
96. par(opar)
    
97. 
    
98. # basic regression diagnostics for states data
    
99. 
    
100. fit <- lm(Murder ~ Population + Illiteracy + Income +
     
101.     Frost, data = states)
     
102. par(mfrow = c(2, 2))
     
103. plot(fit)
     
104. par(opar)
     
105. 
     
106. # Assessing normality
     
107. library(car)
     
108. fit <- lm(Murder ~ Population + Illiteracy + Income +
     
109.     Frost, data = states)
     
110. qqPlot(fit, labels = FALSE, simulate = TRUE, main = "Q-Q Plot")
     
111. 
     
112. 
     
113. # Listing 8.6 Function for plotting studentized residuals
     
114. 
     
115. residplot <- function(fit, nbreaks=10) {
     
116.     z <- rstudent(fit)
     
117.     hist(z, breaks=nbreaks, freq=FALSE,
     
118.     xlab="Studentized Residual",
     
119.     main="Distribution of Errors")
     
120.     rug(jitter(z), col="brown")
     
121.     curve(dnorm(x, mean=mean(z), sd=sd(z)),
     
122.         add=TRUE, col="blue", lwd=2)
     
123.     lines(density(z)$x, density(z)$y,
     
124.         col="red", lwd=2, lty=2)
     
125.     legend("topright",
     
126.         legend = c( "Normal Curve", "Kernel Density Curve"),
     
127.         lty=1:2, col=c("blue","red"), cex=.7)
     
128. }
     
129. 
     
130. residplot(fit)
     
131. 
     
132. #  Durbin Watson test for Autocorrelated Errors
     
133. 
     
134. durbinWatsonTest(fit)
     
135. 
     
136. # Component + Residual Plots
     
137. 
     
138. crPlots(fit, one.page = TRUE, ask = FALSE)
     
139. 
     
140. # Listing 8.7 - Assessing homoscedasticity
     
141. 
     
142. library(car)
     
143. ncvTest(fit)
     
144. spreadLevelPlot(fit)
     
145. 
     
146. # Listing 8.8 - Global test of linear model assumptions
     
147. 
     
148. library(gvlma)
     
149. gvmodel <- gvlma(fit)
     
150. summary(gvmodel)
     
151. 
     
152. # Library 8.9 - Evaluating multi-collinearity
     
153. 
     
154. vif(fit)
     
155. sqrt(vif(fit)) > 2
     
156. 
     
157. # --Section 8.4--
     
158. 
     
159. # Assessing outliers
     
160. 
     
161. library(car)
     
162. outlierTest(fit)
     
163. 
     
164. # Index plot of hat values
     
165. # use the mouse to identify points interactively
     
166. 
     
167. hat.plot <- function(fit){
     
168.     p <- length(coefficients(fit))
     
169.     n <- length(fitted(fit))
     
170.     plot(hatvalues(fit), main = "Index Plot of Hat Values")
     
171.     abline(h = c(2, 3) * p/n, col = "red", lty = 2)
     
172.     identify(1:n, hatvalues(fit), names(hatvalues(fit)))
     
173. }
     
174. 
     
175. hat.plot(fit)
     
176. 
     
177. # Cook's D Plot
     
178. # identify D values > 4/(n-k-1)
     
179. 
     
180. cutoff <- 4/(nrow(states) - length(fit$coefficients) - 2)
     
181. plot(fit, which = 4, cook.levels = cutoff)
     
182. abline(h = cutoff, lty = 2, col = "red")
     
183. 
     
184. # Added variable plots
     
185. # use the mouse to identify points interactively
     
186. 
     
187. avPlots(fit, ask = FALSE, onepage = TRUE, id.method = "identify")
     
188. 
     
189. # Influence Plot
     
190. # use the mouse to identify points interactively
     
191. 
     
192. influencePlot(fit, id.method = "identify", main = "Influence Plot",
     
193.     sub = "Circle size is proportial to Cook's Distance")
     
194. 
     
195. # Listing 8.10 - Box-Cox Transformation to normality
     
196. 
     
197. library(car)
     
198. summary(powerTransform(states$Murder))
     
199. 
     
200. # Box-Tidwell Transformations to Linearity
     
201. 
     
202. library(car)
     
203. boxTidwell(Murder ~ Population + Illiteracy, data = states)
     
204. 
     
205. # Listing 8.11 - Comparing nested models using the anova function
     
206. 
     
207. fit1 <- lm(Murder ~ Population + Illiteracy + Income +
     
208.     Frost, data = states)
     
209. fit2 <- lm(Murder ~ Population + Illiteracy, data = states)
     
210. anova(fit2, fit1)
     
211. 
     
212. # Listing 8.12 - Comparing models with the Akaike Information Criterion
     
213. 
     
214. fit1 <- lm(Murder ~ Population + Illiteracy + Income +
     
215.     Frost, data = states)
     
216. fit2 <- lm(Murder ~ Population + Illiteracy, data = states)
     
217. AIC(fit1, fit2)
     
218. 
     
219. # Listing 8.13 - Backward stepwise selection
     
220. 
     
221. library(MASS)
     
222. fit1 <- lm(Murder ~ Population + Illiteracy + Income +
     
223.     Frost, data = states)
     
224. stepAIC(fit, direction = "backward")
     
225. 
     
226. # Listing 8.14 - All subsets regression
     
227. # use the mouse to place the legend interactively
     
228. # in the second plot
     
229. 
     
230. library(leaps)
     
231. leaps <- regsubsets(Murder ~ Population + Illiteracy +
     
232.     Income + Frost, data = states, nbest = 4)
     
233. plot(leaps, scale = "adjr2")
     
234. 
     
235. library(car)
     
236. subsets(leaps, statistic = "cp",
     
237.     main = "Cp Plot for All Subsets Regression")
     
238. abline(1, 1, lty = 2, col = "red")
     
239. 
     
240. # Listing 8.15 - Function for k-fold cross-validated R-square
     
241. shrinkage <- function(fit, k = 10) {
     
242.     require(bootstrap)
     
243.     # define functions
     
244.     theta.fit <- function(x, y) {
     
245.         lsfit(x, y)
     
246.     }
     
247.     theta.predict <- function(fit, x) {
     
248.         cbind(1, x) %*% fit$coef
     
249.     }
     
250.    
     
251.     # matrix of predictors
     
252.     x <- fit$model[, 2:ncol(fit$model)]
     
253.     # vector of predicted values
     
254.     y <- fit$model[, 1]
     
255.    
     
256.     results <- crossval(x, y, theta.fit, theta.predict, ngroup = k)
     
257.     r2 <- cor(y, fit$fitted.values)^2
     
258.     r2cv <- cor(y, results$cv.fit)^2
     
259.     cat("Original R-square =", r2, "\n")
     
260.     cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
     
261.     cat("Change =", r2 - r2cv, "\n")
     
262. }
     
263. 
     
264. # using shrinkage()
     
265. 
     
266. fit <- lm(Murder ~ Population + Income + Illiteracy +
     
267.     Frost, data = states)
     
268. shrinkage(fit)
     
269. 
     
270. fit2 <- lm(Murder ~ Population + Illiteracy, data = states)
     
271. shrinkage(fit2)
     
272. 
     
273. #  Calculating standardized regression coefficients
     
274. zstates <- as.data.frame(scale(states))
     
275. zfit <- lm(Murder ~ Population + Income + Illiteracy +
     
276.     Frost, data = zstates)
     
277. coef(zfit)
     
278. 
     
279. # Listing 8.16 - relweights() function for calculating relative
     
280. # importance of predictors
     
281. 
     
282. ########################################################################
     
283. # The relweights function determines the relative importance of each   #
     
284. # independent variable to the dependent variable in an OLS regression. #
     
285. # The code is adapted from an SPSS program generously provided by      #
     
286. # Dr. Johnson.                                                         #
     
287. #                                                                      #
     
288. # See Johnson (2000, Multivariate Behavioral Research, 35, 1-19) for   #
     
289. # an explanation of how the relative weights are derived.              #
     
290. ########################################################################
     
291. relweights <- function(fit, ...) {
     
292.     R <- cor(fit$model)
     
293.     nvar <- ncol(R)
     
294.     rxx <- R[2:nvar, 2:nvar]
     
295.     rxy <- R[2:nvar, 1]
     
296.     svd <- eigen(rxx)
     
297.     evec <- svd$vectors
     
298.     ev <- svd$values
     
299.     delta <- diag(sqrt(ev))
     
300.    
     
301.     # correlations between original predictors and new orthogonal variables
     
302.     lambda <- evec %*% delta %*% t(evec)
     
303.     lambdasq <- lambda^2
     
304.    
     
305.     # regression coefficients of Y on orthogonal variables
     
306.     beta <- solve(lambda) %*% rxy
     
307.     rsquare <- colSums(beta^2)
     
308.     rawwgt <- lambdasq %*% beta^2
     
309.     import <- (rawwgt/rsquare) * 100
     
310.     lbls <- names(fit$model[2:nvar])
     
311.     rownames(import) <- lbls
     
312.     colnames(import) <- "Weights"
     
313.    
     
314.     # plot results
     
315.     barplot(t(import), names.arg = lbls, ylab = "% of R-Square",
     
316.         xlab = "Predictor Variables", main = "Relative Importance of Predictor Variables",
     
317.         sub = paste("R-Square = ", round(rsquare, digits = 3)),
     
318.         ...)
     
319.     return(import)
     
320. }
     
321. 
     
322. # using relweights()
     
323. 
     
324. fit <- lm(Murder ~ Population + Illiteracy + Income +
     
325.     Frost, data = states)
     
326. relweights(fit, col = "lightgrey")

复制代码

1. #------------------------------------------------------------------ #
   
2. # R in Action: Chapter 9                                            #
   
3. # requires that the multcomp, gplots, car, HH, effects,             #
   
4. #     rrcov, mvoutlier, MASS packages have been installed           #
   
5. # install.packages(c('multcomp', 'gplots', 'car', 'HH', 'effects',  #
   
6. #     'rrcov', 'mvoutlier', 'MASS'))                                #
   
7. #-------------------------------------------------------------------#
   
8. 
   
9. # pause for each graph
   
10. par(ask = TRUE)
    
11. 
    
12. # save original graphical parameters
    
13. opar <- par(no.readonly = TRUE)
    
14. 
    
15. # Listing 9.1 - One-way ANOVA
    
16. 
    
17. library(multcomp)
    
18. attach(cholesterol)
    
19. table(trt)
    
20. aggregate(response, by = list(trt), FUN = mean)
    
21. aggregate(response, by = list(trt), FUN = sd)
    
22. fit <- aov(response ~ trt)
    
23. summary(fit)
    
24. library(gplots)
    
25. plotmeans(response ~ trt, xlab = "Treatment", ylab = "Response",
    
26.     main = "Mean Plot\nwith 95% CI")
    
27. detach(cholesterol)
    
28. 
    
29. # Listing 9.2 - Tukey HSD pairwise group comparisons
    
30. 
    
31. TukeyHSD(fit)
    
32. par(las = 2)
    
33. par(mar = c(5, 8, 4, 2))
    
34. plot(TukeyHSD(fit))
    
35. par(opar)
    
36. 
    
37. # Multiple comparisons the multcomp package
    
38. 
    
39. library(multcomp)
    
40. par(mar = c(5, 4, 6, 2))
    
41. tuk <- glht(fit, linfct = mcp(trt = "Tukey"))
    
42. plot(cld(tuk, level = 0.05), col = "lightgrey")
    
43. par(opar)
    
44. 
    
45. # Assessing normality
    
46. 
    
47. library(car)
    
48. qqPlot(lm(response ~ trt, data = cholesterol), simulate = TRUE,
    
49.     main = "QQ Plot", labels = FALSE)
    
50. 
    
51. # Assessing homogeneity of variances
    
52. 
    
53. bartlett.test(response ~ trt, data = cholesterol)
    
54. 
    
55. # Assessing outliers
    
56. 
    
57. library(car)
    
58. outlierTest(fit)
    
59. 
    
60. # Listing 9.3 - One-way ANCOVA
    
61. 
    
62. data(litter, package = "multcomp")
    
63. attach(litter)
    
64. table(dose)
    
65. aggregate(weight, by = list(dose), FUN = mean)
    
66. fit <- aov(weight ~ gesttime + dose)
    
67. summary(fit)
    
68. 
    
69. # Obtaining adjusted means
    
70. 
    
71. library(effects)
    
72. effect("dose", fit)
    
73. 
    
74. # Listing 9.4 - Multiple comparisons using user supplied contrasts
    
75. 
    
76. library(multcomp)
    
77. contrast <- rbind(`no drug vs. drug` = c(3, -1, -1, -1))
    
78. summary(glht(fit, linfct = mcp(dose = contrast)))
    
79. 
    
80. # Listing 9.5 - Testing for Homegeneity of Regression Slopes
    
81. 
    
82. library(multcomp)
    
83. fit2 <- aov(weight ~ gesttime * dose)
    
84. summary(fit2)
    
85. 
    
86. # Visualizing a one-way ANCOVA
    
87. 
    
88. library(HH)
    
89. ancova(weight ~ gesttime + dose, data = litter)
    
90. 
    
91. # Listing 9.6 - Two way ANOVA
    
92. 
    
93. attach(ToothGrowth)
    
94. table(supp, dose)
    
95. aggregate(len, by = list(supp, dose), FUN = mean)
    
96. aggregate(len, by = list(supp, dose), FUN = sd)
    
97. fit <- aov(len ~ supp * dose)
    
98. summary(fit)
    
99. 
    
100. # plotting interactions
     
101. 
     
102. interaction.plot(dose, supp, len, type = "b", col = c("red",
     
103.     "blue"), pch = c(16, 18),
     
104.     main = "Interaction between Dose and Supplement Type")
     
105.    
     
106. library(gplots)
     
107. plotmeans(len ~ interaction(supp, dose, sep = " "),
     
108.     connect = list(c(1, 3, 5), c(2, 4, 6)),
     
109.     col = c("red", "darkgreen"),
     
110.     main = "Interaction Plot with 95% CIs",
     
111.     xlab = "Treatment and Dose Combination")
     
112.    
     
113. library(HH)
     
114. interaction2wt(len ~ supp * dose)
     
115. 
     
116. # Listing 9.7 - Repeated measures ANOVA with one between
     
117. # and within groups factor
     
118. 
     
119. w1b1 <- subset(CO2, Treatment == "chilled")
     
120. fit <- aov(uptake ~ (conc * Type) + Error(Plant/(conc)),
     
121.     w1b1)
     
122. summary(fit)
     
123. 
     
124. par(las = 2)
     
125. par(mar = c(10, 4, 4, 2))
     
126. with(w1b1, interaction.plot(conc, Type, uptake, type = "b",
     
127.     col = c("red", "blue"), pch = c(16, 18), main = "Interaction Plot for Plant Type and Concentration"))
     
128. boxplot(uptake ~ Type * conc, data = w1b1, col = (c("gold",
     
129.     "green")), main = "Chilled Quebec and Mississippi Plants",
     
130.     ylab = "Carbon dioxide uptake rate (umol/m^2 sec)")
     
131. par(opar)
     
132. 
     
133. # Listing 9.8 - One-way MANOVA
     
134. 
     
135. library(MASS)
     
136. attach(UScereal)
     
137. y <- cbind(calories, fat, sugars)
     
138. aggregate(y, by = list(shelf), FUN = mean)
     
139. cov(y)
     
140. fit <- manova(y ~ shelf)
     
141. summary(fit)
     
142. summary.aov(fit)
     
143. 
     
144. # Listing 9.9 - Assessing multivariate normality
     
145. # identify points interactively with the mouse
     
146. 
     
147. center <- colMeans(y)
     
148. n <- nrow(y)
     
149. p <- ncol(y)
     
150. cov <- cov(y)
     
151. d <- mahalanobis(y, center, cov)
     
152. coord <- qqplot(qchisq(ppoints(n), df = p), d, main = "QQ Plot Assessing Multivariate Normality",
     
153.     ylab = "Mahalanobis D2")
     
154. abline(a = 0, b = 1)
     
155. identify(coord$x, coord$y, labels = row.names(UScereal))
     
156. 
     
157. # multivariate outliers
     
158. library(mvoutlier)
     
159. outliers <- aq.plot(y)
     
160. outliers
     
161. 
     
162. # Listing 9.10 - Robust one-way MANOVA
     
163. # this may take a while...
     
164. 
     
165. library(rrcov)
     
166. Wilks.test(y, shelf, method = "mcd")
     
167. 
     
168. # Listing 9.11 - A regression approach to the ANOVA
     
169. # problem in section 9.3
     
170. 
     
171. fit.lm <- lm(response ~ trt, data=cholesterol)
     
172. summary(fit.lm)
     
173. 
     
174. contrasts(cholesterol$trt)

复制代码