Murakami's Memorandum

Home Page > My Memo               災害科学ホワイトボード   災害科学コース   理学部

Rメモ

(工事中)


R Tips

高知大学ではRの方がメジャーになりつつあるので,そろそろ準備を。。。いつ終わるのだか?
amazonプラグインは存在しません。
数学の野間口先生が訳された本です。
正誤表など

 Rの準備

  • R version 2.7.2 (2008-08-25)
  • 追加パッケージ: Rcmdr, outliers


  • R version 2.9.2 (2009-08-24)
  • 追加パッケージ:Rcmdr ver.1.5-1
    • ショートカットのプロパティ指定(右クリック)で
C:\Programs\R\R-2.9.2\bin\Rgui.exe --internet2 --sdi

のようにオプションを追加する。1番目はproxy対策で,2番目はRcmdrの新バージョン用の指定。

 Rcmdrの起動

  • コンソールから起動
>library(Rcmdr)
  • メニューから起動
[パッケージ]-[パッケージの読み込み]-[Rcmdr]
  • 自動起動

Rフォルダー中にあるetcフォルダーの中にあるRprofile.siteファイルをエディタなどで開き,最終行にlibrary(Rcmdr)を追加し保存する。

    • ライブラリーのアンロード

Rcmdrを自動起動にしておくと,Rcmdrの窓で「コマンダーの終了」を選び終了してもそのままの状態ではパッケージの更新時に「警告: packages 'car, Rcmdr' are in use and will not be installed」と言われて更新できないことがある。その場合には以下のようにしてパッケージのアンロードをする。
>detach("package:Rcmdr")


 Rcmdrで基礎統計

Rcmdrからデータを入力

  1. [データ]-[新しいデータセット]を選択する。
  2. データセット名を設定する(例えば,Dataset)。
  3. データ入力窓が開くので,var1に変数名を設定する(例えば,xなど)。numericを設定する。
  4. 後は数値データを入力してゆく。
      • PCの状況によるようだが,数値を入力してenterを押すとハングする場合がある。その場合にはいったんRを終了し,再起動する。同じくデータ入力窓を開き,数値を入れたあと下矢印キーを押して次の入力セルに移動するという方法をとる。
  5. 入力が終了したら,データ入力窓を閉じる。

コンソール窓でデータを確認した例。

> Dataset$x
 [1] 3470 2550 2920 2530 3280 2840 2520 3350 3610 3430

データ分布をグラフで見る(Rcmdr)

これらのグラフが簡単に描けるだけでもR+Rcmdrを使う価値はある。

  • 一変数
    • [グラフ]-[ヒストグラム]
    • [グラフ]-[幹歯表示]
    • [グラフ]-[箱ひげ図]
    • [グラフ]-[QQプロット]

  • 二変数
    • [グラフ]-[散布図]

(関連項目)

      • [統計量]-[要約]-[相関行列]    :相関係数
      • [統計量]-[]-[要約]-[相関の検定] :相関係数の検定

五数要約(Rcmdr)

  • [統計量]-[要約]-[アクティブデータセット]

> summary(Dataset)
       x       
 Min.   :2520  ;最小値
 1st Qu.:2622  ;第一四分位
 Median :3100 ;中央値  
 Mean   :3050  ;平均値
 3rd Qu.:3410  ;第三四分位
 Max.   :3610  ;最大値

  • [統計量]-[要約]-[数値による要約]

> numSummary(Dataset[,"x"], statistics=c("mean", "sd", "quantiles"), quantiles=c( 0,.25,.5,.75,1 ))
 mean       sd   0%    25%  50%  75% 100%  n
 3050 426.8229 2520 2622.5 3100 3410 3610 10

標準偏差(sd)が出力されている

平均値の信頼区間(Rcmdr)

とりあえず

  • [統計量]-[平均]-[1標本t検定]

> t.test(Dataset$x, alternative='two.sided', mu=0.0, conf.level=.95)

	One Sample t-test

data:  Dataset$x 
t = 22.5971, df = 9, p-value = 3.086e-09
alternative hypothesis: true mean is not equal to 0 
95 percent confidence interval:        ;95%信頼区間
 2744.669 3355.331           ;
sample estimates:
mean of x 
     3050 

マニュアルで平均・標準偏差

