R_fromOldHtml2_2 - RとLinuxと...

RとLinuxと...


R_fromOldHtml2_2

個人メモを整理中です.不適切な記述が多々あるかと思います.お気づきの際は,ishida-m(この部分を"@"に変更下さい)ias.tokushima-u.ac.jp までご連絡ください.

R_fromOldHtml2

_ バートレットの多群の等分散性 bartlett.test 2006 04 01

永田靖著『統計的多重比較法の基礎』サイエンティスト社

 bartlett.test(p.37$x ~ p.37$name)

_ 群ごとに分散などを求める 2006 04 01

nagata.aov.R 永田靖著『統計的多重比較法の基礎』サイエンティスト社

 tapply(p.37$x, list(p.37$name), var)

_ Dunnett の方法による多重比較 2006 04 01

これは対象群を中心に,それと処理群の平均値の差を調べる

永田 靖 p.40

p.43<- data.frame(x = c(7,9,8,6,9,8,11,10,8,8,
                 8,9,10,8,9,9,10,12,
                 11,12,12,10,11,13,9,10,
                 13,12,12,11,14,12,11,10),
                 name = (rep(c("A1","A2","A3","A4"), 
                 c(10,8,8,8)))
              )

mulcomp ライブラリを使う

 library(multcomp)
 # p.43 のデータ形式のままだと A1 を引く形になる。
 # そこでコントラストを変える
 c.m<- matrix(c(1,-1,0,0,
                 1,0,-1,0,
                 1,0,0,-1), nrow = 3, byrow = T)
 tapply(p.43$x, list(p.43$name), mean)
 # 8.400 - 9.375 # -0.975
 summary(simint(x ~ name, data = p.43, cmatrix = c.m,  
   alternative = "less"))

_ ロジット とは 2006 04 03

このサイトを参考にさせていただきました.以下引用です!(石田).

 対数オッズは,通常「ロジット logit」と呼ばれる.
 log[p/(1‐p)]=βx+α 
 式を書き直すと,
 p=[1+exp{-(βx+α)}]-1
 と表わすことができる.
 2群のロジットの差をとると,
   logit(p1)‐logit(p2)=log[p1(1-p2)/[p2(1-p1)]=β
 となる.
 これは対数オッズ比log(ψ)がβであることを示している.
  すなわち,ψ=exp(β) である.

Faraway p.27, p.35

_ アルファベット順以外のカテゴリ分け関数 gl 2006 04 05

Faraway p.70, p.101

カテゴリの基準をアルファベット順以外で構成する時などに使う

 y <- c(320, 14, 80, 36)
 particle <- gl(2, 1, 2*2, labels = c("no","yes"))
 quality <- gl(2, 2, labels = c("good", "bad"))
 wafer <- data.frame(y, particle, quality)

_ GLM.summary の見方 2006 04 05

Faraway p. 70

 > wafer
     y particle quality
 1 320       no    good
 2  14      yes    good
 3  80       no     bad
 4  36      yes     bad

 Call:
 glm(formula = y ~ particle + quality, family = poisson)

 Deviance Residuals: 
      1       2       3       4  
  1.324  -4.350  -2.370   5.266  

 Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
 (Intercept)   5.6934     0.0572  99.535  <2e-16 ***
 particleyes  -2.0794     0.1500 -13.863  <2e-16 ***
 qualitybad   -1.0575     0.1078  -9.813  <2e-16 ***

 > predict(modl)
        1        2        3        4 
 5.693358 3.613916 4.635807 2.556366

 > model.matrix(modl)
   (Intercept) particleyes qualitybad
 1           1           0          0
 2           1           1          0
 3           1           0          1
 4           1           1          1
 attr(,"assign")
 [1] 0 1 2
 attr(,"contrasts")
 attr(,"contrasts")$particle
 [1] "contr.treatment"

 attr(,"contrasts")$quality
 [1] "contr.treatment"

解釈は以下のように,ここで bad ではなく good が基準となっているのは gl 関数による

 particle no,   quality good      5.693358
 particle yes,  quality good      5.6934 - 2.0794(x1) 

# この行に対応する Pr(z)は<2e-16 は particle no との差が有意なのを表す

 particle no,   quality bad       5.6934 - 1.0575(x2) /

# この行に対応する Pr(z)は<2e-16 は quality good との差が有意なのを表す

 particle yes,  quality bad       
 5.6934 - 2.0794(x1) - 1.0575(x2)

 > predict(modl,type = "response")
         1         2         3         4 
 296.88889  37.11111 103.11111  12.88889 
 ---

_ outer 同時確率を出す 2006 04 05

Faraway p.72

 > wafer
     y particle quality
 1 320       no    good
 2  14      yes    good
 3  80       no     bad
 4  36      yes     bad

 > outer(qp,pp)
        particle
 quality        no        yes
    good 0.6597531 0.08246914
    bad  0.2291358 0.02864198

_ deviance を Pearson X^2 で測る, すなわち chisq の別の方法 2006 04 05

Faraway p.73, p.76

 > wafer
     y particle quality
 1 320       no    good
 2  14      yes    good
 3  80       no     bad
 4  36      yes     bad

xtabs を使う.

 (ov<- xtabs(y ~ quality + particle))
 prop.test(ov) 

_ 特異値分解 2006 04 05

Faraway p.77

関数 svd の対象は residuals

_ plot 関数への asp 引数 2006 04 05

このサイトを参考にしました.

 Note that if asp is a finite positive value then 
 the window is set up 
 so that one data unit in the x direction is equal 
 in length to asp * one 
 data unit in the y direction. 

_ pchisq への lower = F 引数 2006 04 06

Faraway p.80

 pchisq(deviance(mods), df.residual(mods), lower =F)
 # [1] 0.003762852
 pchisq(deviance(mods), df.residual(mods))
 # [1] 0.9962371

_ margin.table の使い方 2006 04 06

 < ct
         left
 right    best second third worst
   best   1520    266   124    66
   second  234   1512   432    78
   third   117    362  1772   205
   worst    36     82   179   492

 margin.table(ct, 1)
 right
   best second  third  worst 
   1976   2256   2456    789 

_ テーブルをデータフレームに変換 2006 04 06

Faraway p.89

 data(nes96)
 # xtabs  使う
 xtabs(~ PID + educ, nes96)
 (partyed as.data.frame.table(xtabs( ~ PID + educ, nes96))) 

_ xtabs の使い方 2006 04 06

 xtabs(target.values ~ categorical.var1 + categorical.var2, 
      data=data.name)

あるいは

 > xtabs( ~ Q1 + Q2, data = mydata[,2:3])
    Q2
 Q1  1 2 3 4 5
   1 3 2 3 2 0
   2 0 0 2 0 0
   3 0 0 1 1 1
   4 0 0 0 0 1

_ cut 関数の使い方 2006 04 06

Faraway p.98

cutinc さらに ここ も参照

_ log による summary 出力を antilog (exp) 化して解釈する 2006 04 07

Faraway p.105, p. 34 に簡易的な見方

 > summary(binmodw)

 Call:
 glm(formula = cbind(CNS, NoCNS) ~ Water + Work, 
    family = binomial,  data = cns)
 Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
 (Intercept)   -4.4325803  0.0897889 -49.367 < 2e-16 ***
 Water         -0.0032644  0.0009684  -3.371 0.000749 ***
 WorkNonManual -0.3390577  0.0970943  -3.492 0.000479 ***

ここで WorkNonManual? の Coefficient -0.3390577 は次のように解釈できる。 log ベースで基準より -0.339 少ないとは,antilog に戻すと,

 exp(-0.3390577) 0.7124413

であり,これは 基準 * 0.7124413, すなわち基準より約 29%少ないと解釈できる。

_ P-value と null distribution 2006 04 08

P-value と null distribution

以下はこのサイトからの引用です! 石田

帰無仮説どおりのランダムサンプリングによって得られたデータの分布を、 null distribution と呼ぶ。例えば、 T-検定を使用するような状況で t-統計量の null distribution は t-分布に従うし、chi^2-検定を使用するような状況で 統計量の null distribution は chi^2-分布に従う。

ところが、どんな検定でも検定統計量のP-valueを計算すると、 そのnull distribution は[1,0]の一様分布となる。 逆に、P-value はそうなるように定義されていると解釈してもよい。 つまり、 「P-value とは、その値が小さければ小さいほど仮説の有意性が高いと解釈できるような指標であって、帰無仮説(null hypothesis)が真であるときに、その値の出現頻度が[1,0]で一様分布するもの」である。さらに他の言葉で言えば、「P-value は有意性指標を正規化したもの」である。 どんな帰無仮説に基いてどんな検定統計量を用いるどんな検定でも、 P-valueが同じである仮説同士は同じぐらい有意 であると言うことができる。

_ Wald test と,その問題点 2006 04 08

Faraway p.30, p.122 - 123

_ データを中心化 (center) する 2006 04 15

Faraway p.186

_ iid 独立同型分布 2006 04 15

iid は独立同型分布(independently and identically distributed)

_ expand.grid 変数の値を変化させる(ただし他変数は固定) 2006 04 18

Faraway p.273

 expand.grid(temp = seq(-3, 3, 0.1), ibh = 0, ibt = 0)

ibh, ibt の値は平均(標準化された平均値)に固定し,temp を[-3,3] の間で変化

_ diagnostic.plot の見方 2006 04 19

Wood, GAM, p.5 下 - p.6

_ 任意のヒストグラムの描き方 2006 04 27

例えば平均値が 64.68, 標準偏差が 13.97 のヒストグラムであれば

 x <- c(64.58-(13.97*3) : 64.58+(13.97*3))
 y <- dnorm(x, 64.68, 13.97)
 plot(x,y, type="l")
 lines(c(64.68, 64.68), c(max(y), min(y)))
 lines(c(64.68 - 13.97, 64.68 - 13.97), c(0, .0169))
 lines(c(64.68 + 13.97, 64.68 + 13.97), c(0, .0169))
 text(66, -0.0005, "64.68")
 lines(c(64.68 - 13.97, 64.68), c(.0169, .0169), 
    lty = 3)
 text(58, .0159, "13.97")
 lines(c(min(x), max(x)), c(0,0))

_ lm 関数出力の F 値 2006 04 28

Wood, p.30

推測されたパラメータが,切片しか含まないモデルから生成された確率。

_ lm 関数出力の R^2 がマイナスのケース 2006 04 28

Wood, GLM, p.33

モデルに不必要な説明変数が多く存在する

_ A % in % B による条件判定とデータ抽出 2006 04 28

Wood, GLM, p.35

 sperm.comp1$m.vol <- 
     sperm.comp2$m.vol
             [sperm.comp2$pair %in% sperm.comp1$subject]

他に Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.165

 offset の意味 2006 04 28  
 Wood, GLM, p.45
 モデルの含める際,そのパラメータが 1 と仮定する。

_ layout 関数によるプロット座標の分割 2006 05 04

 attach(iris)
 layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE), 
  c(2,1), c(1,2), TRUE)
 layout.show(4)
 for(i in 1:4){
   boxplot(iris[, i] ~ Species, main = colnames(iris)[i])
 }

