R_Lattice_Graphics

4506 days ago by takepwave

# Rのユーティリティ関数を読み込む attach(DATA+'RUtil.py') 
       

1.1章

英国における1997年のAレベルテストの化学の試験結果を使ってLatticeのデモをする。

最初に頻度表の表示
#r('install.packages("mlmRev")') r('library(mlmRev)') 
        
 [1] "mlmRev"    "lme4"      "Matrix"    "lattice"   "stats"    
"graphics"  "grDevices" "utils"    
 [9] "datasets"  "methods"   "base"     
 [1] "mlmRev"    "lme4"      "Matrix"    "lattice"   "stats"     "graphics"  "grDevices" "utils"    
 [9] "datasets"  "methods"   "base"
r('data(Chem97, package="mlmRev")') r('xtabs( ~ score, data = Chem97)') 
        
score
   0    2    4    6    8   10 
3688 3627 4619 5739 6668 6681 
score
   0    2    4    6    8   10 
3688 3627 4619 5739 6668 6681

1.1.1  グループ別ヒストグラム

最初にLatticeのパッケージを読み込む
#r('install.packages("lattice")') r('library(lattice)') 
        
 [1] "mlmRev"    "lme4"      "Matrix"    "lattice"   "stats"    
"graphics"  "grDevices" "utils"    
 [9] "datasets"  "methods"   "base"     
 [1] "mlmRev"    "lme4"      "Matrix"    "lattice"   "stats"     "graphics"  "grDevices" "utils"    
 [9] "datasets"  "methods"   "base"

以下の例では、scoreが平均gcsescoreに左右されるかを知るために、socreの値毎にgcsescoreのヒストグラムを表示する。
fig1_1 = preGraph("fig1_1.pdf") # Sageではそのままhistogramを使って表示させるとうまく表示できないので、 # 一度変数にセットして、それをplot関数でプロットする方法を取る r('tp1<-histogram(~ gcsescore | factor(score), data = Chem97)') r('plot(tp1)') postGraph(fig1_1) 
       

~ gcsescoreは、gcsescoreで回帰することを意味し、

| factore(score)は、条件変数で、変数を因子に変換と説明があるがこの部分がよくわからない。

2.1.3の説明によると条件変数は、ほとんどの場合カテゴリカル変数で、因子(factor)であると説明されている。
graph = preGraph("fig1_1a.pdf") r('tp1<-histogram(~ gcsescore , data = Chem97)') r('plot(tp1)') postGraph(graph) 
       

1.1.3 カーネル密度プロット

ヒストグラムの代わりにカーネル密度推定をプロットするのが、densityplot関数です。
graph = preGraph("fig1_2.pdf") r('tp1<-densityplot(~ gcsescore | factor(score), data = Chem97, plot.points = FALSE, ref = TRUE)') r('plot(tp1)') postGraph(graph) 
       

1.2 重ね合わせ

groupsにグループ化する変数をセットして、カーネル密度推定を重ね合わせることができます。
graph = preGraph("fig1_3.pdf") r('tp2<-densityplot(~ gcsescore , data = Chem97, groups = score, plot.points = FALSE, ref = TRUE, auto.key = list(columns = 3))') r('plot(tp2)') postGraph(graph) 
       

1.4 複数の図を１枚に出力する

splitオプションで２つのプロットを１つのグラフに出力する例を示します。
graph = preGraph("fig1_4.pdf") r('plot(tp1, split = c(1, 1, 1, 2))') r('plot(tp2, split = c(1, 2, 1, 2), newpage = FALSE)') postGraph(graph) 
       
#r('install.packages("MEMSS")') r('library(MEMSS)') 
        
 [1] "MEMSS"     "mlmRev"    "lme4"      "Matrix"    "lattice"   "stats"
"graphics"  "grDevices"
 [9] "utils"     "datasets"  "methods"   "base"     
 [1] "MEMSS"     "mlmRev"    "lme4"      "Matrix"    "lattice"   "stats"     "graphics"  "grDevices"
 [9] "utils"     "datasets"  "methods"   "base"
r('data(Oats, package= "MEMSS")') 
        
[1] "Oats"
[1] "Oats"
graph = preGraph("fig2_1.pdf") r('tp1.oats <- xyplot(yield ~ nitro | Variety + Block, data = Oats, type = "o")') r('plot(tp1.oats)') postGraph(graph) 
       
r('dim(tp1.oats)') 
        
[1] 3 6
[1] 3 6
r('dimnames(tp1.oats)') 
        
$Variety
[1] "Golden Rain" "Marvellous"  "Victory"    

$Block
[1] "I"   "II"  "III" "IV"  "V"   "VI" 
$Variety
[1] "Golden Rain" "Marvellous"  "Victory"    

$Block
[1] "I"   "II"  "III" "IV"  "V"   "VI"
r('xtabs(~ Variety + BVlock, data = Oats)') 
        
