- 追加された行はこの色です。
- 削除された行はこの色です。
-[[英語学研究2008-2]] | [[英語学研究2009-2]] | [[英語学研究2010-後期]] | [[英語学研究2011]] | [[Rのメモ]] | [[SLWRMemo]] | [[Gries SFLWR>http://www.linguistics.ucsb.edu/faculty/stgries/research/sflwr/sflwr.html]]
#contents
**ヘルプ参照 [#z5ba9ed0]
> help.start() # 全体のヘルプを出す場合
> help(rnorm) # function のヘルプを参照する場合、ただ help()に入れればよい
**GUIの言語を変えたい [#pcb26f21]
-Rの起動オプション内で指定するのが楽な模様です。
-起動オプションを付与するには、Rを起動するショートカットを右クリックして、「プロパティ」から「リンク先(T):」の最後にコマンドを書き込みます。
> LANGUAGE=en #英語で起動したい場合
> "C:\Program Files\R\R-3.0.0\bin\x64\Rgui.exe" LANGUAGE=en #書き込むと左のようなショートカットになります。
**Function と argument [#g5e0231d]
-log(x) log() が function, x が argument
-plot(x,y) plot() が function, x, y が argument
**Vector ベクトル [#gffc451f]
-Numerical vector: 数値のベクトル
-Character vector: 文字列のベクトル
> wordlist <- c("able", "ability", "abuse")
> wordlist
[1] "able" "ability" "abuse"
**パッケージの読み込み [#kea71742]
-コマンドラインでライブラリを読み込む場合：library(xxxx) と打つ
> library(survival)
要求されたパッケージ splines をロード中です
次のパッケージを付け加えます: 'survival'
The following object(s) are masked from package:ISwR :
lung
-library を閉じる場合：
> detach("package:survival")
-利用可能なlibrary を参照する場合：
> library()
**テキストファイルからの読み込み [#o74de443]
-（例）c:\temp にある lemma.num を bnclemma というデータフレームにする：
> bnclemma <- read.table("c:/temp/lemma.num", header=T)
header=T で１行目をヘッダとして読み込む
**テキストファイルへの書き出し [#k6a104af]
>write.table (データフレーム名, file="ファイル名")
**データエディタの起動 [#d9e3bbd5]
-作成したデータフレームをスプレッドシートで編集したい場合：
> lm <- edit(bnclemma)
--編集したデータフレームは lm という名前（何でもよい）になって保存される
--bnelemma という元のフレームはそのまま
-上書きで編集するには fix() とする
> fix(bnclemma)
-まったく新しいファイルを書くには
>dd <- data.frame()
>fix(dd)
**ワークスペースの履歴・保存 [#i65b2a7d]
-ls() 関数
-今までの実行したワークスペースのコマンド一覧を出す
> ls()
[1] "bmi" "d" "energy" "exp.lean" "exp.obese" "fpain"
[7] "h" "height" "hh" "intake.post" "intake.pre" "l"
[13] "mylist" "oops" "pain" "thue2" "thuesen" "weight"
[19] "x" "y" "ylim"
-ワークスペースをファイルに保存する
>save.image()
**データフレームの部分加工 [#n225bbb0]
>icci2 <-icci[c(1,2)]
#icci のデータフレームの第1、2コラムを取り出して icci2 を作る
**データフレームのサブセットの取り出し [#h0326953]
> icci_A1 <- subset(icci, CEFR == "A1")
#subset()関数で、データ icci の CEFRの変数の値が A1 のものだけを icci_A1 に格納
**データフレームのマージ [#g73594ce]
>icci <- merge(icci1, icci2)
# これでお互いのファイルに共通する変数名を参照しながらマージする
# もしそれぞれにない個体があった場合に、それをすべて含めたければ all=TRUE とする
>icci <- merge(icci1, icci2, all=TRUE)
**データフレームのマージ（行の結合） [#aad41245]
> total <- rbind(A1_sample, A2_sample, B1_sample)
# rbind で2つ以上のデータフレームをくっつける。カラム変数は同じでないとダメ。
**データフレームからランダムに X 件を抽出 [#z177b7b2]
> A1_sample <- icci_A1[sample(1:nrow(icci_A1), 1000, replace=FALSE),]
# icci_A1 の行データ（1から個体数すべて）のなかから1000件をランダムにサンプルして A1_sample に格納
**分布 [#w2f350a6]
-Density or point probability 密度確率・
--the probability density function 確率密度関数
--the probability distribution function 確率分布関数
-確率分布関数＝確率密度関数を積分したもの
-確率密度関数＝確率分布を微分したもの
--統計データが十分にたくさんあれば、
--ヒストグラム ＝ 確率密度関数になり
--累積頻度分布 ＝ 確率分布関数になる。
***正規曲線の書き方１ [#a43c63aa]
> x <- seq(-4,4,0.1)
-4 から 4 まで 0.1ずつの感覚で sequence を作って x に格納
> plot(x, dnorm(x),type="l")
X軸に x を Y 軸に x にそって確率密度の正規曲線を生成
type="l" は line かな
> plot(x, dnorm(x),type="h")
同じプロットを h にすると histogram のようにそれぞれの x 軸に縦棒が伸びる形式になる
***正規曲線の書き方２ [#dbcc2dff]
> curve(dnorm(x), from=-4, to=4)
これでも同じカーブを描くことができる
***累積分布関数 [#wd8f1286]
-正規分布のどの位置にあるか、その位置までの cumulative distribution を出すことで、その値の正規分布での外れ具合（有意確率）を見る
--平均１３２、標準偏差が１３の正規分布の中で、１６０という数字よりも小さい値に入る
---正規分布上の値のとる確率は何パーセントか、を pnorm が計算
---それを 1 から引くことで、160以上の値をとる確率を得ている
> 1-pnorm(160,mean=132,sd=13)
[1] 0.01562612
> 1-pbinom(15,size=20,prob=.5)
[1] 0.005908966
**記述統計 [#kaed79a2]
-サンプルデータとしてランダムに５０個の正規分布を成す数値を発生させる
> x <- rnorm(50)
> x
[1] 1.11403983 -1.15165895 0.80461867 1.71557037 0.03169234 0.18838381 -0.80416122
[8] -0.10692887 -0.96407973 1.16309160 0.61545475 0.09261578 -0.76893852 -0.78379193
[15] 0.12446743 0.77148385 1.19391658 -0.78081663 -0.83306264 0.35587091 -0.76966258
[22] -0.15209363 -0.49532705 -0.52829549 0.64414085 1.59186102 -0.07530280 -0.35851963
[29] 1.55049255 0.21063906 0.22529564 0.34462737 0.19350551 -0.65114617 -1.49958077
[36] 1.15710520 -0.69795906 -1.12325445 -2.08717916 0.01182514 0.13297142 1.40239832
[43] 1.04982555 0.46580345 0.55212419 -0.53356749 -0.31117847 0.64688955 -1.44102550
[50] -0.06345214
-個別の記述統計の出し方
--平均
> mean(x)
[1] 0.02739456
-欠損値があって飛ばしたい場合の処理
> library(ISwR) ＃ISwR というライブラリを読み込む（データのため）
> data(juul) ＃data juul を読み込む
> juul
age menarche sex igf1 tanner testvol
1 NA NA NA 90 NA NA
2 NA NA NA 88 NA NA
3 NA NA NA 164 NA NA
4 NA NA NA 166 NA NA
5 NA NA NA 131 NA NA
--juul のデータの igf1 を参照したいが、いちいち juul$igf1 と書くのは面倒くさいので attach 関数で
> attach(juul) ＃こうすると juul$…を省略できる
> mean(igf1)
[1] NA ＃NA は欠損値があるので平均値を出せない、とエラーメッセージを返している
> mean(igf1,na.rm=T) ＃"na.rm=T" not available, remove の意味 欠損値を除く、という指定
[1] 340.168
--標準偏差
> sd(x)
[1] 0.876461
--分散
> var(x)
[1] 0.768184
--中央値（メディアン）
> median(x)
[1] 0.06215406
--分位数（quantile）
> quantile(x)
0% 25% 50% 75% 100%
-2.08717916 -0.68625583 0.06215406 0.63696932 1.71557037
---デフォルトではこのように正規分布の下から 25,50, 75,100% の位置にある値を出す
---もし quantile を特定の観察したい割合で切りたい場合には以下のようにする
> pvec <- seq(0,1,0.1)
pvec という変数に 0 から 1 までを 0.1 区切りに格納
> pvec
[1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
> quantile(x,pvec)
0% 10% 20% 30% 40% 50% 60%
-2.08717916 -0.97999720 -0.77189339 -0.52987709 -0.12499478 0.06215406 0.20035893
70% 80% 90% 100%
0.49169967 0.77811082 1.16617410 1.71557037
---quantile に２つ目の argument で何％で切るかを 0.00-1.00 範囲で指定する →この例は 0.1 = 10% ずつ
**記述統計を一発で出す方法 [#fc2b1dfa]
-summary 関数は 最小値、第１分位点、メジアン、平均、第３分位点、最大値、欠損数を出す。
> summary(igf1)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
25.0 202.2 313.5 340.2 462.8 915.0 321.0
-データフレーム全体の summary も可能
> summary(juul)
age menarche sex igf1 tanner
Min. : 0.170 Min. : 1.000 Min. :1.000 Min. : 25.0 Min. : 1.000
1st Qu.: 9.053 1st Qu.: 1.000 1st Qu.:1.000 1st Qu.:202.2 1st Qu.: 1.000
Median :12.560 Median : 1.000 Median :2.000 Median :313.5 Median : 2.000
Mean :15.095 Mean : 1.476 Mean :1.534 Mean :340.2 Mean : 2.640
3rd Qu.:16.855 3rd Qu.: 2.000 3rd Qu.:2.000 3rd Qu.:462.8 3rd Qu.: 5.000
Max. :83.000 Max. : 2.000 Max. :2.000 Max. :915.0 Max. : 5.000
NA's : 5.000 NA's :635.000 NA's :5.000 NA's :321.0 NA's :240.000
testvol
Min. : 1.000
1st Qu.: 1.000
Median : 3.000
Mean : 7.896
3rd Qu.: 15.000
Max. : 30.000
NA's :859.000
**ヒストグラムを書く [#xb4b861d]
> x <- rnorm(100) # normal density of pseudo-random numbers
> x
[1] -0.23852918 0.80007985 2.59488947 1.92574080 0.04232929 -1.06697286 -1.28712104
[8] -0.07618785 -1.59706302 0.82031697 -1.20817887 1.46262997 -1.29586205 -0.76564449
[15] 0.75213738 -0.06597683 0.23352029 -0.16270841 -1.64600217 0.75705638 0.46523458
[22] -0.12607458 0.09170498 -0.10259290 -0.45916474 0.25041968 -0.51401120 -0.23879437
[29] 1.45689028 0.64155579 -0.33764720 -1.25446597 0.70214707 0.44216508 -0.04646707
[36] -0.26736349 -0.97858879 1.45124404 -1.22135682 0.08365315 -0.26954490 -0.22825442
[43] -0.45558029 1.61702489 -1.51530741 1.12869612 1.03321808 1.67374952 -0.46476407
[50] -0.73521276 -0.12948352 -0.09378984 0.93129898 0.77798970 0.54310765 0.52949773
[57] 0.25044045 1.01388932 0.82988086 -1.59370935 -0.60922676 1.53202118 0.19678465
[64] -0.38521575 1.30834830 1.20246317 -0.03461022 -0.60884186 0.89151459 0.48934389
[71] 0.99960792 0.87056316 1.17321910 0.02667156 0.52161582 -0.09861786 -0.63170346
[78] -0.85243182 -0.87548459 0.20321376 -1.46930772 1.31395373 -1.09770929 -0.19700751
[85] -0.68914995 1.06938790 1.29973824 0.65570765 0.71533448 -2.32653615 0.06533285
[92] -0.25530870 0.46701640 1.27102185 1.73760250 1.31901418 1.08818204 -1.94065162
[99] -0.32908069 -0.49896537
> hist(x, freq=F) # write histogram of x, freq=F means "not based on frequency"
> curve(dnorm(x),add=T) # write the normal curve graph, add means "overplot"
-縦軸の目盛りが切れてしまった場合の処理（少し複雑）
> h <- hist(x, plot=F)
> ylim <- range(0, h$density, dnorm(0))
> hist(x, freq=F, ylim=ylim)
> curve(dnorm(x), add=T)
**変数:数値と文字列の操作 [#a703d170]
-たとえば energy というデータは expend, stature という２変数をもち、expend は数値データ、stature はカテゴリーデータ
> data(energy)
> energy
expend stature
1 9.21 obese
2 7.53 lean
3 7.48 lean
4 8.08 lean
5 8.09 lean
6 10.15 lean
7 8.40 lean
8 10.88 lean
9 6.13 lean
10 7.90 lean
11 11.51 obese
12 12.79 obese
13 7.05 lean
14 11.85 obese
15 9.97 obese
16 7.48 lean
17 8.79 obese
18 9.69 obese
19 9.68 obese
20 7.58 lean
21 9.19 obese
22 8.11 lean
-lean の１行だけとりだす
> exp.lean <- energy$expend[energy$stature=="lean"]
-obese の行だけとりだす
> exp.obese <- energy$expend[energy$stature=="obese"]
> exp.lean
[1] 7.53 7.48 8.08 8.09 10.15 8.40 10.88 6.13 7.90 7.05 7.48 7.58 8.11
> exp.obese
[1] 9.21 11.51 12.79 11.85 9.97 8.79 9.69 9.68 9.19
-同様のことが以下の split()関数でもできる
> l <-split(energy$expend, energy$stature)
>
> l
$lean
[1] 7.53 7.48 8.08 8.09 10.15 8.40 10.88 6.13 7.90 7.05 7.48 7.58 8.11
$obese
[1] 9.21 11.51 12.79 11.85 9.97 8.79 9.69 9.68 9.19
**データ・フレームの一部取り出し [#h2e2b342]
（例）verbs というデータから AnimacyOfTheme が animate のものだけを取り出す
> verbs[verbs$AnimacyOfTheme == "animate",]
RealizationOfRec Verb AnimacyOfRec AnimacyOfTheme LengthOfTheme
58 NP give animate animate 1.0986123
100 NP give animate animate 2.8903718
143 NP give inanimate animate 2.6390573
390 NP lend animate animate 0.6931472
506 NP give animate animate 1.9459101
736 PP trade animate animate 1.6094379
> subset(verbs, AnimacyOfTheme == "animate")
RealizationOfRec Verb AnimacyOfRec AnimacyOfTheme LengthOfTheme
58 NP give animate animate 1.0986123
100 NP give animate animate 2.8903718
143 NP give inanimate animate 2.6390573
390 NP lend animate animate 0.6931472
506 NP give animate animate 1.9459101
736 PP trade animate animate 1.6094379
**標準化得点を出す方法 [#rb5c6194]
> x=rnorm(10,3,0.1)
> x
[1] 2.987822 3.087673 2.972487 2.901589 3.056770 3.079720 2.994168 2.985317 2.869461 3.120720
> x-mean(x)
[1] -0.01775045 0.08210025 -0.03308569 -0.10398376 0.05119739 0.07414762 -0.01140460 -0.02025609 -0.13611205 0.11514739
> y = x-mean(x)
> y/sd(x)
[1] -0.2184789 1.0105193 -0.4072305 -1.2798693 0.6301557 0.9126354 -0.1403719 -0.2493192 -1.6753157 1.4172752
> scale(x)
[,1]
[1,] -0.2184789
[2,] 1.0105193
[3,] -0.4072305
[4,] -1.2798693
[5,] 0.6301557
[6,] 0.9126354
[7,] -0.1403719
[8,] -0.2493192
[9,] -1.6753157
[10,] 1.4172752
attr(,"scaled:center")
[1] 3.005573
attr(,"scaled:scale")
[1] 0.08124561
**さまざまなプロット [#o8ba5c7c]
-棒グラフ
--ratings 中の１つの変数 Length に関して crosstab を作り、棒グラフにする、グラフのラベルや色を指定する
>barplot(xtabs( ~ ratings$Length), xlab="word length", col="grey")
-ヒストグラム
--ratings 中の１つの変数 Length に関して、ヒストグラムを作成
--この truehist というコマンドは MASS というパッケージ中の関数
>truehist(ratings$Length, xlab="word length", col="grey")
-一般的な描画関数 plot()
>plot(h)
>plotsort(lexdec$RT), ylab = "log RT") # RT という変数でソートしたものをプロットする
-箱ひげ図を書く
>boxplot(lexdec$RT) # ふつうに RT という変数を箱ひげ図にしたもの
>boxplot(exp(lexdec$RT)) # RT は対数なのでそれをふつうの数字に戻してプロットする
-すべての変数間のプロットを一発で出す → scatterplot matrix 関数は pairs()
>pairs(rating[ , -c(1, 6:8, 10:14)])
-カテゴリーデータをモザイク図にする
>mosaicplot(verbs.xtabs, main = "dative")
--上位階層にしたいものを先に xtab で処理しておく
-上位階層を変数でくくって箱髭図を書く
**複数の表を並べて作成する場合 [#cb1be9cd]
-par()関数を使うと、以後作成する表を１つ１つ並べて表示してくれる
>par(mfrow = c(3,2))
--2列3行で表を並べてくれる
**表を画像にはき出す場合 [#rcdf67f6]
-jpeg で barplot.jpeg というファイルに書き出す
>jpeg("barplot.jpeg", width = 400, height=420)
>truehist ... # ここで作表を行う
>dev.off() # 画像モードを抜ける
**単語リストの操作 [#v2bfd808]
-テキスト・ファイルを読み込む
（例）c:\temp にある lemma.num を bnclemma というデータフレームにする：
> bnclemma <- read.table("c:/temp/lemma.num", header=T)
header=T で１行目をヘッダとして読み込む
-データフレームを編集する
--作成したデータフレームをスプレッドシートで編集したい場合：
> lm <- edit(bnclemma)
---編集したデータフレームは lm という名前（何でもよい）になって保存される
---lemmalist として lm というデータフレームを作成した
-単語リストの抽出
--データフレームを行列として参照する
---ランク１０位の単語（[,3]はwordのカラムを示す）を参照
> lm[10,3]
[1] to
---ランク１０位の品詞（[,4]）を参照
> lm[10,4]
[1] prep
Levels: a adv conj det infinitive-marker interjection modal n prep pron v
※このようにアウトプットの最後に、その変数に何レベルあるかが列挙される
---lm 中の word のカラムの最初の４個を取り出す
> lm$word[1:4]
[1] the be of and
5464 Levels: a abandon abbey ability able abnormal abolish abolition abortion about ... zone
---lm 中の pos のカラムの１番目を取り出す
> lm$pos[1]
[1] det
Levels: a adv conj det infinitive-marker interjection modal n prep pron v
---lm 中のwordの項目数をカウントする → length()を使用
> length(lm$word)
[1] 6318
---rank が１００位以内の word を取り出す
> lm$word[lm$rank < 100]
[1] the be of and a in to have it to for
[12] i that you he on with do at by not this
[23] but from they his that she or which as we an
[34] say will would can if their go what there all get
[45] her make who as out up see know time take them
[56] some could so him year into its then think my come
[67] than more about now last your me no other give just
[78] should these people also well any only new very when may
[89] way look like use her such how because when as good
---rank が１００位以内で、pos が名詞のものを取り出す
> lm$word[lm$rank < 100 & lm$pos == "n"]
[1] time year people way
※ [ ] で条件設定し、& で２つの条件を組み合わせている
※ logical operator としては、&, |(=or), !(=not) が使える
---rank が１００位以内で、pos が名詞または形容詞のものを取り出す
> lm$word[lm$rank < 100 & (lm$pos == "n"|lm$pos == "a")]
[1] time year last other people new way good
※ logical operator | が使われていることに注意
---１０位以内のリストをすべてとりだす
> lm[lm$rank <10,]
rank freq word pos
1 1 6187267 the det
2 2 4239632 be v
3 3 3093444 of prep
4 4 2687863 and conj
5 5 2186369 a det
6 6 1924315 in prep
7 7 1620850 to infinitive-marker
8 8 1375636 have v
9 9 1090186 it pron
※１０位まで取り出すならば、<11 とする
※[]の中で条件を付けるが、最後の ,] とカンマをつけるところが重要。
---１００位以内で品詞が動詞である行だけをとりだす
> lm[(lm$rank <101 & lm$pos == "v"),]
rank freq word pos
2 2 4239632 be v
8 8 1375636 have v
18 18 559596 do v
34 34 333518 say v
40 40 249540 go v
44 44 220940 get v
46 46 217268 make v
51 51 191661 see v
52 52 185534 know v
54 54 179220 take v
64 64 153881 think v
66 66 151871 come v
76 76 131417 give v
90 90 111058 look v
92 92 108820 use v
100 100 98899 find v
-１００位以内の単語の頻度プロットを作る
> lm100 <- lm[(lm$rank <101),] # lm100 というデータフレームに１００以内を格納
> lm100
rank freq word pos
1 1 6187267 the det
2 2 4239632 be v
3 3 3093444 of prep
4 4 2687863 and conj
5 5 2186369 a det
6 6 1924315 in prep
7 7 1620850 to infinitive-marker
8 8 1375636 have v
9 9 1090186 it pron
10 10 1039323 to prep
11 11 887877 for prep
12 12 884599 i pron
13 13 760399 that conj
14 14 695498 you pron
15 15 681255 he pron
16 16 680739 on prep
17 17 675027 with prep
18 18 559596 do v
19 19 534162 at prep
20 20 517171 by prep
※出力はこんな感じ
プロットの関数は plot()
> plot(lm100$rank, log(lm100$freq), type="n")
※ランキングを横軸、頻度の自然対数を縦軸にとる。
※type="n" はプロットの記号を消す、という意味。後で単語を載せるためここでは出力しない。
> text(lm100$rank, log(lm100$freq), as.character(lm100$word), cex=0.5)
※この text() という関数が交差する点に as.character() で指定した文字列を載せる機能
※cex=で文字サイズを指定する
--頻度リスト＆プロット：応用 [#c0377bc4]
---１００語以内の動詞のプロットを作る
> vlm <- lm100[lm100$pos == "v",] #動詞を１００位以内から抽出 vlm に格納
> vlm
rank freq word pos
2 2 4239632 be v
8 8 1375636 have v
18 18 559596 do v
34 34 333518 say v
40 40 249540 go v
44 44 220940 get v
46 46 217268 make v
51 51 191661 see v
52 52 185534 know v
54 54 179220 take v
64 64 153881 think v
66 66 151871 come v
76 76 131417 give v
90 90 111058 look v
92 92 108820 use v
100 100 98899 find v
> plot(vlm$rank, log(vlm$freq), type="n") #横軸ランク、縦軸に頻度の対数をとる、type="n" でプロットの記号を消しておく
> text(vlm$rank, log(vlm$freq), as.character(vlm$word), cex=0.9)
※vlm$rank と vlm$freq の交差する点に vlm$word から取ってきた文字列を配置するというコマンド
※cex=0.9 で少し大きめに
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