Murrell, p.78

# layout を使うなら
 layout(rbind(c(1,2),c(3,4)))
# layout.show(4)
# layout(matrix(c(2,0,1,3), 2, 2, 
      byrow = TRUE), c(2,1), c(1,2), TRUE)
 for(i in 1:4){
   boxplot(iris[, i] ~ Species,  
      notch = TRUE, main = colnames(iris)[i])
 }

Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.31 上

Wood, p.36

 layout(matrix(c(2,0,1,3), 2, 2, byrow = TRUE), 
     c(2,1), c(1,2), TRUE)

はプロットの横軸に関しては,左が全体の 2/3 右が 1/3 のサイズとし 縦軸に関しては,上が全体の 1/3 下が 2/3 とする設定である。

 > matrix(c(2,0,1,3), 2, 2, byrow = TRUE)
      [,1] [,2]
 [1,]    2    0
 [2,]    1    3
  pairs と panel の使い方 2006 05 06 

Everitt & Hothorn, p.68

 pairs(means[,-1], panel = function(x, y)
 {text(x,y, abbreviate(levels(skulls$epoch)))})

_ substitute() の使い方 2006 05 08

substitute のヘルプより

 f1<- function(x, y = x)             { x<- x + 1; y }
 s1<- function(x, y = substitute(x)) { x<- x + 1; y }
 s2<- function(x, y) { if(missing(y)) 
 y<- substitute(x); x<- x + 1; y }

 a<- 10
 f1(a)# 11
 s1(a)# 11
 s2(a)# a

 > x<- 1
 > substitute(x+1) 
 x + 1
 > substitute(x+1, list(x = 2))
 # x+1 の変数を指定された環境(ここではリスト)で置き換える
 2 + 1

他に Crawley S-PLUS p. 831

_ bubble plot の効果的な使い方 2006 05 12

Everitt & Hothorn, p.98