[1] "Oats"
[1] "Oats"
r('summary(tp1.oats)') 
        
Call:
xyplot(yield ~ nitro | Variety + Block, data = Oats, type = "o")

Number of observations:
             Block
Variety       I II III IV V VI
  Golden Rain 4  4   4  4 4  4
  Marvellous  4  4   4  4 4  4
  Victory     4  4   4  4 4  4
Call:
xyplot(yield ~ nitro | Variety + Block, data = Oats, type = "o")

Number of observations:
             Block
Variety       I II III IV V VI
  Golden Rain 4  4   4  4 4  4
  Marvellous  4  4   4  4 4  4
  Victory     4  4   4  4 4  4
# パネルは配列として扱える # 下から１行目の図 graph = preGraph("fig2_2.pdf") r('plot(tp1.oats[,1])') postGraph(graph) # 下から４行１列目の図 graph = preGraph("fig2_2a.pdf") r('plot(tp1.oats[1,4])') postGraph(graph) 
       
graph = preGraph("fig2_3.pdf") r('plot(update(tp1.oats, aspect="xy"))') postGraph(graph) 
       

layout引数の最初に0を指定すると、２番目の要素がページあたりの総パネル数と解釈される。
graph = preGraph("fig2_4.pdf") r('plot(update(tp1.oats, aspect="xy", layout = c(0, 18)))') postGraph(graph) 
       

2.3 グループ化された画像

groupsを指定すれば、自動でグループ化された画像が作成される。
graph = preGraph("fig2_6.pdf") r('tp1<- dotplot(variety ~ yield | site, barley, layout = c(1, 6), aspect = c(0.7), groups = year, auto.key = list(space = "right"))') r('plot(tp1)') postGraph(graph) 
       

2.4.1 凡例の補足
graph = preGraph("fig2_7.pdf") r('key.variety <- list(space = "right", text = list(levels(Oats$Variety)), points = list(ch = 1:3, col = "black"))') r('tp1<- xyplot(yield ~ nitro | Block, data = Oats, aspect = "xy", type = "o", groups = Variety, key = key.variety, lty = 1, pch = 1:3, col.line = "darkgrey", col.symbol = "black", xlab = "Nitrogen concentration (cwt/acre)", ylab = "Yield (bushels/acre)", main = "Yield of three varieties of oats", sub = "A 3 x 4 split-plot experiment with 6 blocks")') r('plot(tp1)') postGraph(graph) 
       

2.5.1 縮尺と軸

タイタニックのデータで等級、性別、年齢のカテゴリ毎にクロス集計の結果をプロット
graph = preGraph("fig2_8.pdf") r('tp1<- barchart(Class ~ Freq | Sex + Age, data = as.data.frame(Titanic), groups = Survived, stack = TRUE, layout = c(4, 1), auto.key = list(title = "Surviced", columns = 2))') r('plot(tp1)') postGraph(graph) 
       

スケールをパネル毎に変更するには、scalesを指定する。
graph = preGraph("fig2_9.pdf") r('tp1<- barchart(Class ~ Freq | Sex + Age, data = as.data.frame(Titanic), groups = Survived, stack = TRUE, layout = c(4, 1), auto.key = list(title = "Surviced", columns = 2), scales = list( x = "free"))') r('plot(tp1)') postGraph(graph) 
       

panelは関数として与えられるので、これに独自の関数を定義すると表示形式を変更することができる。

図2.10では背景に基準線を追加している。
graph = preGraph("fig2_10.pdf") r('bc.taitanic<- barchart(Class ~ Freq | Sex + Age, data = as.data.frame(Titanic), groups = Survived, stack = TRUE, layout = c(4, 1), auto.key = list(title = "Surviced", columns = 2), scales = list( x = "free"))') r('plot(update(bc.taitanic, panel = function(...) { panel.grid(h = 0, v=-1); panel.barchart(...) }))') postGraph(graph) 
       

3.1 密度プロット

間欠泉のfaithfulのデータで噴出継続時間の分布を表示する。
graph = preGraph("fig3_1.pdf") r('tp1<- densityplot(~ eruptions, data = faithful)') r('plot(tp1)') postGraph(graph) 
       

kernelとして矩形カーネル（rect）、バンド幅bwを指定する。

以下の図は、平均移動ヒスとブラムとも呼ばれる。
graph = preGraph("fig3_2.pdf") r('tp1<- densityplot(~ eruptions, data = faithful, kernel = "rect", bw = 0.2, plot.points = "rug", n = 200)') r('plot(tp1)') postGraph(graph) 
       
#r('install.packages("latticeExtra")') r('library("latticeExtra")') 
        
 [1] "latticeExtra" "RColorBrewer" "MASS"         "mlmRev"       "lme4" 
"Matrix"      
 [7] "lattice"      "stats"        "graphics"     "grDevices"    "utils"
