ヘルプ参照

> help.start()   # 全体のヘルプを出す場合
> help(rnorm)     # function のヘルプを参照する場合、ただ help()に入れればよい

GUIの言語を変えたい

  • Rの起動オプション内で指定するのが楽な模様です。
  • 起動オプションを付与するには、Rを起動するショートカットを右クリックして、「プロパティ」から「リンク先(T):」の最後にコマンドを書き込みます。
> LANGUAGE=en #英語で起動したい場合
> "C:\Program Files\R\R-3.0.0\bin\x64\Rgui.exe" LANGUAGE=en #書き込むと左のようなショートカットになります。

Function と argument

  • log(x) log() が function, x が argument
  • plot(x,y) plot() が function, x, y が argument

Vector ベクトル

  • Numerical vector: 数値のベクトル
  • Character vector: 文字列のベクトル
> wordlist <- c("able", "ability", "abuse")
> wordlist
[1] "able"    "ability" "abuse"  

パッケージの読み込み

  • コマンドラインでライブラリを読み込む場合:library(xxxx) と打つ
> library(survival)
要求されたパッケージ splines をロード中です 

次のパッケージを付け加えます: 'survival' 

       The following object(s) are masked from package:ISwR :

        lung 
  • library を閉じる場合:
> detach("package:survival")
  • 利用可能なlibrary を参照する場合:
> library()

テキストファイルからの読み込み

  • (例)c:\temp にある lemma.num を bnclemma というデータフレームにする:
> bnclemma <- read.table("c:/temp/lemma.num", header=T)

	header=T で1行目をヘッダとして読み込む

テキストファイルへの書き出し

>write.table (データフレーム名, file="ファイル名")

データエディタの起動

  • 作成したデータフレームをスプレッドシートで編集したい場合:
> lm <- edit(bnclemma)
  • 編集したデータフレームは lm という名前(何でもよい)になって保存される
  • bnelemma という元のフレームはそのまま
  • 上書きで編集するには fix() とする
> fix(bnclemma)
  • まったく新しいファイルを書くには
>dd <- data.frame()
>fix(dd)

ワークスペースの履歴・保存

  • 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()

データフレームの部分加工

>icci2 <-icci[c(1,2)]
#icci のデータフレームの第1、2コラムを取り出して icci2 を作る

データフレームのサブセットの取り出し

> icci_A1 <- subset(icci, CEFR == "A1")
#subset()関数で、データ icci の CEFRの変数の値が A1 のものだけを icci_A1 に格納

データフレームのマージ

>icci <- merge(icci1, icci2)
# これでお互いのファイルに共通する変数名を参照しながらマージする
# もしそれぞれにない個体があった場合に、それをすべて含めたければ all=TRUE とする
>icci <- merge(icci1, icci2, all=TRUE)

データフレームのマージ(行の結合)

> total <- rbind(A1_sample, A2_sample, B1_sample)
# rbind で2つ以上のデータフレームをくっつける。カラム変数は同じでないとダメ。

データフレームからランダムに X 件を抽出

> A1_sample <- icci_A1[sample(1:nrow(icci_A1), 1000, replace=FALSE),]
# icci_A1 の行データ(1から個体数すべて)のなかから1000件をランダムにサンプルして A1_sample に格納

分布

  • Density or point probability 密度確率・
    • the probability density function 確率密度関数
    • the probability distribution function 確率分布関数
  • 確率分布関数=確率密度関数を積分したもの
  • 確率密度関数=確率分布を微分したもの
    • 統計データが十分にたくさんあれば、
    • ヒストグラム = 確率密度関数になり
    • 累積頻度分布 = 確率分布関数になる。

正規曲線の書き方1

> 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 軸に縦棒が伸びる形式になる
	
	

正規曲線の書き方2

> curve(dnorm(x), from=-4, to=4)
これでも同じカーブを描くことができる

累積分布関数

  • 正規分布のどの位置にあるか、その位置までの cumulative distribution を出すことで、その値の正規分布での外れ具合(有意確率)を見る
  • 平均132、標準偏差が13の正規分布の中で、160という数字よりも小さい値に入る
    • 正規分布上の値のとる確率は何パーセントか、を pnorm が計算
    • それを 1 から引くことで、160以上の値をとる確率を得ている
> 1-pnorm(160,mean=132,sd=13)
[1] 0.01562612
> 1-pbinom(15,size=20,prob=.5)
[1] 0.005908966

