74. その他の分析方法の紹介


詳しくは (その関数のパッケージを library() で呼び出してから) それぞれの関数のヘルプを御覧頂きたい.

因子分析

データ v2 はデータ v1 にノイズを加えたデータで,v4 と v3 ,v6 と v5 も同様な関係となっている.

 
 v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6); v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)
 v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6); v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)
 v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5); v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)
 m1 <- cbind(v1,v2,v3,v4,v5,v6)
 cor(m1)

          v1        v2        v3        v4        v5        v6
v1 1.0000000 0.9393083 0.5128866 0.4320310 0.4664948 0.4086076
v2 0.9393083 1.0000000 0.4124441 0.4084281 0.4363925 0.4326113
v3 0.5128866 0.4124441 1.0000000 0.8770750 0.5128866 0.4320310
v4 0.4320310 0.4084281 0.8770750 1.0000000 0.4320310 0.4323259
v5 0.4664948 0.4363925 0.5128866 0.4320310 1.0000000 0.9473451
v6 0.4086076 0.4326113 0.4320310 0.4323259 0.9473451 1.0000000

 factanal(m1, factors=3)              # varimax is the default

Call:
factanal(x = m1, factors = 3)
Uniquenesses:
   v1    v2    v3    v4    v5    v6 
0.005 0.101 0.005 0.224 0.084 0.005 
Loadings:                             # 固有ベクトルそのもの
   Factor1 Factor2 Factor3
v1 0.944   0.182   0.267  
v2 0.905   0.235   0.159  
v3 0.236   0.210   0.946  
v4 0.180   0.242   0.828  
v5 0.242   0.881   0.286  
v6 0.193   0.959   0.196  
 
               Factor1 Factor2 Factor3
SS loadings      1.893   1.886   1.797
Proportion Var   0.316   0.314   0.300
Cumulative Var   0.316   0.630   0.929
 
The degrees of freedom for the model is 0 and the fit was 0.4755 

よって,v2 と v1 ,v4 と v3 ,v6 と v5 が上手く選別できていることが分かる.

正準相関分析

正準相関分析を行なう場合は関数 cancor() を用いる.

 
 library(mva)
 data(LifeCycleSavings)
 pop <- LifeCycleSavings[, 2:3]
 oec <- LifeCycleSavings[, -(2:3)]
 cancor(pop, oec)                               # 一般的な教科書で知られている出力とは異なる

$cor                                            # 変数群の固有値 (正準相関係数)
[1] 0.8247966 0.3652762
$xcoef                                          # 変数 x の推定係数
              [,1]        [,2]
pop15 -0.009110856 -0.03622206
pop75  0.048647514 -0.26031158
$ycoef                                          # 変数 y の推定係数
             [,1]          [,2]          [,3]
sr   0.0084710221  3.337936e-02 -5.157130e-03
dpi  0.0001307398 -7.588232e-05  4.543705e-06
ddpi 0.0041706000 -1.226790e-02  5.188324e-02
$xcenter                                        # 変数 x を調節するのに使った値
  pop15   pop75 
35.0896  2.2930 
$ycenter                                        # 変数 y を調節するのに使った値
       sr       dpi      ddpi 
   9.6710 1106.7584    3.7576 

判別分析

線形を仮定した判別分析を行なう場合はパッケージ MASS の中の関数 lda() を用いる.関連のある関数として predict.lda() というものもある.

 
 library(MASS)
 data(iris3)
 Iris <- data.frame(rbind(iris3[,,1], iris3[,,2], iris3[,,3]), 
                    Sp = rep(c("s","c","v"), rep(50,3)))
 train <- sample(1:150, 75)
 table(Iris$Sp[train])       # your answer may differ

 c  s  v                     #  c  s  v
25 25 25                     # 22 23 30

 z <- lda(Sp ~ ., Iris, prior = c(1,1,1)/3, subset = train)
 predict(z, Iris[-train, ])$class

 [1] s s s s s s s s s s s s s s s s s s s s s s s s s c c c c c c c c c c c c c
[39] c c v c c c c c c c c c v v v v v v v v v v v v v v v v c v v v v v v v v
Levels: c s v

二次の判別分析を行なう場合はパッケージ MASS の中の関数 qda() を用いる.関連のある関数として predict.qda() というものがある.

 
 data(iris3)
 tr <- sample(1:50, 25)
 train <- rbind(iris3[tr,,1], iris3[tr,,2], iris3[tr,,3])
 test <- rbind(iris3[-tr,,1], iris3[-tr,,2], iris3[-tr,,3])
 cl <- factor(c(rep("s",25), rep("c",25), rep("v",25)))
 z <- qda(train, cl)
 predict(z,test)$class

 [1] s s s s s s s s s s s s s s s s s s s s s s s s s c c c c c c c c c c c c c
