Paul Johnson 氏の R Tips 集訳 1/2

後半 P_Johnson_tips_2 へ移動

2003/08/09 現在

注:これは P. Johnson さんが r-help の記事から役にたちそうな箇所を抜き出したものです。オリジナルは http://pj.freefaculty.org/R/Rtips.html。様々な話題があり、参考になります。(なお元そのものが r-help 記事の抄約ですから、著作権問題はクリアされていると考えますが?)
注意:あまりに長くて表示に時間がかかるので二つに分割しました。


(私はこれを "StatsRus" と名付けるのが気がきいていると思ったのですが、玩具店の法律顧問が電話してきて、「ご存知のように ...」)

もしあなたが R の Tips 集が必要なら、これがそうです。これは R ドキュメントの代用品ではなく、電子メイルで読んだことを覚えておくために使っているものです。

Brian D. Ripley の言葉を注意しましょう、 "One enquiry came to me yesterday which suggested that some users of binary distributions do not know that R comes with two Guides in the doc/manual directory plus an FAQ and the help pages in book form. I hope those are distributed with all the binary distributions (they are not made nor installed by default). Windows-specific versions are available in PDF in rw1001d1.zip and rw1001d2.zip"

FaqManager ソフトに対しその著者の Stas Bekman に謝辞を捧げます。


1 データの入出力

1.1 生の数値を R に取り込む

(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)

これは本当に簡単。ファイル "myData"で、空白区切りの数値があり、第1行目にはカラム名があるとすると、Rのコマンドは以下のようになり、

myDataFrame<-read.table("myData",header=TRUE)

このコマンドは素晴らしく動く。

"?read.table" とタイプすれば、他の区切り記号を使っているファイルのインポート方法が示される。

タブ区切りのデータを持っていて、空白を欠損値を示すのに使っている場合は,以下のようにする。

> myDataFrame<-read.table("myData",sep="\t",na.strings=" ",header=TRUE)

(R1.3 の新しいフィーチャー) データが gzipで圧縮されたファイル myData.gz である場合。このようにする。

myDataFrame <- read.table(gzfile("myData.gz"), header=T)

別々のファイルの列を読み込む際には、以下のようにして、これらを組み合わせて、一つのデータフレームを作成する。

variable1 <- scan("file1")
variable2 <- scan("file2")
mydata <- cbind(variable1, variable2)

もしくは、以下を使えば、同様なことができる。

mydata <- data.frame(variable1, variable2)

R オブジェクトのデータフレームを、以下のコマンドで保存することもできる。

write.table(mydata, file="filename3")

1.2 データアクセスに関する基本的注意

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

データフレーム "x" の列を列番号でアクセスするには、例えば、x[,1], とすれば、最初の列を取得できる。列の名前を知っている場合は、例えば "pjVar1" とすれば、 x$pjVar1 と同じものになる。
注:原文にも抜けているのだが,x["pjVar1"] は x$pjVar1 と同じになるといいたいのだろう。
注の注:「第1列の名前が"pjVar1"であることを知っていれば」,x[,1] と書くのと x$pjVar1 とが同じになるといいたいのだ。

もし,そのデータフレームに対して attach(x) を行うと,単に pjVar1 と参照すれば,R はデータフレームの中でそれを見つけるだろう。このことで私はいくつかのトラブルを体験した。

回帰分析において,data=x オプションを含めることにより,あなたは pjVar1 を変数名として使える。それは一番クリーンなプログラムコードになる。

リストの要素を得るためには x[[1]] とすればよい。

例えば,もしデータフレームが "nebdata" ならば,

  nebdata[1, 2]

あるいは

  nebdata[nebdata$Volt=="ABal", 2]

あるいは,データフレームにアタッチすることで,

  attach(nebdata)
  nebdata[Volt=="ABal", 2]
  detach(nebdata)

(from Diego Kuonen)

1.3 新しい Data Import/Export マニュアルのチェック

(13/08/2001 Paul Johnson <pauljohn@ukans.edu>)

R-1.2 において,R チームは R にデータを読み込んだり R からデータを書き出すためのマニュアルを公開した。ヘルプが必要なら,まず見るべきものである。現在それは http://lib.stat.cmu.edu/R/CRAN/doc/manuals/R-data.pdf にある。

1.4 データを R と他のプログラム (Excelなど)の間で交換する

(04/05/2001 Paul Johnson <pauljohn@ukans.edu>)

私はこの覚え書きを新しい R のデータインポート・データエクスポートマニュアルが出る前に書いた。今はそのマニュアルを読めばよい(前項参照のこと)。

Excel -> R

エクセルで,例えば table.txt などという名前のテキストファイルを,「空白区切りのテキストファイル(.prn)として保存」を使って,作りなさい。次に,R は,data.matrix(read.table("table.txt")) により,行と列のラベル(名前)を持った二次元配列を戻り値として返す。

R -> Excel

テキストファイルに書き出すときには,sink("table") とする('table' は R オブジェクトの名前)。Excel で,そのファイルを開き,インポート・ウィザードの手順に従う。 (from Griffith Feeney)

R-Excel インターフェースは試したかい? http://lib.stat.cmu.edu/R/CRAN/contrib/extra/index.html これを見るとよい。 (from Charles Raux)

私の学生達に教えようと思っていることを書いておく(インポートのみ。エクスポートは write.table を使う)。

今のところ,最も簡単なやり方は,えとね,データをテキストファイルとしてエクスポーし,それを R でデータフレームとして読むためには,read.table,read.csv または read.csv2 のどれかを使う。

現在の最も普遍的な三つの PC プログラムにおけるテクニックは以下のようなものである。若干のバリエーションを含むタブ区切りが一般的である。

SPSS:データエディタの File, Export メニューにより,タブ区切りファイルとしてエクスポートする(説明を簡単にするため C: ドライブのトップディレクトリの mydata.txt としよう)。変数名からなるヘッダー行を含むことを確実にするために,read.csv2("C:/mydata.txt",sep="\t") により読み込む。

このやり方はカンマが小数点として使われているかどうか locale を参照する。そのやり方を避けたいなら,read.csv を使え。

Excel: File, Save as... を使い,タブ区切り(.txt)粗選択するか,クリップボードを使う。それ以降のデータの扱いは SPSS の場合と同じ。ヘッダー行は元のスプレッドシートにあるときにのみ,存在する。

SAS: カンマ区切りかタブ区切りでエクスポートできる。SAS は locale に関係なく小数点には "." を使う。ファイルを読むには,read.csv("C:/mydata.csv") とか read.csv("C:/mydata.csv",sep="\t") を使う。

(from Peter Dalgaard)

1.5 データフレームのマージ

(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)

experiment<-data.frame(times=c(0,0,10,10,20,20,30,30),expval=c(1,1,2,2,3,3,4,4))
simul<-data.frame(times=c(0,10,20,30),simul=c(3,4,5,6))

以下のように、マージされたデータセットが欲しい。

>   times expval simul
> 1     0      1     3
> 2     0      1     3
> 3    10      2     4
> 4    10      2     4
> 5    20      3     5
> 6    20      3     5
> 7    30      4     6
> 8    30      4     6

答え

merge(experiment, simul)

が,貴方の望む全てだろう。 (from Brian D. Ripley)

あるいは,以下を検討する

 exp.sim<-data.frame(experiment,simul=
               simul$simul[match(experiment$times,simul$times)]) 

(from Jim Lemon)

このようなデータフレームを持っている。

state count1 percent1
CA    19    0.34
TX    22    0.35
FL    11    0.24
OR    34    0.42
GA    52    0.62
MN    12    0.17
NC    19    0.34

state count2 percent2
FL    22     0.35
MN    22     0.35
CA    11     0.24 
TX    52     0.62

以下が欲しい

state count1 percent1 count2 percent2
CA    19    0.34      11     0.24
TX    22    0.35      52     0.62
FL    11    0.24      22     0.35
OR    34    0.42       0     0
GA    52    0.62       0     0
MN    12    0.17      22     0.35
NC    19    0.34       0     0

( from Yu-Ling Wu )

Ben Bolker のレスは以下のとおり

state1 <- c("CA","TX","FL","OR","GA","MN","NC")
count1 <- c(19,22,11,34,52,12,19)
percent1 <- c(0.34,0.35,0.24,0.42,0.62,0.17,0.34)
state2 <- c("FL","MN","CA","TX") 
count2 <- c(22,22,11,52) 
percent2 <-   c(0.35,0.35,0.24,0.62)
data1 <- data.frame(state1,count1,percent1) 
data2 <-  data.frame(state2,count2,percent2)
datac <- data1
m <- match(data1$state1,data2$state2,0) 
datac$count2 <- ifelse(m==0,0,data2$count2[m]) 
datac$percent2 <- ifelse(m==0,0,data2$percent2[m])
(日本語訳注!)上記のmatchを利用した結果は,以下のようになるので誤りです。(2004/1/30)
state count1 percent1 count2 percent2
CA    19    0.34      11     0.24
TX    22    0.35      52     0.62
FL    11    0.24      22     0.35
OR    34    0.42       0     0
GA    52    0.62       0     0
MN    12    0.17      52 <-- 0.62 <-- ここ二箇所が変です。 
NC    19    0.34       0     0

(訳注その2)次のようにしたらうまくいくけど,どうでしょうか。(2004/05/20)

m <- match(data2$state2,data1$state1,0) 
datac[,"percent2"] <- datac[,"count2"] <- 0
datac[m, "count2"] <- count2
datac[m, "percent2"] <- percent2

If you didn't want to keep all the rows in both data sets (but just the shared rows) you could use

merge(data1,data2,by=1)

1.6 一時に一つの行を加える

(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)

Question: I would like to create an (empty) data frame with "headings" for every column (column titles) and then put data row-by-row into this data frame (one row for every computation I will be doing), i.e.

no. time temp pressure <---the headings
1   0     100   80     <---first result
2   10    110   87     <---2nd result .....

Answer: Depends if the cols are all numeric: if they are a matrix would be better.

If you know the number of results in advance, say, N, do this

df <- data.frame(time=numeric(N), temp=numeric(N), pressure=numeric(N))
df[1, ] <- c(0, 100, 80)
df[2, ] <- c(10, 110, 87)
...

または

m <- matrix(nrow=N, ncol=3)
colnames(m) <- c("time", "temp", "pressure")
m[1, ]  <- c(0, 100, 80)
m[2, ] <- c(10, 110, 87)

The matrix form is better size it only needs to access one vector (a matrix is a vector with attributes), not three.

If you don't know the final size you can use rbind to add a row at a time, but that is substantially less efficient as lots of re-allocation is needed. It's better to guess the size, fill in and then rbind on a lot more rows if the guess was too small.(from Brian Ripley) [top]

1.7 データフレームに対するもう一つのマージ法

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

(from Stephen Arthur)

例のデータ: File1

C A T

File2

1 2 34 56
2 3 45 67
3 4 56 78

これは以下を生成:

> C A T 1 2 34 56
> C A T 2 3 45 67
> C A T 3 4 56 78

以下の関数は動作する:

repcbind <- function(x,y) {
 nx <- nrow(x)
 ny <- nrow(y)
 if (nx<ny)
   x <- apply(x,2,rep,length=ny)
 else if (ny<nx)
   y <- apply(y,2,rep,length=nx)
 cbind(x,y)
}

(from Ben Bolker)

1.8 NULL オブジェクトがあるかどうか検査する

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

R がロードする際に、R はロードするものがないとき警告は出さず、これを NULL として記録する。従って、これらのものをチェックする必要がある。以下のコマンドを使ってチェックする。

is.null(list$component)

上のコマンドは list という名のものにある component と名づけられたものをチェックする。

存在しないデータフレームのカラムを "[" を使ってアクセスするとエラーが出るので、その代わりに以下を実行する。

