2012年4月28日土曜日

R: ヒストグラムの一部に色をつける

追記:書き直しました(2013/05/25)


ヒストグラムの一部に色をつける関数を書いた。
あまりエレガントなやり口とはいえないが……


使い方)
R のコンソールに以下をコピーして貼り付ける↓
 histcol <-function(data,part,haba, iro="#0080ff", honsu=20){
  xv0 =seq(part[1],part[2],haba)
  xv =sort(rep(xv0,2))
  xv =xv[-c(1,length(xv))]
  yv=numeric(0)
  for(k in  1:(length(xv0)-1)){
    yv  =c(yv, rep(length(data[data >= xv0[k] & data < xv0[k+1]]),2))
  }
polygon(c(xv, rev(xv)) ,c(rep(0,length(yv)), rev(yv)), col=iro ,density=honsu, border = NA)
}
引数 part に色をつける範囲を x 座標で入力、引数 haba にヒストグラムの幅を入力する。
ヒストグラムの幅はヒストグラムを描いてみて目で見て確認する。
iro は線の色、honsu は線の本数。


使用例)

x <- rnorm(100)
#png("survreg.png")
hist(x)
histcol(data=x,part=c(-1,1),haba=0.5)
#dev.off()
 


2012年4月21日土曜日

R: survreg の分布一覧

このページ(とこの本)では、パラメトリック生存モデルを当てはめる関数 survreg(survival パッケージ)が紹介されている。
そして指定できる確率分布として、ワイブル分布(デフォルト)、指数分布(exponential)、ガウス分布(gaussian) 、対数正規分布 (log-normal) 、ロジスティック分布(logistic)、対数ロジスティック分布(log-logistic)が挙げられている。
しかし実際には、以下の通り、極値分布(extreme)、レイリー分布(rayleigh)、t分布も使える。
> library(survival)
> names(survreg.distributions)
 [1] "extreme"     "logistic"    "gaussian"    "weibull"     "exponential"
 [6] "rayleigh"    "loggaussian" "lognormal"   "loglogistic" "t"     
これらの分布を使う必要がある場面は、たぶんあまり多くはないが、以下に使い方をメモして置く。

例)"Gehan" データにパラメトリックモデルを当てはめる。

パッケージ Fadist を読み込むと対数ロジスティック分布などが使える。
library(survival)
library(FAdist)
モデル当てはめ。
sfs =survfit(Surv(time, cens)~treat,data=gehan)

srs =survreg(Surv(time, cens) ~ treat, dist="exponential", data = gehan)
srs_r =survreg(Surv(time, cens) ~ treat, dist="rayleigh", data = gehan)
srs_w =survreg(Surv(time, cens) ~ treat, dist="weibull", data = gehan)
srs_g =survreg(Surv(time, cens) ~ treat, dist="gaussian", data = gehan)
srs_l =survreg(Surv(time, cens) ~ treat, dist="logistic", data = gehan)
srs_lg =survreg(Surv(time, cens) ~ treat, dist="loggaussian", data = gehan)
srs_ll =survreg(Surv(time, cens) ~ treat, dist="loglogistic", data = gehan)
srs_ext =survreg(Surv(time, cens) ~ treat, dist="extreme", data = gehan)
survreg の "extreme(極値)" はグンベル分布のことらしい。FAdist にもグンベル分布は入っているのだが、パラメータの置き方がよく分からなかったので、ここでは自分で定義することにする。
my_ext = function(t,a=1,l=0){
  1-exp(-exp((t-l)/a))
  }
プロット。
#png("survreg.png")
par(mfrow= c(2,4))

plot(sfs, main="exponential",lty=1:2)
curve(1-pexp(x, rate=1/exp(coef(srs)[1])) , add=TRUE, lty=1)
curve(1-pexp(x, rate=1/exp(coef(srs)[1] + coef(srs)[2])) , add=TRUE, lty=2)

plot(sfs, main="rayleigh",lty=1:2)
curve(1-pweibull(x,shape=1/0.5, scale=exp(coef(srs_r)[1]) ), add=TRUE, lty=1)
curve(1-pweibull(x, shape=1/0.5, scale=exp(coef(srs_r)[1] + coef(srs_r)[2])) , add=TRUE, lty=2)