_ multinom 関数の使い方 2006 05 17

 rmultinom(n, size, prob)

 rmultinom(10, size = 15, prob=c(0.9,0.9))
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    7    6    6    8    8    7    6    8    4     8
 [2,]    8    9    9    7    7    8    9    7   11     7

10 は試行の回数, 15 は玉の最大個数 prob は箱の個数と,それぞれに入る玉の個数,ここでは箱は2個。

_ glm 出力から Null deviance Residual deviance を抽出 2006 07 19

Faraway p.41

 modl<- glm(cbind(dead, alive) ~ conc, 
          family = binomial, bliss)
 summary(modl)

 modl$dev
 modl$null

とする.

_ ポアソン分布 2006 07 21

工学のためのデータサイエンス, p.47 - 48

 mase.47<- c(rep(0,55), rep(1,72), rep(2,40), rep(3,9), 
     rep(4,5), rep(5,1))
 mean(mase.47) #= (72 + 40 *2 + 3 * 9 + 4 *5 + 5) / 181
 var(mase.47)

村上征勝先生の「工業統計学」, p.44

 ppois(0,2)
 ppois(1,2) - ppois(0,2)

_ table 関数で作成したカウントのカテゴリ頻度 0 の調整 2006 07 25

 # 高頻度では、データ区間が飛び飛びになっているので、これを補正する
 y<- min(buntyo):max(buntyo)
 y<- 0:max(buntyo)
 bun.df<- data.frame(cate = y, freq = c(rep(0, length(y))))
 z<- 0
 bun.df[1,] = c(0,0)
 for(i in 1:nrow(bun.df)){
   if(bun.orig.df[i - z,]$buntyo == bun.df[i,]$cate){
     bun.df[i,]$freq <-
       bun.orig.df[bun.orig.df$buntyo == 
        bun.df[i,]$cate, ]$Freq}
   else{
     z<- z + 1;
     next;
   }
 }
 ##    

_ カイ二乗検定での期待度数 2006 07 27

青木先生の説明を参照

http://aoki2.si.gunma-u.ac.jp/lecture/Cross/warning.html

以下のいずれかに該当する分割表に対して χ2 検定を行ったり, カイ二乗値から導き出される関連性の指標を用いる場合には注意が必要である (Cochran, W. G.: Some methods for strengthening the common χ2 tests. Biometrics, 10, 417-451, 1954.)。

_ ks.test を使った適合度の検定 2006 07 27

 x <- rnbinom(30, size = 5, prob = 0.5)
 ks.test(x,"pnbinom", size = 5, prob = 0.5)

 # 間瀬先生著 p.51 下
 n<- 4+9+16+13+9+7+5+4+3
 x<- rnbinom(n, size = 10.59, prob = 0.7609)
 ks.test(x, "pnbinom",  size = 10.59, prob = 0.7609)
 table(x)
 ks.test(table(x), "pnbinom",  size = 10.59, prob = 0.7609)
 ks.test(observed, "pnbinom", size = my.k,  
         prob = (my.k / (my.k + mean(x))))

ここで size は 2項分布のパラメータ k (= size) prob は prob<- size/(size + mu) で計算される.

 ks.test(rnorm(100), "pnorm")
 ks.test(rpois(100, 5), "ppois", lambda = 5)

_ glm を使った適合度の検定 2006 07 29

null model を当てはめてみる

 # データは負の二項分布に従っているか
 x<- rnbinom(100, size = 2, prob=0.5)
 x.glm.nb<- glm.nb(x ~ 1)
 x.glm<- glm(x ~ 1, family = negative.binomial(1.57))
 1 - pchisq( x.glm$deviance, x.glm$df.residual)
 [1] 0.2081144

ポアソン分布ではどうか

 x.glm.poi<- glm(x ~ 1, family = poisson)
 1 - pchisq( x.glm.poi$deviance, x.glm.poi$df.residual)
 [1] 4.327188e-09 # poisson による当てはめは悪い

他にも実験

 # ポアソン分布を負の二項分布で当てはめる
 x<-  rpois(100, 3)
 x.glm.nb<- glm.nb(x ~ 1)
 summary( x.glm.nb)
 1 - pchisq( x.glm.nb$deviance, x.glm.nb$df.residual)
 # [1] 0.9046653

正規分布を負の二項分布で当てはめる

 x<-  rnorm(100, mean = 10, sd = 2)
 x.glm.nb<- glm.nb(x ~ 1)
 summary( x.glm.nb)
 1 - pchisq( x.glm.nb$deviance, x.glm.nb$df.residual)
 # [1] 1 当てはまってしまう

ここも

Crawley S-PLUS p.214, p.511, p.517, p.539 も

Faraway p.120 (さらに 122 ページまで説明続く).

このような記事もあった また次のような記事もあった

  From: Dan Kehler (kehler@mathstat.dal.ca)
 Date: Sat 29 May 2004 - 02:54:48 EST