"datasets"    
[13] "methods"      "base"        
 [1] "latticeExtra" "RColorBrewer" "MASS"         "mlmRev"       "lme4"         "Matrix"      
 [7] "lattice"      "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
[13] "methods"      "base"
r('data(gvhd10)') 
        
[1] "gvhd10"
[1] "gvhd10"
# プロット graph = preGraph("fig3_3.pdf") r('tp1<- densityplot(~ log(FSC.H) | Days, data = gvhd10, plot.points = FALSE, ref = TRUE, layout = c(2, 4))') r('plot(tp1)') postGraph(graph) 
       

3.3 ヒストグラム

Fig3.3をヒストグラムで表示したもの、データ数が多いのでピン区間の数nintを50にセットしている。
graph = preGraph("fig3_4.pdf") r('tp1<- histogram(~ log(FSC.H) | Days, data = gvhd10, xlab = "log Forward Scatter", type = "density", nint = 50, layout = c(2, 4))') r('plot(tp1)') postGraph(graph) 
       

3.4 正規Q-Qプロット

分位点-分位点（Q-Q）プロット（よくわからないが、大局的な特徴を判断するツールの一つで、観測データの分位点をプロットしたもので、直線からずれていたら、理論的な分布と適合していないと判断できるらしい）
graph = preGraph("fig3_5.pdf") r('tp1<- qqmath(~gcsescore | factor(score), data = Chem97, f.value = ppoints(100))') r('plot(tp1)') postGraph(graph) 
       
graph = preGraph("fig3_6.pdf") r('tp1<- qqmath(~gcsescore | gender, data = Chem97, groups = score, aspect = "xy", f.value = ppoints(100), auto.key = list(space = "right"), xlab = "Standar NOrmal Quantiles", ylab = "Average GCSE Score")') r('plot(tp1)') postGraph(graph) 
       

 

3.4.1 正規性とボックス・コックス

 

3.4.1 正規性とボックス・コックス

ボックス・コックス変換は、尺度と位置の移動を伴うベキ乗変換で以下の式で与えられる。
$f_{\lambda} (x) $

 

 
r('library("MASS")') r('Chem97.pos <- subset(Chem97, gcsescore > 0)') graph = preGraph("fig3_7a.pdf") r('with(Chem97.pos, boxcox(gcsescore ~ score * gender, lambda = seq(0, 4, 1/10)))') postGraph(graph) 
       

3.5 経験累積分布関数

ECDFは、累積分布関数Fの最尤推定値であり、ecdfplot関数で作成できる。
graph = preGraph("fig3_8.pdf") r('tp1 <- ecdfplot(~ gcsescore | factor(score), data = Chem97, groups = gender, auto.key = list(columns = 2), subset = gcsescore > 0, xlab = "Average GSSE Score")') r('plot(tp1)') postGraph(graph) 
       

3.6 2標本Q-Qプロット

二つの観測データの集合を比較するためにもQ-Qプロットが有効である。

女性の方が男性よりも平均が高いことが分かる。（分散が小さいことはどうすればこの図から分かるのだろうか？）
graph = preGraph("fig3_10.pdf") r('tp1 <- qq(gender ~ gcsescore | factor(score), data = Chem97, f.value = ppoints(100), aspect = 1)') r('plot(tp1)') postGraph(graph) 
       

3.7 箱ひげ図

性別で条件付けした最終スコアごとの平均GCSEスコアの箱ひげ図、女性の方が平均点が高く、分散も小さいことがこの図では分かる。
graph = preGraph("fig3_11.pdf") r('tp1 <- bwplot(factor(score) ~ gcsescore | gender, data = Chem97, xlab = "Average GCSE Score")') r('plot(tp1)') postGraph(graph) 
       
graph = preGraph("fig3_13.pdf") r('tp1 <- bwplot(Days ~ log(FSC.H), data = gvhd10, xlab = "log(Forward Scatter)", ylab = "Days Past transplant")') r('plot(tp1)') postGraph(graph) 
       

3.7.1 バイオリン・プロット

分布が単峰ではない場合に、箱ひげ図よりもバイオリンプロットの方が、分布が分かりやすい。
graph = preGraph("fig3_14.pdf") r('tp1 <- bwplot(Days ~ log(FSC.H), data = gvhd10, panel = panel.violin, xlab = "log(Forward Scatter)", ylab = "Days Past transplant")') r('plot(tp1)') postGraph(graph) 
       
graph = preGraph("fig3_15.pdf") r('tp1 <- stripplot(factor(mag) ~ depth, quakes)') r('plot(tp1)') postGraph(graph) 
       
# 軸を逆にプロット, alphaで半透明にセット, ジッタ指定で重なって部分も分かるようにしている。 graph = preGraph("fig3_16.pdf") r('tp1 <- stripplot(depth ~ factor(mag) , quakes, jitter.data = TRUE, alpha = 0.6, xlab = "Magnitude (Richter)", ylab = "Depth (km)")') r('plot(tp1)') postGraph(graph)