> x<-c(1,2,3,4,5)
> mean(x)              ;平均値
[1] 3
> sd(x)        ;標準偏差
[1] 1.581139

t分布のパーセント点をマニュアルで計算

>qt(1-a/2,df)      ;a:有意水準    ,df:自由度
  • 平均値の95%信頼区間を求めた例
> z = c(23.5,16.6,25.4,19.1,19.3,22.4,20.9,24.9)
> ave = mean(z)       ;平均値
> u = sd(z)                ;標準偏差
> t = qt(1-0.05/2,8-1)         ;t分布のパーセント点
> t
[1] 2.364624
> upper=ave+t*u/sqrt(8)     ;95%信頼区間上限
> lower=ave-t*u/sqrt(8)       ;95%信頼区間下限
> upper
[1] 24.09023
> lower
[1] 18.93477





 クリップボードのデータをR(Windows R+Rcmdr)に渡す


 Excelシート上のデータをクリップボード経由でRに渡したい時のTips.

  1. Excelシート上のデータを選択しコピー(Ctl+C)
  2. R+Rcmdrを立ち上げる
  3. 「データ」->「データのインポート」->「テキストファイルまたはクリップボードから...」
  4. 窓が開くので,データセット名を設定し,「クリップボードからデータを読み込む」ボタンを設定,「フィールド区切り記号」で空白を設定。

 Shapiro-Wilk normality test (Rcmdr)

  • [統計量]-[要約]-[シャピロウ−ウィルクの正規性の検定]

マニュアルでやるには

> x<-c(3470,2550,2920,2530,3280,2840,2520,3350,3610,5430)
> shapiro.test(x)

        Shapiro-Wilk normality test

data:  x 
W = 0.7819, p-value = 0.008743

標本 x1, ..., xnが正規母集団からサンプリングされたものであるという帰無仮説を検定する検定方法。
帰無仮説:正規分布である
対立仮説:正規分布でない
p値が0.05より小さいので(帰無仮説を棄却し),「有意水準5%で正規分布でない」可能性がある。

 Grubbs Test

> x<-c(3470,2550,2920,2530,3280,2840,2520,3350,3610,5430)
> library(outliers)
> grubbs.test(x)

        Grubbs test for one outlier

data:  x 
G = 2.5155, U = 0.2188, p-value = 0.003454
alternative hypothesis: highest value 5430 is an outlier 

標本 x1,..., xn において一つの値が外れ値(異常値)かどうかを検定する方法。
帰無仮説:5430は外れ値でない(標本は同じ正規母集団からの標本である)
対立仮説:5430は外れ値である(標本は同じ正規母集団からの標本ではない)
p値が0.05より小さいので(帰無仮説を棄却し),「有意水準5%で外れ値である」とみなせる。

 回帰直線(Rcmdr)

  1. 散布図を描く [グラフ]-[散布図]
  2. 回帰直線をあてはめる [統計量]-[モデルへの適合]-[線形回帰]
  3. モデル残差のQQプロット [モデル]-[グラフ]-[残差QQプロット]

入力データ

> Dataset
     x      y
1  0.0  0.012
2  0.1  0.121
3  0.2 -0.097
4  0.3 -0.061
5  0.4 -0.080
6  0.5  0.037
7  0.6  0.196
8  0.7  0.077
9  0.8  0.343
10 0.9  0.448
11 1.0  0.434
> 

Rcmdr 線形回帰の実行例(Rcmdr の下の窓に表示された結果)

> RegModel.1 <- lm(y~x, data=Dataset)
> summary(RegModel.1)
Call:
lm(formula = y ~ x, data = Dataset)
Residuals:
     Min       1Q   Median       3Q      Max 
-0.16191 -0.09391  0.01791  0.09559  0.18336 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.11045    0.07195  -1.535  0.15913   ;切片
x            0.48091    0.12162   3.954  0.00333 **  ;傾き
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.1276 on 9 degrees of freedom
Multiple R-squared: 0.6347,	Adjusted R-squared: 0.5941
F-statistic: 15.64 on 1 and 9 DF,  p-value: 0.003334 
   ※Excelで回帰直線を実行した時の決定係数R^2は上記のMultiple R-squared

マニュアルで回帰直線をあてはめるには,この画面にある lm 関数を使用する。

 重み付き回帰直線

lm 関数を使うと重み付きの回帰直線を求めることができる。