* Next message: Richard Valliant: "[R] dotchart questions"
* Previous message: Ted Harding: 
   "RE: [R] Generate a sequence of random integer values" 

 Message-id:<Pine.GSO.3.96.1040528133631.18220B-100000@chase>
 Using R 1.8.1, and the negative binomial glm 
 implemented in MASS,
 the default when using anova and a chi-square test is 
 to divide the   deviance by the estimated dispersion. 
 Using my UNIX version  of S-plus (v  3.4), 
 and the same MASS functions, the deviances are *not* 
 divided by the
 estimated dispersion.

 Firstly, I'm wondering if anyone can enlighten about 
 the correct procedure   (I thought the F-test was more 
 appropriate  when dispersion is estimated)?

 Secondly, after a bit of muddling with the negative 
 binomial pdf, 
 I concluded that, like for the Poisson, phi is actually 1. 
 This result is   borne out by simulations. Is this correct?

 # an example in R 1.81 with library(MASS)
  y<-rnegbin(n=100,mu=1,theta=1)
  x<-1:length(y) 
 # 石田 追加,これは便宜上の説明変数を作ったに過ぎない
  model<-glm(y~x,family=neg.bin(1))
  summary(model)$dispersion
  [1] 1.288926

  anova(model,test='Chisq")
 #...
      Df Deviance Resid. Df Resid. Dev P(>|Chi|)
 NULL 99 102.038
 x 1 0.185 98 101.853 0.705

 # But the "real" chi-square probability is

   1-pchisq(0.185,1)
 [1] 0.6671111

 1-pchisq(0.185/1.288926, 1) # 石田 追加
 [1] 0.7047963               # 石田 追加

 Thanks in advance,

 Dan

これは,フルモデルにたいして,ヌルモデルを比べ,観測値の変動をヌルモデルだけで どれだけ説明できるかを示したものである. フルモデルは,本来,deviance がゼロのはずであるので,ヌルモデルとフルモデルの デヴィアンスの差とは,ヌルモデルのデヴィアンスのことではないのか.

また

  x<-1:length(y)

は単変量であって,これを使って回帰を行ってもフルモデルを設定したことには ならないのではないか.そもそもフルモデルは自由度 0 のはずである. その意味で,ヌルモデルそのものがある分布に適しているかどうかを調べるには 先の方法でも十分かも知れない. 下の場合,もともと負の二項分布にのみ従った乱数で,他に変動要因が無い. 尤も,この場合,誤差は負の二項分布の誤差に従っているということになるのか?

 y<- 1:length(x)
 x.glm.nb2<- glm.nb(x ~ y)
 anova(x.glm.nb2, test = "Chi")
 これは,以下に等しい
 pchisq( x.glm.nb$null.deviance -x.glm.nb$deviance , 
         x.glm.nb$df.null - x.glm.nb$df.residual,
         lower = F)
 anova(x.glm.nb, x.glm.nb2, test = "Chi")

以下の数値とは異なる数値が出る

 x.glm.nb2$null.deviance # NULL モデルの deviance
 x.glm.nb$df.null        # NULL モデルの自由度
 pchisq( x.glm.nb2$null.deviance , x.glm.nb$df.null, 
    lower = F)

_ residual の種類,data$residuals と residual関数は working 2006 07 31

Faraway p. 124

response, pearson, working とあるが, 解析結果の $residuals で出力されるのは working その他の residual を得るには関数 residual を使う.

_ 警告メッセージを抑制する 2006 08 01

suppressWarnings 関数の引数として,命令を実行する.

_ 分布の適合度を一般化線形モデルで行う 2006 08 02

Faraway p.29 -30. またカイ自乗分布は,平均 = 自由度,分散 = 2*自由度に従う.

_ 飽和モデルの自由度は 0 2006 08 02

Faraway p.33

_ predicted と fitted の違い 2006 08 02

Faraway p.36

_ expression の使い方 2006 08 18

描画を行う

 si.op1<- sichel.opti$par[1]
 si.op2<- sichel.opti$par[2]

 text.main<- paste("Sichel's Compound Poisson Distribution\n", text.main)
 barplot(both, col = rep(c("white", "green"), length(kekka)), 
   names.arg = cate.label, ylab = "Frequency", 
   xlab = "Cases", 
   main = text.main )

 legend(legend.x - 25, legend.y, c("Observed", "Expected"), 
    fill = c("white", "green"))
 # legend(legend.x + 15, legend.y - 8, 
    c(expression(paste(alpha, sichel.opti$par[1], beta,
       sichel.opti$par[2], chi^2, " P (freq > 1) =" ))),
    chi.score1))
 legend(legend.x, legend.y - 8, 
        c(expression(paste(alpha, " =  ")), si.op1,
          expression(paste(beta, " = "))),  si.op2, 
          expression(paste(chi^2, " P (freq > 1) =" ))),
          chi.score1))

_ dnorm を使った期待値の計算 2006 08 19

正確に行うには 青木先生のサイト

http://aoki2.si.gunma-u.ac.jp/lecture/GoodnessOfFitness/normaldist.html http://aoki2.si.gunma-u.ac.jp/R/fit_normal.html

にあるように,各階級の級中心の確率を求め,その階級の期待度数は,上の級の 確率との差にデータ数を乗じることで求まる.

Crawley, Statistics using R, p.56

 means<- numeric(10000)
 for(i in 1:10000){
   means[i]<- mean(runif(5) * 10)
 }

 # sum(means)# 49946.35
 length(means)
 hist.breaks<- hist(means, ylim=c(0,1600))$breaks
 hist.breaks
 # 1  から 10 の範囲 10個を 2 倍して,0.5, 1, ... としている.
 # 従って,length(means) の半分を dnorm の計算結果に乗じる
  hist(means, ylim=c(0,1600))
 yv lines(xv, yv)
 ##################### これを形式的に行う
 length(hist.breaks)

 length(dnorm(hist.breaks, mean = mean(means), 
     sd = sd(means)))

 (dnorm.test<- sum(dnorm(hist.breaks, mean = mean(means), 
   sd = sd(means))) * length(means))

 sum(means)# mean(means) * length(means)
 # 従って 
 # length(dnorm(hist.breaks, mean = mean(means), 
  sd = sd(means)))
 sum(means)#
 sum(dnorm.test<- sum(dnorm(hist.breaks, 
   mean = mean(means), 
     sd = sd(means))) * length(means))
 # sum(dnorm.test <- sum(dnorm(hist.breaks, 
       mean = mean(means), 
      sd = sd(means))) * sum(means))
 test.yv <- dnorm(hist.breaks, mean = mean(means), 
      sd = sd(means)) * length(means)  / 
                 sum (dnorm(hist.breaks, 
                 mean = mean(means), sd = sd(means)) * 
                 length(means)) / sum(means)#

 test.yv<- dnorm(hist.breaks, mean = mean(means), 
            sd = sd(means)) * length(means)  *  
            sum (dnorm(hist.breaks, mean = mean(means), 
            sd = sd(means)) * length(means)) / sum(means)#
 hist(means, ylim=c(0,1600))
 lines(hist.breaks, test.yv)

# 実測値
 # log.bread<- seq(min(log.V1) - 0.5, max(log.V1) 
     + 0.5, 0.25)
 log.bread<- seq(min(log.V1), max(log.V1)+1, 0.25)

 log.V1.ob<- as.vector(table(cut(log.V1, 
     breaks = log.bread)))
 length(log.V1.ob)
 sum(log.V1.ob)
 # 期待値
 log.V1.dnorm<- dnorm(log.bread[-(length(log.bread))],  
  mean = mean(log.V1), sd = sd(log.V1))

 log.V1.exp<- log.V1.dnorm * length(log.V1) / 
      (sum(log.V1.dnorm * length(log.V1)) / sum(log.V1.ob) )
 sum(log.V1.exp)

_ merge 関数でデータフレームに要素を追加 2006 09 16

  test<- data.frame(x = 1, y = 2)
 > test
   x y
 1 1 2

 > merge(test, data.frame(x = 3, y = 4), all = T)
   x y
 1 1 2
 2 3 4

