![[Chubo's Web Site]](./pimage/title.gif "[Chubo's Web Site]")
|
- The Chubo's World -
目次 Contents■トップへ■自己紹介 ■研究内容 ■研究業績 ■単語検索
実習関係■ネチケットについて■実習 ■メモ・雑文 ■最近のマイブーム ■リンク
■更新履歴 ■連絡先 古いHP■旧ホームページはこちらへ■目次 Contents ■私の作ったPDSのダウンロード 最新の15件
2010-01-28
2009-12-16
2009-12-15
2009-12-10
2009-11-09
2009-11-08
■Wiki形式でのページを2005年1月26日に開始しました． ■リンクはフリーですが，メールでお知らせいただければ幸いです． ■ご意見・ご感想・その他のお問い合わせがございましたら mailto:chubo@hiroshima-u.ac.jp までどうぞ． Online:
1
|
Rで統計 - R入門とりあえず使ってみたいけど，ようわからん人のために作りました．まだ作りかけです．
News
Rとは - Rってなあに
心の準備 - Rのインストール（for Windows）
使ってみよう - とりあえずやってみるべしRの実行Rのショートカットをダブルクリックして起動する．コンソール画面が現れ，以下のメッセージが表示された後，コマンドプロンプトが表示される．以下，コマンドプロンプトを>の記号で表す．
R : Copyright 2005, The R Foundation for Statistical Computing Version 2.0.1 Patched (2005-03-23), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > 英語だけど，一応説明は1回くらい読みましょう．
Rの基本 - これくらいは知っておこう>という記号は，プロンプトです．実際のコマンド入力の際には入力してはいけません．勝手に表示されます．以下の説明でも入れてある場合もありますが，実際の実行には不要です．
#コメントシャープがある場合は，それ以降はコメントです．これも入力不要です．
ls()Rがおぼえていること（オブジェクト）表示します．
ヘルプ - 困ったときには・・・
help()でヘルプが別窓で表示されます．でも英語．
help(command)でcommandと言う名前のコマンドのヘルプが表示されます（commandは適当に置き換えてください）．例えば， help(help)でhelpのヘルプが表示されます．
help(solve)または， ?solveなども試してみてください．基本的な文法などが表示されます． 他にも， help.start()で，html版のヘルプが表示されます．
Rの終了終わる場合はq()または quit()で終了．
簡単な計算例> 1+1 #足し算 [1] 2 > x <- 1 #xへ１を代入 > y <- 2 #yへ2を代入 > x+y #xとyの和を表示 [1] 3 > result <- (x+y)*10 #xとyの和を10倍したものをresultへ代入 > result #resultを表示 [1] 30 > result*10 #resultを10倍した値を表示 [1] 300 > sum(1,2,3,4,5,6,7,8,9,10) #1から10までの和を表示 [1] 55 > sum(1:10) #1から10までの和を表示 [1] 55 >
デモ - こんなことができるんですどんなことができるか，デモを見るといいです．demo()でいろいろなデモのオプションが表示されます． 例えば， demo(image) や
demo(graphics) でデモを見ることができます．グラフなどは別窓に表示されて，メインの窓に"Hit <Return> to see next plot:" と言われるので，適当に進んでください．
少しだけRしてみよう準備30個のデータを用意します．これはアカマツの葉の長さ（単位はcm）です．x <- c(7.4, 7.2, 6.5, 7.0, 7.6, 6.8, 7.8, 7.3, 7.7, 7.5, 7.1, 7.7, 7.1, 6.8, 7.3, 6.9, 6.7, 7.1, 7.2, 7.6, 7.6, 7.5, 7.8, 7.5, 7.2, 7.4, 7.0, 7.1, 7.3, 7.6)データがxというベクトル変数に格納されます．もちろん，ここをカットして，R上でペーストするだけで実行されます．
xとすると変数xの内容が表示されます． ついでに，アカマツの葉の生重量（単位はg）も別の変数に入れてください． y <-c(0.056, 0.059, 0.046, 0.052, 0.059, 0.056, 0.063, 0.058, 0.063, 0.047, 0.044, 0.048, 0.050, 0.045, 0.061, 0.053, 0.046, 0.050, 0.055, 0.055, 0.058, 0.052, 0.059, 0.056, 0.043, 0.047, 0.046, 0.057, 0.054, 0.063)
基本中の基本からそれでは，いろいろな基本統計量を計算してみましょう．> sum(x) #総和 [1] 218.3 > mean(x) #平均 [1] 7.276667 > median(x) #中央値 [1] 7.3 > var(x) #不偏分散 [1] 0.1149540 > sd(x) #標準偏差 [1] 0.3390487 > max(x) #最大値 [1] 7.8 > min(x) #最小値 [1] 6.5 > range(x) [1] 6.5 7.8 #範囲（最小値と最大値） という具合に，簡単に基本統計量を算出してくれます． さらに，2変数の間の基本統計量も算出できます． > var(x) #xの不偏分散 [1] 0.1149540 > var(y) #yの不偏分散 [1] 3.720575e-05 > var(x,y) #xとyの不偏共分散 [1] 0.001143333 > cor(x,y) #相関係数 [1] 0.552848なども計算できます．
関数の定義Rでは新しい関数の定義もできます．統計量の歪度と尖度を定義してみます．waido <- function(x) mean((x-mean(x))^3)/(sd(x)^3) #歪度 sendo <- function(x) mean((x-mean(x))^4)/(sd(x)^4) #尖度 実際に実行すると， > waido(x) [1] -0.3294951 > sendo(x) [1] 2.218392ちなみに，データが正規分布に近い場合は，それぞれ歪度が0，尖度が3に近い値になる．
ヒストグラムと箱ひげ図Rでは，ヒストグラムや箱ひげ図も簡単に描けます．リターンキーを押すと，別窓に表示されます．
> hist(x) #ヒストグラム
Hit <Return> to see next plot:
> boxplot(x, names=c("x")) #箱ひげ図
Hit <Return> to see next plot:
変数から変数への値の代入例えば，変数の単位をgからmgに変更したい場合，簡単に変数に代入できます．> ymg <- y*1000 > y [1] 0.056 0.059 0.046 0.052 0.059 [6] 0.056 0.063 0.058 0.063 0.047 [11] 0.044 0.048 0.050 0.045 0.061 [16] 0.053 0.046 0.050 0.055 0.055 [21] 0.058 0.052 0.059 0.056 0.043 [26] 0.047 0.046 0.057 0.054 0.063 > ymg [1] 56 59 46 52 59 56 63 58 63 47 44 [12] 48 50 45 61 53 46 50 55 55 58 52 [23] 59 56 43 47 46 57 54 63
t検定してみよう - 1標本t検定Rではいろいろな検定を行うため，関数があらかじめ用意去れいます． t検定を例に，使い方を見てみましょう．取り合え時，1標本のt検定です．母平均が7であるかどうか検定して，さらに95%の信頼区間を計算します．
> t.test(x, mu=7)
One Sample t-test
data: x
t = 4.4695, df = 29, p-value = 0.0001105
alternative hypothesis: true mean is not equal to 7
95 percent confidence interval:
7.150064 7.403270
sample estimates:
mean of x
7.276667
p-valueが0.05よりも小さいので優位差があります．つまり，帰無仮説「標本平均値と母平均値7には差がない」が棄却され，対立仮説「母平均が7ではない」が採用でき，結局母平均は7とは等しくないということになります．
それでは，母平均が7.25であるかの検定をするとどうなるでしょう．
> t.test(x, mu=7.25)
One Sample t-test
data: x
t = 0.4308, df = 29, p-value =
0.6698
alternative hypothesis: true mean is not equal to 7.25
95 percent confidence interval:
7.150064 7.403270
sample estimates:
mean of x
7.276667
p-valueが0.05よりも大きいので優位差がなく，帰無仮説は保留されます．結局，母平均の値が7.25ではないとは言えないことになります．簡潔に言えば，まあ母平均はそんな値でしょうということです．
t検定してみよう - 2標本t検定2つのグループの平均が等しいかどうかを検定してみます．例は良くありませんが，先から使っているxとyを例にして計算します．とりあえず，等分散を仮定しません．
> t.test(x,y)
Welch Two Sample t-test
data: x and y
t = 116.6713, df = 29.019, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
7.09668 7.34992
sample estimates:
mean of x mean of y
7.27666667 0.05336667
p-valueが0.05よりも明らかに小さいので帰無仮説が棄却され，yの平均が，xの平均よりも小さいことが読み取れます．等分散を仮定する場合は，var.equal=Tを入れることで計算できます．
> t.test(x, y, var.equal=T)
Two Sample t-test
data: x and y
t = 116.6713, df = 58, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
7.09937 7.34723
sample estimates:
mean of x mean of y
7.27666667 0.05336667
ちなみに， t.test(x, y, var.equal=F)とすると，ウェルチの検定（Welch Two Sample t-test）になります．
データをプロットしてみようデータをプロットしてみましょう．"Hit <Return> to see next plot: "と表示されたら，これまでと同様にリターンキーを押してください．
plot(x,y) 別窓でプロットされたはずです．
回帰してみようプロットされたら，回帰したくなるのが人情です．
> result <- lm(y ~ x)
> abline(result)
> summary(result)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.0096041 -0.0030987 0.0004202 0.0041459 0.0074013
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.019007 0.020637 -0.921 0.36490
x 0.009946 0.002833 3.511 0.00153 **
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 0.005173 on 28 degrees of freedom
Multiple R-Squared: 0.3056, Adjusted R-squared: 0.2808
F-statistic: 12.32 on 1 and 28 DF, p-value: 0.001533
データの入出力エクセルからの出力これまでデータシートの作成にエクセルMS-Excelを利用している場合が少なくないと思います．これまでにあるデータをエクセルで利用したい場合，CSVファイルを経由してRへの取り込みができます．
エクセルからの入力CSV形式で保存されたファイルがあれば，Rに容易にデータを取り込める．
x <- read.csv("data.csv")
エクセルから，クリップボードを使ったデータ取り込みエクセルのデータを取り込む関数です．青木先生のサイトにあります．
excel.w <- function(nc)
{
matrix(scan("clipboard"), byrow=TRUE, ncol=nc)
}
使い方は，エクセルで使いたいデータをコピーして，関数を実行して変数に代入するだけです． 例えばここでは30個のデータを変数xに入れます．
x <- excel.w(30)
もう少しプロットしてみよう正規分布をプロットする標準正規分布をplotを使ってプロットしてみましょう． まず，正規分布の密度関数を定義します．
gauss.density <- function(x) {
1/sqrt(2*pi)*exp(-x^2/2)
}
さっそく，プロットしてみましょう． plot(gauss.density, -3, 3) 表示する範囲も変更できます． plot(gauss.density, -4, 4) これでもプロットできます． curve(dnorm, -4, 4)
お決まりどこの世界もそうですが，はじめてだと何のことだかわからない暗黙の了解というのがRの世界にもあります．例えばこれ．> data(iris)これは，多変量解析を少しでもかじったことのある人ならば知っているであろう，かの有名なフィッシャーのアヤメの分類データを読み込むコマンドです．データirisを準備しなくても，このコマンド一発で読み込んでくれます．次のコマンドで内容を確認できます．IrisData
> 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
.
.
.
146 6.7 3.0 5.2 2.3 virginica
147 6.3 2.5 5.0 1.9 virginica
148 6.5 3.0 5.2 2.0 virginica
149 6.2 3.4 5.4 2.3 virginica
150 5.9 3.0 5.1 1.8 virginica
散布図やグラフを描いてみようIris Dataを例に，散布図やグラフを描いてみます．
散布図の例デモにもありますが，おきまりのやつを．data(iris)をしておいて， pairs(iris[1:4], main = "Iris Data", font.main = 2, pch = 19)とすると，黒丸でプロットしてくれます．pchの値19というのが黒丸指定です． "Hit <Return> to see next plot: "と出るので，リターンキーを押してください． 色を付けてみましょう．
pairs(iris[1:4], main = "Iris Data", font.main = 2, pch = 21,
bg = c("red", "green3", "blue")[unclass(iris$Species)])
で色が付きます．色を変えることもできます．
pairs(iris[1:4], main = "Iris Data", font.main = 2, pch = 21,
bg = c("springgreen3", "red1", "steelblue3")[unclass(iris$Species)])
ちなみに使える色は， colors()で表示されます．
データを重ねてプロットI. setosaとI. virginicaの花弁とがく片のデータを重ねてプロットします．
data(iris)
table(iris$Species)
iS <- iris$Species == "setosa"
iV <- iris$Species == "virginica"
matplot(c(1,8), c(0,5), type="n",
xlab="Length", ylab="Width",
main="Petal and Sepal Dimensions in Iris")
matpoints(iris[iS,c(1,3)], iris[iS,c(2,4)], pch="sS",
col=c("red", "green3"))
matpoints(iris[iV,c(1,3)], iris[iV,c(2,4)], pch="vV",
col=c("red", "green3"))
大文字が花弁（Petal）で，小文字ががく片（Sepal）です．ちなみに単子葉植物なので，基本的に花弁が3枚，がく片が3枚あり，アヤメ属Irisの場合，内側にある花弁よりも外側にあるがく片の方が大きくなって目立ちます．日本にも分布するI. setosaの場合，花弁が小さく，がく片が大きくなり，両者の差が非常に大きいのですが， I. virginicaの場合，花弁とがく片の差がI. setosaに比べて小さいことが分かります．
主成分分析青木先生のページの主成分分析より
pca <- function(dat)
{
dat <- dat[complete.cases(dat),]
nr <- nrow(dat)
nc <- ncol(dat)
heikin <- apply(dat, 2, mean)
sd <- sqrt(bunsan <- apply(dat, 2, var))
eval <- (result <- eigen(r <-cor(dat)))$values
evec <- result$vectors
cum.contr <- cumsum(contr <- eval/nc*100)
fl <- t(sqrt(eval)*t(evec))
fs <- scale(dat)%*%evec*sqrt(nr/(nr-1))
# add names
names(heikin) <- names(bunsan) <- names(sd) <- rownames(r) <- colnames(r) <- rownames(fl) <- colnames(dat)
names(eval) <- names(contr) <- names(cum.contr) <- colnames(fl) <- colnames(fs) <- paste("PC", 1:nc)
invisible(list(mean=heikin, variance=bunsan, standard.deviation=sd, r=r, factor.loadings=fl, eval=eval, contribution=contr, cum.contribution=cum.contr, fs=fs))
}
Fisherのirisのデータでの計算例 data(iris) result <- pca(iris[1:4]) result #第1主成分と第2主成分の値でプロットする． plot(result$fs[,1:2], pch=as.integer(iris$Species)) #他の組み合わせでプロットしたい場合 pc1 <- result$fs[,1] pc2 <- result$fs[,2] pc3 <- result$fs[,3] plot(pc1, pc3, pch=as.integer(iris$Species))
発展オブジェクトの消去objectという名前のオブジェクトを消去します．
rm(object)
パッケージ（ライブラリ）についてRではいろいろな関数などをまとめてライブラリにしています．普通これをパッケージと言ったりします．今入っているパッケージを確認したい場合は， library()とすると別窓に一覧が表示されます．
Packages in library 'C:/PROGRA~1/R/RW2001~1/library':
base The R Base Package
boot Bootstrap R (S-Plus) Functions (Canty)
class Functions for Classification
cluster Functions for clustering (by Rousseeuw et al.)
datasets The R Datasets Package
foreign Read Data Stored by Minitab, S, SAS, SPSS,
Stata, Systat, dBase, ...
graphics The R Graphics Package
grDevices The R Graphics Devices and Support for Colours
and Fonts
grid The Grid Graphics Package
KernSmooth? Functions for kernel smoothing for Wand &
Jones (1995)
lattice Lattice Graphics
MASS Main Package of Venables and Ripley's MASS
methods Formal Methods and Classes
mgcv GAMs with GCV smoothness estimation and GAMMs
by REML/PQL
nlme Linear and nonlinear mixed effects models
nnet Feed-forward Neural Networks and Multinomial
Log-Linear Models
rpart Recursive Partitioning
spatial Functions for Kriging and Point Pattern
Analysis
splines Regression Spline Functions and Classes
stats The R Stats Package
stats4 Statistical Functions using S4 Classes
survival Survival analysis, including penalised
likelihood.
tcltk Tcl/Tk Interface
tools Tools for Package Development
utils The R Utils Package
それぞれのパッケージの中身が知りたい場合は， library(help="base") などとすれば，
Information on Package 'base'
Description:
Package: base
Version: 2.0.1
Priority: base
Title: The R Base Package
Author: R Development Core Team and contributors worldwide
Maintainer: R Core Team <R-core@r-project.org>
Description: Base R functions
License: GPL Version 2 or later.
Built: R 2.0.1; ; 2005-03-24 03:27:30; windows
Index:
.Device Lists of Open Graphics Devices
.Internal Call an Internal Function
.Last.value Value of Last Evaluated Expression
.Library Search Paths for Packages
.Machine Numerical Characteristics of the Machine
.
.
.
こんな感じで表示されます．
参考文献・インターネットリソース
Counter: 4331,
today: 1,
yesterday: 24
|
Last-modified: 2007-01-04(木) 09:11:52