[39] c c c c c v c c c c c c v v v v v v v v v v v v v v v v v v c v v v v v v
Levels: c s v

主成分分析

主成分分析を行なう場合は関数 prcomp() を用いる.

 
 data(USArrests)
 prcomp(USArrests)

Standard deviations:
[1] 83.732400 14.212402  6.489426  2.482790
Rotation:
                 PC1         PC2         PC3         PC4
Murder   -0.04170432 -0.04482166  0.07989066  0.99492173
Assault  -0.99522128 -0.05876003 -0.06756974 -0.03893830
UrbanPop -0.04633575  0.97685748 -0.20054629  0.05816914
Rape     -0.07515550  0.20071807  0.97408059 -0.07232502

 prcomp(USArrests, scale = TRUE)

Standard deviations:
[1] 1.5748783 0.9948694 0.5971291 0.4164494
Rotation:
                PC1        PC2        PC3         PC4
Murder   -0.5358995 -0.4181809  0.3412327 -0.64922780
Assault  -0.5831836 -0.1879856  0.2681484  0.74340748
UrbanPop -0.2781909  0.8728062  0.3780158 -0.13387773
Rape     -0.5434321  0.1673186 -0.8177779 -0.08902432

 plot(prcomp(USArrests))
 summary(prcomp(USArrests, scale = TRUE))

Importance of components:
                        PC1   PC2    PC3    PC4
Standard deviation     1.57 0.995 0.5971 0.4164
Proportion of Variance 0.62 0.247 0.0891 0.0434
Cumulative Proportion  0.62 0.868 0.9566 1.0000

クラスター分析

クラスター分析を行なう場合は関数 hclust() を用いる.

 
 data(USArrests)
 hc <- hclust(dist(USArrests), "ave")
 plot(hc)
 plot(hc, hang = -1) 
 hc <- hclust(dist(USArrests)^2, "cen") 
 memb <- cutree(hc, k = 10); cent <- NULL 
 for(k in 1:10) cent <- rbind(cent, colMeans(USArrests[memb == k, , drop = FALSE])) 
 hc1 <- hclust(dist(cent)^2, method = "cen" , members = table(memb)) 
 opar <- par(mfrow = c(1, 2)) 
 plot(hc, labels = FALSE, hang = -1, main = "Original Tree") 
 plot(hc1, labels = FALSE, hang = -1, main = "Re-start from 10 clusters") 
 par(opar) 

 

 

k-means アルゴリズムによるクラスター分析を行う場合は以下のようにする.

 
 library(mva)
 x <- rbind(matrix(rnorm(100, sd = 0.3), ncol = 2),
            matrix(rnorm(100, mean = 1, sd = 0.3), ncol = 2))  # 2D example
 cl <- kmeans(x, 2, 20)
 plot(x, col = cl$cluster)
 points(cl$centers, col = 1:2, pch = 8)

回帰ツリー(CART ; Classification And Regression Trees)

回帰ツリーを用いて分析を行なう場合は関数 rpart() を用いる.

 
 library(rpart)
 result <- rpart(Species ~ ., data=iris)
 summary(result)  # 結果の要約
 par(xpd=T)       # ラベルが欠けないための処理
 plot(result)     # 回帰ツリーを描く
 text(result)     # ラベルを追記

時系列解析

R には時系列解析のための関数が多数用意されている.詳しくは『Rによる統計解析の基礎』 (中澤 港 著,ピアソン・エデュケージョン) ,『THE R BOOK』 岡田 昌史 他 著 (九天社) を参照されたい.

関数

機能

ts()

時系列オブジェクトを生成する.as.ts(),is.ts() も用意されている.

ts.union()

時系列オブジェクトを合併する.

ts.intersection()

時系列オブジェクトの共通部分を求める.

window()

時系列オブジェクトの一部分を切りとる.

diff()

要素の差分を取る.

diffinv()

要素の差分をどんどん足していく.

filter()

要素の線形フィルタリングを行う.

lag()

時間軸を過去にずらす.

acf()

時系列オブジェクトの自己共分散と自己相関係数を求める.

aggregate()

データを部分集合に分割し(引数 nfrequency で指定),それぞれの要約統計量を計算する.

ar()

時系列オブジェクトへの AR モデルの当てはめを行う.関連する関数として ar.burg(),ar.yw(),ar.mle(),ar.ols() がある.

arima()

時系列オブジェクトへの ARIMA モデルの当てはめを行う.関連する関数として arima.sim() がある.

ARMAacf()