記述統計

  • サンプルデータとしてランダムに50個の正規分布を成す数値を発生させる
> 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 に2つ目の argument で何%で切るかを 0.00-1.00 範囲で指定する →この例は 0.1 = 10% ずつ

記述統計を一発で出す方法

  • summary 関数は 最小値、第1分位点、メジアン、平均、第3分位点、最大値、欠損数を出す。

 > 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

ヒストグラムを書く

> 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)

変数:数値と文字列の操作

  • たとえば energy というデータは expend, stature という2変数をもち、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 の1行だけとりだす
> 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

データ・フレームの一部取り出し

(例)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

標準化得点を出す方法

> 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

さまざまなプロット

  • 棒グラフ
    • ratings 中の1つの変数 Length に関して crosstab を作り、棒グラフにする、グラフのラベルや色を指定する
>barplot(xtabs( ~ ratings$Length), xlab="word length", col="grey")
  • ヒストグラム
    • ratings 中の1つの変数 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 で処理しておく
  • 上位階層を変数でくくって箱髭図を書く

複数の表を並べて作成する場合

  • par()関数を使うと、以後作成する表を1つ1つ並べて表示してくれる
>par(mfrow = c(3,2))     
  • 2列3行で表を並べてくれる

表を画像にはき出す場合

  • jpeg で barplot.jpeg というファイルに書き出す
>jpeg("barplot.jpeg", width = 400, height=420)
>truehist ...   # ここで作表を行う
>dev.off()      # 画像モードを抜ける

単語リストの操作

  • テキスト・ファイルを読み込む
(例)c:\temp にある lemma.num を bnclemma というデータフレームにする:

> bnclemma <- read.table("c:/temp/lemma.num", header=T)

	header=T で1行目をヘッダとして読み込む
  • データフレームを編集する
  • 作成したデータフレームをスプレッドシートで編集したい場合:
> lm <- edit(bnclemma)
  • 編集したデータフレームは lm という名前(何でもよい)になって保存される
  • lemmalist として lm というデータフレームを作成した
  • 単語リストの抽出
  • データフレームを行列として参照する
  • ランク10位の単語([,3]はwordのカラムを示す)を参照
> lm[10,3]
[1] to
  • ランク10位の品詞([,4])を参照
> lm[10,4]
[1] prep
Levels: a adv conj det infinitive-marker interjection modal n prep pron v

※このようにアウトプットの最後に、その変数に何レベルあるかが列挙される
  • lm 中の word のカラムの最初の4個を取り出す
> lm$word[1:4]
[1] the be  of  and
5464 Levels: a abandon abbey ability able abnormal abolish abolition abortion about ... zone
  • lm 中の pos のカラムの1番目を取り出す
> 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 が100位以内の 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 が100位以内で、pos が名詞のものを取り出す
> lm$word[lm$rank < 100 & lm$pos == "n"]
[1] time   year   people way   

※ [ ] で条件設定し、& で2つの条件を組み合わせている
※ logical operator としては、&, |(=or), !(=not) が使える
  • rank が100位以内で、pos が名詞または形容詞のものを取り出す
> lm$word[lm$rank < 100 & (lm$pos == "n"|lm$pos == "a")]
[1] time   year   last   other  people new    way    good  

※ logical operator | が使われていることに注意
  • 10位以内のリストをすべてとりだす
> 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 
※10位まで取り出すならば、<11 とする
※[]の中で条件を付けるが、最後の ,] とカンマをつけるところが重要。
  • 100位以内で品詞が動詞である行だけをとりだす
> 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
  • 100位以内の単語の頻度プロットを作る
> lm100 <- lm[(lm$rank <101),]		# lm100 というデータフレームに100以内を格納
> 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]
  • 100語以内の動詞のプロットを作る
> vlm <- lm100[lm100$pos == "v",]		#動詞を100位以内から抽出 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 で少し大きめに 

‐サブコーパス間のランク差の統計情報を得る

# 複数の頻度表を、見出し語・ランク・頻度・コーパスの4コラムで整理し brownclone というデータフレームに格納した場合。numerical summary の関数で見出し語ごとの各コーパスのランクの記述統計でテーブルを作成、それを CSV に書き込むという操作
sum <- numSummary(brownclone[,"Rank"], groups=brownclone$lemma, 
 statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.25,.5,.75,1))
write.csv(sum$table, file="numsummary_brownclone.csv")

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2013-12-18 (水) 10:11:02 (1404d)