ここでは,sがyの標準偏差を表すものとする。weight=1/(s*s)で,正規方程式の重みとなる。

> Dataset
   x   y   s weight
1 -2 2.1 0.2     25
2  0 2.4 0.2     25
3  2 2.5 0.5      4
4  4 3.5 0.1    100

重み付きの回帰直線を求める例

> RegModel.1 <- lm(y~x,data=Dataset,weight=Dataset$weight)    ;重み付き回帰直線を計算
> summary(RegModel.1)                        ;結果を表示

Call:
lm(formula = y ~ x, data = Dataset, weights = Dataset$weight)

Residuals:
      1       2       3       4 
 0.4585 -0.4923 -0.9772  0.2123 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  2.49846    0.09569  26.109  0.00146 **
x            0.24508    0.02867   8.549  0.01341 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.8522 on 2 degrees of freedom
Multiple R-squared: 0.9734,	Adjusted R-squared:  0.96 
F-statistic: 73.08 on 1 and 2 DF,  p-value: 0.01341 

重みを考慮しない場合の例

> RegModel.2 <- lm(y~x, data=Dataset)
> summary(RegModel.2)

Call:
lm(formula = y ~ x, data = Dataset)

Residuals:
    1     2     3     4 
 0.12 -0.01 -0.34  0.23 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  2.41000    0.16568  14.546  0.00469 **
x            0.21500    0.06764   3.179  0.08635 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.3025 on 2 degrees of freedom
Multiple R-squared: 0.8348,	Adjusted R-squared: 0.7521 
F-statistic:  10.1 on 1 and 2 DF,  p-value: 0.08635 

なれるまでは難しそうですが作図例

> ucl <- Dataset$y+Dataset$s                                      #y+s
> lcl <- Dataset$y-Dataset$s                                      #y-s
> plot(Dataset$x,Dataset$y, ylim=range(c(lcl,ucl)))               #散布図
> arrows(Dataset$x,ucl,Dataset$x,lcl,length=.05,angle=90,code=3)  #エラーバー
> RegModel.1 <- lm(y~x,data=Dataset)                              #回帰式を計算
> RegModel.2 <- lm(y~x,data=Dataset,weight=Dataset$weight)        #重み付き回帰式を計算
> abline(RegModel.1)                                              #回帰式をプロット
> abline(RegModel.2,col=4)                                        #重み付き回帰式をプロット


 多項式のあてはめ(Rcmdr)

  1. [統計量]-[モデルへの適合]-[線形モデル]
  2. 変数としてx,yを設定している場合に,線形モデル窓の式の設定で
原点を通る直線: y ~ -1 + x
一次式(直線): y ~ x
二次式: y ~ x + I(x^2)
三次式: y ~ x + I(x^2) +I(x^3)

やはりlm関数を使用している。

 モデルのAIC(赤池情報量基準)を計算(Rcmdr)

 あるデータセットを説明するモデル(例えば多項式)を作る場合,どのモデルを選べばよいのかの一つの基準としてAICがある。AICが最少になるモデルを最適モデルとして選ぶ。

  1. [統計量]-[モデルへの適合]-[線形モデル] or [線形回帰] でモデルを決定する。
  2. [モデル]-[赤池情報量基準]
  3. 複数のモデルのAICを計算し,最少値になるモデルを選択する。

 重回帰式を求める(Rcmdr)

  1. [統計量]-[モデルへの適合]-[線形モデル]
  2. 目的変数,説明変数を設定する。式が決まっていて係数を決めるだけならここまで。
  3. 複数モデルを比較するなら,[モデル]-[赤池情報量基準]でAICを計算して比較する。

入力データ

> DatasetA
   緯度X1 経度X2 高度X3 温度Y