ARMA モデルの自己相関係数を計算する.関連する関数として ARMAtoMA() がある.

Box.test()

時系列オブジェクトについて独立性の検定を行う.

ccf()

二つの一次元時系列間の相関係数共分散を求める.

decompose()

時系列を,移動平均により季節,傾向,不規則成分に分解する.似たような関数に stl() がある.

pacf()

偏自己相関係数(partial autocorrelations)を計算する.

PP.test()

時系列オブジェクトが単位根をもつかどうかの検定を行う.

spectrum()

スペクトル密度関数を推定する.関連する関数として spec.ar() や spec.pgram(),spec.taper() がある.

cpgram()

累積ピリオドグラムをプロットする.

lag.plot()

時系列のラグプロットを行う.

monthplot()

時系列の季節成分をプロットする.

ts.plot()

複数の時系列を一つの画面にプロットする.

関数 ts.plot() は任意個の時系列ベクトルを与えてそれらを一度にプロットさせることが出来る.また,色や線の種類,点の記号などを時系列ベクトル毎に変更することも出来る.

 
 library(ts)
 data(UKLungDeaths)
 ts.plot(ldeaths, mdeaths, fdeaths,
         gpars=list(xlab="year", ylab="deaths", lwd=c(1:3), col=c(2:4)))

plot() に時系列ベクトルを与える場合は,一つの時系列だけしか与えることが出来ないことに注意.

生存時間解析

急性骨隋白血病データ aml の内訳は

となっている.カプラン・マイヤー法によるメディアン生存関数は関数 survfit() で求めることが出来,その結果をプロットすることでカプラン・マイヤープロットが得られる.

 
 library(survival)

要求されたパッケージ splines をロード中です

 MYDATA <- aml
 MYDATA$time2 <- Surv(MYDATA$time, MYDATA$status)   # time2 の打ち切りには + がつく
 ( result <- survfit(time2 ~ x, data=MYDATA) )      # メディアン生存関数

Call: survfit(formula = time2 ~ x, data = MYDATA)
                 n events median 0.95LCL 0.95UCL
x=Maintained    11      7     31      18     Inf
x=Nonmaintained 12     11     23       8     Inf

 summary(result)                                    # 詳しい要約を表示(結果は省略)
 plot(result, col=c(1,2))                           # カプラン・マイヤープロット
 legend(100, 0.8, c("Maintained","Nonmaintained"), col=c(1:5), lwd=1, bg="gray")

関数 Surv(観察開始時間,観察終了時間,打ち切りフラグ) で「0:右側打ち切り,1:イベント,2:左側打ち切り,3:区間打ち切り」という打ち切りフラグを生成することが出来る.

ちなみに,ログランク検定は関数 survdiff() で,コックス回帰は関数 coxph() で実行できる.

 
 coxph(time2 ~ x + strata(DEMOG) , data=MYDATA)     # 層によってベースラインハザードが異なると仮定

ベースライン(時間0)で重症度が異なる場合,関数 strata() でベースライン値を層別因子として(共変量として)分析することが出来る.詳しくは『Rによる統計解析の基礎』 (中澤 港 著,ピアソン・エデュケージョン) ,『THE R BOOK』 岡田 昌史 他 著 (九天社) を参照されたい.

ノンパラメトリック回帰

詳しくは『みんなのためのノンパラメトリック回帰 (上・下)』 (竹澤 邦夫 著,吉岡書店) を参照されたい.まず,ナダラヤ・ワトソン推定量を用いた平滑化は以下のようにすれば良い.

 
 library(modreg)
 data(cars)
 with(cars, {
        plot(speed, dist)
        lines(ksmooth(speed, dist, "normal", bandwidth=2), col=2)
        lines(ksmooth(speed, dist, "normal", bandwidth=5), col=3)
      })

平滑化スプラインによる平滑化を行う場合は以下のようにする.

 
 data(cars)
 attach(cars)
 plot(speed, dist, main = "data(cars)  &  smoothing splines")
 ( cars.spl <- smooth.spline(speed, dist) )

Call:
smooth.spline(x = speed, y = dist)
 
Smoothing Parameter  spar= 0.7801305  lambda= 0.1112206 (11 iterations) 
Equivalent Degrees of Freedom (Df): 2.635278 
Penalized Criterion: 4337.638 
GCV: 244.1044 

 lines(cars.spl, col = "blue")
 lines(smooth.spline(speed, dist, df=10), lty=2, col = "red")
 legend(5,120,c(paste("default [C.V.] => df =",round(cars.spl$df,1)),
        "s( * , df = 10)"), col = c("blue","red"), lty = 1:2, bg='bisque')