[读者文库]Mastering Scientific Computing with R

Nicolle 发表于 2015-5-9 08:37

我已经上传过了

1. Analysis of variance
   
2. 
   
3. Analysis of variance (Anova) is used to fit data to a linear model when all explanatory variables are categorical. Each of these categorical explanatory variables is known as a factor, and each factor can have two or more levels. When a single factor is present with three or more levels, we use a one-way Anova to analyze the data. If we had a single factor with two levels, we would use a student's t-test to analyze the data. When there are two or more factors, we would use a two-way or three-way Anova. You can easily perform an Anova in R using the aov() function.
   
4. 
   
5. In the first example, we will look at the effect of the dose of drugA on the level of fatigue reported by the 20 patients as follows:
   
6. 
   
7. > patient.fatigue <- read.table(header=TRUE, text='
   
8.    patients fatigue drugA_dose
   
9. 1         1     low        0.2
   
10. 2         2     low        0.2
    
11. 3         3     med        0.2
    
12. 4         4     med        0.2
    
13. 5         5     med        0.2
    
14. 6         6     low        0.4
    
15. 7         7     low        0.4
    
16. 8         8     low        0.4
    
17. 9         9     med        0.4
    
18. 10       10     med        0.4
    
19. 11       11     med        0.8
    
20. 12       12    high        0.8
    
21. 13       13     med        0.8
    
22. 14       14     med        0.8
    
23. 15       15    high        0.8
    
24. 16       16    high        1.2
    
25. 17       17    high        1.2
    
26. 18       18    high        1.2
    
27. 19       19    high        1.2
    
28. 20       20     med        1.2 ')
    
29. 
    
30. >attach(patient.fatigue)
    
31. > aov(drugA_dose ~ fatigue)

复制代码

1. To use an alternative link function when fitting our data to a GLM model, we can change the link argument as follows:
   
2. 
   
3. glm(y, family=binomial(link=probit))
   
4. 
   
5. Now, let's go over a detailed example of how to fit your data to a GLM in R. In the first example, we will look at the effect of the dose of a compound on the death of 20 male and female mice. Let's take a look at this in the following lines of code:
   
6. 
   
7. > cmp1.ld <- read.table(header=TRUE, text='
   
8.    lethaldose sex numdead numalive
   
9. 1           0   M       1       19
   
10. 2           1   M       3       17
    
11. 3           2   M       9       11
    
12. 4           3   M      14        6
    
13. 5           4   M      17        3
    
14. 6           5   M      20        0
    
15. 7           0   F       0       20
    
16. 8           1   F       2       18
    
17. 9           2   F       2       18
    
18. 10          3   F       3       17
    
19. 11          4   F       4       16
    
20. 12          5   F       6       14
    
21. ')
    
22. > attach(cmp1.ld)
    
23. We can plot the data to take a look at the relationship between the dose of the compound and the proportions of deaths by mouse gender. Let's take a look at this in the following lines of code:
    
24. 
    
25. > proportion_dead <- numdead/20
    
26. > plot(proportion_dead ~ lethaldose, pch=as.character(sex))
    
27. The result is shown in the following plot:
    
28. 
    
29. Generalized linear models
    
30. Now we will fit our data to a GLM model. First, we combine the number of dead and alive mice into a counts matrix as follows:
    
31. 
    
32. > counts <- cbind(numdead, numalive)
    
33. > cmp1.ld.model <- glm( counts ~ sex * lethaldose, family=binomial)
    
34. We can summarize the results for the GLM model with the summary() function as follows:
    
35. 
    
36. > summary(cmp1.ld.model)

复制代码

1. Generalized additive models
   
2. 
   
3. Generalized additive models (GAMs) are a non-parametric extension of GLMs in which the linear predictor depends linearly on unknown smooth functions of some predictor variables. GAMs are typically used to let the data "speak for themselves" since you don't need to specify the functional form of the relationship, the response, and the continuous explanatory variables beforehand. To fit your data to a GAM, you will need to obtain the gam() function from the mgcv package. It is similar to the glm() function except that you add s() to each explanatory variable you wish to add smooths. For example, if you wish to describe the relationship between y and to smooth three continuous explanatory variables x, z, and w, you would enter model <- gam(y ~ s(x) + s(z) + s(w)). Now let's go through a detailed example in R.
   
4. 
   
5. In this example, we will explore the relationship between the total number of pregnancies and a variety of measurements taken from 300 mice.
   
6. 
   
7. Let's first simulate a dataset for 300 mice with pregnancies, glucose, pressure, insulin, and weight as explanatory variables. Since the sample() and rnorm() functions generate pseudorandom numbers, you will get different results when you run the following code in R. Let's simulate our explanatory variables with the following lines of code:
   
8. 
   
9. > pregnancies <- sample(0:25, 300,replace=T)
   
10. > glucose <- sample(65:200, 300,replace=T)
    
11. > pressure <-  sample(50:120, 300,replace=T)
    
12. > insulinD <- abs(rnorm(150, 450, 100))
    
13. > insulinN <- abs(rnorm(150, 65, 75))
    
14. > insulin <- c(insulinD, insulinN)
    
15. > weight <- sample(20:70, 300,replace=T)
    
16. Now let's use the gam() function part of the mgcv package to explore the relationship between diabetes and pregnancies, glucose, pressure, insulin, and weight as follows:
    
17. 
    
18. > library("mgcv")
    
19. > mouse.data.gam <- gam(pregnancies ~ s(glucose) + s(pressure) + s(insulin) + s(weight))
    
20. 
    
21. > summary(mouse.data.gam)

复制代码

1. Linear discriminant analysis
   
2. 
   
3. Linear discriminant analysis (LDA) is used to find the linear combinations of explanatory variables that give the best possible separation between the groups in our dataset. More specifically, LDA is used as a linear supervised classifier to form a classification model from the data provided. For example, we can use LDA to classify fish based on their length, weight, and speed underwater. Let's simulate a dataset based on the Ontario warm water fish species Bluegill, Bowfin, Carp, Goldeye, and Largemouth Bass.
   
4. 
   
5. > set.seed(459)
   
6. > Bluegill.length <- sample(seq(15, 22.5, by=0.5), 50, replace=T)
   
7. > Bluegill.weight <- sample(seq(0.2, 0.8, by=0.05), 50, replace=T)
   
8. > Bowfin.length <- sample(seq(46, 61, by=0.5), 50, replace=T)
   
9. > Bowfin.weight <- sample(seq(1.36, 3.2, by=0.5), 50, replace=T)
   
10. > Carp.length <- sample(seq(30, 75, by=1), 50, replace=T)
    
11. > Carp.weight <- sample(seq(0.2, 3.5, by=0.1), 50, replace=T)
    
12. > Goldeye.length <- sample(seq(25, 38, by=0.5), 50, replace=T)
    
13. > Goldeye.weight <- sample(seq(0.4, 0.54, by=0.01), 50, replace=T)
    
14. > Largemouth_Bass.length <- sample(seq(22, 55, by=0.5), 50, replace=T)
    
15. > Largemouth_Bass.weight <- sample(seq(0.68, 1.8, by=0.01), 50, replace=T)
    
16. 
    
17. > weight <-c(Bluegill.weight, Bowfin.weight, Carp.weight, Goldeye.weight, Largemouth_Bass.weight)
    
18. 
    
19. > length <-c(Bluegill.length, Bowfin.length, Carp.length, Goldeye.length, Largemouth_Bass.length)
    
20. 
    
21. > speed <- rnorm(50*5, 7.2, sd=1.8)
    
22. 
    
23. > fish <- c(rep("Bluegill", 50), rep("Bowfin", 50), rep("Carp", 50), rep("Goldeye", 50), rep("Largemouth_Bass", 50))
    
24. > fish.data <- data.frame(length, weight, speed, fish)
    
25. 
    
26. > str(fish.data)

复制代码

1. Principal component analysis
   
2. 
   
3. Principal component analysis (PCA) is another exploratory method you can use to separate your samples into groups. PCA converts a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components. PCA is widely used as a dimension reduction technique to help visualize data. PCA is different from LDA because it doesn't rely on class information to explore the relationship between the variable values and the sample group numbers. For example, let's perform a PCA to explore our simulated fish.data dataset. Before performing PCA, it is important to remember that the magnitude of the variables and any skews in the data will influence the resulting principal components. So, we need to scale and transform our data.
   
4. 
   
5. First, we recommend you to log transform the data (if necessary). Then, run PCA using the prcomp() function as follows:
   
6. 
   
7. > fish.data.mx <- as.matrix(fish.data[, 1:3])
   
8. > fish.data.log <- log(fish.data.mx)
   
9. > fish.log.pca <- prcomp(fish.data.log, scale=T, center=T)
   
10. > summary(fish.log.pca)
    
11. 
    
12. Instead of the summary, you might want to see the standard deviation of each principal component, and their rotation (or loadings), which are the coefficients of the linear combinations of the continuous variables. To get this information, we can use the print() function as follows:
    
13. > print(fish.log.pca)
    
14. 
    
15. > plot(fish.log.pca, ylim=c(0, 2)) #plot not shown
    
16. 
    
17. > library("ggplot2")
    
18. > loadings <- as.data.frame(fish.log.pca$rotation)
    
19. 
    
20. > # Add a column with the name of each variable to the loadings data frame.
    
21. > loadings$variables <- colnames(fish.data[,1:3])
    
22. 
    
23. > # Plot figure with qplot()
    
24. > q <- qplot(x = variables, y = PC1, data = loadings, geom = "bar", stat="identity")
    
25. 
    
26. > # Adjust axis label sizes
    
27. > q + theme(axis.title = element_text(face="bold", size=20), axis.text  = element_text(size=18))

复制代码

1. Clustering
   
2. 
   
3. An alternative approach to PCA is k-means (unsupervised) clustering, which partitions the data into k clusters in which each observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster. We can perform k-means clustering with the kmeans() function and plot the results with plot3d() as follows:
   
4. 
   
5. > set.seed(44)
   
6. > cl <- kmeans(fish.data[,1:3],5)
   
7. > fish.data$cluster <- as.factor(cl$cluster)
   
8. > plot3d(fish.log.pca$x[,1:3], col=fish.data$cluster, main="k-means clusters")
   
9. 
   
10. > with(fish.data, table(cluster, fish))
    
11.       
    
12. > di <- dist(fish.data[,1:3], method="euclidean")
    
13. > tree <- hclust(di, method="ward")
    
14. 
    
15. > fish.data$hcluster <- as.factor((cutree(tree, k=5)-2) %% 3 +1)
    
16. > plot(tree, xlab="", cex=0.2)
    
17. > rect.hclust(tree, k=5, border="red")
    
18. > with(fish.data, table(hcluster, fish))

复制代码

1. Structural Equation Modeling
   
2. 
   
3. We went over the matrix representation of SEM models in theory. Now, we will show how this can be practically done in R using the RAM matrix specifications mentioned earlier. We will use the political democracy dataset available in the lavaan package for this example. Notably, this dataset is ordered such that the first eight columns are the y variables and the last three are the x variables. Let's take a look at the following example:
   
4. 
   
5. library(lavaan)
   
6. data(PoliticalDemocracy)
   
7. We will create the covariance matrix of this dataset, as follows:
   
8. 
   
9. pd.cov <- cov(PoliticalDemocracy)
   
10. Now, we will create each of our matrices A, S, F, and I. In a real SEM example, we would iteratively create new A and S matrices in an attempt to find the best values to match the implied to the observed covariance matrix. Think of this as starting values for the matrices. In this case, we will choose starting values that pretty closely match a final solution, so that we do not have to iterate through dozens of times to find a good solution. Remember that we have 11 observed and three unobserved variables for a total of 14 variables.
    
11. 
    
12. First, let's take a look at the A matrix of paths, which is a matrix of 14*14:
    
13. 
    
14. mat.A <- matrix(
    
15.   c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,  0,
    
16.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,  0,
    
17.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,  0,
    
18.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,  0,
    
19.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  1,
    
20.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  1,
    
21.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  1,
    
22.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  1,
    
23.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,  0,
    
24.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,  0,
    
25.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,  0,
    
26.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  0,
    
27.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,1.5,0,  0,
    
28.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0.5,0.5,0
    
29.   ), nrow = 14, byrow = TRUE
    
30. )
    
31. Then, let's take a look at the S matrix of residual variances or covariances, which is also a matrix of 14*14:
    
32. 
    
33. mat.S <- matrix(
    
34.   c(2, 0, 0, 0,.5, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    
35.     0, 7, 0, 1, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0,
    
36.     0, 0, 5, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
    
37.     0, 1, 0, 3, 0, 0, 0,.5, 0, 0, 0, 0, 0, 0,
    
38.    .5, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    
39.     0, 2, 0, 0, 0, 5, 0, 1, 0, 0, 0, 0, 0, 0,
    
40.     0, 0, 1, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0,
    
41.     0, 0, 0,.5, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0,
    
42.     0, 0, 0, 0, 0, 0, 0, 0,.1, 0, 0, 0, 0, 0,
    
43.     0, 0, 0, 0, 0, 0, 0, 0, 0,.1, 0, 0, 0, 0,
    
44.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,.5, 0, 0, 0,
    
45.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,.5, 0, 0,
    
46.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,
    
47.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,.2
    
48.   ), nrow = 14, byrow = TRUE
    
49. )
    
50. Next, let's take a look at the filter matrix, which selects the observed variables is a matrix of 11*14:
    
51. 
    
52. mat.F <- matrix(
    
53.   c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    
54.     0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    
55.     0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    
56.     0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    
57.     0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    
58.     0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
    
59.     0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
    
60.     0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
    
61.     0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
    
62.     0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
    
63.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0
    
64.   ), nrow = 11, byrow = TRUE
    
65. )
    
