irisの解析

有名なヤマメ（だっけ？）の解析。

データ構造の確認

> head(iris)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

> str(iris)

'data.frame':   150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

基本統計量

> summary(iris)

  Sepal.Length   Sepal.Width    Petal.Length   Petal.Width 
 Min.   :4.30   Min.   :2.00   Min.   :1.00   Min.   :0.1  
 1st Qu.:5.10   1st Qu.:2.80   1st Qu.:1.60   1st Qu.:0.3  
 Median :5.80   Median :3.00   Median :4.35   Median :1.3  
 Mean   :5.84   Mean   :3.06   Mean   :3.76   Mean   :1.2  
 3rd Qu.:6.40   3rd Qu.:3.30   3rd Qu.:5.10   3rd Qu.:1.8  
 Max.   :7.90   Max.   :4.40   Max.   :6.90   Max.   :2.5  
       Species  
 setosa    :50  
 versicolor:50  
 virginica :50  

散布図

> plot(iris[, 1:4], col = iris$Species)

ggplot2でプロット

縦長(long-format)データにします。

> library(reshape2)
> iris2 <- melt(iris)

Using Species as id variables

> head(iris2)

  Species     variable value
1  setosa Sepal.Length   5.1
2  setosa Sepal.Length   4.9
3  setosa Sepal.Length   4.7
4  setosa Sepal.Length   4.6
5  setosa Sepal.Length   5.0
6  setosa Sepal.Length   5.4

ボックスプロット

> library(ggplot2)

Find out what's changed in ggplot2 with
news(Version == "0.9.1", package = "ggplot2")

> ggplot(iris2, aes(variable, value, colour = Species)) + 
+     geom_boxplot()

> ggplot(iris2, aes(variable, value, fill = Species)) + 
+     stat_summary(fun.y = mean, geom = "bar", position = "dodge")

t検定で対比較

> library(plyr)
> dlply(iris2, .(variable), function(x) pairwise.t.test(x$value, 
+     x$Species))

$Sepal.Length

    Pairwise comparisons using t tests with pooled SD 

data:  x$value and x$Species 

           setosa  versicolor
versicolor 1.8e-15 -         
virginica  < 2e-16 2.8e-09   

P value adjustment method: holm 

$Sepal.Width

    Pairwise comparisons using t tests with pooled SD 

data:  x$value and x$Species 

           setosa  versicolor
versicolor < 2e-16 -         
virginica  9.1e-10 0.0031    

P value adjustment method: holm 

$Petal.Length

    Pairwise comparisons using t tests with pooled SD 

data:  x$value and x$Species 

           setosa versicolor
versicolor <2e-16 -         
virginica  <2e-16 <2e-16    

P value adjustment method: holm 

$Petal.Width

    Pairwise comparisons using t tests with pooled SD 

data:  x$value and x$Species 

           setosa versicolor
versicolor <2e-16 -         
virginica  <2e-16 <2e-16    

P value adjustment method: holm 

attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")
      variable
1 Sepal.Length
2  Sepal.Width
3 Petal.Length
4  Petal.Width