>data(trees)
>trees$aardvark
NULL
>trees[,"aardvark"]
Error in [.data.frame(trees, , "aardvark") : subscript out of bounds 

(from Thomas Lumley)

1.9 乱数を発生する

(22/08/2000 Paul Johnson <pauljohn@ukans.edu>)

runif() --uniform random on [0,1) 
rnorm() --normal random

multivariate normal: Given covar matrix:

0.25,    0.20
0.20,    0.25

Brian Ripley said: "mvrnorm in package MASS (part of the VR bundle).

mvrnorm(2, c(0,0), matrix(c(0.25, 0.20, 0.20, 0.25), 2,2))

and Peter Dalgaard observed "a less general solution for this particular case would be

rnorm(1,sd=sqrt(0.20)) + rnorm(2,sd=sqrt(0.05))

Wait. You want randomly drawn integers? Here:

If you mean sampling without replacement:

>sample(1:10,3,replace=FALSE)

If you mean with replacement:

>sample(1:10,3,replace=TRUE)

(from Bill Simpson)

1.10 平均/分散を固定して乱数を発生

(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)

任意のツールで乱数を生成する場合、指定したちょうどの平均を持つサンプルが得られない。平均0のジェネレータの生成するサンプルは平均値がばらばらになる。これって、正しいよね?

そこで、平均0のサンプルをどのようにして強制的に作成するのか。以下の2段階のアプローチがある。

> x <- rnorm(100, mean = 5, sd = 2) # ここで mean, sd を指定しても無意味
> x <- (x - mean(x))/sqrt(var(x))
> mean(x)
[1] 1.385177e-16
> var(x)
[1] 1

日本語訳を見ての注:ここまでの結果を得るには scale 関数を使うべし

> x <- scale( rnorm(100) ) # これで x は 平均0,分散1に正規化される
> mean(x); sd(x)
[1] 1.297573e-17 # set.seed していないので上とは結果が違うが気にしない
[1] 1

そして平均5、標準偏差2のサンプルができる。

# 以下は同じ
> x <- x*2 + 5
> mean(x)
[1] 5
> var(x)
[1] 4

(from Torsten.Hothorn)

1.11 重みつきの複製

(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)

x<-c(10,40,50,100)  # income vector for instance
w<-c(2,1,3,2)       # the weight for each observation in x with the same
rep(x,w)
[1]  10  10  40  50  50  50 100 100

(from P. Malewski)

1.12 重みつきのデータを表すようにデータセットを拡大する

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

Also, most of the functions that have weights have frequency weights rather than probability weights: that is, setting a weight equal to 2 has exactly the same effect as replicating the observation.

If you have frequency weights you may have to expand your dataset.

 expanded.data<-as.data.frame(lapply(compressed.data,
                  function(x) rep(x,compressed.data$weights)))

should do it.

(from Thomas Lumley)

1.13 分割表をデータフレームに変換する

(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)

もし8次元のクロス表があるなら、表にセルの8つのファクターと1つのカラムのデータフレームが欲しい。

R1.2 はこれを扱うのに as.data.frame.table() 関数を導入している。

この処理はマニュアルでもできる。ここに関数がある(expand.grid を中心にしたシンプルなラッパーである)。

dfify <- function(arr, value.name = "value", dn.names =
                           names(dimnames(arr))) {
Version <- "$Id: dfify.sfun,v 1.1 1995/10/09 16:06:12 d3a061 Exp $" 
dn <- dimnames(arr <- as.array(arr)) 
if(is.null(dn)) stop("Can't data-frame-ify an array without dimnames") 
names(dn) <- dn.names 
ans <- cbind(expand.grid(dn), as.vector(arr)) 
names(ans)[ncol(ans)] <- value.name ans }

関数名は "data-frame-i-fy" を縮めたものである。

質問の例では、多重配列が適切な dimnames を持っていると推定されるので、以下のようにすればよい。

my.data.frame <- dfify(my.array, value.name="frequency")

(from Todd Taylor)

1.14 書き込み: アスキーファイルのデータ

(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)

たとえば 28 x 28 データマトリックスを生成するコマンドを持っている場合。どのようにして、マトリックスを txt ファイルを出力できるか(マトリックスをコピー/ペーストするより)?

write.table(mat,file="filename")

Pesl Thomas は「"tab"という名前のデータフレームの内容をファイルをディスクに出力したい」と言った。

> tab
or 1 2 3 4 5 6 1 57 80 25 23 46 23 2 106 35 59 8 51 4

これが使える

write.table(unclass(tab))

これも使える

 write(t(tab),ncol=NCOL(tab))

なぜなら、 write() は行と列を転置する。(from Thomas Lumley)

注目してもらいたいのは、MASS ライブラリには、write.matrix 関数がある。こちらの方が、数値マトリックスを書き出す必要がある場合には速い。大仕事には打って付けである。

2 データフレームを使った作業: 記録、合体

2.1 変数をデータフレーム(又はリスト)に加える

(02/06/2003 Paul Johnson <pauljohn@ku.edu>)

In r-help, I asked " I keep finding myself in a situation where I want to calculate a variable name and then use it on the left hand side of an assignment. For example

> iteration <- 1
> varName <- paste("run",iteration,sep="")
> myList$parse(text=varName) <- aColumn

I got many useful answers that opened my eyes!

Brian Ripley said: For a data frame you could use

mydf[paste("run",iteration,sep="")] <- aColumn

and for a list or a data frame

Robject[[paste("run",iteration,sep="")]] <- aColumn

And Thomas Lumley added: " If you wanted to do something of this sort for which the above didn't workyou could also learn about substitute():

 eval(substitute(myList$newColumn<-aColumn),
        list(newColumn=as.name(varName)))

as an alternative to parse().

2.2 変数名を「その場で」作る

(10/04/2001 Paul Johnson <pauljohn@ukans.edu>)

変数名とそれにセットする値を一緒に貼り付ける。assign を使う。以下のように。

> assign(paste("file", 1, "max", sep=""), 1)
> ls()
[1] "file1max"

(Brian Ripley, June 18, 2001)

2.3 一つの列を記録する、値を別の列に出力する

??Recode one column, output values into another column (20/06/2001 Paul Johnson <pauljohn@ukans.edu>)

Aleksey Naumov がこの質問をした。例えば、ここに簡単なデータフレームがある。

 var1
1    1
2    2
3    3

どのようにして、上の "var1" をグループ化した値 'var2' を新規作成の列としてデータに追加するのか?たとえば以下のようにすればよい。

  var1 var2
1    1    4
2    2    4
3    3    5

data$var2 <- c(4,4,5)[data$var1]

または

data <- transform(data, var2=c(4,4,5)[var1])

(from Peter Dalgaard)

I admit the previous seems complicated, and I'm inclined to take the easy road. Some cases are particularly simple because of the way arrays are processed. Suppose you create a variable, and then want to reset some values to missing. Go like this:

 x <- rnorm(10000) 
 x[x > 1.5] <- NA

And if you don't want to replace the original variable, create a new one first ( xnew <- x ) and then do that same thing to xnew.

You can put other variables inside the brackets, so if you want x to equal 234 if y equals 1, then

x[y == 1] <- 234

So, to do the v1 v2 example above, I think it is clearest to do:

> v1 <- c(1,2,3)

#now initialize v2
> v2 <- rep( -9, length(v1))  
#now recode v2
> v2[v1==1] <- 4
> v2[v1==2]<-4
> v2[v1==3]<-5
> v2
[1] 4 4 5

Note that R's "ifelse" command can work too:

x<-ifelse(x>1.5,NA,x)

I asked Mark Myatt (author of the ReX guide mentioned below) for more examples:

For example, suppose I get a bunch of variables coded on a scale

  1= no
  6  = yes
  8 = tied
  9 = missing
 10 = not applicable.

Recode that into a new variable name with 0=no, 1=yes, and all else NA, for example.

It seems like the replace() function would do it for single values but you end up with empty levels in factors but that can be fixed by re- factoring the variable. Here is a basic recode() function:

recode <- function(var, old, new)
{
x <- replace(var, var == old, new)
if(is.factor(x)) factor(x) else x
}

For the above example:

test <- c(1,1,2,1,1,8,1,2,1,10,1,8,2,1,9,1,2,9,10,1)
test
test <- recode(test, 1, 0)
test <- recode(test, 2, 1)
test <- recode(test, 8, NA)
test <- recode(test, 9, NA)
test <- recode(test, 10, NA)
test

Although it is probably easier to use replace():

test <- c(1,1,2,1,1,8,1,2,1,10,1,8,2,1,9,1,2,9,10,1)
test
test <- replace(test, test == 8 | test == 9 | test == 10, NA)
test <- replace(test, test == 1, 0)
test <- replace(test, test == 2, 1)
test

I suppose a better function would take from and to lists as arguments:

recode <- function(var, from, to)
{
x <- as.vector(var)
for (i in 1:length(from))
{
  x <- replace(x, x == from[i], to[i])
}
if(is.factor(var)) factor(x) else x
} 

For the example:

test <- c(1,1,2,1,1,8,1,2,1,10,1,8,2,1,9,1,2,9,10,1)
test
test <- recode(test, c(1,2,8:10), c(0,1))
test

and it still works with single values.

Suppose somebody gives me a scale from 1 to 100, and I want to collapse it into 10 groups, how do I go about it?

Mark says: Use cut() for this. This cuts into 10 groups:

test <- trunc(runif(1000,1,100))
groups <- cut(test,seq(0,100,10))
table(test, groups)

To get ten groups without knowing the minimum and maximum value you can use pretty():

groups <- cut(test,pretty(test,10))
table(test, groups)

You can specify the cut-points:

groups <- cut(test,c(0,20,40,60,80,100))
table(test, groups)

And they don't need to be even groups:

groups <- cut(test,c(0,30,50,75,100))
table(test, groups)

Mark added, "I think I will add this sort of thing to the REX pack."

2.4 指示(ダミー)変数をつくり出す

(20/06/2001 Paul Johnson <pauljohn@ukans.edu>)

以下に2例ある。

c is a column, you want dummy variable, one for each valid value. First, make it a factor, then use model.matrix():

> x <- c(2,2,5,3,6,5,NA)
> xf <- factor(x,levels=2:6)
> model.matrix(~xf-1)
  xf2 xf3 xf4 xf5 xf6
1   1   0   0   0   0
2   1   0   0   0   0
3   0   0   0   1   0
4   0   1   0   0   0
5   0   0   0   0   1
6   0   0   0   1   0
attr(,"assign")
[1] 1 1 1 1 1

(from Peter Dalgaard )

質問;5つのカテゴリーのアル変数があり、各カテゴリに対するダミー変数を作りたい。

答え;行の添え字または model.matrix を使う。

ff <- factor(sample(letters[1:5], 25, replace=TRUE))
diag(nlevels(ff))[ff,]

または

model.matrix(~ ff - 1)

(from Brian D. Ripley)

2.5 時系列回帰のための変数の遅延値をつくり出す

(02/06/2003 Paul Johnson <pauljohn@ku.edu>)

Peter Dalgaard explained, "the simple way is to create a new variable which shifts the response, i.e.

yshft <- c(y[-1], NA) # pad with missing
summary(lm(yshft ~ x + y))

Alternatively, lag the regressors:

N <- length(x)
xlag <- c(NA, x[1:(N-1)])
ylag <- c(NA, y[1:(N-1)])
summary(lm(y ~ xlag + ylag))

2.6 関連する観測値を持たない因子水準をデータセットから取り除くには?

(08/01/2002 Paul Johnson <pauljohn@ukans.edu>)

もっともいいのが、水準を落とすことである。たとえば、以下のようになる、

problem.factor <- problem.factor[,drop=TRUE] (from Brian D. Ripley)

2.7 データフレームから観測値を選択する、部分集合を取り出す

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

ある変数が欠損値である行を省きたいなら、John Fox の言うように以下を実行する。

x2 <- na.omit(x)

これにより、欠損値を含む行すべてが除かれたデータフレーム x がコピーされる。関数 na.exclude も利用できる。欠損値についてより詳細な情報が欲しければ、ヘルプ : ?na.exclude を参照されたい。

欠損値の除外については、Peter Dalgaard は以下のコマンドを好んでいる。

subset(x,complete.cases(x)) or x[complete.cases(x),]

x != "NA" を追加するよりも、 "is.na(x)” の方が望ましい。

Suppose you want observations with c=1 in df1. This makes a new data frame you want.

df2 <- df1[df1$c == 1,]

and note that indexing is pretty central to using S (the language), so it is worth learning all the ways to use it. (from Brian Ripley)

Or use "match" select values from the column "d" by taking the ones that match the values of another column, as in

> d <- t(array(1:20,dim=c(2,10)))
> i <- c(13,5,19)
> d[match(i,d[,1]), 2]
[1] 14  6 20

(from Peter Wolf)

This gets lines from the dataframe that meet the test like so:

> d[d[,1] %in% i,]
     [,1] [,2]
[1,]    5    6
[2,]   13   14
[3,]   19   20

(from R. Woodrow Setzer, Jr.)

Till Baumgaertel wanted to select observations for men over age 40, and sex was coded either m or M. Here are two working commands:

maleOver40 <- subset(d, sex %in% c("m","M") & age > 40)
maleOver40 <- d[(d$sex == "m" | d$sex == "M") & d$age > 40,]

2.8 最初の観測値をそれぞれから抹消する

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

Given data like:

1  ABK     19910711     11.1867461      0.0000000       108
2  ABK     19910712     11.5298979     11.1867461      111
6 CSCO   19910102      0.1553819      0.0000000         106
7 CSCO   19910103      0.1527778      0.1458333         166

remove the first observation for each value of the "sym" variable (the one coded ABK,CSCO, etc). . If you just need to remove rows 1,5, and 13, do:

newhilodata <- hilodata[-c(1,6,13),]

To solve the more general problem of omitting the first in each group, assuming "sym" is a factor, try something like

newhilodata <- subset(hilodata, diff(c(0,as.integer(sym))) != 0)

(actually, the as.integer is unnecessary because the c() will unclass the factor automagically)

(from Peter Dalgaard)

Alternatively, you could use the match function because it returns the first match. Suppose jm is the data set. Then:

> match(unique(jm$sym), jm$sym)
 [1]  1  6 13
> jm <- jm[ -match(unique(jm$sym), jm$sym), ]

(from Douglas Bates)

As Robert pointed out to me privately: duplicated() does the trick

 subset(hilodata, duplicated(sym))

has got to be the simplest variant.

(from Peter Dalgaard)

2.9 データのランダム標本を選ぶ

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

sample(N, n, replace=F)

and

seq(N)[rank(runif(N)) <= n]

is another general solution.

(from Brian D. Ripley)

2.10 モデルの変数選択:subset 関数を忘れるな

(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)

直接行その他を削除することでデータを管理できるが、subset() を使えばデータをまったくいじらずに同じことができる。もっと多くのことを知りたければ ?select を実行すればよい。Subset はまた lm のような統計関数の多くでオプションになっている。

Peter Dalgaard が「組み込み」データセット airquality を使ったこの例をくれた。

data(airquality)
names(airquality)
lm(Ozone~.,data=subset(airquality,select=Ozone:Month))
lm(Ozone~.,data=subset(airquality,select=c(Ozone:Wind,Month)))
lm(Ozone~.-Temp,data=subset(airquality,select=Ozone:Month))

lm の RHS の "." は「すべての変数」を意味しており、rhs の subset コマンドはデータセットより別の変数を選び出す。"x1:x2" は x1 と x2 の間に含まれる変数を意味している。

2.11 すべての数値変数を処理し、文字変数を無視するには?

??Process all numeric variables, ignore character variables? (11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

data<-read.table("lastline.txt",header=T,as.is = TRUE) 
indices<-1:dim(data)[2] 
indices<-na.omit(ifelse(indices*sapply(data,is.numeric),indices,NA)) 
mean<-sapply(data[,indices],mean) sd<-sapply(data[,indices],sd)

2.12 二つ以上の変数に関しソート

(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)

Can someone tell me how can I sort a list that contains duplicates (name) but keeping the duplicates together when sorting the values.

name  M
1234  8
1234  8.3
4321  9
4321  8.1

I also would like to set a cut-off, so that anything below a certain values will not be sorted.(from Kong, Chuang Fong)

I take it that the cutoff is on the value of M. OK, suppose it is the value of `coff'.

sort.ind <- order(name, pmax(coff, M))  #  sorting index
name <- name[sort.ind]
M <- M[sort.ind]

Notice how using pmax() for "parallel maximum" you can implement the cutoff by raising all values below the mark up to the mark thus putting them all into the same bin as far as sorting is concerned.

If your two variables are in a data frame you can combine the last two steps into one, of course.

sort.ind <- order(dat$name, pmax(coff, dat$M))
dat <- dat[sort.ind, ]

In fact it's not long before you are doing it all in one step:

dat <- dat[order(dat$name, pmax(coff, dat$M)), ]

(from Bill Venables)

I want the ability to sort a data frame lexicographically according to several variables

Here's how:

spsheet[order(name,age,zip),] 

(from Peter Dalgaard)

2.13 ある因子により定義されるサブグループ内で順位を取る

(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)

by() のヘルプを読みなさい

> by(x[2], x$group, rank)
x$group: A
[1] 4.0 1.5 1.5 3.0
------------------------------------------------------------ 
x$group: B
[1] 3 2 1

> c(by(x[2], x$group, rank), recursive=T) 
A1 A2 A3 A4 B1 B2 B3 4.0 1.5 1.5 3.0 3.0 2.0 1.0

(from Brian D. Ripley)

2.14 データフレームから欠損値を除外するために na.omit を使う。

(22/08/2000 Paul Johnson <pauljohn@ukans.edu>)

欠損値が "insure" という名前のデータセットより除去されるようにする。

thisdf <- na.omit(insure[,c(1, 19:39)])
body.m <- lm(BI.PPrem ~ ., data=thisdf, na.action=na.fail)

2.15 各行から一つの集約値を計算

(16/08/2000 Paul Johnson <pauljohn@ukans.edu>)

質問:200行で,各行に32個の実数データを含むテキストファイルから値を読み込んで,32個の数値の平均値を求めたい。そして,1〜200を要素として持つ x ベクトルと,200個の平均値を要素とする y ベクトルを作りたい。

Peter Dalgaard は以下を提案した。

y <- apply(as.matrix(read.table("myfile")), 1, mean)
x <- seq(along=y)

(ファイル形式によっては,read.table に,いくつかのオプションを追加する必要があろう)

2.16 各因子に対する集約値を保持する新しいデータフレームをつくり出す。

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

これでどう

Z<-aggregate(X, f, sum)

(X の変数はすべて合計できるものとする)

または、(もし)Xが因子も含んでいる場合、要約統計の計算に必要な変数を選択することになる。そこで、以下を使った。

Z <- data.frame(f=levels(f),x1=as.vector(tapply(x1,f,sum)))

(from Wolfgang Koller)

2.17 データフレームを選択的に集約

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

以下のようなデータフレームがあって

10 20  23 44 33
10 20  33 23 67

つぎのような結果を望むなら(訳注:何を表しているかはプログラムを読むしかない(^_^))

10 20   56 67 100

以下をやってみなさい。

dat<-read.table("data.dat",header=TRUE)
aggregate(dat[,-(1:2)], by=list(std=dat$std, cf=dat$cf), sum)

最初の2列は [,-(1:2)] で除外され,by オプションで出力のために取っておかれている。

2.18 実数から各桁数値を一時に一つずつ取り出す

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

小数点以下1桁目の数値を取り出し出現度数をプロットする。次いで,小数点以下2桁目の数値についても同じことをする。

a<- log(1:1000)
d1<-floor(10*(a-floor(a)))  # first decimal
par(mfrow=c(2,2))
hist(d1,breaks=c(-1:9))
table(d1)
d2<-floor(10*(10*a-floor(10*a))) # second decimal
hist(d2,breaks=c(-1:9))
table(d2)

(from Yudi Pawitan)

x <- 1:1000
ndig <- 6
(ii <- as.integer(10^(ndig-1) * log(x)))[1:7] 
(ci <- formatC(ii, flag="0", wid= ndig))[1:7] 
cm <- t(sapply(ci, function(cc) strsplit(cc,NULL)[[1]]))
cm [1:7,]
apply(cm, 2, table) #--> Nice tables

# The plots :

par(mfrow= c(3,2), lab = c(10,10,7)) 
for(i in 1:ndig) 
 hist(as.integer(cm[,i]), breaks = -.5 + 0:10, 
         main = paste("Distribution of ", i,"-th digit"))

(from Martin Maechler)

2.19 リスト中の各副行列から一つの項目を取り出す

(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)

Let Z denote the list of matrices. All matrices have the same order. Suppose you need to take element [1,2] from each.

lapply(Z, function(x) x[1,2])

should do this, giving a list. Use sapply if you want a vector. (Brian Ripley)

2.20 データセット中の値を示すベクトルを得る

(10/04/2001 Paul Johnson <pauljohn@ukans.edu>)

xlevels<-sort(unique(x))

Which you can use in a contour plot, as in

contour(xlevels,ylevels,zvals)

2.21 Rの命令を表す文字列の値を計算する

(13/08/2000 Paul Johnson <pauljohn@ukans.edu>)

Use the eval(parse()) pairing:

String2Eval <- "A valid R statement"
       eval(parse(text = String2Eval)) 

(from Mark Myatt)

Also check out substitute(), as.name() et al. for other methods of manipulating expressions and function calls (from Peter Dalgaard)

Or

eval(parse(text="ls()"))

Or

eval(parse(text = "x[3] <- 5"))

2.22 "which" はあるテストをパスする観測値の添字値を抜き出すことができる

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

大きなベクトルのデータを箱ひげ図で分析して、外れ値を見つけるには、以下を試してみる

which(x==boxplot(x,range=1)$out)

2.23 行列/データフレーム中の一意的な行を見出す

(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)

Jason Liao writes:

"I have 10000 integer triplets stored in A[1:10000, 1:3]. I would like to find the unique triplets among the 10000 ones with possible duplications."

Peter Dalgaard answers, "As of 1.4.0 (!):

unique(as.data.frame(A))"

2.24 数値変数だけを選ぶ

(12/11/2002 Paul Johnson <pauljohn@ku.edu>)

sapply(dataframe, is.numeric)

または

which(sapply(data.frame, is.numeric))

Thomas Lumley 11/11/2002

3 行列とベクトルの演算

3.1 単位行列をつくり出すには?

(16/08/2000 Paul Johnson <pauljohn@ukans.edu>)

diag(n)

Or, go the long way:

n<-c(5)
I <- matrix(0,nrow=n,ncol=n)
I[row(I)==col(I)] <- 1

E.D.Isaia

3.2 行列 "m" を一つの長いベクトルに変換

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

dim(m)<-NULL, or c(m) (from Peter Dalgaard )

3.3 変わった数列 (1 2 3 4 1 2 3 1 2 1) をつくり出す

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

なぜこれが有用なのかは分かりませんが、いくつかの列と行列を見せてくれる。

結局は Brian Ripley の方法を使うことになった。この方法は最初に手に入れ、動いたし、以下のようにね。

A <- matrix(1, n-1, n-1) 
rA <- row(A) 
rA[rA + col(A) <= n] 

しかしながら、 その後、Andy Royle の御蔭で、もっと簡単ですばらしい解決法を見つけた。

sequence((n-1):1) 
[1] 1 2 3 4 1 2 3 1 2 1 

Karen Kotschy

3.4 すべての n 個めの項目を選ぶ

(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)

非常に長いベクトルよりn個目の要素をすべて取り出す、これは vec かな?

vec[seq(n, length(vec), n)]  #(from Brian Ripley)

seq(1,11628,length=1000) はインデックスとして使える 1:11628 より1000 個の等間隔の数を与える。(from Thomas Lumley)

My example:

vec <- rnorm(1999)
newvec <- vec[1, length(vec), 200]
> newvec
 [1]  0.2685562  1.8336569  0.1371113  0.2204333 -1.2798172  0.3337282
 [7] -0.2366385  0.5060078  0.9680530  1.2725744

This shows the items in vec at indexes (1, 201, 401, ..., 1801)

3.5 ベクトル中の 1.5 に最も近い値の添字を見出す

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

n <- 1000
x <- sort(rnorm(n))
x0 <- 1.5
dx <- abs(x-x0)
which(dx==min(dx))

(from Jan Schelling)

which(abs(x - 1.5) == min(abs(x - 1.5)))

(from Uwe Ligges)

3.6 ベクトル中の非負値項目の添字を見つける

(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)

which(x!=0) (from Uwe Ligges)
rfind <- function(x) seq(along=x)[x != 0]  #(from Brian D. Ripley)

スピードとメモリーの効率性に関しては、以下がある。

as.logical(x)

これは以下よりもよく

x!=0

そして

seq(along=x)[as.logical(x)]

これは以下よりもよく

which(as.logical(x))

このため

which(x!=0)

は最も短いが

rfind <- function(x)seq(along=x)[as.logical(x)]

以下はコンピュータではもっとも効率的であるようである。

(from Jens Oehlschl?gel-Akiyoshi)

3.7 欠損値の添字を見出す

(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)

Pes というベクトルが600の観測値を持つとする。欠損値の添え字を取り出すために,以下のようにしてはいけない

(1:600)[is.na(Pes)]

正しいやり方は,次の通り。

seq(along=Pes)[is.na(Pes)]

この場合にはどちらでも良いことではあるが,添え字の長さが0ということもあるので,ライブラリやパッケージを作るときにはこのような点に注意していないと,作者が思っている以上に困った問題が生じる。

(from Brian Ripley)

しかし,私は以下の解を与える

which(is.na(Pes))

is the one I still really recommend; it does deal with 0-length objects, and it keeps names when there are some, and it has an `arr.ind = FALSE' argument to return array indices instead of vector indices when so desired. (from Martin Maechler)

3.8 ベクトル中の最大値を持つ項の添字を見出す

(16/08/2000 Paul Johnson <pauljohn@ukans.edu>)

T[which(A==max(A,na.rm=TRUE))]

(pj: note na.rm deprecated? I think it is now na.omit)

3.9 行列中の値を変更する

(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)

> tmat <- matrix(rep(0,3*3),ncol=3)
> tmat
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0
> tmat[tmat==0]<-1 
> tmat
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1 #(from Jan Goebel)

la がデータフレームなら、データを強制的にマトリックス形式にする必要がある。以下を使えばよい。

la <- as.matrix(la)
la[la==0] <- 1

日本語訳を見てのコメント
as.matrix は不要になっているようだけど?

> la <- data.frame(x=c(1,3,2,0,4), y=c(0,1,0,0,3))
> la
  x y
1 1 0
2 3 1
3 2 0
4 0 0
5 4 3
> class(la)
[1] "data.frame"
> la[la==0] <- 1
> la
  x y
1 1 1
2 3 1
3 2 1
4 1 1
5 4 3

以下を試してみよう

la <- ifelse(la == 0, 1, la)

(from Martyn Plummer)

日本語訳を見てのコメント
la がデータフレームのとき,何もせずいきなり la[la==0] <- 1 とやったら悲惨な結果になるよね

3.10 行列から特定の行を取り除く

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

> x <- matrix(1:10,,2)
> x[x[,1]%in%c(2,3),]
> x[!x[,1]%in%c(2,3),]

(from Peter Malewski)

mat[!(mat$first %in% 713:715),]

(from Peter Dalgaard )

3.11 ある基準を満たす項目の総数を数える

(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)

以前の質問で述べた which() の結果に"length()" を以下のように適用する
length(which(v<7))

または

sum( v<7,na.rm=TRUE)

3.12 相関行列から偏相関係数を計算する

(08/12/2000 Paul Johnson <pauljohn@ukans.edu>)

多変量変数間の相関係数から,二変数以外の変数を制御して二変数間の偏相関係数を計算しなければならないのだけど。

実際問題として,これはきわめて簡明直裁。R が相関係数行列であるとき,

       Rinv<-solve(R)
       D<-diag(1/sqrt(diag(Rinv)))
       P<- -D %*% Rinv %*% D

P の非対角要素は他の変数の影響を除いた(コントロールした)偏相関係数である(なぜこのようなことをしようと思うのかが別の質問であるが)。

(from John Fox )

一般的に,分散・共分散行列の逆行列を求めて,対角成分が1になるように再スケールする。非対角成分は他の変数を制御した偏相関係数の符号を替えたものになっている。

pcor2 <- function(x){
  conc <- solve(var(x))
  resid.sd <- 1/sqrt(diag(conc))
  pcc <- - sweep(sweep(conc, 1, resid.sd, "*"), 2, resid.sd, "*")
  return(pcc)
}
pcor2(cbind(x1,x2,x3))

see J. Whittaker's book "Graphical models in applied multivariate statistics" from Martyn Plummer)

以下は,私が今使っているもので,それぞれの偏相関係数の有意検定(H0: coeff=0)も同時に行う。

f.parcor <-  function (x, test = F, p = 0.05)
{
   nvar <- ncol(x)
   ndata <- nrow(x)
   conc <- solve(cor(x))
   resid.sd <- 1/sqrt(diag(conc))
   pcc <- -sweep(sweep(conc, 1, resid.sd, "*"), 2, resid.sd,
       "*")
   colnames(pcc) <- rownames(pcc) <- colnames(x)
   if (test) {
       t.df <- ndata - nvar
       t <- pcc/sqrt((1 - pcc^2)/t.df)
       pcc <- list(coefs = pcc, significant = t > qt(1 - (p/2),
           df = t.df))
   }
   return(pcc)
}

(from Kaspar Pflugshaupt)

3.13 正準変数を計算

(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)

同じ観察対象について,データ行列 A と B があるとする。

cancor <- function(A, B) 
{
  Ap <- prcomp(scale(A, T, F), retx=T)
  Apc <- Ap$x %*% diag(1/Ap$sdev)
  Bp <- prcomp(scale(B, T, F), retx=T)
  Bpc <- Bp$x %*% diag(1/Bp$sdev)
  Sigma <- crossprod(Apc, Bpc)/(nrow(A) - 1)
  s <- svd(Sigma, ncol(A), ncol(B))
  return(list(cor=s$d, canvar.x=Apc %*% s$u, canvar.y=Bpc %*% s$v))
}

この関数には,トリックがある。正準変量の平均は0で,分散は1である(S-PLUSとは異なる)(from Brian D. Ripley)

3.14 多次元行列(R 配列)をつくり出す

(20/06/2001 ? <?>)

Brian Ripley の話では

my.array<-array(0,dim=c(10,5,6,8))

上のようにすれば、 10 x 5 x 6 x 8 の 4次元の配列が与えられる。

または以下のような方法もある。

array.test <- array(1:64,c(4,4,4))
array.test[1,1,1]
1
array.test[4,4,4]
64

3.15 多数の行列を結合する

(20/06/2001 ? <?>)

100個の n×2 行列のリストがあるとする。どうやったら,これを大きなN×2行列にできるか。

list は行列のリストである。

nr <- sapply(list, nrow)
cs <- cumsum(nr)
st <- c(0, cs[-length(cs)]) + 1
res <- matrix(, sum(nr), 2)
for(i in seq(along=nr)) res[st[i]:cs[i],] <- list[[i]]

(from Brian Ripley)

3.16 特定の論理に基づき "近傍" 行列をつくり出す

(20/06/2001 ? <?>)

あるセルがある位置に近傍があるかどうかの有無を0と1で示したマトリックスが欲しい。

N <- 3
x <- matrix(1:(N^2),nrow=N,ncol=N)
rowdiff <- function(y,z,mat)abs(row(mat)[y]-row(mat)[z])
coldiff <- function(y,z,mat)abs(col(mat)[y]-col(mat)[z])
rook.case <- function(y,z,mat){coldiff(y,z,mat)+rowdiff(y,z,mat)==1}
bishop.case <- function(y,z,mat){coldiff(y,z,mat)==1 & rowdiff(y,z,mat)==1}
queen.case <- function(y,z,mat){rook.case(y,z,mat) | bishop.case(y,z,mat)}
matrix(as.numeric(
            sapply(x,function(y)sapply(x,rook.case,y,mat=x))),ncol=N^2,nrow=N^2)
matrix(as.numeric(
            sapply(x,function(y)sapply(x,bishop.case,y,mat=x))),ncol=N^2,nrow=N^2)

matrix(as.numeric(

            sapply(x,function(y)sapply(x,queen.case,y,mat=x))),ncol=N^2,nrow=N^2)

(from Ben Bolker)

3.17 ある "key" により二つの数値列をマッチング

(20/06/2001 ? <?>)

The question was:

I have a matrix of predictions from an proportional odds model (using the polr function in MASS), so the columns are the probabilities of the responses, and the rows are the data points. I have another column with the observed responses, and I want to extract the probabilities for the observed responses.

小さな例として,以下の x, y

> x <- matrix(c(1,2,3,4,5,6),2,3)
> y <- c(1,3)

and I want to extract the numbers in x[1,1] and x[2,3] (the columns being indexed from y), what do I do?

Is

x[cbind(seq(along=y), y)]

what you had in mind? The key is definitely matrix indexing. (from Brian Ripley)

3.18 上・下三角行列をつくり出す

(20/06/2001 ? <?>)

これにはたくさんのやり方がある。r-help に投稿されるたびに,異なった解法が示される。

以下のような行列を考えよう。

a^0  0    0
a^1  a^0  0
a^2  a^1  a^0

なぜこれが必要かはわからないが,引き出された多数の解の相違を見てみよう。

 ma <- matrix(rep(c(a^(0:n), 0), n+1), nrow=n+1, ncol=n+1)
 ma[upper.tri(ma)] <- 0

(from Uwe Ligges)

> n <- 3
> a <- 2
> out <- diag(n+1)
> out[lower.tri(out)] <-
       a^apply(matrix(1:n,ncol=1),1,function(x)c(rep(0,x),1:(n-x+1)))[lower.tri(out)]
> out

(from Woodrow Setzer)

tmpmat <- function(a,n) {
  x <- matrix(a,nrow=n,ncol=n)
  x <- x^pmax(0,row(x)-col(x))
  x[row(x)<- 0
  x
}

(from Ben Bolker)

これらを見て,私も作ってみた

> out <- matrix(0, n+1, n+1)
> out[lower.tri(out, diag=FALSE)] <-  a ^ sequence(n:1)

(青木繁伸)

上三角行列がほしいなら,以下に示す ifelse 関数の使い方はすばらしいと思う。

fun <- function(x,y) { ifelse(y>x, x+y, 0) }

(そしてこれを outer() とともに使うか) または

m <- outer(x,y,FUN="+")
m[lower.tri(m, diag=T)] <- 0

(from Kaspar Pflugshaupt)

3.19 X の逆行列を計算

solve(A)は A の逆行列を返す

3.20 行列添字の面白い使い方

(20/06/2001 ? <?>)

Here's a similar problem.

Mehdi Ghafariyan said "I have two vectors A=c(5,2,2,3,3,2) and B=c(2,3,4,5,6,1,3,2,4,3,1,5,1,4,6,1,4) and I want to make the following matrix using the information I have from the above vectors.

     0 1 1 1 1 1 
     1 0 1 0 0 0 
     0 1 0 1 0 0 
     1 0 1 0 1 0
     1 0 0 1 0 1 
     1 0 0 1 0 0 

so the first vector says that I have 6 elements therefore I have to make a 6 by 6 matrix and then I have to read 5 elements from the second vector , and put 1 in [1,j] where j=2,3,4,5,6 and put zero elsewhere( i.e. in [1,1]) and so on. Any idea how this can be done in R ?

Use matrix indices:

a<-c(5,2,2,3,3,2) 
b<-c(2,3,4,5,6,1,3,2,4,3,1,5,1,4,6,1,4)
n<-length(a)
M<-matrix(0,n,n)
M[cbind(rep(1:n,a),b)]<-1

(from Peter Dalgaard )

3.21 固有値の例

(20/06/2001 ? <?>)

Tapio Nummi asked about double-checking results of spectral decomposition

In what follows, m is this matrix:

 0.4015427  0.08903581 -0.2304132
 0.08903581 1.60844812  0.9061157
-0.23041322 0.9061157   2.9692562

Brian Ripley posted:

> sm <- eigen(m, sym=TRUE)
> sm
$values
[1] 3.4311626 1.1970027 0.3510817
$vectors 
[,1] [,2] [,3]
[1,] -0.05508142 -0.2204659 0.9738382 
[2,] 0.44231784  -0.8797867 -0.1741557 
[3,] 0.89516533 0.4211533 0.1459759

> V <- sm$vectors 
> t(V) %*% V 
[,1] [,2] [,3] 
[1,] 1.000000e+00 -1.665335e-16 -5.551115e-17 
[2,] -1.665335e-16 1.000000e+00 2.428613e-16 
[3,] -5.551115e-17 2.428613e-16 1.000000e+00
 
> V %*% diag(sm$values) %*% t(V) 
[,1] [,2] [,3]
 [1,] 0.40154270 0.08903581 -0.2304132 
[2,] 0.08903581 1.60844812 0.9061157 
[3,] -0.23041320 0.90611570 2.9692562

4 関数 tapply の適用等

4.1 一つの関数から複数の値を返す

(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)

Below is a summary on returning variable from subfunction, Thanks to Brian Ripley and Douglas Bates:

To access the variables within a function, return a list or data structure assigned to a variable in the calling function. For more than one varaible, assign them to a list variable and separate the components in the calling function using listname$variablename.

"test" <- function() {
   "hello" <-
     function()    
       {
         x <- 1:3
         y <- 23:34
         z <- cbind(x,y)
         return(x,y,z)              
       }
   c <- hello()
   return(c$z,c$x, c$y)
 }

Sam McClatchie

4.2 有意性検定のリストから "p" 値を掴み出す

(22/08/2000 Paul Johnson <pauljohn@ukans.edu>)

chisq1M11 が htest オブジェクトのリストだとする。そして,その中から P 値を取り出したい。Kjetil Kjernsmo は以下のプログラムが動くことを確認した。

> apply(cbind(1:1000), 1, function(i) chisq1M11[[i]]$p.value)

Peter Dalgaard はもっとエレガントなやり方を示した。

sapply(chisq1M11,function(x)x$p.value)

このどちらも,"function" と呼ばれる関数を作り,それぞれの項目に適用する。 (訳注:本当は無名の関数を作っているのだが)

4.3 ifelse の用法

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

y が 0 から 1 の範囲にあるときにのみある計算をしたいときには

ifelse(y==0,0,y*log(y))

訳注:上のは,y が 0 であるときとそれ以外で別の計算をする例だ

(from Ben Bolker)

4.4 apply を用いて各セル毎の確率の行列を作る

(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)

ガンマ関数のスケールパラメータと形状パラメータの様々な値の組み合わせでガンマ関数の密度関数を計算したいとしよう。scales と shapes というベクトルを作り expand.grid によりグリッドを作る。次いで,関数を書き,適用する(この例は,Jim Robison-Cox の好意による)

gridS <- expand.grid(scales, shapes)
survLilel <- function(ss) sum(dgamma(Survival,ss[1],ss[2]))
Likel <- apply(gridS,1,survLilel)

実際には私は以下のように使う

sc <- 8:12; sh <- 7:12
args <- expand.grid(scale=sc, shape=sh)
matrix(apply(args, 1, function(x) sum(dgamma(Survival, scale=x[1], shape=x[2], 
   log=T))), length(sc), dimnames=list(scale=sc, shape=sh))

(from Brian Ripley)

4.5 外積 outer

(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)

outer(x,y,f) does just one call to f with arguments created by stacking x and y together in the right way, so f has to be vectorised. (from Thomas Lumley)

4.6 あるものが公式か関数か検査する

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

公式(モデル式)は "formula" クラスを持つので,inherits(obj, "formula") が一番良さそうだ。(from Prof Brian D Ripley)

もし,Y~X+Z+W ではなくて ~X+Z+W であることが分かっているならば,以下が使える。

inherits(obj,"formula") && (length(obj)==2)

(from Thomas Lumley)

4.7 変数のベクトルに関する最適化

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

The function being optimized has to be a function of a single argument. If alf and bet are both scalars you can combine them into a vector and use

opt.func <- function(arg) 
  t(Y-(X[,1] * arg[1] + X[,2] * arg[2])^delta) %*% covariance.matrix.inverse %*%
           (Y-(X[,1] * arg[1] + X[,2] * arg[2])^delta)

(from Douglas Bates)

4.8 S+ のような slice.index

(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)

slice.index<-function(x,i){k<-dim(x)[i];sweep(x,i,1:k,function(x,y)y)}

(from Peter Dalgaard)

5 グラフ

5.1 グラフ描きの前に par で特性を変更

(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)

par() はグラフの品質に影響を与える多目的コマンド。例を以下に含める。

1. 以降のプロットを一つの絵にまとめる:

par(mfrow=c(3,2))

プロットの 3x2 のマトリックスを作成する。

2. 余白を次の例のように変更

par(mfrow=c(2,1))
par(mar=c(10, 5, 4, 2))
x <- 1:10
plot(x)
box("figure", lty="dashed")
par(mar=c(5, 5, 10, 2))
plot(x)
box("figure", lty="dashed")

5.2 グラフ出力を保存

(15/08/2001 Paul Johnson <pauljohn@ukans.edu>)

グラフ保存には2つの方法がある。グラフの保存前に、出力形式を決めることができる。プロットの前に、スクリーン上に表示する。Unixでは、デフォルトデバイス、x11を使用して、モニターまたは「スクリーン」を参照する。MS Windows では、デフォルトのスクリーンデバイスは「windows」と呼んでいる。Windows の R は、 x11() をデバイスとして使えば、同じデバイスを使用することになる。

There are other devices representing kinds of file output. You can choose among device types, postscript, png, jpeg, bitmap, etc. (windows users can choose windows metafile). Type

>?device

Rがシステムで見つけるデバイスのリストが見れる。

各デバイスについて、より多くを見つけられる。

?x11

または

?png

また MS Windowsでは

?windows

(等々)

Before you start using a device, you can configure it. The png, jpeg, and bitmap are "picture" formats, they take a "snapshot", and if you want rescale them, it can get ugly with image processing software. Postscript, pdf, xfig, (and windows metafile) are scalable vector graphics and can be edited somewhat more gracefully. But you can't put them on a web page (or, at least, they won't show anything :)).

You start a device with a configuration command, like

png()

for the default png (a gif-like) format.

After you create a device, check what you have on:

> dev.list()

And, if you have more than one device active, you tell which one to use with dev.set(n), for the n device. Suppose we want the third device:

> dev.set(3)

Now here is the big thing. The devices are configured differently. The x11, postscript, and pdf devices are configured in inches, that's inches of "screen space" or the equivalent. The other devices, like png and jpeg, are configured in pixels. The number of pixels per inch as shown on your monitor will depend on your computer. You can also set the point size. Peter Dalgaard told me that on postscript devices, a point is 1/72nd of an inch. If you expect to shrink an image for final presentation, set a big pointsize, like 20 or so, in order to avoid having really small type.

Suppose we want png output. We create an instance of a png device and adjust features for that device, as in:

> png(filename="mypicture.png", width=480, height=640, pointsize=12)

Note, if no filename is here, it will print to the default printer. You can then run your graph, and when it is finished you have to close the device:

> dev.off()

Some print devices can accept multiple plots, and dev.off() is needed to tell them when to stop recording. On a single-plot device, if you try to run another graph before closing this one, you'll get an error.

The postscript device has options to output multiple graphs, one after the other. If you want to set options to the postscript device, but don't actually want to open it at the moment, use ps.options.

Here's an example of how to make a jpeg:

jpeg(filename="plot.jpg")
plot(sin, 2*pi)
dev.off()

For me, the problem with this approach is that I don't usually know what I want to print until I see it. If I am using the jpeg device, there is no screen output, no picture to look at. So I have to make the plot, see what I want, turn on an output device and run the plot all over in order to save a file. It seems complicated, anyway.

So what is the alternative?

If you already have a graph on the screen, and did not prepare a device, use dev.copy() or dev.print(). This may not work great if the size of your display does not match the target device.

dev.copy() opens the device you specify. It is as if (but not exactly like) you ran the png() command before you made the graph:

dev.copy(device=png, file="foo", width=500, height=300)
dev.off()

The dev.copy() function takes a device type as the first argument, and then it can take any arguments that would ordinarily be intended for that device, and then it can take other options as well. Be careful, some devices want measurements in inches, while others want them in pixels. Note that, if you do not specify a device, dev.copy() expects a which= option to be specified telling it which pre-existing device to use.

If you forget the last line, it won't work. It's like a write command.

If you use dev.copy or dev.print, you may run into the problem that your graph elements have to be resized and they don't fit together the way you expect. The default x11 size is 7in x 7in, but postscript size is 1/2 inch smaller than the usable papersize. That mismatch means that either you should set the size of the graphics window on the monitor display device to match the eventual output device, or you fiddle the dev.copy() or dev.print() command to make sure the sizes are correct. Recently I was puzzled over this and Peter Dalgaard said you can force the sizes to match, as in

1. Force the pdf output to match the size of the monitor:

dev.copy(pdf, file="foo", height=7, width=7)

or
2. change the display size before creating the graphs so that its size matches the intended device you might use for output:

x11(height=6, width=6)

There are some specialized functions that can do this in a single step, as in dev.copy2eps() or dev2bitmap(). dev.copy2eps creates an eps file, and I'm not sure how it is different from the eps-compatible output created by the postscript() device. dev.copy2eps *is* an indirect way of calling dev.copy(device=postscript,...). It takes the dimensions from the current plot, and sets a couple of options. The same is true of dev.print().

The dev.print() function, as far as I can see, is basically the same as dev.copy(), except it has two appealing features. First, the graph is rescaled according to paper dimensions, and the fonts are rescaled accordingly. Due to the problem mentioned above, not everything gets rescaled perfectly, however, so take care. Second, it automatically turns off the device after it has printed/saved its result.

dev.print(device=postscript,file="yourFileName.ps")

Before you run the graphs, you can do use the postscript() command to start a PostScript device that dumps the graphs you make onto disk.

If you want pdf output, as of R 1.3 there is a pdf device. Just use pdf() like postscript() or whatever. If you have an old R (say it ain't so...), look into the bitmap device: bitmap() does the postscript + epstopdf bit for you if you select type = "pdfwrite". (from Brian Ripley)

If you just need to test your computer setup, Bill Simpson offered this. Here is a plot printing example

x<-1:10
y<-x
plot(x,y)
dev.print(width=5,height=5, horizontal=FALSE)

A default jpeg is 480x480, but you can change that:

jpeg(filename="plot.jpg" width = 460, height = 480, pointsize = 12, quality = 85)

The format of png use is the same.

As of R1.1, the dev.print() and dev.copy2eps() will work when called from a function, for example:

> ps <- function(file="Rplot.eps", width=7, height=7, ...) {
       dev.copy2eps(file=file, width=width, height=height, ...)
  }
> data(cars)
> plot(cars)
> ps()

That didn't work before, but it does work now.

In r-help's list, I've been one of the countless new users to be stumped by the problem of copying a screen device to a file format. At the current time, on Unix, all I can say is that either you resize the display device before doing the graphs so it matches (or is smaller than) the output device, or be careful that the output command specifies width and height (and probably pointsize) to match the screen device. Even then it may not be perfect.

The mismatch of "size" between devices even comes up when you want to print out an plot. This command will print to a printer:

> dev.print(height=6, width=6, horizontal=FALSE)

You might want to include pointsize=20 or whatever so the text is in proper proportion to the rest of your plot.

One user observed, "Unfortunately this will also make the hash marks too big and put a big gap between the axis labels and the axis title...", and in response Brian Ripley observed: "The problem is your use of dev.print here: the ticks change but not the text size. dev.copy does not use the new pointsize: try

x11(width=3, height=3, pointsize=8)
x11(width=6, height=6, pointsize=16)
dev.set(2)
plot(1:10)
dev.copy()

Re-scaling works as expected for new plots but not re-played plots. Plotting directly on a bigger device is that answer: plots then scale exactly, except for perhaps default line widths and other things where rasterization effects come into play. In short, if you want postscript, use postscript() directly.

The rescaling of a windows() device works differently, and does rescale the fonts. (I don't have anything more to say on that now)

5.3 プロット出力を別個のファイルに保存するには

(10/04/2001 Paul Johnson <pauljohn@ukans.edu>)

postscript(file="test3.%d.ps",onefile=FALSE) 

ファイル名を test3.1.ps, test3.2.ps などと付ける。(from Thomas Lumley)

5.4 用紙サイズをコントロールする

(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)

options(papersize="letter")

whereas ps.options is only for the postscript device

5.5 R のグラフを文章に統合する: Tex と EPS

(20/06/2001 Paul Johnson <pauljohn@ukans.edu>)

Advice seems to be to output R graphs in a scalable vector format, eps on Linux or eps or metafile on windows. If you are using Tex for document prep, on linux use "tetex" and on windows it has an implementation called http://www.tile.net/listserv/fptex.html and another called Miktex www.miktex.org I've installed Miktex and its pretty good.

Keith Reckdahl's Using EPS in LaTeX documents is freely available on the CTAN sites. (Brian Ripley)

For a description of EPSF and various other vector and bitmap formats, see the Encyclopedia of Graphics File Formats from O'Reilly at http://www.oreilly.com/catalog/gffcd/ (from Joel West)

I personally swear by the LaTeX Companion and the LaTeX Graphics companion, both published by Addison Wesley. (from James Marca)

Thanks to Prof. Ripley for the hint, maybe someone other is interested, so I will shortly outline what "to use pstricks with LaTeX" means:

-get PSfrag (e.g. www.dante.de)
- get graphics (same url)
- create your ps-graphic in R and define "tags" e.g.
       
       postscript("myfile.ps")
       plot(..., xlab="rho", ylab="wrho")
       dev.off()
       system("ps2epsi myfile.ps")
add
       ?usepackage{graphics}
       ?usepackage{psfrag}

to your LaTeX-File

and now we reach the climax:

       ?begin{figure}
       ?begin{center}
       ?psfrag{rho}{$?varrho$}         <-- this replaces rho with ?varrho
       ?psfrag{wrho}{$W(?varrho)$}     <-- this replaces wrho with W(?varrho)
       ?includegraphics{myfile.epsi}
       ?caption{wonderful R plotting ?label{fig}}
       ?end{center}
       ?end{figure}

run latex and dvips and enjoy :-) (from Torsten Hothorn)

If you need to create a preview inside an eps file, you can do it with GSView (available free).

5.6 "Snapshot" グラフとそれらの閲覧

(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)

Ordinarily, the graphics device throws away the old graph to make room for the new one.

Have a look to ?recordPlot to see how to keep snapshots. (Doo myGraph<-recordPlot() to save snap, then replayPlot(myGraph) to see it again);

On Windows you can choose "history"-"recording" in the menu and scroll through the graphs using the PageUp/Down keys. (from Uwe Ligges)

5.7 scatterplot の行列の上三角部分を選択する

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

必要なのはプロットする点を指定するだけ:

## upper triangle of the correct size:
temp <- upper.tri(mmtop94.2, diag = TRUE) 
z <- mmtop94.2[temp]  # values
y <- nrow(mmtop94.2) + 1 - row(mmtop94.2)[temp]  # row indices
x <- col(mmtop94.2)[temp]  # col indices
scatterplot3d(x, y, z, type="h") 

(from Uwe Ligges)

5.8 密度関数をプロットする (例えば正規分布)

(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)

x<-seq(-3,3,.1) plot(x,dnorm(x,0,1), type="l")

(Bill Simpson)

5.9 エラーバー付きでプロット

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

以下にやり方を示す。

x<-c(1,2,3,4,5)
y<-c(1.1, 2.3, 3.0, 3.9, 5.1)
ucl<-c(1.3, 2.4, 3.5, 4.1, 5.3)
lcl<-c(.9, 1.8, 2.7, 3.8, 5.0)
plot(x,y, ylim=range(c(lcl,ucl)))
arrows(x,ucl,x,lcl,length=.05,angle=90,code=3) 

または

segments(x,ucl,x,lcl)

(from Bill Simpson)

5.10 推定密度関数付きのヒストグラム

(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)

同じスケールで密度推定とヒストグラムが欲しければ、以下のようにしてみればよい。

> IQR <- diff(summary(data)[c(5,2)])
> dest <- density(data, width = 2*IQR)  # or some smaller width, maybe,
> hist(data, xlim = range(dest$x), xlab = "x", ylab = "density",
         probability = TRUE)    # <<<--- this is the vital argument
> lines(dest, lty=2)

(from Bill Venables)

diffはLagged Differencesを返す関数ですので、上記ではIQRが負になりエラーが出ます。

IQR <- diff(summary(data)[c(2,5)])

とすべきです。又、欠損値(NA)が含まれているとエラーが出ますので、

dest <- density(data, width = 2*IQR,na=T)

としてNAを排除してください。

5.11 いくつもの線プロットを重ね描きするには?

(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)

前の出力がかき消されないようにするためには,以下のように par("new"=TRUE) を使う。

tmp1<-plot(acquaint~T,type='l', ylim=c(0,1),ylab="average
   proportion",xlab="PERIOD",lty=1,pch=1,main="")
par("new"=TRUE)
tmp2<-plot(harmony~T,type='l', ylim=c(0,1),ylab="average
proportion",xlab="PERIOD",lty=2,pch=1,main="")
par("new"=TRUE)
tmp3<-plot(identical~T,type='l', ylim=c(0,1),ylab="average
proportion",xlab="PERIOD",lty=3,pch=1,main="")

par() は,高水準のプロット命令を重ねるものであるが,abline()のような低水準プロット命令は既に描画されているものに付け加えられるので par() は必要ではない。

Here's another idea: form a new data.frame and pass it through as y:

plot(cbind(ts1.ts,ts2.ts),xlab="time", ylab="", plot.type="single")

以下のようにすれば,もう少しまし。

plot(cbind(ts1.ts,ts2.ts),plot.type="single",
      col=c("yellow","red"),lty=c("solid","dotted")) #colours and patterns

'c(ts1.ts,ts2.ts)' と 'cbind(ts1.ts,ts2.ts)' を比較するとためになるだろう。

(from Guido Masarotto)

時系列データについては特別な関数がある。

ts.plot(ts1.ts, ts2.ts) # same as above

plot.ts と ts.plot のヘルプを見よ。 (from Brian D. Ripley)

以下は,ヒストグラムの例である。

> hist(rnorm(100, mean = 20, sd =12), xlim=range(0,100), ylim=range(0,50))
> par(new = TRUE)
> hist(rnorm(100, mean = 88, sd = 2), xlim=range(0,100), ylim=range(0,50))

横軸ラベルがめちゃくちゃになる。これを避けるには,xlab="" オプションを加える。 (訳注:タイトルもめちゃくちゃになるので,main="" も加える。しかし,二番目の xlab と mai は反映されないので,根本的な解決ではない)

散布図を書くときにうまい方法がある。 (蛇足の訳注:xlim, ylim で xlim=c(min(x1,x2),max(x1,x2)) 等としなくても range を使えばよいと言うこと)

> x1<-1:10
> x2<-2:12
> y1<-x1+rnorm(length(x1))
> y2<-.1*x2+rnorm(length(x2))
> plot(x1,y1,xlim=range(x1,x2), ylim=range(y1,y2), pch=1)
> points(x2,y2,pch=2)

(from Bill Simpson)

5.12 グラフを行列状に配置する

(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)

ページに複数のグラフを配置する標準的な方法は,以下のように par(mfrow=) または par(mfcol=) オプションを使うことである。

par(mfrow=c(2,3))

そうすれば,次の 6 つのグラフは2行3列に配置される

The layout option can give a powerful set of controls for that kind of thing as well, e.g. (from Paul Murrell):

layout(matrix(c(1, 0, 2), ncol=1),
          heights=c(1, lcm(5), 1))
plot(x)
box("figure", lty="dashed")
plot(x)
box("figure", lty="dashed")

I think you have a data with a like:

xx <- data.frame(y=rnorm(100), x1=as.factor(round(runif(100,1,4))),  
          x2=as.factor(round(runif(100,1,4))) ) attach(xx) by(y,list(x1,x2),plot) 
by(y,list(x1,x2),function(x)print(x)) xn <- as.integer(x1) 
## either use coplot or  
par(mfrow=) approach: coplot(y~xn|x2) 
##or you use a loop: 
par(mfrow=c(2,2)) for( i in unique(as.integer(x1))) 
  plot.default(x2 [as.integer(x1)== i] ,  y[as.integer(x1)==i ] ,
 main=paste("Code:",i) ) 
##of course there might be ways to use tapply, lapply etc.

(from Peter Malewski. Per Peter Dalgaard's advice, I've replaced usage of codes with as.integer. pj)

And to label the whole page, use the "mtext" function

5.13 線と棒グラフを結合する

(07/12/2000 Paul Johnson <pauljohn@ukans.edu>)

David James was kind enough to help me out and enlighten me to the fact that the x-scales used by barplot are independent of those from my data. He also sent me a function which I include below (with a minor modification). This does exactly what I was looking for. Thanks!

xbarplot <- function(y, col=par("col"), border = par("fg"), gap=gap, ...)
{
  ny <- length(y)
  x <- seq(0,ny)+0.5
  n <- length(x)
  space <- gap * mean(diff(x))
  old <- par(xaxt = "n")
  on.exit(par(old))
  plot(range(x, na.rm=T), range(y, na.rm=T), bty="n",
     xlab="", ylab = "", type = "n", ...)
  rect(x[-n]+space/2, rep(0,length(y)), 
      x[-1]-space/2, y, col = col, border = border)
}

(from Dirk Eddelbuettel)

5.14 回帰散布図: グラフに当てはめ回帰直線を加える

(17/08/2001 Paul Johnson <pauljohn@ukans.edu>)

suppose you have

reg <-lm( y ~ x)
plot(x,y)

then do

abline(reg)

5.15 散布図中のプロット文字を制御する

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

Plots show dots:

plot(rnorm(100),rnorm(100),pch=".") (from Peter Malewski)

Plots with numbers indicating values of a third variable:

text(tx, ty, label = ch, col = "red", bg = "yellow", cex = 3)

If you specify pch, only the first character is taken as the symbol for your point.

Matt Wiener suggested creating a color vector and then using the value of another value to control both the lable and the color:

> col.vec <- c("black", "red", "blue")
> text(chem.predict[,1:2], labels = km$cluster, color =
         col.vec[km$clusters])

5.16 散布図: プロット文字を制御する (men vs women, 等)

(11/11/2002 Paul Johnson <pauljohn@ku.edu>)

pch は,プロット記号を選択するために整数を用いることを覚えておこう。

x <- 1:10
y <- 1:10
res <- -4:5
plot(x, y, pch = ifelse(res < 0, 20, 1))

もし論理式が真であれば,記号20が,さもなくば記号1が使われる。色々な数を入れて,実験してみよ。

Then note you can create a vector of integers to control symbols, as in

g <- 15:24plot (x,y, pch= g)

(From post by Uwe Ligges (11/7/2002)

Here's another cool example from Roger Bivand (11/7/2002),

> x <- rnorm(25)
> y <- rnorm(25)
> z <- rnorm(25)
> pchs <- c(1, 20)
> pchs[(z < 0)+1]
[1] 20 20 20 20 1 1 20 20 1 20 20 20 20 20 1 20 1 ...
> plot(x, y, pch=pchs[(z < 0)+1])
> text(x, y, labels=round(z, 2), pos=1, offset=0.5)

Roger adds, "This "cuts" the z vector at zero, using the convenient slight-of-hand that TRUE and FALSE map to integers 1 and 0, and thus gives the pch argument in plot() or points() a vector of values of indices of the pchs vector. More generally, use cut() to break a numeric vector into class intervals (possibly within ordered())."

5.17 サイズと色を修正した散布図

(12/11/2002 Paul Johnson <pauljohn@ku.edu>)

test<-c(2,6,4,7,5,6,8,3,7,2) 
plot(test,type="n",main="Color/size test plot",ylab="Size scale",xlab="Colorscale") 
colsiz<-function (yvec) { 
  ymin <- min(yvec) 
 cexdiv <- max(yvec)/2 for (i in 1:length(yvec)) {
    nextcex <- (yvec[i] - ymin)/cexdiv + 1   
 par(cex = nextcex) 
    points(i, yvec[i], type = "p", col = i) 
 }
}
colsiz(test)

Note that the size here ranges from 1->2. You can change that by fiddling withthe calculation of 'cexdiv'. (from Jim Lemon)

Paul Murrell posted this nice example (11/7/2002):

x <- rnorm(50)
y <- rnorm(50)
z <- rnorm(50)
pch <- rep(1, 50)
pch[z < 0] <- 20
cex <- (abs(z)/max(abs(z))) * 4.9 + 0.1
plot(x, y, pch=pch, cex=cex)

5.18 散布図:三番目の変数に関してサイズを調整する

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

Have a look at:

?symbols

e.g.:

symbols(1:10,1:10,circles=1:10, inches=0.1)

(from Uwe Ligges)

5.19 散布図: 点を結ぶ線を滑らかにする

(02/06/2003 Paul Johnson <pauljohn@ku.edu>)

Rado Bonk がいくつかの点をなめらかに結ぶ曲線を求めている。一つの答えは「spline 関数はいかがですか?」例は以下のようになる

x <- 1:5
y <- c(1,3,4, 2.5,2)
plot(x, y)
sp <- spline(x, y, n = 50)
lines(sp)

Roger Peng

5.20 回帰散布図: 推定値をプロットに加える

(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)

I don't know how to get the estimate of the line printed, but I have seen one example of how to get the R-square:

I used this to indicate the R-squared value on a regression plot:

text(5, 0.6, as.expression(substitute(R^2 == r, list(r=round(R2.the,3)))))

where R2.the was computed a few lines above. (from Emmanuel Paradis)

Or

that <- 1
plot(1:10)
title(substitute(hat(theta) == that, list(that=that)))

(from Brian Ripley)

5.21 軸の制御: チックマーク、チックマークをなくす、数、等

(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)

チックマークは表示するが、チックマークの数値ラベルは表示しない

 plot(1:10, xaxt="n")
 axis(1, labels=FALSE)

X軸およびY軸のタイトルを省くが、各チックの数値ラベルは残す。

plot(x,y,ann=FALSE)

もっと多くのグラフの制御をしたいなら、以下のようなことができる。

plot(x,y,ann=FALSE,axes=FALSE)
box()  
axis(side=1,labels=TRUE)
axis(side=2,labels=TRUE)
mtext("Y", side=2, line=2)
mtext("X", side=1, line=2)

? axis と入力してみよ (from Bill Simpson)

Control range of axes through the graph command itself. This shows the "image" function:

> x <- 1:50  
> y <- 1:50
> z <- outer(x,y)
> image(z)
> image(z, xlim=c(0.2,0.5))
> image(z, xlim=c(0.2,0.5), ylim=c(0,0.5))

The same works with contour(). (from Kaspar Pflugshaupt)

5.22 軸: ラベルの回転

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

一度で、par(las=) は見かけ上これをする。入力は以下の通り。

par(las=2) or par(las=3)

上では軸のラベルを90度回転する。

More recently, Paul Murrell has observed:

"The drawing of xlab and ylab ignores the par(las) setting.

...You can override this "sensible" default behaviour if it annoys you. For example, you can stop plot() drawing xlab and ylab by something like plot(ann=F) or plot(xlab="", ylab="") AND you can draw them yourself using mtext(), which does listen to par(las).

A couple of examples of what I mean ...

   par(mfrow=c(2,2), mar=c(5.1, 4.1, 4.1, 2.1), las=0)
   plot(0:10, 0:10, type="n")
   text(5, 5, "The default axes")
   box("figure", lty="dashed")
   par(las=1)
   plot(0:10, 0:10, type="n")
   text(5, 5, "The user says horizontal text please?n?nbut R knows better !")
   box("figure", lty="dashed")
   par(las=1, mar=c(5.1, 6.1, 4.1, 2.1))
   plot(0:10, 0:10, type="n", ann=F)
   mtext("0:10", side=2, line=3)
   mtext("0:10", side=1, line=3)
   text(5, 5, "The user fights back !")

(posted 6/11/2001)

5.23 軸: 軸に整形した日付を加える

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

CRANで入手できる"date" パッケージを試してみる。

作図例は以下のとおり。

library(date)
plot.date(.....) 
## if x-values are dates, you can also use plot(...) after
library(date)

(from Uwe Ligges)

5.24 軸: プロット中の軸を逆転する

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

R1.1.1 にはバグがあるが、後のバージョンでは、Ross Ihaka は次のようにいっている。

x <- rnorm(100)
y <- rnorm(100)
plot(x, y, xlim=rev(range(x)), ylim=rev(range(y)))

and get the result you expect.

Ben Bolker said, at this time (R1.1.1), you have to do something like

x <- rnorm(100)
y <- rnorm(100)
plot(max(x)-x,max(y)-y,axes=FALSE)
xlabvals <- pretty(x,5)
ylabvals <- pretty(y,5)
axis(side=1,at=max(x)-xlabvals,labels=xlabvals)
axis(side=2,at=max(y)-ylabvals,labels=ylabvals)

I think it's time to put this one in the FAQ ... (although, again, I don't know if I would have found it just searching the archives if I hadn't remembered that it was there in the first place).

Another person said try (if you don't want axes):

plot(x, -y, axes=FALSE)
axis(1)
axis(2)

5.25 軸: 日付によるラベル付き軸

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

以下のようにデータセットをプロットする必要がある。

> x <- c("05/27/00","06/03/00","08/22/00",...)
> y <- c(10,20,15,...)

"date"パッケージを試すべきだ。このパッケージは CRAN で入手できる。以下に例を示す。

library(date)
x <- as.date(c("05/27/2000","06/03/2000","08/22/2000"))
y <- c(10,20,15)
plot(x,y)

Uwe Ligges

これが chron を利用したものだ。

library(chron)
x <- dates(c("05/27/00","06/03/00","08/22/00"))
y <- c(10,20,15)
plot(x, y)

これは貴方が合衆国にいて最初の文字列が 2000年5月27日、すなわち,ISO 8601 の日付フォーマットで 2000-05-27 を意味しているのであれば動くだろう。(R 1.2.0 はこれらと合衆国およびヨーロッパの簡略表現を容易に読み込める。)(Brian D. Ripley より)

5.26 軸: 軸ラベルの上付き添字

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

plot(c(1:10), ylab= expression(x^"++")) (from Peter Dalgaard)

5.27 軸: 位置の補正

(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)

質問:プロットするとx軸はグラフの下に,y軸はグラフの左に描かれるが,原点で交差するx軸とy軸を描きたい。

以下のようにすればよいでしょう

plot(...., xaxt="n", yaxt="n")
axis(1, pos=0)
axis(2, pos=0)

Uwe Ligges

訳注:このままだと,x,y 軸のラベルは元の位置に書かれてしまう。

plot(...., xaxt="n", yaxt="n", xlab="", ylab="")

として,

axis(1, at=ほげほげ, labels=モーモー, pos=0)

等とすべし

5.28 散布図に "エラー矢印" を加える

(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)

arrows(..., code=3, angle=90, ...) を使えば、まったく簡単にできる、たとえば以下のように。

       x <- rnorm(10)
       y <- rnorm(10)
       se.x <- runif(10)
       se.y <- runif(10)
       plot(x, y, pch=22)
       arrows(x, y-se.y, x, y+se.y, code=3, angle=90, length=0.1)
       arrows(x-se.x, y, x+se.x, y, code=3, angle=90, length=0.1)

最初の arrows() は y のエラーバーで、2番目のは x のエラーバーである。 'code=3' で矢印の両側の頭を描き、'angle=' で矢印の主軸に対する矢印の頭の角度を指定する。また、通常のパラメータ (col, lwd, ...)を追加できる。

Emmanuel Paradis

5.29 時系列: 一つのグラフに複数の線をプロットする

(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)

同じグラフで2つを描画するのは以下のようにする。

plot(x)
lines(y)

または

points(y)

または

matplot(cbind(x,y),type="l")

グラフを上下に分割して描画するには以下のようにすればよい。

par(mfrow=c(2,1))
plot(x)
plot(y)

色の設定、範囲、軸のラベルなどもできる。

(from Jon Baron 12/29/2001)

5.30 時系列: 当てはめ値と実際のデータをプロット

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

> library(ts) 
> data(LakeHuron) 
> md <- arima0(LakeHuron, order=c(2,0,0), xreg=1:98)  plot(LakeHuron) 
> lines(LakeHuron-md$resid,col="red")

(from Adrian Trapletti)

5.31 プロットにテキストを挿入

(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)

次の文字列をプロットに入れたい。"The average temperature is 23.1 degrees."

text(x,y,paste("The average temperature is ", variable, "degrees"))

(from Thomas Lumley)

5.32 はみ出る変数をプロット

(07/12/2000 Paul Johnson <pauljohn@ukans.edu>)

Graphs blow up if a variable is unbounded, so if you are plotting tan(x), for example:

You probably want to plot something like

pmax(-10,pmin(10,tan(x))) 

or

ifelse(abs(tan(x))>10,NA,tan(x))

(from Thomas Lumley)

5.33 動的に生成された内容/数式マークアップを持つラベル

(16/08/2000 Paul Johnson <pauljohn@ukans.edu>)

Don Wingate wrote "More generally, what I need to do is dynamically create a vector of math expressions (using the powers of string manipulation), which can then be used as the textual values for a plot function."

Uve Ligges answered: What about something like:

  numbers <- 1:2  # the dynamical part
  my.names <- NULL
  for(i in numbers) 
     my.names <- c(my.names, 
        eval(as.expression(substitute(expression(x[i]^2), list(i=i)))))
  dotplot(1:2, labels = my.names)

5.34 プロットのタイトルに数式等の複雑なものを使う

(11/11/2002 Paul Johnson <pauljohn@ku.edu>)

substitute() を使う。なにかこのようになる。

e<-substitute(expression(paste("Power plot of ", x ^ alpha, " for ",
        alpha == ch.a)), list(ch.a=formatC(alpha,wid=1)))
plot(x, x^alpha, main=e)

(from Peter Dalgaard)

以下のようにしているとすると

> that <- 1
> plot(1:10)

そして、タイトルに "that" を指定する必要がある。

expression(paste(hat(theta),'= ',that)))

Generally, to get an expression either use parse on your pasted character string or substitute on an expression. The neatest is

title(substitute(hat(theta) == that, list(that=that)))

(note it is == not =)

(from Roger Kroenker, from Brian Ripley)

Want a variable value in the title, as in "Data at the 45o North:

lat<-45
plot(1,main=substitute("Data at "*lat*degree*" North",list(lat=lat)))

(from Thomas Lumley)

You can use different way: First export plot for R as an (e)ps file and than use LaTeX psfrag package to add text and formulae. (from Jan Krupa)

I want a Greek letter subscripted in the title, as in

[Greek] gamma [subscript]1 = [Roman] Threshold 1

Try:

 plot(1,1,main=expression(gamma[1]=="Threshold 1"))

(from Thomas Lumley)

Try

help(plotmath)

5.35 三番目の変数の欠損状態を示す色分けした点を散布図に使う

(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)

plot(x,y,color=ifelse(is.na(z),"red","black"))  #(from Peter Dalgaard)
plot(x[!is.na(z)],y[!is.na(z)],xlim=range(x),ylim=range(y)) 
points(x[is.na(z)],y[is.na(z)],col=2) (from Ross Darnell)

5.36 lattice: その他の例

(12/11/2002 Paul Johnson <pauljohn@ku.edu>)

The lattice library came along recently, I've not explored it much. But here is a usage example, emaphsizing the aspect setting from Paul Murrell:

library(lattice)
# 10 rows
x <- matrix(rnorm(100), ncol=10)
lp1 <- levelplot(x, colorkey=list(space="bottom"),  aspect=10/10)
print(lp1)
# 5 rows -- each row taller than in previous plot
x <- matrix(rnorm(50), ncol=5)
lp2 <- levelplot(x, colorkey=list(space="bottom"),aspect=5/10)
print(lp2)

This one from Deepayan Sarkar (deepayan@stat.wisc.edu) is especially neat.

library(lattice)
xyplot(rnorm(100) ~ rnorm(100) | gl(2, 50),
      strip = function(factor.levels, ...) {
strip.default(factor.levels = expression(1 * degree, 2 * degree),
 ...)
})

5.37 三次元の散布図

(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)

以下をタイプすると

?persp
?contour

ベースパッケージにあるこれらの特徴を参照できる。

CRAN には、scatterplot3d パッケージがある。

以下も参照のこと:

http://www.statistik.uni-dortmund.de/leute/ligges.htm (from Uwe Ligges)

以下もまた参照のこと:

Get Xgobi (an X windows graphics library) and the R contrib package "xgobi". The web page for Xgobi is http://lib.stat.cmu.edu/general/XGobi/ This does 3-D graphs, rotations, etc. Very nice. Mainly for Unix/X users.

5.38 値を反映した線種を用いた三次元等高線図

(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)

正の値が実線で、負の値が点線のコンターを描画したい場合,

種明かしをしてみると levels で間隔を指定して、 add= 引数を加えてみればよい.

       x <- seq(-1, 1, length=21)
       f <- function(x,y) (x^2 + y^2) / 2 - 0.5
       z <- outer(x, x, f)
       contour(z, levels = seq(-0.5, -0.1, by = 0.1), lty = "dotted")
       contour(z, levels = 0, lty = "dashed", add = TRUE)
       contour(z, levels = seq(0.1, 1, by = 0.1), add = TRUE)

(from Ross Ihaka)

訳を見て:あまりたいした変更ではないが以下のようなものも

       x <- seq(-1, 1, length=21)
       f <- function(x,y) (x^2 + y^2) / 2 - 0.5
       z <- outer(x, x, f)
       k <-seq(-0.5,1, by = 0.1)
       contour(z, levels = k , lty = c("dotted", "dashed", "solid")[sign(k)+2])

このようにしておけば,例えば色も変えられる

       contour(z, levels = k, lty = c("dotted", "dashed", "solid")[sign(k)+2], col=c("red", "black", "blue")[sign(k)+2])

5.39 グラフのアニメーション

(13/08/2000 Paul Johnson <pauljohn@ukans.edu>)

ImageMagick (www.imagemagick.org)を使えば本当に容易にできる。ImageMagick はほとんどのOSにインストールできる。わたしはこのシンプルな手順を試してみたが、ImageMagick は見事に動いた。

Rで、15個のヒストグラムを作成する。

> for(i in 1:5) {
+ jpeg(paste("fig", i, ".jpg", sep = ""))
+ hist(rnorm(100))
+ dev.off()
+ }

それから、コマンドライン(わたしはLinuxで試したが、他のプラットフォームでも同様のはずだ)で以下のようにタイプする。

% convert -delay 50 -loop 50 fig*.jpg animated.gif

これは animated.gif を作成した。これは私の一連のファイルの素晴らしいアニメーションである。 -delay および-loopを使ってアニメーションのタイミングを制御できる(Matthew R. Nelson, Ph.D. より)。

5.40 グラフのユーザー領域の背景色を余白と異なったものにする

(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)

"par" を使って、 boundaries/axes の描画位置を見つけ出し、"rect" を実行して画像を描画する ... demo(graphics) の6番目の作図が同様なことをする。

例えば,

## code for coloring background
x <- runif(10) 
y <- runif(10)
plot(x,y,type="n") 
u <- par("usr") 
rect(u[1],u[3],u[2],u[4],col="magenta") 
points(x,y)

(from Ben Bolker)

(pj: see next because it includes a usage of do.call)

I think one is stuck with things like

plot(2:10,type="n")
do.call("rect",c(as.list(par("usr")[c(1,3,2,4)]),list(col="pink")))
points(2:10)

or for the middle line, somewhat simpler

bb <- par("usr")
rect(bb[1],bb[3],bb[2],bb[4], col="pink")

(from Peter Dalgaard )

5.41 知ってて良かったと思う気の効いたグラフ知識(その他)

(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)

1. example(rug). 'nuff said. 2. jitter

5.42 実際に動くグラフコード幾つか (その他)

(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)

点プロット、直線上書き、凡例、PS 出力

x <- c(1,2,3,4)
y1 <- c(1,2,3,4)
y2 <- c(2,2,2,2)
Fred <- c(1,2)
postscript(file="d:/Bob/Papers/IFM/try2.ps") 
plot(x,y1, type="l") 
lines(x,y2,lty="33")  
legend(1,4, c("y1","y2"), lty=Fred) 
graphics.off()

(from Anon.)

サイン・コサイン関数のグラフ。layout 関数の使用例

layout.my <- function (m, n) {
    par(mfrow = c(m, n))
}
x <- 0:12566 / 1000 # Range from 0 to 4*pi
layout.my( 1, 2 ) 
plot( sin(x), type = "l", xaxs = "i", yaxs = "i", axes = F, 
         xlab = "x", ylab = "sin(x)", main = "Y = sin(x), x = [ 0, 720 ]" ) 
axis( 2, at = seq( -1, 1, by=1 ),las = 2 ) box(lty="dotted") abline( h = 0, lwd = 1 )
plot( cos(x), type = "l", xaxs = "i", yaxs = "i", axes = F, 
         xlab = "x", ylab = "cos(x)", main = "Y = cos(x), x = [ 0, 720 ]" ) 
axis( 2, at = seq( -1, 1, by=1 ),las = 2 ) 
box(lty="dotted") 
abline( h = 0, lwd = 1 )

(from Ko-Kang Wang)

正多角形のプロット

Below is a function that generates regular polygons, filled or borders, of n sides (n>8 => circle), with "diameter" prescribed in mm, for use alone or with apply.

ngon <- function (xydc, n=4, type=1) 
  # draw or fill regular polygon 
  # xydc a four element vector with 
  # x and y of center, d diameter in mm, and c color 
  # n number of sides of polygon, n>8 => circle 
  # if n odd, vertex at (0,y), else midpoint of side 
  # type=1 => interior filled, type=2 => edge 
  # type=3 => both 
  { u <- par("usr") 
    p <- par("pin") 
    d <- as.numeric(xydc[3]) 
    inch <- d/25.4 
    rad <- inch*((u[2]-u[1])/p[1])/2 
    ys <- inch*((u[4]-u[3])/p[2])/2/rad 
    if (n > 8) n <- d*4 + 1 
    th <- pi*2/n 
    costh <- cos (th) 
    sinth <- sin (th) 
    x <- y <- rep (0,n+1) 
    if (n %% 2) { x[1] <- 0 y[1] <- rad } 
       else { x[1] <- -rad*sin(th/2) y[1] <- rad*cos(th/2) }  
    for (i in 2:(n+1)) { 
        xl <- x[i-1] 
        yl <- y[i-1] 
        x[i] <- xl*costh - yl*sinth 
        y[i] <- xl*sinth + yl*costh 
    } 
    x <- x + as.numeric(xydc[1]) 
    y <- y*ys + as.numeric(xydc[2]) 
    if (type %% 2) polygon (x,y,col=xydc[4],border=0) 
    if (type %/% 2) lines (x,y,col=xydc[4]) invisible() }

(from Denis White)

自前の軸

plot(1:1000, rnorm(1000), axes=FALSE) # plot with no axes 
axis(1, at = seq(0, 1000, by = 100)) # custom x axis 
axis(2) # default y axis 
box() # box around the plot

(from Ross Ihaka)

プロットの中心に軸を持つ三次関数を描く

 #  Draw the basic curve (limits were established by trial and error).
 #  The axes are turned off blank axis labels used.
 curve(x * (2 * x - 3) * (2 * x + 3), from = -1.85, to = 1.85, 
   xlim = c(-2, 2), ylim = c(-10, 10), axes = FALSE, xlab = "", ylab = "")
 # Draw the axes at the specified positions (crossing at (0, 0)). 
 # The ticks positions are specified (those at (0,0) are omitted). 
 # Note the use of las=1 to produce horizontal tick labels on the 
 # vertical axis.
 axis(1, pos = 0, at = c(-2, -1, 1, 2)) 
 axis(2, pos = 0, at = c(-10, -5, 5, 10), las = 1)
 # Use the mathematical annotation facility to label the plot.
 title(main = expression(italic(y==x * (2 * x - 3) * (2 * x + 3))))

(from Ross Ihaka)

3次元曲面を等高線プロットの様に塗分ける、例えば大きな z に対しては黒、小さな z に対しては白

# Create a simple surface f(x,y) = x^2 - y^2
 nx <- 21
 ny <- 21
 x <- seq(-1, 1, length = nx)
 y <- seq(-1, 1, length = ny)
 z <- outer(x, y, function(x,y) x^2 - y^2)

多面体の各面の隅の値を平均し [0,1] 間の値にスケール化

# We will use this to select a gray for colouring the facet.
hgt <- 0.25 * (z[-nx,-ny] + z[-1,-ny] + z[-nx,-1] + z[-1,-1]) 
hgt <- (hgt - min(hgt))/ (max(hgt) - min(hgt))

各面の色を指定した曲面を描く

persp(x, y, z, col = gray(1 - hgt), theta = 35) 
persp(x, y, z, col = cm.colors(10)[floor(9 * hgt + 1)], theta = 35)

(from Ross Ihaka)

polygon 関数使用例

y <- c(1.84147098480790, 1.90929742682568, 1.14112000805987,
         0.243197504692072, -0.958924274663133,0.720584501801074,
        1.65698659871879, 1.98935824662338, 1.41211848524176,
        0.45597888911063)
plot(y,ylim=c(0.0000000001,2),type="l",log="y")
x <- 1:10
temp <- which(y<0) 
polygon(x[-temp], y[-temp], col="red")

Uwe Ligges

matplot による並置出力

(Peter Dalgaard, June 11, 2001)
p <- 1:3 
u <- matrix(c(1,1,1,2,2,2,3,3,3),3,) 
r <- p*u
X11(height=3.5,width=7) 
par(mfcol=c(1,3),cex=1.5,mex=0.6, mar=c(5,4,4,1)+.1) 
matplot(p,r,type="b",main="A",col="black") 
matplot(log(p),log(r),type="b",main="B",col="black") 
r <- p+u 
matplot(p,r,type="b",main="C",col="black")

6章以降は P_Johnson_tips_2 へ移動


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-03-25 (土) 11:19:16