plot(sfs,  main="Weibull")
curve(1-pweibull(x,shape=1/srs_w$scale, scale=exp(coef(srs_w)[1])) ,add=TRUE, lty=1)
curve(1-pweibull(x,shape=1/srs_w$scale, scale=exp(coef(srs_w)[1] + coef(srs_w)[2])) ,add=TRUE, lty=2)

plot(sfs,lty=1:2,  main="Gaussian")
curve(1-pnorm(x,sd=srs_g$scale, mean=exp(coef(srs_lg)[1])) ,lty=1, add=TRUE)
curve(1-pnorm(x,sd=srs_g$scale, mean=exp(coef(srs_lg)[1] + coef(srs_lg)[2])) ,lty=2, add=TRUE)

plot(sfs,lty=1:2,  main="log-Gaussian")
curve(1-plnorm(x,sdlog=srs_lg$scale, meanlog=(coef(srs_lg)[1])) ,lty=1, add=TRUE)
curve(1-plnorm(x,sdlog=srs_lg$scale, meanlog=(coef(srs_lg)[1] + coef(srs_lg)[2])) ,lty=2, add=TRUE)
 
plot(sfs,lty=1:2,  main="logistic")
curve(1-plogis(x,scale=srs_l$scale, location=(coef(srs_l)[1])),lty=1, add=TRUE)
curve(1-plogis(x,scale=srs_l$scale, location=(coef(srs_l)[1] + coef(srs_l)[2])),lty=2, add=TRUE)
 
plot(sfs,lty=1:2,  main="log-logistic")
curve(1-pllog(x,shape=srs_ll$scale, scale=(coef(srs_ll)[1])),lty=1, add=TRUE)
curve(1-pllog(x,shape=srs_ll$scale, scale=(coef(srs_ll)[1] + coef(srs_ll)[2])),lty=2, add=TRUE)

plot(sfs,lty=1:2,  main="extreme")
curve(1-my_ext(x,a=srs_ext$scale, l=(coef(srs_ext)[1])),lty=1, add=TRUE)
curve(1-my_ext(x,a=srs_ext$scale, l=(coef(srs_ext)[1] + coef(srs_ext)[2])),lty=2, add=TRUE)
#dev.off()
t分布については、よく分からなかった。

2012年4月13日金曜日

R: グラフの一部を塗りつぶす

いつもながら R-Tips からの孫引き

例)
標準正規分布の上側 5% 部分に色を付けたい時は以下のようにする。
 plot(dnorm, -4, 4)
 xvals <- seq(qnorm(0.95), 4, length=10)   # 10等分
 dvals <- dnorm(xvals)  # 対応するグラフの高さ
 polygon(c(xvals,rev(xvals)),c(rep(0,10),rev(dvals)),col="gray") 

polygon の引数 density に正の整数を与えると、塗りつぶしが斜線になる。
density=1 だと 1 本、density=2 だと 2 本の斜線が引かれる。
angle で斜線の傾きを指定できる(省略可。デフォルトだと 45°)。
 polygon(c(xvals,rev(xvals)),  c(rep(0,10),rev(dvals)),col="magenta",angle=60,density=20)


2012年4月3日火曜日

R: 空のリストを作る

追記(1/18/2013)
以下ぜんぶなしで、やっぱり単に l <-list() でやりましょう

追記(10/7)
最近はこれより lapply(1:4, function(i){numeric()}) でやることが多い

かんたんに空のリストを作る方法


例)
長さ4のリストを作る
> l <-mapply(rep,1:4,0)
> l
[[1]]
integer(0)

[[2]]
integer(0)

[[3]]
integer(0)

[[4]]
integer(0)

文字列や行列も代入できる
> l[[1]]="hello"
> l[[2]]= matrix(1:6, nrow=2, ncol=3) 
> l
[[1]]
[1] "hello"

[[2]]
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

[[3]]
integer(0)

[[4]]
integer(0)