_ リストの中身を全て取り出す方法 2006 09 16

 ################ リストの中身を全て取り出す
 x.test<- 1:5
 y.test<- 6:8
 z.test<- 9:12
 xyz <- list(x.test, y.test, z.test)

 for(i in 1:length(xyz)){
   for(j in 1:length(xyz[[i]])){
     print (xyz[[i]][j])
   }
 }

_ if()else{} でelse で始まる文は認められない 2006 09 18

 if(bun.p >= 4){
   ex.freq<- 5
 }else{
   ex.freq<- 1
 }

は良いが

 if(bun.p >= 4){
   ex.freq<- 5
 }else{
   ex.freq<- 1
 }

はエラーとなる.

_ exists 関数でオブジェクトを確認 2006 09 18

 if(!exists("myobject"))

_ 自作プログラム,観測データとその期待値の区間を調整して,頻度をはかる 2006 09 18

 ~/daigaku/GakubuKeihi/bun.cut.R より
 ########## 頻度表の区間をまとめる
 # オリジナル
 #
 bun.df0<- bun.df[-1,]# ゼロ頻度をけずる
 #
 length(bun.df0$freq)
 # bun.p<- 3

 ########### 指定された区間でデータをまとめる
 # bun.p が 1 ならば,元データと同じこと
 bun.cut<- data.frame(cate = 1, freq = 0)
 xx<- 0
 for(i in 1:length(bun.df0$freq)){
   xx<- xx + bun.df0[i,2]
   z1<- ceiling(i / bun.p)
   z2<- i %% bun.p
   if(z2 == 0) {
     bun.cut[z1,1]<- z1
     bun.cut[z1,2]<- xx
     xx<- 0
   }
   if((z2 != 0) && (i == length(bun.df0$freq))){
     bun.cut[z1,1]<- z1
     bun.cut[z1,2]<- xx
   } 
 }
 length(bun.cut$freq)

 ##########################
 # 期待値 の頻度を指定された区間でまとめる
 #
 theo.cut<- data.frame(cate = 1, freq = 0)
 xx<- 0
 for(i in 1:length(theo.bun.0)){
   xx<- xx + theo.bun.0[i]
   z1<- ceiling(i / bun.p)
   z2<- i %% bun.p
   if(z2 == 0) {
     theo.cut[z1,1]<- z1
     theo.cut[z1,2]<- xx
     xx<- 0
   }
   if((z2 != 0) && (i == length(bun.df0$freq))){
     theo.cut[z1,1]<- z1
     theo.cut[z1,2]<- xx
   } 
 }
 ####################################
 # 期待値が ex.freq 未満のセルをまとめる
 # ex.freq<- 5
 theo.under<- data.frame(cate = 0, freq = 0)
 un<- 1
 xx<- 0
 yy<- numeric(0)
 zz<- list(0)
 for(i in 1: length(theo.cut$freq)){
   if(theo.cut[i,2]< ex.freq){
     yy<- c(yy,i)# ex.freq 以下の行番号を記録
     # print(yy)
     xx<- xx + theo.cut[i,2]
     if(xx >= ex.freq){
       zz[[un]]<- yy #マージが実行された行を記録
       yy<- numeric(0)
       theo.under[un,1]<- i #マージ実行時のカテゴリを記録
       theo.under[un,2]<- xx ##マージによる融合合計頻度
       un<- un + 1
       xx<- 0
     }
   }
   else {
     if((xx != 0) && xx< ex.freq){
            # 頻度は5以上だが,前のセルが残っている
       zz[[un]]<- c(yy,i) #マージが実行された行を記録
       yy<- numeric(0)
       theo.under[un,1]<- i #マージ実行時のカテゴリを記録
       theo.under[un,2]<- xx + theo.cut[i,2] 
              #マージによる融合合計頻度
       un<- un + 1
       xx<- 0
     }
   }
  ###### 末尾に達した
   if(i ==  length(theo.cut$freq) ){
     if(xx > 0){# そして最後の行のセルが記録されていない
       if(un == 1) {# ここまで指定頻度以下のセルが無かった場合
         theo.under[un,1]<- i
         theo.under[un,2]<- xx
         zz[[un]]<- i
       }else{
       # 5 以下の頻度であるが,
       # そのまま出力する#最後は5以下でもよしとする
         theo.under[un,1]<- i
         theo.under[un,2]<- xx
         zz[[un]]<- yy
       }
     }
   }
 }
 ###併合結果をマージする
 if(zz[[1]][1] != 0){#併合が行われているならば
   zz.vec<- unlist(zz)
   theo.target<- theo.cut[-zz.vec,]
   bun.target<- bun.cut[-zz.vec,]
   if(nrow(theo.target) == 0){#全部が併合対象になっている場合
       theo.target<- theo.cut
       bun.target<- bun.cut
   }
   #######
   for(i in 1:length(zz)){
     theo.target<- merge(theo.target, theo.under[i,], all = T)
     bun.target<- merge(bun.target, 
         data.frame(cate = theo.under[i,1],
	  freq = sum(bun.cut[zz[[i]],2])), all = T)
   }
 }else{
   theo.target<- theo.cut
   bun.target<- bun.cut
 }
 ##### 以下三つの出力 (文の数) は等しくなければならない
 sum(theo.target$freq)
 sum(bun.target$freq)
 length(buntyo$V1)

_ 0-負の二項分布のパラメータ推定 2006 09 29

 source("research/statistics/wyshak.R")

あらかじめ次の計算をしておくこと なお nb0.table ではゼロ頻度のカテゴリも登録しておくこと

 (nb0.mean<- mean(nb0.data)) # mean(nb0.table) ではない
 (nb0.sum<- sum(nb0.table))  # sum(nb0.data)   ではない 
 (nb0.var<- var(nb0.data))
 (nb0.ratio<- nb0.table[1] / sum(nb0.table))#length(nb0.data)

 (w0<- nb0.mean / nb0.var * ( 1- nb0.ratio))
 (k0<- (w0 *  nb0.mean  - nb0.ratio) / (1 - w0))