1   45.42 141.68    2.8  -8.0
2   43.77 142.37  111.9 -13.6
3   43.05 141.33   17.2  -9.5
4   40.82 140.78    3.0  -5.4
5   39.70 141.17  155.2  -6.7
6   38.27 140.90   38.9  -3.2
7   36.55 136.65   26.1  -0.1
8   36.67 138.20  418.2  -5.5
9   36.15 137.25  560.2  -7.6
10  36.33 138.55  999.1 -10.0
11  35.17 136.97   51.1  -0.9
12  35.52 137.83  481.8  -4.7
13  35.68 139.77    5.3  -0.4
14  35.48 134.23    7.1   0.5
15  35.02 135.73   41.4  -0.6
16  34.37 132.43   29.3   0.2
17  33.58 130.38    2.5   1.5
18  31.57 130.55    4.3   2.0
19  33.55 133.53    1.9   0.1
20  26.23 127.68   34.9  13.5

 20点の観測点のデータセット(緯度,経度,高度,気温)を使い,気温を目的変数とする重回帰式を作成する時にどの変数を説明変数とするのが最適かをAICを用いて判断する。以下に一部のモデルの計算結果を示すが,最終的には「気温~緯度+高度」のモデルのAICが最少となる。
 ※このような解析をする前には,散布図行列を描き変数間の相関をあらかじめ確認しておく必要がある。緯度と経度は本来独立な変数であるが,日本列島の形状から観測点の緯度と経度の間には相関がみられるので独立な変数として扱うのは適切ではない。

  • Model1 T(緯度,経度,高度)
> ML.1 <- lm(温度Y  ~ 緯度X1  + 経度X2  + 高度X3 , data=DatasetA)
> summary(ML.1)

Call:
lm(formula = 温度Y ~ 緯度X1 + 経度X2 + 高度X3, data = DatasetA)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8296 -0.6381 -0.1699  0.7527  3.6469 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 38.324695  23.946132   1.600    0.129      ;説明変数として必要?
緯度X1      -1.171541   0.215599  -5.434 5.52e-05 ***   ;説明変数として必要!
経度X2       0.023061   0.226127   0.102    0.920         ;必要?
高度X3      -0.009830   0.001679  -5.855 2.44e-05 ***   ;必要!
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.754 on 16 degrees of freedom
Multiple R-squared: 0.9246,	Adjusted R-squared: 0.9105 
F-statistic: 65.42 on 3 and 16 DF,  p-value: 3.362e-09 

> AIC(ML.1)
[1] 84.77708

  • Model2 T(緯度,経度)
> ML.2 <- lm(温度Y ~ 緯度X1 + 経度X2, data=DatasetA)
> summary(ML.2)

Call:
lm(formula = 温度Y ~ 緯度X1 + 経度X2, data = DatasetA)

Residuals:
   Min     1Q Median     3Q    Max 
-6.433 -1.892  0.381  2.493  4.365 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  94.2581    37.7618   2.496   0.0231 * ;説明変数として必要
緯度X1       -0.6935     0.3432  -2.021   0.0593 .  ;必要?
経度X2       -0.5242     0.3541  -1.480   0.1570    ;必要?
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 3.017 on 17 degrees of freedom
Multiple R-squared: 0.7631,	Adjusted R-squared: 0.7353 
F-statistic: 27.38 on 2 and 17 DF,  p-value: 4.824e-06 

> AIC(ML.2)
[1] 105.6765

  • Model3 T(緯度,高度)
> ML.3 <- lm(温度Y ~ 緯度X1 + 高度X3 , data=DatasetA)
> summary(ML.3)

Call:
lm(formula = 温度Y ~ 緯度X1 + 高度X3, data = DatasetA)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8425 -0.6485 -0.1034  0.7275  3.5930 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 40.742078   3.296653  12.359 6.40e-10 *** ;必要!
緯度X1      -1.151646   0.089081 -12.928 3.19e-10 *** ;必要!
高度X3      -0.009759   0.001484  -6.578 4.70e-06 ***  ;必要!
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.702 on 17 degrees of freedom
Multiple R-squared: 0.9246,	Adjusted R-squared: 0.9157 
F-statistic: 104.2 on 2 and 17 DF,  p-value: 2.878e-10 

> AIC(ML.3)
[1] 82.79008

※説明変数として意味があるかどうかは十分検討が必要であるが,有意水準5%で判断するのであればlm関数が出力する係数の後ろの記号が目安になる。

*** , ** , *  :説明変数として必要
. , なし  :説明変数としては不要

 LaTeX or html形式での計算結果の出力

>library(xtable)
>RegM.1 <- lm(y~x)
>a.tbl <- xtable(RegM.1) 
>print(a.tbl)                                   # LaTeX code output
>print(a.tbl,type=html)                   # HTML code output
RegWeight.jpeg

最終更新時間:2009年09月01日 19時12分51秒


HomePage > FrontPage