66. Finally the 14 x 14 identity matrix:
    
67. mat.I <- diag(rep(1, 14), nrow = 14)
    
68. We can write a function that will perform all of the matrix operations of the McArdle McDonald equation (variables have a 0 in the name), as follows:
    
69. 
    
70. RAM.implied.covariance <- function(A.0, S.0, F.0, I.0) {
    
71.   implied.covariance <- F.0 %*% solve(I.0-A.0) %*% S.0 %*% t(solve(I.0-A.0)) %*% t(F.0)
    
72.   return(implied.covariance)
    
73. }
    
74. And, finally, we can estimate an implied covariance matrix based on our starting A and S matrices, which is shown in the following matrix rounded to 2 decimal places:
    
75. 
    
76. > round(RAM.implied.covariance(mat.A, mat.S, mat.F, mat.I), 2)
    
77.       [,1]  [,2]  [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
    
78. [1,] 7.12  5.12  5.12 5.12 3.44 2.94 2.94 2.94 0.75  1.50  1.50
    
79. [2,] 5.12 12.12  5.12 6.12 2.94 4.94 2.94 2.94 0.75  1.50  1.50
    
80. [3,] 5.12  5.12 10.12 5.12 2.94 2.94 3.94 2.94 0.75  1.50  1.50
    
81. [4,] 5.12  6.12  5.12 8.12 2.94 2.94 2.94 3.44 0.75  1.50  1.50
    
82. [5,] 3.44  2.94  2.94 2.94 3.98 1.98 1.98 1.98 0.63  1.25  1.25
    
83. [6,] 2.94  4.94  2.94 2.94 1.98 6.98 1.98 2.98 0.63  1.25  1.25
    
84. [7,] 2.94  2.94  3.94 2.94 1.98 1.98 4.98 1.98 0.63  1.25  1.25
    
85. [8,] 2.94  2.94  2.94 3.44 1.98 2.98 1.98 4.98 0.63  1.25  1.25
    
86. [9,] 0.75  0.75  0.75 0.75 0.63 0.63 0.63 0.63 0.60  1.00  1.00
    
87. [10,] 1.50  1.50  1.50 1.50 1.25 1.25 1.25 1.25 1.00  2.10  2.00
    
88. [11,] 1.50  1.50  1.50 1.50 1.25 1.25 1.25 1.25 1.00  2.00  2.50
    
89. Compare the preceding matrix to the observed covariance matrix as follows:
    
90. 
    
91. > round(pd.cov, 2)
    
92.      y1    y2    y3    y4   y5    y6    y7    y8   x1   x2   x3
    
93. y1 6.88  6.25  5.84  6.09 5.06  5.75  5.81  5.67 0.73 1.27 0.91
    
94. y2 6.25 15.58  5.84  9.51 5.60  9.39  7.54  7.76 0.62 1.49 1.17
    
95. y3 5.84  5.84 10.76  6.69 4.94  4.73  7.01  5.64 0.79 1.55 1.04
    
96. y4 6.09  9.51  6.69 11.22 5.70  7.44  7.49  8.01 1.15 2.24 1.84
    
97. y5 5.06  5.60  4.94  5.70 6.83  4.98  5.82  5.34 1.08 2.06 1.58
    
98. y6 5.75  9.39  4.73  7.44 4.98 11.38  6.75  8.25 0.85 1.81 1.57
    
99. y7 5.81  7.54  7.01  7.49 5.82  6.75 10.80  7.59 0.94 2.00 1.63
    
100. y8 5.67  7.76  5.64  8.01 5.34  8.25  7.59 10.53 1.10 2.23 1.69
     
101. x1 0.73  0.62  0.79  1.15 1.08  0.85  0.94  1.10 0.54 0.99 0.82
     
102. x2 1.27  1.49  1.55  2.24 2.06  1.81  2.00  2.23 0.99 2.28 1.81
     
103. x3 0.91  1.17  1.04  1.84 1.58  1.57  1.63  1.69 0.82 1.81 1.98
     
104. As we can see, it is pretty close. Can we improve it? Yes, but it will take some trial and error with new values in the A (and subsequently the S) matrix. This is what optimizers do in statistical software, and the R packages that we described here do just that. Notably, there are actually an infinite number of possible solutions for many SEM models, so some constraining of values is often needed. In this case, we would want to constrain a few particular path values. In this case, constraining one path going from each of the latent variables to a manifest variable will do the trick.

复制代码