author = {G. Wyshak}, title = {A Program for Estimating the Parameters of the Truncated Negative Binomial Distribution}, journal = {Journal of the Royal Statistical Society. Series C, Applied statistics},

 ## start

 ## Brass のデータ
 ## start

 nb0.data<- c(rep(1,49), rep(2,56), rep(3,73), rep(4,41), 
       rep(5,43), rep(6,23), rep(7,18), rep(8,18), rep(9,7), 
       rep(10,7), rep(11,3), rep(12,2))

 (nb0.mean<- mean(nb0.data)) # mean(nb0.table) ではない
 (nb0.sum<- sum(nb0.table))  # sum(nb0.data)   ではない 
 (nb0.var<- var(nb0.data))
 (nb0.ratio<- nb0.table[1] / sum(nb0.table)) 
                #length(nb0.data)
 (w0<- nb0.mean / nb0.var * ( 1- nb0.ratio))
 (k0<- (w0 *  nb0.mean  - nb0.ratio) / (1 - w0))
 ################### end
 ## Sampford のデータ
 ## start
 nb0.data<-  c(rep(1,11), rep(2,6),rep(3,4), 
          rep(4,5), 6, rep(8,2), 9, 11, 13)
 (nb0.table<- table(nb0.data))

 # 0 頻度のカテゴリを登録する
 #
 nb0.df<- data.frame(nb0.table)
 y<- 1:max(nb0.data)

 nb0.df2<- data.frame(cate = y, freq = c(rep(0, length(y))))
 z<- 0
 # nb0.df2[1,] = c(0,0)
 for(i in 1:nrow(nb0.df2)){
   if(nb0.df2[i,]$cate == nb0.df[i-z,]$nb0.data){
     nb0.df2[i,]$freq<- nb0.df[i-z,]$Freq
   }
   else{
     z<- z + 1;
     next;
   }
 }
 (nb0.mean<- mean(nb0.data)) # mean(nb0.table) ではない
 (nb0.sum<- sum(nb0.df2$freq))  # sum(nb0.data)   ではない 
 (nb0.var<- var(nb0.data))
 (nb0.ratio<- nb0.df2[1,]$freq / sum(nb0.df2$freq))
                  #length(nb0.data)

 (w0<- nb0.mean / nb0.var * ( 1- nb0.ratio))
 (k0<- (w0 *  nb0.mean  - nb0.ratio) / (1 - w0))

 nb0.table<- nb0.df2$freq
############### end

 ## start
 wyshak.wk     
 nw <- (nb0.sum *z[2] * log(z[1])) - 
         (nb0.sum * log(1-z[1]^z[2]))
     r<- length(nb0.table)
         # 頻度表 途中のゼロ頻度を省略してはいけない

     nx<- 0
     for(j in 1:r){
       nx<- nx + j * nb0.table[j]# * log(1 - z[1])
     }
     ny<- 0
     for(j in 1:r){
       ny<- ny + nb0.table[j] * log(factorial(j))
     }
     nz<- 0
     for(j in 1:r){
       nz2<- 0
       for(j2 in 1:j){
         nz2<- nz2 + log(z[2] + j2 - 1)
       }
       nz<- nz + nb0.table[j] * nz2
     }
     return( nw + nx * log(1 - z[1])  - ny + nz )
   }
 #
 optim(c(w0, k0), wyshak.wk, control = list(fnscale = -1))
 # 最大化する
 system.time(optim(c(w0, k0), wyshak.wk, 
       control = list(fnscale = -1)))

_ optimize の使いかた 2006 11 08

渡部『ベイズ統計学入門』 p.72

モデル分布がベルヌーイ分布である母集団から n = 3 の標本を抽出したとき,その観測値は

 c(1,0,1)

であった.母数の最尤推定値を求めよ.

 f <- function(x) x^2 * (1-x)
 optimize(f,c(0.001, 0.999),maximum=T)

_ fitdistr の使いかた 2006 11 08

RjpWikiを参照:http:www.okada.jp.org/RWiki/index.php?R%A4%CE%C5%FD%B7%D7%B2%F2%C0%CF%B4%D8%BF%F4Tips

 library(MASS)
 mydt<- function(x, m, s, df) dt((x - m)/s, df)/s 
 # 密度関数を独自定義
 fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, 
      lower = c(-Inf, 0))
          m             s    # 2パラメータの推定値と推定標準偏差
    -0.01069635    1.04409435
   ( 0.07222562) ( 0.05434249)
 set.seed(123) # 乱数種
 x4<- rnegbin(500, mu = 5, theta = 4)
        # 負の二項分布に従う乱数500個発生
 fitdistr(x4, "Negative Binomial")  
       # 負の二項分布を当てはめ
     size         mu    # 2パラメータの推定値と推定標準偏差
  4.2178713   4.9439803
  0.5048242) (0.1465538)
  Warning messages: # 関連警告

_ poisson 分布で予測する 2006 11 08

 x<- c(0,0,0,0,0,2)
 library(MASS)
 fitdistr(x, "Poisson")
 ##     lambda  
 ##   0.3333333 
 ##  (0.2357023)
 # 0 人現れる確率
  exp(-0.3333) * 0.3333^0 / factorial(0)
 #[1] 0.7165552
 # 1 人現れる確率
  0.3333/(0+1) * 0.7165552
 # [1] 0.2388278
 # 2 人現れる確率
  0.3333/(1+2) * 0.2388278

_ 多元分割表への対数線形モデル 2006 11 23

質的情報の多変量解析 (松田 紀之 著)# p.121

 matsuda.121<- array(c(32,86, 11, 35, 61,73,41,70),  
   dim = c(2,2,2), 
   dimnames = list(height=c("over","lower"), 
   diam = c("lower","heigher"), 
   species = c("sagrei", "distichus")))

 class( matsuda.121 )<- "table"
 matsuda.121 
# p.129
 matsuda.121.loglm <- 
    loglm(~ height *  diam * species, data = matsuda.121 )
 stepAIC(matsuda.121.loglm )

 Step:  AIC= 14.03 
  ~height + diam + species + height:species + diam:species 

                  Df    AIC
 <none>              14.026
 - height:species  1 22.431
 - diam:species    1 24.632
 Call:
 loglm(formula = ~height + diam + species + 
         height:species + diam:species, 
     data = matsuda.121, evaluate = FALSE)

 Statistics:
                       X^2 df  P(> X^2)
 Likelihood Ratio 2.025647  2 0.3631921
 Pearson          2.017364  2 0.3646994

同じことだが

 matsuda.121.loglin <-
   loglin( matsuda.121, list(c(1,3), c(2,3)) )
 matsuda.121.loglin 

c(1,2) を外してよいか.つまり heigh と diam は独立だとみなし,無視してよいか 実際に外した結果を,飽和モデルと比較する 0.3631921 なので,外してよい ↓

  1 - pchisq(matsuda.121.loglin$lrt, 
   matsuda.121.loglin$df)
  1 - pchisq(matsuda.121.loglin$pearson, 
   matsuda.121.loglin$df)#

よって species さえ分かれば,height や diam の間に特に関係は無い.

_ ホームディレクトリへのパッケージのインストール 2006 12 02

次のように,デフォルト以外のライブラリの場所を確認

 Sys.getenv("R_LIBS")

必要があれば

 Sys.setenv(R_LIBS="C:/R")
 install.packages("MyPackages", lib = "C:/R") 
    # 第二引数はオプション
 library("MyPackages")

_ 日本語フォントの利用 2006 12 05

RjpWikiを参照:http:www.okada.jp.org/RWiki/index.php?%C6%FC%CB%DC%B8%EC%B2%BD%B7%C7%BC%A8%C8%C4 より

フォントの確認方法
  postscriptFonts()
  postscriptFonts()$Japan1

_ Fonts の設定 X11fonts

R のデフォルトは

options()$X11fonts
[1] "-adobe-helvetica-%s-%s-*-*-%d-*-*-*-*-*-*-*"
[2] "-adobe-symbol-medium-r-*-*-%d-*-*-*-*-*-*-*"
 Fedora core5 では以下のXLFD 情報などがきれい.
options(X11fonts = 
  c("-shinonome-gothic-%s-%s-normal--%d-*-*-*-*-*-*-*",
   "-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*-0"))

 options(X11fonts = 
 c("-shinonome-gothic-%s-%s-normal--%d-*-*-*-*-*-*-*-*",
   "-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
options(X11fonts = 
  c("-shinonome-mincho-%s-%s-*-*-%d-*-*-*-*-*-*-*",
    "-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))

ただし,座標ラベルが指定されていない場合,無意味な文字が出力される.

 x<- c("あ","い","う","え","お")
 y<- c("ア","イ","ウ","エ","オ")
 plot(1:5, 1:5, xlab = "これはテスト", type = "n", ylab = "")
 text(1:5, 1:5, x, cex = 0.7)
 legend("topright" , legend = c("日本語ひらがら"))

RWikiを参照

_ データフレームから,列名を指定して抽出 2006 12 07

例えばEveritt & Hothorn, A Handbook of Statistical Analyses Using R p.67 で

 data(skulls, package = "HSAUR")
 means<- aggregate(skulls[, c("mb", "bh", "bl", "nh")], 
                 list(epoch = skulls$epoch), mean)

# 他に

 cov(iris[, -5])    # 共分散行列.var(iris[,-5] でも同じ
 var(iris[, colnames(iris) != "Species"])  
               # 共分散行列.条件指定を変更
 cor(iris[sapply(iris, is.numeric)])       
               # 相関行列  # 条件指定を変更

_ Contrasts Matrix の意味 2006 12 08

Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.75 下

_ 外れ値を表示する方法 2006 12 08

Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.78

 data("clouds", package = "HSAUR")
 layout(matrix(1:2, nrow = 2))
 bxpseeding       xlab = "Seeding")
 bxpsecho       xlab = "Echo motion")
 rownames(clouds)[clouds$rainfall %in% 
    c(bxpseeding$out, bxpsecho$out)]

_ subset の使いかた 2006 12 08

_ legend の使いかた 2006 12 08

Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.51

 cmh_test( classification ~ treatment, data = Lanza, 
     scores = list(classification = c(0, 1, 6, 17, 30)), 
     subset  = Lanza$study == "I")
 cmh_test( classification ~ treatment, data = Lanza, 
     scores = list(classification = c(0, 1, 6, 17, 30)), 
     subset  = Lanza$study == "II")

# p.82

 psymb<- as.numeric(clouds$seeding)
 plot(rainfall ~ cloudcover , data = clouds, pch = psymb)
 abline(lm(rainfall ~ cloudcover, data = clouds, 
       subset  = seeding == "no"))
 abline(lm(rainfall ~ cloudcover, data = clouds, 
       subset  = seeding == "yes"), lty = 2)
 legend("topright" , legend = c("No seeding", "Seeding"),
 pch = 1:2, lty = 1:2, bty = "n")

_ quasilikelihood による overdispersion の考慮 2006 12 13

Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.106

また overdispersion が生じるのは,しばしばデータが独立でないため

_ mixture model 2006 12 15

Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.117

_ grep 関数の使いかた 2006 12 22

Everitt & Hothorn, A Handbook of Statistical Analyses Using R p. 166

例えば以下のコードは

 BtheB[, grep("bdi", names(BtheB))]

BtheB データフレームで,列名に "bdi" がついている列番号を抽出している.

さらに以下は

 subset(BtheB, treatment == "TAU")
              [, grep("bdi", names(BtheB))]

BtheB データで,列 treatment が "TAU" であるサブセットから, "bdi" を名前として含む列を取り出している.

_ reshape 関数の使いかた 2006 12 22

Everitt & Hothorn, A Handbook of Statistical Analyses Using R chapter 10 p. 165

 data("BtheB", package = "HSAUR")
 BtheB$subject<- factor(rownames(BtheB))
 nobs<- nrow(BtheB)
 BtheB.long<- reshape(BtheB, idvar = "subject",
   varying = c("bdi.2m", "bdi.4m", "bdi.6m", "bdi.8m"), 
    direction = "long")

# 石村貞夫「分散分析のはなし」の例.データは p.77

# データフレームの作成

 # kaeru<- read.csv("africanFlog.csv",header=F)
 # 始めに横長形式のデータとして与えられていたとする.
 frogs<- data.frame(stage = c("A1", "A2", "A3", "A4", "A5"),
                  sample1 = c(12.2, 22.2, 20.8, 26.4, 24.5),
                  sample2 = c(18.8, 20.5, 19.5, 32.6, 21.2),
                  sample3 = c(18.2, 14.6, 26.3, 31.3, 22.4))
 frog.res<- reshape(frogs, idvar = "stage",
    varying = list(names(frogs[2:4])),v.names = "data", 
       direction = "long")

# ここで作成される time 変数は,各数値の順番を表す

# (時間ごとにデータを取っている場合などに重要)

 frog2<-  frog.res[order(frog.res$stage),] # 並び替え
 frog2

_ 判別分析での適合率 2006 12 26

以下は Rjp Wiki よりの引用です! 石田

判別分析では、モデルに含まれるパラメーター数が多くなればなるほど、あてはめ誤差(残差平方和)は小さくなるということがあります。そこで、単に残差平方和の大小を比較するだけでなく、パラメータ数も考慮したモデルの選択を行います。その際の基準となるものがAICです。wleライブラリ中のmle.aic関数を使ってモデルを選択します。

library(MASS)
library(wle)
result<-mle.aic(lda(class~.,data=DATA))
summary(result,num=20)

mle.cvでやった場合にも判別率は出力されませんでした。AIC値の代わりにcv値が出力されていますが、どういう値なのかよくわかりません。 上記に続いて、判別関数を求める手順、交差妥当性(クロス・バリデーション)のチェックをやって判別率を比べる手順を示します。

(z<-lda(class~.,data=DATA))
apply(z$means%*%Z$scaling,2,mean)

この結果、それぞれの変数の係数と定数が出力されます。

predict(z)$x
predict(z)$class
predict(z)$posterior

と入力すればこの判別関数を用いた判別得点、各個体が判別されたグループのラベル、各個体がどのグループに判別されているかに関する事後確率(0から1)が出力されます。判別率を計算するためのクロス表の出力は次のようにします。

 table(data[,],predict(z)$class)

交差妥当性のチェックを行ったクロス表は

 DATA.CV<-lda(class~.,data=DATA,CV=T)
 (lda.tab<-table(DATA[,],data=DATA.CV$class))

で出力できます。判別率、誤判別率はそれぞれ

 sum(lda.tab[row(lda.tab)==col(lda.tab)])/sum(lda.tab)
 sum(lda.tab[row(lda.tab)!=col(lda.tab)])/sum(lda.tab)

交差確認(cross validation)はデータを2分し,一方のサンプルで判別式をつくりその判別式で残りのサンプルを判別して, 判別式の有効性を見る方法です。データを2分する以外に3,4,5,10分する方法もありますが、どれを用いるかは明確な基準がなくデータサイズに依存するそうです。 この特殊なケースとして1 つ取っておき法(reave-one-out cross-validation)があります。 データから1つの個体を取り除いて判別分析を行い、取り除かれた個体で判別モデルの評価を行うという作業を全ての個体に繰り返して行います。 「R」でこれが出来るということなので紹介します。ここでは「R」の中にあるirisデータを例に使います。

  library(MASS)
  iris.CV<-lda(Species~.,data=iris,CV=T)
  (lda.tab<-table(iris[,5],iris.CV$class))

デフォルトではCV=FALSEとなっているのでCV=Tにすると1つ取っておき法(reave-one- out cross-validation)ができます。 ちなみに判別関数lda(Linear Discriminant Analysis),qda(Quadratic Discriminant Analysis)には1つ取っておき法(reave-one-out cross-validation)以外の交差確認機能は用意されていないそうです。 詳しい説明は,同志社大学Jin先生のページhttp://www1.doshisha.ac.jp/~mjin/R/の「Rと判別分析」を見ると書いてあります。

1つ取っておき法とjackknife法はほとんど同義語として使われるようです。 しかし,jackknife法よりも良さそうな方法,bootstrap (ブートストラップ)法があります。 データからm個の標本を復元抽出(sampling wtih replacement)し,それらの標本から繰り返しパラメータ推定値を求める方法です。 これを線型モデル(回帰分析)で使った例が,ヴェナブルズとリプリーの「S-PLUSによる統計解析」シュプリンガー・フェアラーク東京(2001)の211〜214ぺージにあります。

bootstrap法をつかって判別分析の交差確認を実行するには, Rのipredライブラリのerrorest()関数が使えます。ipredはimproved predictorsの略。

library(ipred)
mypredict.lda<- function(object, newdata) 
    predict(object, newdata=newdata)$class
errorest(class ~ ., data=soccer, model=lda, 
    estimator="boot", predict=mypredict.lda)

mypredict.ldaはclass labelだけを返させるのに必要な手続きらしい。

_ 判別分析での変数選択 2006 12 26

  • 線型判別分析
    • klaRライブラリ中のstepclass関数を使う。mehodに線型判別分析関数のldaを与える。ステップワイズの方向は変数増減がデフォルトとなっている。
 library(klaR); stepclass(Y ~ ., data=DATA, method="lda")
  • ロジスティック回帰
    • MASSライブラリ中のstepAIC関数を使ってステップワイズ変数選択(変数増減法がデフォルト)をおこなう。引数には一般化線型モデルの関数glmの応答分布の関数に二項分布(binomial)の結果を与える(logitがリンク関数のデフォルトとなっている)。
library(MASS); stepAIC(glm(Y ~ ., data=DATA,
        family=binomial))
  

_ データフレームと文字列 2006 12 27

 test<- data.frame(v1 = c("A","B","C"), 
   v2 = c(1,2,3))
 test$v1[test$v1 == test$v1[1]]<- "AA"

これは

 Warning message:
 無効な因子水準です。NA が発生しました in: 
  `[<-.factor`(`*tmp*`, iseq, value = "AA") 

となる.そこで,

 test$v1<- as.character(test$v1)
 is.factor(test$v1); FALSE
 test$v1[test$v1 == test$v1[1]]<- "AA"

とすれば,交換はできる

 test<- data.frame(v1 = c("X","Y","Z","X","Y"), 
       v2 = c(1,2,3,4,5))
 test$v1<- as.character(test$v1)
 test.lab<- c("X","Y","Z")
 for(i in 1:length(test.lab)){
  test$v1[test$v1 == test.lab[i]]<- LETTERS[i]
 }

_ parse eval の使いかた 2006 12 27

もしも列名などで,ループ処理したいときは,以下のようにする.

 test<- data.frame(v1 = c("A","B","C"), v2 = c(1,2,3))
 sum( eval (parse(text =  paste("test$v", 2, sep="") ) ))
 [1] 6

ちなみに,以下は動作しない

 xxx<- "v2"
 sum( eval (parse(text =  "test$xxx") ))

ただし,一般には

 sum(text[,i])

として列を指定してループ処理する.

_ mtext の使いかた 2006 12 28

「Rの基礎とプログラミング技法」(シュプリンガージャパン)169ページのグラフィックス

ligges.R

 par(oma=rep(3, 4), bg="white")
 plot(c(0, 1), c(0, 1), type="n", ann=FALSE, axes=FALSE)
  rect(-1, -1, 2, 2, col="white")
 box("figure")
 par(xpd=FALSE)
 rect(-1, -1, 2, 2, col="white")
 box("plot", lty="dashed")
 text(.5, .5, "plot region", cex = 1.6)
 mtext("figure region", side=3, line=2, adj = 1, cex = 1.4)

 for (i in 1:4)
     mtext(paste("inner margin", i), side=i, 
    line=1, outer= FALSE)

 for (i in 1:4)
     mtext(paste("outer margin", i), side=i, 
    line=1, outer=TRUE)

 mtext("device region", side=3, line=2, 
     outer=TRUE, adj = 1, cex = 1.2)

 box("outer", col="black")
 
Link: Rの備忘録(1823d) R_fromOldHtml2(4099d)
Last-modified: 2008-02-04 (月) 18:22:58 (3967d)