P_Johnson_tips_1
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
COLOR(red){SIZE(25){Paul Johnson 氏の R Tips 集訳 1/2}}
後半 [[P_Johnson_tips_2]] へ移動~
2003/08/09 現在~
注：これは P. Johnson さんが r-help の記事から役にたちそ...
注意：あまりに長くて表示に時間がかかるので二つに分割しま...
----
(私はこれを "StatsRus" と名付けるのが気がきいていると思っ...
もしあなたが R の Tips 集が必要なら、これがそうです。これ...
Brian D. Ripley の言葉を注意しましょう、
"One enquiry came to me yesterday which suggested that so...
distributions do not know that R comes with two Guides in...
directory plus an FAQ and the help pages in book form. I...
distributed with all the binary distributions (they are n...
installed by default). Windows-specific versions are ava...
rw1001d1.zip and rw1001d2.zip"
FaqManager ソフトに対しその著者の Stas Bekman に謝辞を捧...
#contents
~
* 1 データの入出力 [#k0b83b2c]
**1.1 生の数値を R に取り込む [#t84e9115]
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
これは本当に簡単。ファイル "myData"で、空白区切りの数値が...
//This is truly easy. Suppose you've got numbers in a spa...
myDataFrame<-read.table("myData",header=TRUE)
このコマンドは素晴らしく動く。
//command works great.
"?read.table" とタイプすれば、他の区切り記号を使っている...
//If you type "?read.table" it tells about importing file...
タブ区切りのデータを持っていて、空白を欠損値を示すのに使...
//Suppose you have tab delimited data with blank spaces t...
> myDataFrame<-read.table("myData",sep="\t",na.strings="...
(R1.3 の新しいフィーチャー) データが gzipで圧縮されたフ...
//(new R1.3 feature) Suppose your data is in a gzipped fi...
myDataFrame <- read.table(gzfile("myData.gz"), header=T)
別々のファイルの列を読み込む際には、以下のようにして、こ...
//If you read columns in from separate files, combine int...
variable1 <- scan("file1")
variable2 <- scan("file2")
mydata <- cbind(variable1, variable2)
もしくは、以下を使えば、同様なことができる。
//or use the equivalent:
mydata <- data.frame(variable1, variable2)
R オブジェクトのデータフレームを、以下のコマンドで保存す...
//Optionally save dataframe in R object file with:
write.table(mydata, file="filename3")
**1.2 データアクセスに関する基本的注意 [#y50a7128]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
データフレーム "x" の列を列番号でアクセスするには、例えば...
//To access the columns of a data frame "x" with the colu...
注：原文にも抜けているのだが，x["pjVar1"] は x$pjVar1 と...
注の注：「第1列の名前が"pjVar1"であることを知っていれば」...
もし，そのデータフレームに対して attach(x) を行うと，単に...
//If you attach that data frame with attach(x) then you c...
回帰分析において，data=x オプションを含めることにより，あ...
//In a regression, include the option data=x and then you...
リストの要素を得るためには x[[1]] と...
//もとの和訳ではそこを考慮していなかったので，意味不明の...
//Grab an element in a list as x[[1]]
例えば，もしデータフレームが "nebdata" ならば，
//For instance, if a data frame is "nebdata" then
nebdata[1, 2]
あるいは
//or
nebdata[nebdata$Volt=="ABal", 2]
あるいは，データフレームにアタッチすることで，
//or by attaching the data frame
attach(nebdata)
nebdata[Volt=="ABal", 2]
detach(nebdata)
(from Diego Kuonen)
**1.3 新しい Data Import/Export マニュアルのチェック [#r9...
(13/08/2001 Paul Johnson <pauljohn@ukans.edu>)
R-1.2 において，R チームは R にデータを読み込んだり R か...
//With R-1.2, the R team released a manual about how to g...
**1.4 データを R と他のプログラム (Excelなど)の間で交換す...
(04/05/2001 Paul Johnson <pauljohn@ukans.edu>)
私はこの覚え書きを新しい R のデータインポート・データエク...
//I made these notes before the new R Data Import/Export ...
Excel -> R
エクセルで，例えば table.txt などという名前のテキストファ...
//Use 'save as space delimited text (.prn)' in Excel to c...
R -> Excel
テキストファイルに書き出すときには，sink("table") とする...
//Use sink("table"), where 'table' is the name of the R o...
(from Griffith Feeney)
R-Excel インターフェースは試したかい？
http://lib.stat.cmu.edu/R/CRAN/contrib/extra/index.html ...
(from Charles Raux)
//Did you try R-Excel interface ? see http://lib.stat.cmu...
私の学生達に教えようと思っていることを書いておく（インポ...
//Here's what I'm planning to teach my students (Import o...
今のところ，最も簡単なやり方は，えとね，データをテキスト...
//Currently the easiest way is to to ... export data as a...
現在の最も普遍的な三つの PC プログラムにおけるテクニック...
//最小公分母というのは，日本ではあまり使わないな
//The current techniques for the three most common PC pro...
SPSS：データエディタの File, Export メニューにより，タブ...
//SPSS: Use |File, Export| menu in the data editor and ex...
このやり方はカンマが小数点として使われているかどうか loca...
//This refers to locales where the comma is used as decim...
Excel: File, Save as... を使い，タブ区切り（.txt）粗選択...
//Excel: Use |File,Save as...| and choose tab-sep (|.txt|...
SAS: カンマ区切りかタブ区切りでエクスポートできる。SAS は...
//SAS: It is possible to export in a comma-separated or t...
(from Peter Dalgaard)
**1.5 データフレームのマージ [#y72d434c]
(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)
experiment<-data.frame(times=c(0,0,10,10,20,20,30,30),ex...
simul<-data.frame(times=c(0,10,20,30),simul=c(3,4,5,6))
以下のように、マージされたデータセットが欲しい。
//I want a merged datatset like:
> 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
答え
//Answers
merge(experiment, simul)
が，貴方の望む全てだろう。
//does all the work for you.
(from Brian D. Ripley)
あるいは，以下を検討する
//Or consider:
exp.sim<-data.frame(experiment,simul=
simul$simul[match(experiment$times,simul$...
(from Jim Lemon)
このようなデータフレームを持っている。
//I have dataframes like this:
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
以下が欲しい
//And I want
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 のレスは以下のとおり
//In response, Ben Bolker said
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を利用した結果は，以下のようにな...
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
（訳注その２）次のようにしたらうまくいくけど，どうでしょ...
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...
merge(data1,data2,by=1)
**1.6 一時に一つの行を加える [#q885c512]
(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)
Question: I would like to create an (empty) data frame wi...
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 ...
If you know the number of results in advance, say, N, do ...
df <- data.frame(time=numeric(N), temp=numeric(N), press...
df[1, ] <- c(0, 100, 80)
df[2, ] <- c(10, 110, 87)
...
または
//or
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 on...
If you don't know the final size you can use rbind to add...
[top]
**1.7 データフレームに対するもう一つのマージ法 [#v32711f4]
(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 オブジェクトがあるかどうか検査する [#p66901e7]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
R がロードする際に、R はロードするものがないとき警告は出...
//If you load things, R does not warn you when they are n...
is.null(list$component)
上のコマンドは list という名のものにある component と名づ...
//to check a thing named component in a thing named list.
存在しないデータフレームのカラムを "[" を使ってアクセスす...
//Accessing non-existent dataframe columns with "[" does ...
>data(trees)
>trees$aardvark
NULL
>trees[,"aardvark"]
Error in [.data.frame(trees, , "aardvark") : subscript o...
(from Thomas Lumley)
**1.9 乱数を発生する [#c21fab27]
(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 ...
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 ...
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 平均/分散を固定して乱数を発生 [#bbb2a147]
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
任意のツールで乱数を生成する場合、指定したちょうどの平均...
//If you generate random numbers with a given tool, you d...
そこで、平均０のサンプルをどのようにして強制的に作成する...
//So, how to force the sample to have a mean of 0? Take a...
> 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
そして平均５、標準偏差２のサンプルができる。
//and now create your sample with mean 5 and sd 2:
# 以下は同じ
> x <- x*2 + 5
> mean(x)
[1] 5
> var(x)
[1] 4
(from Torsten.Hothorn)
**1.11 重みつきの複製 [#z9bcb072]
(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...
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 freque...
If you have frequency weights you may have to expand your...
expanded.data<-as.data.frame(lapply(compressed.data,
function(x) rep(x,compressed.data$weig...
should do it.
(from Thomas Lumley)
**1.13 分割表をデータフレームに変換する [#s9c0f947]
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
もし８次元のクロス表があるなら、表にセルの８つのファクタ...
//Given a 8 dimensional crosstab, you want a data frame w...
R1.2 はこれを扱うのに as.data.frame.table() 関数を導入し...
//R1.2 introduces a function as.data.frame.table() to han...
この処理はマニュアルでもできる。ここに関数がある（expand....
//This can also be done manually. Here's a function(it's ...
dfify <- function(arr, value.name = "value", dn.names =
names(dimnames(arr))) {
Version <- "$Id: dfify.sfun,v 1.1 1995/10/09 16:06:12 d3...
dn <- dimnames(arr <- as.array(arr))
if(is.null(dn)) stop("Can't data-frame-ify an array with...
names(dn) <- dn.names
ans <- cbind(expand.grid(dn), as.vector(arr))
names(ans)[ncol(ans)] <- value.name ans }
関数名は "data-frame-i-fy" を縮めたものである。
//The name is short for "data-frame-i-fy".
質問の例では、多重配列が適切な dimnames を持っていると推...
//For your example, assuming your multi-way array has pro...
my.data.frame <- dfify(my.array, value.name="frequency")
(from Todd Taylor)
**1.14 書き込み: アスキーファイルのデータ [#o0598ccd]
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
たとえば 28 x 28 データマトリックスを生成するコマンドを持...
//Say I have a command that produced a 28 x 28 data matri...
write.table(mat,file="filename")
Pesl Thomas は「"tab"という名前のデータフレームの内容をフ...
//Pesl Thomas said, "I want to output a file to my disc w...
> tab
or 1 2 3 4 5 6 1 57 80 25 23 46 23 2 106 35 59 8 51 4
これが使える
//You can use
write.table(unclass(tab))
これも使える
//You could also do
write(t(tab),ncol=NCOL(tab))
なぜなら、 write() は行と列を転置する。(from Thomas Lumley)
//since write() transposes rows and columns. (from Thomas...
注目してもらいたいのは、MASS ライブラリには、write.matrix...
//Note MASS library has a function write.matrix which is ...
*2 データフレームを使った作業: 記録、合体 [#d7e9a175]
**2.1 変数をデータフレーム(又はリスト)に加える [#l5ca6160]
(02/06/2003 Paul Johnson <pauljohn@ku.edu>)
In r-help, I asked " I keep finding myself in a situation...
> 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 ...
eval(substitute(myList$newColumn<-aColumn),
list(newColumn=as.name(varName)))
as an alternative to parse().
-(Thomas Lumley)
**2.2 変数名を「その場で」作る [#oe2e9c8b]
(10/04/2001 Paul Johnson <pauljohn@ukans.edu>)
変数名とそれにセットする値を一緒に貼り付ける。assign を使...
//Paste together a variable name, set it to a value. Use ...
> assign(paste("file", 1, "max", sep=""), 1)
> ls()
[1] "file1max"
(Brian Ripley, June 18, 2001)
**2.3 一つの列を記録する、値を別の列に出力する [#ofcc96a7]
??Recode one column, output values into another column
(20/06/2001 Paul Johnson <pauljohn@ukans.edu>)
Aleksey Naumov がこの質問をした。例えば、ここに簡単なデー...
//Aleksey Naumov asked this question. For example, here i...
var1
1 1
2 2
3 3
どのようにして、上の "var1" をグループ化した値 'var2' を...
//How do I add a new column to data - 'var2' in which I g...
var1 var2
1 1 4
2 2 4
3 3 5
data$var2 <- c(4,4,5)[data$var1]
または
//or
data <- transform(data, var2=c(4,4,5)[var1])
(from Peter Dalgaard)
I admit the previous seems complicated, and I'm inclined ...
x <- rnorm(10000)
x[x > 1.5] <- NA
And if you don't want to replace the original variable, c...
You can put other variables inside the brackets, so if yo...
x[y == 1] <- 234
So, to do the v1 v2 example above, I think it is clearest...
> 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 bel...
For example, suppose I get a bunch of variables coded on ...
1= no
6 = yes
8 = tied
9 = missing
10 = not applicable.
Recode that into a new variable name with 0=no, 1=yes, an...
It seems like the replace() function would do it for sing...
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...
test <- replace(test, test == 1, 0)
test <- replace(test, test == 2, 1)
test
I suppose a better function would take from and to lists ...
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 wa...
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...
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...
**2.4 指示(ダミー)変数をつくり出す [#lb43f1a5]
//Create indicator (dummy) variables
(20/06/2001 Paul Johnson <pauljohn@ukans.edu>)
以下に2例ある。
c is a column, you want dummy variable, one for each vali...
> 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 )
質問；５つのカテゴリーのアル変数があり、各カテゴリに対す...
//Question: I have a variable with 5 categories and I wan...
答え；行の添え字または model.matrix を使う。
//Answer: Use row indexing or model.matrix.
ff <- factor(sample(letters[1:5], 25, replace=TRUE))
diag(nlevels(ff))[ff,]
または
//or
model.matrix(~ ff - 1)
(from Brian D. Ripley)
**2.5 時系列回帰のための変数の遅延値をつくり出す [#vf0da2...
//Create lagged values of variables for time series regre...
(02/06/2003 Paul Johnson <pauljohn@ku.edu>)
Peter Dalgaard explained, "the simple way is to create a ...
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 関連する観測値を持たない因子水準をデータセットから...
//How to drop factor levels for datasets that don't have ...
(08/01/2002 Paul Johnson <pauljohn@ukans.edu>)
もっともいいのが、水準を落とすことである。たとえば、以下...
//The best way to drop levels, BTW, is
problem.factor <- problem.factor[,drop=TRUE] (from Brian...
**2.7 データフレームから観測値を選択する、部分集合を取り...
//Select/subset observations out of a dataframe
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
ある変数が欠損値である行を省きたいなら、John Fox の言うよ...
//Want to drop out observations that are mising on a vari...
x2 <- na.omit(x)
これにより、欠損値を含む行すべてが除かれたデータフレーム ...
//produces a copy of the data frame x with all rows that ...
欠損値の除外については、Peter Dalgaard は以下のコマンドを...
//For exclusion of missing, Peter Dalgaard likes
subset(x,complete.cases(x)) or x[complete.cases(x),]
//"is.na(x)”を追加することよりも、 x != "NA" の方が望ま...
x != "NA" を追加するよりも、 "is.na(x)” の方が望ましい。
//adding "is.na(x) is preferable to x != "NA"
Suppose you want observations with c=1 in df1. This makes...
df2 <- df1[df1$c == 1,]
and note that indexing is pretty central to using S (the ...
Or use "match" select values from the column "d" by takin...
> 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 lik...
> 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 ov...
maleOver40 <- subset(d, sex %in% c("m","M") & age > 40)
maleOver40 <- d[(d$sex == "m" | d$sex == "M") & d$age > ...
**2.8 最初の観測値をそれぞれから抹消する [#k5d732bf]
//Delete first observation for each
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
Given data like:
1 ABK 19910711 11.1867461 0.0000000 ...
2 ABK 19910712 11.5298979 11.1867461 111
6 CSCO 19910102 0.1553819 0.0000000 ...
7 CSCO 19910103 0.1527778 0.1458333 ...
remove the first observation for each value of the "sym" ...
newhilodata <- hilodata[-c(1,6,13),]
To solve the more general problem of omitting the first i...
newhilodata <- subset(hilodata, diff(c(0,as.integer(sym)...
(actually, the as.integer is unnecessary because the c() ...
(from Peter Dalgaard)
Alternatively, you could use the match function because i...
> 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 ...
subset(hilodata, duplicated(sym))
has got to be the simplest variant.
(from Peter Dalgaard)
**2.9 データのランダム標本を選ぶ [#r4bf29d7]
//Select a random sample of data
(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 関数を忘れるな [#j1cad700]
//Selecting Variables for Models: Don't forget the subset...
(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)
直接行その他を削除することでデータを管理できるが、subset(...
//You can manage data directly by deleting lines or so fo...
Peter Dalgaard が「組み込み」データセット airquality を使...
//Peter Dalgaard gave this example, using the "builtin" d...
data(airquality)
names(airquality)
lm(Ozone~.,data=subset(airquality,select=Ozone:Month))
lm(Ozone~.,data=subset(airquality,select=c(Ozone:Wind,Mo...
lm(Ozone~.-Temp,data=subset(airquality,select=Ozone:Mont...
lm の RHS の "." は「すべての変数」を意味しており、rhs の...
//The "." on the RHS of lm means "all variables" and the ...
**2.11 すべての数値変数を処理し、文字変数を無視するには？...
??Process all numeric variables, ignore character variabl...
(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),...
mean<-sapply(data[,indices],mean) sd<-sapply(data[,indic...
**2.12 二つ以上の変数に関しソート [#x55a44eb]
//Sorting by more than one variable
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
Can someone tell me how can I sort a list that contains d...
name M
1234 8
1234 8.3
4321 9
4321 8.1
I also would like to set a cut-off, so that anything belo...
I take it that the cutoff is on the value of M. OK, suppo...
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 im...
If your two variables are in a data frame you can combine...
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 ...
dat <- dat[order(dat$name, pmax(coff, dat$M)), ]
(from Bill Venables)
I want the ability to sort a data frame lexicographically...
Here's how:
spsheet[order(name,age,zip),]
(from Peter Dalgaard)
**2.13 ある因子により定義されるサブグループ内で順位を取る...
//Rank within subgroups defined by a factor
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
by() のヘルプを読みなさい
//Read the help for 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 ...
//Use na.omit to screen missings from a data frame:
(22/08/2000 Paul Johnson <pauljohn@ukans.edu>)
欠損値が "insure" という名前のデータセットより除去される...
//To make sure missings are omitted from a dataset called...
thisdf <- na.omit(insure[,c(1, 19:39)])
body.m <- lm(BI.PPrem ~ ., data=thisdf, na.action=na.fail)
**2.15 各行から一つの集約値を計算 [#d9eb624c]
//Aggregate values, one for each line
(16/08/2000 Paul Johnson <pauljohn@ukans.edu>)
質問:200行で，各行に32個の実数データを含むテキストファイ...
//Question: I want to read values from a text file - 200 ...
Peter Dalgaard は以下を提案した。
//Peter Dalgaard says do thi:
y <- apply(as.matrix(read.table("myfile")), 1, mean)
x <- seq(along=y)
（ファイル形式によっては，read.table に，いくつかのオプシ...
//(possibly adding a couple of options to read.table, dep...
**2.16 各因子に対する集約値を保持する新しいデータフレーム...
//Create new data frame to hold aggregate values for each...
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
これでどう
//How about this:
Z<-aggregate(X, f, sum)
(X の変数はすべて合計できるものとする)
//(assumes all X variables can be summed)
または、（もし）Xが因子も含んでいる場合、要約統計の計算に...
//Or: [If] X contains also factors. I have to select vari...
Z <- data.frame(f=levels(f),x1=as.vector(tapply(x1,f,sum...
(from Wolfgang Koller)
**2.17 データフレームを選択的に集約 [#d43804cb]
//Selectively aggregate a data frame
(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),...
//note the first two columns are excluded by [,(1:2)] and...
最初の2列は [,-(1:2)] で除外され，by オプションで出力のた...
**2.18 実数から各桁数値を一時に一つずつ取り出す [#yc94b1e9]
//Rip digits out of real numbers one at a time
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
//I want to "take out" the first decimal place of each ou...
小数点以下1桁目の数値を取り出し出現度数をプロットする。次...
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 リスト中の各副行列から一つの項目を取り出す [#t2307...
//Grab an item from each of several matrices in a List
(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)
Let Z denote the list of matrices. All matrices have the ...
lapply(Z, function(x) x[1,2])
should do this, giving a list. Use sapply if you want a v...
**2.20 データセット中の値を示すベクトルを得る [#y017ddba]
//Get vector showing values in a dataset
(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の命令を表す文字列の値を計算する [#t38d404f]
//Calculate the value of a string representing an R comma...
(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 m...
Or
eval(parse(text="ls()"))
Or
eval(parse(text = "x[3] <- 5"))
**2.22 "which" はあるテストをパスする観測値の添字値を抜き...
//"Which" can grab the index values of cases satisfying a...
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
大きなベクトルのデータを箱ひげ図で分析して、外れ値を見つ...
//To analyse large vectors of data using boxplot to find ...
which(x==boxplot(x,range=1)$out)
**2.23 行列/データフレーム中の一意的な行を見出す [#k57449...
//Find unique lines in a matrix/data frame
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
Jason Liao writes:
"I have 10000 integer triplets stored in A[1:10000, 1:3]....
Peter Dalgaard answers, "As of 1.4.0 (!):
unique(as.data.frame(A))"
**2.24 数値変数だけを選ぶ [#h40cdbff]
//Select only numerical variables
(12/11/2002 Paul Johnson <pauljohn@ku.edu>)
sapply(dataframe, is.numeric)
または
//or
which(sapply(data.frame, is.numeric))
Thomas Lumley 11/11/2002
* 3 行列とベクトルの演算 [#yb7940fe]
**3.1 単位行列をつくり出すには？ [#pa2417c8]
//How to create an identity matrix?
(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" を一つの長いベクトルに変換 [#z97db96f]
//Convert matrix "m" to one long vector
(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) をつくり出す [#e...
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
なぜこれが有用なのかは分かりませんが、いくつかの列と行列...
//I don't know why this is useful, but it shows some row ...
結局は Brian Ripley の方法を使うことになった。この方法は...
//I ended up using Brian Ripley's method, as I got it fir...
A <- matrix(1, n-1, n-1)
rA <- row(A)
rA[rA + col(A) <= n]
しかしながら、 その後、Andy Royle の御蔭で、もっと簡単で...
//However, thanks to Andy Royle I have since discovered t...
// sequence(n-1:1) もとからこうなっていたのだろうが，良...
sequence((n-1):1)
[1] 1 2 3 4 1 2 3 1 2 1
Karen Kotschy
**3.4 すべての n 個めの項目を選ぶ [#o8322b1e]
//Select every n'th item
(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)
非常に長いベクトルよりｎ個目の要素をすべて取り出す、これ...
//extract every nth element from a very long vector, vec?
vec[seq(n, length(vec), n)] #(from Brian Ripley)
seq(1,11628,length=1000) はインデックスとして使える 1:116...
//seq(1,11628,length=1000) will give 1000 evenly spaced n...
My example:
vec <- rnorm(1999)
newvec <- vec[1, length(vec), 200]
> newvec
[1] 0.2685562 1.8336569 0.1371113 0.2204333 -1.2798...
[7] -0.2366385 0.5060078 0.9680530 1.2725744
This shows the items in vec at indexes (1, 201, 401, ...,...
**3.5 ベクトル中の 1.5 に最も近い値の添字を見出す [#b9a08...
//Find index of a value nearest to 1.5 in a vector
(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 ベクトル中の非負値項目の添字を見つける [#ha85c609]
//Find index of nonnegative items in vector
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
which(x!=0) (from Uwe Ligges)
rfind <- function(x) seq(along=x)[x != 0] #(from Brian ...
スピードとメモリーの効率性に関しては、以下がある。
//Concerning speed and memory efficiency I find
as.logical(x)
これは以下よりもよく
//is better than
x!=0
そして
//and
seq(along=x)[as.logical(x)]
これは以下よりもよく
//is better than
which(as.logical(x))
このため
//thus
which(x!=0)
は最も短いが
//is shortest and
rfind <- function(x)seq(along=x)[as.logical(x)]
以下はコンピュータではもっとも効率的であるようである。
//seems to be computationally most efficient
(from Jens Oehlschl?gel-Akiyoshi)
**3.7 欠損値の添字を見出す [#j42f9f57]
//Find index of missing values
(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)
//Suppose the vector "Pes" has 600 observations. Don't do...
Pes というベクトルが600の観測値を持つとする。欠損値の添え...
(1:600)[is.na(Pes)]
//The `approved' method is
正しいやり方は，次の通り。
seq(along=Pes)[is.na(Pes)]
//In this case it does not matter as the subscript is of ...
この場合にはどちらでも良いことではあるが，添え字の長さが0...
(from Brian Ripley)
//However, the solution I gave
しかし，私は以下の解を与える
which(is.na(Pes))
is the one I still really recommend; it does deal with 0-...
**3.8 ベクトル中の最大値を持つ項の添字を見出す [#hdd3d6af]
//Find index of largest item in vector
(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 行列中の値を変更する [#ad45fc77]
//Replace values in a matrix
(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 がデータフレームなら、データを強制的にマトリックス形式...
//If la is a data frame, you have to coerce the data into...
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
以下を試してみよう
//Try this:
la <- ifelse(la == 0, 1, la)
(from Martyn Plummer)
日本語訳を見てのコメント~
la がデータフレームのとき，何もせずいきなり la[la==0] <- ...
**3.10 行列から特定の行を取り除く [#rdc52cd6]
//Delete particular rows from matrix
(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 ある基準を満たす項目の総数を数える [#a8a81810]
//Count number of items meeting a criterion
(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)
以前の質問で述べた which() の結果に"length()" を以下のよ...
//Apply "length()" to results of which() described in pre...
length(which(v<7))
または
//or
sum( v<7,na.rm=TRUE)
**3.12 相関行列から偏相関係数を計算する [#pcb94fcf]
//Compute partial correlation coefficients from correlati...
(08/12/2000 Paul Johnson <pauljohn@ukans.edu>)
多変量変数間の相関係数から，二変数以外の変数を制御して二...
//I need to compute partial correlation coefficients betw...
実際問題として，これはきわめて簡明直裁。R が相関係数行列...
//Actually, this is quite straightforward. Suppose that R...
Rinv<-solve(R)
D<-diag(1/sqrt(diag(Rinv)))
P<- -D %*% Rinv %*% D
P の非対角要素は他の変数の影響を除いた（コントロールした...
//The off-diagonal elements of P are then the partial cor...
(from John Fox )
一般的に，分散・共分散行列の逆行列を求めて，対角成分が1に...
//In general you invert the variance-covariance matrix an...
pcor2 <- function(x){
conc <- solve(var(x))
resid.sd <- 1/sqrt(diag(conc))
pcc <- - sweep(sweep(conc, 1, resid.sd, "*"), 2, resid...
return(pcc)
}
pcor2(cbind(x1,x2,x3))
see J. Whittaker's book "Graphical models in applied mult...
以下は，私が今使っているもので，それぞれの偏相関係数の有...
//This is the version I'm using now, together with a test...
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...
"*")
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 -...
df = t.df))
}
return(pcc)
}
(from Kaspar Pflugshaupt)
**3.13 正準変数を計算 [#y693ff4f]
//Compute Canonical variates
(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)
同じ観察対象について，データ行列 A と B があるとする。
//Suppose you have data matrices A and B on the same obse...
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=Bp...
}
この関数には，トリックがある。正準変量の平均は0で，分散は...
//should do the trick. The canonical variates are zero-me...
**3.14 多次元行列(R 配列)をつくり出す [#nf518cfb]
//Create a multidimensional matrix (R array)
(20/06/2001 ? <?>)
//Brian Ripley said:
Brian Ripley の話では
my.array<-array(0,dim=c(10,5,6,8))
//will give you a 4-dimensional 10 x 5 x 6 x 8 array.
上のようにすれば、 10 x 5 x 6 x 8 の 4次元の配列が与えら...
//Or
または以下のような方法もある。
array.test <- array(1:64,c(4,4,4))
array.test[1,1,1]
1
array.test[4,4,4]
64
**3.15 多数の行列を結合する [#vf71c9d3]
//Combine a lot of matrices
(20/06/2001 ? <?>)
100個の n×2 行列のリストがあるとする。どうやったら，これ...
//Given a list with 100s of (nx2) matrices, how do you co...
list は行列のリストである。
//list is the list of matrices.
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 特定の論理に基づき "近傍" 行列をつくり出す [#kd19f...
//Create "neighbor" matrices according to specific logics...
(20/06/2001 ? <?>)
あるセルがある位置に近傍があるかどうかの有無を０と１で示...
//Want a matrix of 0s and 1s indicating whether a cell ha...
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(...
bishop.case <- function(y,z,mat){coldiff(y,z,mat)==1 & r...
queen.case <- function(y,z,mat){rook.case(y,z,mat) | bis...
matrix(as.numeric(
sapply(x,function(y)sapply(x,rook.case,y,mat...
matrix(as.numeric(
sapply(x,function(y)sapply(x,bishop.case,y,m...
matrix(as.numeric(
sapply(x,function(y)sapply(x,queen.case,y,ma...
(from Ben Bolker)
**3.17 ある "key" により二つの数値列をマッチング [#s3f21f...
//"Matching" two columns of numbers by a "key"
(20/06/2001 ? <?>)
The question was:
I have a matrix of predictions from an proportional odds ...
小さな例として，以下の x, y
//As a toy example, if I have
> 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] (t...
Is
x[cbind(seq(along=y), y)]
what you had in mind? The key is definitely matrix indexi...
**3.18 上・下三角行列をつくり出す [#idb327b7]
//Create Upper or Lower Triangular matrix
(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...
> 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 の逆行列を計算 [#j6878b82]
//Calculate inverse of X(20/06/2001 ? <?>)
solve(A)は A の逆行列を返す
**3.20 行列添字の面白い使い方 [#r4d2915e]
//Interesting use of Matrix Indices
(20/06/2001 ? <?>)
Here's a similar problem.
Mehdi Ghafariyan said "I have two vectors A=c(5,2,2,3,3,2...
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...
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 固有値の例 [#p3eb8b69]
//Eigenvalues example
(20/06/2001 ? <?>)
Tapio Nummi asked about double-checking results of spectr...
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 の適用等 [#wedb7571]
**4.1 一つの関数から複数の値を返す [#dfe54df3]
//Return multiple values from a function
(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)
Below is a summary on returning variable from subfunction...
To access the variables within a function, return a list ...
"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" 値を掴み出す [#v0bd0df1]
//values out of a list of significance tests
(22/08/2000 Paul Johnson <pauljohn@ukans.edu>)
//Suppose chisq1M11 is a list of htest objects, and you w...
chisq1M11 が htest オブジェクトのリストだとする。そして，...
> apply(cbind(1:1000), 1, function(i) chisq1M11[[i]]$p.v...
//And Peter Dalgaard showed a more elegant approach:
Peter Dalgaard はもっとエレガントなやり方を示した。
sapply(chisq1M11,function(x)x$p.value)
このどちらも，"function" と呼ばれる関数を作り，それぞれの...
（訳注：本当は無名の関数を作っているのだが）
//In each of these, a simple R function called "function"...
**4.3 ifelse の用法 [#zca930f7]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
//If you want to calculate something, but only if y is *b...
y が 0 から 1 の範囲にあるときにのみある計算をしたいとき...
ifelse(y==0,0,y*log(y))
訳注：上のは，y が 0 であるときとそれ以外で別の計算をする...
(from Ben Bolker)
**4.4 apply を用いて各セル毎の確率の行列を作る [#je988350]
//Apply to create matrix of probabilities, one for each "...
(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)
ガンマ関数のスケールパラメータと形状パラメータの様々な値...
//Suppose you want to calculate the gamma density for var...
gridS <- expand.grid(scales, shapes)
survLilel <- function(ss) sum(dgamma(Survival,ss[1],ss[2...
Likel <- apply(gridS,1,survLilel)
実際には私は以下のように使う
//Actually, I would use
sc <- 8:12; sh <- 7:12
args <- expand.grid(scale=sc, shape=sh)
matrix(apply(args, 1, function(x) sum(dgamma(Survival, s...
log=T))), length(sc), dimnames=list(scale=sc, shape=s...
(from Brian Ripley)
**4.5 外積 outer [#j2426a97]
(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)
outer(x,y,f) does just one call to f with arguments creat...
**4.6 あるものが公式か関数か検査する [#vd03b825]
//Check if something is a formula/function
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
公式（モデル式）は "formula" クラスを持つので，inherits(o...
//Formulae have class "formula", so inherits(obj, "formul...
もし，Y~X+Z+W ではなくて ~X+Z+W であることが分かっている...
//And if you want to ensure that it is ~X+Z+W rather than...
inherits(obj,"formula") && (length(obj)==2)
(from Thomas Lumley)
**4.7 変数のベクトルに関する最適化 [#h99431a3]
//Optimize with a vector of variables
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
The function being optimized has to be a function of a si...
opt.func <- function(arg)
t(Y-(X[,1] * arg[1] + X[,2] * arg[2])^delta) %*% covar...
(Y-(X[,1] * arg[1] + X[,2] * arg[2])^delta)
(from Douglas Bates)
**4.8 S+ のような slice.index [#t3cf7951]
(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)
slice.index<-function(x,i){k<-dim(x)[i];sweep(x,i,1:k,fu...
(from Peter Dalgaard)
*5 グラフ [#l8916f99]
**5.1 グラフ描きの前に par で特性を変更 [#m4ae1348]
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
par() はグラフの品質に影響を与える多目的コマンド。例を以...
//par() is a multipurpose command to affect qualities of ...
***1. 以降のプロットを一つの絵にまとめる: [#s88a8fc8]
par(mfrow=c(3,2))
プロットの 3x2 のマトリックスを作成する。
//creates a 3x2 matrix of plots.
***2. 余白を次の例のように変更 [#ka7a90a6]
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 グラフ出力を保存 [#c7347e49]
(15/08/2001 Paul Johnson <pauljohn@ukans.edu>)
グラフ保存には２つの方法がある。グラフの保存前に、出力形...
//There are two ways to save and/or print graphs. Before ...
There are other devices representing kinds of file output...
>?device
Rがシステムで見つけるデバイスのリストが見れる。
//to see a list of devices R finds on your system.
各デバイスについて、より多くを見つけられる。
//For each device, you can find out more about it.
?x11
または
//or
?png
また MS Windowsでは
//or, in MS Windows
?windows
(等々)
//(and so forth)
Before you start using a device, you can configure it. Th...
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 wh...
> dev.set(3)
Now here is the big thing. The devices are configured dif...
Suppose we want png output. We create an instance of a pn...
> png(filename="mypicture.png", width=480, height=640, p...
Note, if no filename is here, it will print to the defaul...
> dev.off()
Some print devices can accept multiple plots, and dev.off...
The postscript device has options to output multiple grap...
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 us...
So what is the alternative?
If you already have a graph on the screen, and did not pr...
dev.copy() opens the device you specify. It is as if (but...
dev.copy(device=png, file="foo", width=500, height=300)
dev.off()
The dev.copy() function takes a device type as the first ...
If you forget the last line, it won't work. It's like a w...
If you use dev.copy or dev.print, you may run into the pr...
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 ...
x11(height=6, width=6)
There are some specialized functions that can do this in ...
The dev.print() function, as far as I can see, is basical...
dev.print(device=postscript,file="yourFileName.ps")
Before you run the graphs, you can do use the postscript(...
If you want pdf output, as of R 1.3 there is a pdf device...
If you just need to test your computer setup, Bill Simpso...
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, poin...
The format of png use is the same.
As of R1.1, the dev.print() and dev.copy2eps() will work ...
> ps <- function(file="Rplot.eps", width=7, height=7, .....
dev.copy2eps(file=file, width=width, height=heigh...
}
> 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 user...
The mismatch of "size" between devices even comes up when...
> dev.print(height=6, width=6, horizontal=FALSE)
You might want to include pointsize=20 or whatever so the...
One user observed, "Unfortunately this will also make the...
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-pla...
The rescaling of a windows() device works differently, an...
**5.3 プロット出力を別個のファイルに保存するには [#z501e9...
(10/04/2001 Paul Johnson <pauljohn@ukans.edu>)
postscript(file="test3.%d.ps",onefile=FALSE)
ファイル名を test3.1.ps, test3.2.ps などと付ける。(from T...
//gives files called test3.1.ps, test3.2.ps and so on. (f...
**5.4 用紙サイズをコントロールする [#o2fe077f]
(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 [#p68d2627]
(20/06/2001 Paul Johnson <pauljohn@ukans.edu>)
Advice seems to be to output R graphs in a scalable vecto...
Keith Reckdahl's Using EPS in LaTeX documents is freely a...
For a description of EPSF and various other vector and bi...
I personally swear by the LaTeX Companion and the LaTeX G...
Thanks to Prof. Ripley for the hint, maybe someone other ...
-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...
?psfrag{wrho}{$W(?varrho)$} <-- this replaces...
?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 c...
**5.6 "Snapshot" グラフとそれらの閲覧 [#ma9b3933]
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
Ordinarily, the graphics device throws away the old graph...
Have a look to ?recordPlot to see how to keep snapshots. ...
On Windows you can choose "history"-"recording" in the me...
**5.7 scatterplot の行列の上三角部分を選択する [#s7b37d6b]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
必要なのはプロットする点を指定するだけ:
//You only need to specify the points to be plotted:
## 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 i...
x <- col(mmtop94.2)[temp] # col indices
scatterplot3d(x, y, z, type="h")
(from Uwe Ligges)
**5.8 密度関数をプロットする (例えば正規分布) [#bff44199]
(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 エラーバー付きでプロット [#ufb394cd]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
以下にやり方を示す。
//Here is how to do it.
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)
//or
または
segments(x,ucl,x,lcl)
(from Bill Simpson)
**5.10 推定密度関数付きのヒストグラム [#n8f63224]
(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)
同じスケールで密度推定とヒストグラムが欲しければ、以下の...
//If you want a density estimate and a histogram on the s...
> IQR <- diff(summary(data)[c(5,2)])
> dest <- density(data, width = 2*IQR) # or some smalle...
> hist(data, xlim = range(dest$x), xlab = "x", ylab = "d...
probability = TRUE) # <<<--- this is the vit...
> 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 いくつもの線プロットを重ね描きするには？ [#rc63ca52]
(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)
//Use par("new"=TRUE) to stop the previous output from be...
前の出力がかき消されないようにするためには，以下のように ...
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="")
//Note the par() is used to overlay "high" level plotting...
par() は，高水準のプロット命令を重ねるものであるが，ablin...
Here's another idea: form a new data.frame and pass it th...
plot(cbind(ts1.ts,ts2.ts),xlab="time", ylab="", plot.typ...
//or better something like:
以下のようにすれば，もう少しまし。
plot(cbind(ts1.ts,ts2.ts),plot.type="single",
col=c("yellow","red"),lty=c("solid","dotted")) #co...
//It can also be helpful to contrast 'c(ts1.ts,ts2.ts)' w...
'c(ts1.ts,ts2.ts)' と 'cbind(ts1.ts,ts2.ts)' を比較すると...
(from Guido Masarotto)
//For time series data, there is a special function for t...
時系列データについては特別な関数がある。
ts.plot(ts1.ts, ts2.ts) # same as above
//See help on plot.ts and ts.plot (from Brian D. Ripley)
plot.ts と ts.plot のヘルプを見よ。 (from Brian D. Ripley)
//Here's an example for histograms:
以下は，ヒストグラムの例である。
> hist(rnorm(100, mean = 20, sd =12), xlim=range(0,100),...
> par(new = TRUE)
> hist(rnorm(100, mean = 88, sd = 2), xlim=range(0,100),...
//The label at the bottom is all messed up. You can put t...
横軸ラベルがめちゃくちゃになる。これを避けるには，xlab=""...
（訳注：タイトルもめちゃくちゃになるので，main="" も加え...
//Here's another example that works great for scatterplots
散布図を書くときにうまい方法がある。
（蛇足の訳注：xlim, ylim で xlim=c(min(x1,x2),max(x1,x2))...
> 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 グラフを行列状に配置する [#f3268d7d]
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
//The standard way to array several images across a page ...
ページに複数のグラフを配置する標準的な方法は，以下のよう...
par(mfrow=c(2,3))
//And then the next 6 plot/image commands will be laid ou...
そうすれば，次の 6 つのグラフは2行3列に配置される
The layout option can give a powerful set of controls for...
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(...
x2=as.factor(round(runif(100,1,4))) ) attach(x...
by(y,list(x1,x2),function(x)print(x)) xn <- as.integer(x...
## 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(x...
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 r...
And to label the whole page, use the "mtext" function
**5.13 線と棒グラフを結合する [#d1da03b8]
(07/12/2000 Paul Johnson <pauljohn@ukans.edu>)
David James was kind enough to help me out and enlighten ...
xbarplot <- function(y, col=par("col"), border = par("fg...
{
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 回帰散布図: グラフに当てはめ回帰直線を加える [#v5c...
(17/08/2001 Paul Johnson <pauljohn@ukans.edu>)
suppose you have
reg <-lm( y ~ x)
plot(x,y)
then do
abline(reg)
**5.15 散布図中のプロット文字を制御する [#q6e48599]
(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...
If you specify pch, only the first character is taken as ...
Matt Wiener suggested creating a color vector and then us...
> 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 は，プロット記号を選択するために整数を用いることを覚...
//Remember pch uses integers to select plotting characters.
x <- 1:10
y <- 1:10
res <- -4:5
plot(x, y, pch = ifelse(res < 0, 20, 1))
もし論理式が真であれば，記号20が，さもなくば記号1が使われ...
//If true, use character 20. Else use character 1. Experi...
Then note you can create a vector of integers to control ...
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 ...
**5.17 サイズと色を修正した散布図 [#d297ed49]
(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="Siz...
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 ...
(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 散布図：三番目の変数に関してサイズを調整する [#h03...
(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 散布図: 点を結ぶ線を滑らかにする [#ua9136e8]
(02/06/2003 Paul Johnson <pauljohn@ku.edu>)
//Rado Bonk asked how to smooth a line connecting some po...
Rado Bonk がいくつかの点をなめらかに結ぶ曲線を求めている...
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 回帰散布図: 推定値をプロットに加える [#a01ad6ac]
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
I don't know how to get the estimate of the line printed,...
I used this to indicate the R-squared value on a regressi...
text(5, 0.6, as.expression(substitute(R^2 == r, list(r=r...
where R2.the was computed a few lines above. (from Emmanu...
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>)
チックマークは表示するが、チックマークの数値ラベルは表示...
//Ticks, but no numbers:
plot(1:10, xaxt="n")
axis(1, labels=FALSE)
X軸およびY軸のタイトルを省くが、各チックの数値ラベルは残...
//This leaves off the x and y titles. You still have the ...
plot(x,y,ann=FALSE)
もっと多くのグラフの制御をしたいなら、以下のようなことが...
//If you want a lot of control you can do
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)
//do ?axis
? axis と入力してみよ (from Bill Simpson)
Control range of axes through the graph command itself. T...
> 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 軸: ラベルの回転 [#d9e30912]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
一度で、par(las=) は見かけ上これをする。入力は以下の通り。
//At one time, apparently par(las=) would do this, as in
par(las=2) or par(las=3)
上では軸のラベルを90度回転する。
//to make the 90 degrees to axes labels.
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 ...
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?nb...
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 軸: 軸に整形した日付を加える [#mf83c8fb]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
CRANで入手できる"date" パッケージを試してみる。
//Try package "date", available at CRAN.
作図例は以下のとおり。
//Example for plotting:
library(date)
plot.date(.....)
## if x-values are dates, you can also use plot(...) after
library(date)
(from Uwe Ligges)
**5.24 軸: プロット中の軸を逆転する [#w0080bf7]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
R1.1.1 にはバグがあるが、後のバージョンでは、Ross Ihaka ...
//In R1.1.1, there is a bug, but in later versions Ross I...
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 so...
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 ... (althoug...
Another person said try (if you don't want axes):
plot(x, -y, axes=FALSE)
axis(1)
axis(2)
**5.25 軸: 日付によるラベル付き軸 [#e3bab1b9]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
以下のようにデータセットをプロットする必要がある。
//I have to plot the data set like the following.
> x <- c("05/27/00","06/03/00","08/22/00",...)
> y <- c(10,20,15,...)
"date"パッケージを試すべきだ。このパッケージは CRAN で入...
//You should try the package "date", available at CRAN. E...
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 を利用したものだ。
//Here is one using 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日、す...
//which will work if you were in the USA and your first s...
**5.26 軸: 軸ラベルの上付き添字 [#c283699d]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
plot(c(1:10), ylab= expression(x^"++")) (from Peter Dalg...
**5.27 軸: 位置の補正 [#b93bfb1b]
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
質問：プロットするとx軸はグラフの下に，y軸はグラフの左に...
//Question: When I plot them default x axis takes place o...
以下のようにすればよいでしょう
//Maybe this helps:
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 散布図に "エラー矢印" を加える [#x4c9351f]
(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)
arrows(..., code=3, angle=90, ...) を使えば、まったく簡単...
//I think using arrows(..., code=3, angle=90, ...) is qui...
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, le...
arrows(x-se.x, y, x+se.x, y, code=3, angle=90, le...
最初の arrows() は y のエラーバーで、2番目のは x のエラ...
//The first arrows() draws the error bars for y, and the ...
Emmanuel Paradis
**5.29 時系列: 一つのグラフに複数の線をプロットする [#d80...
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
同じグラフで２つを描画するのは以下のようにする。
//To do the two in the same graph:
plot(x)
lines(y)
または
//or
points(y)
または
//or
matplot(cbind(x,y),type="l")
グラフを上下に分割して描画するには以下のようにすればよい。
//To do separate graphs one above the other:
par(mfrow=c(2,1))
plot(x)
plot(y)
色の設定、範囲、軸のラベルなどもできる。
//You can do other things like set colors, ranges, axis l...
(from Jon Baron 12/29/2001)
**5.30 時系列: 当てはめ値と実際のデータをプロット [#hc154...
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
> library(ts)
> data(LakeHuron)
> md <- arima0(LakeHuron, order=c(2,0,0), xreg=1:98) pl...
> lines(LakeHuron-md$resid,col="red")
(from Adrian Trapletti)
**5.31 プロットにテキストを挿入 [#x4c67d7b]
(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)
次の文字列をプロットに入れたい。"The average temperature ...
//I want this in a plot: "The average temperature is 23.1...
text(x,y,paste("The average temperature is ", variable, ...
(from Thomas Lumley)
**5.32 はみ出る変数をプロット [#r765e81f]
(07/12/2000 Paul Johnson <pauljohn@ukans.edu>)
Graphs blow up if a variable is unbounded, so if you are ...
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 d...
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)...
dotplot(1:2, labels = my.names)
**5.34 プロットのタイトルに数式等の複雑なものを使う [#l2f...
(11/11/2002 Paul Johnson <pauljohn@ku.edu>)
substitute() を使う。なにかこのようになる。
//Use substitute(). Something like
e<-substitute(expression(paste("Power plot of ", x ^ alp...
alpha == ch.a)), list(ch.a=formatC(alpha,wid=1)))
plot(x, x^alpha, main=e)
(from Peter Dalgaard)
以下のようにしているとすると
//Suppose you have
> that <- 1
> plot(1:10)
そして、タイトルに "that" を指定する必要がある。
//and you need to specify "that" in the title.
expression(paste(hat(theta),'= ',that)))
Generally, to get an expression either use parse on your ...
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 45...
lat<-45
plot(1,main=substitute("Data at "*lat*degree*" North",li...
(from Thomas Lumley)
You can use different way: First export plot for R as an ...
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 P...
plot(x[!is.na(z)],y[!is.na(z)],xlim=range(x),ylim=range(...
points(x[is.na(z)],y[is.na(z)],col=2) (from Ross Darnell)
**5.36 lattice: その他の例 [#bc8fdfd6]
//misc examples
(12/11/2002 Paul Johnson <pauljohn@ku.edu>)
The lattice library came along recently, I've not explore...
library(lattice)
# 10 rows
x <- matrix(rnorm(100), ncol=10)
lp1 <- levelplot(x, colorkey=list(space="bottom"), aspe...
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...
print(lp2)
This one from Deepayan Sarkar (deepayan@stat.wisc.edu) is...
library(lattice)
xyplot(rnorm(100) ~ rnorm(100) | gl(2, 50),
strip = function(factor.levels, ...) {
strip.default(factor.levels = expression(1 * degree, 2 *...
...)
})
**5.37 三次元の散布図 [#x867c42a]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
以下をタイプすると
//Type
?persp
?contour
ベースパッケージにあるこれらの特徴を参照できる。
//to see about features in the base package.
CRAN には、scatterplot3d パッケージがある。
//There is a scatterplot3d package on CRAN.
以下も参照のこと：
//See also:
http://www.statistik.uni-dortmund.de/leute/ligges.htm (f...
以下もまた参照のこと：
//See also:
Get Xgobi (an X windows graphics library) and the R contr...
**5.38 値を反映した線種を用いた三次元等高線図 [#gd556439]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
正の値が実線で、負の値が点線のコンターを描画したい場合,
//Suppose you want a contour plot solid contour lines for...
種明かしをしてみると levels で間隔を指定して、 add= 引数...
//The trick is to specify the levels and to use the 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), lt...
contour(z, levels = 0, lty = "dashed", add = TRUE)
contour(z, levels = seq(0.1, 1, by = 0.1), add = ...
(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...
このようにしておけば，例えば色も変えられる
contour(z, levels = k, lty = c("dotted", "dashed"...
**5.39 グラフのアニメーション [#d937e3ff]
(13/08/2000 Paul Johnson <pauljohn@ukans.edu>)
ImageMagick (www.imagemagick.org)を使えば本当に容易にでき...
//It is quite easy to do with ImageMagick (www.imagemagic...
Rで、１５個のヒストグラムを作成する。
//In R, create 15 histograms:
> for(i in 1:5) {
+ jpeg(paste("fig", i, ".jpg", sep = ""))
+ hist(rnorm(100))
+ dev.off()
+ }
それから、コマンドライン（わたしはLinuxで試したが、他のプ...
//Then from the command line (I tried it using Linux, tho...
% convert -delay 50 -loop 50 fig*.jpg animated.gif
これは animated.gif を作成した。これは私の一連のファイル...
//This created animated.gif, a nice animation of my seque...
**5.40 グラフのユーザー領域の背景色を余白と異なったものに...
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
"par" を使って、 boundaries/axes の描画位置を見つけ出し、...
//Use "par" to find out where the boundaries/axes are, an...
//For example:
例えば，
## 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(co...
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 実際に動くグラフコード幾つか (その他) [#x24d401b]
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
*** 点プロット、直線上書き、凡例、PS 出力 [#ldad56be]
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 関数の使用例 [#h...
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 =...
xlab = "x", ylab = "sin(x)", main = "Y = sin(x)...
axis( 2, at = seq( -1, 1, by=1 ),las = 2 ) box(lty="dott...
plot( cos(x), type = "l", xaxs = "i", yaxs = "i", axes =...
xlab = "x", ylab = "cos(x)", main = "Y = cos(x)...
axis( 2, at = seq( -1, 1, by=1 ),las = 2 )
box(lty="dotted")
abline( h = 0, lwd = 1 )
(from Ko-Kang Wang)
***正多角形のプロット [#hb2462cc]
Below is a function that generates regular polygons, fill...
borders, of n sides (n>8 => circle), with "diameter" pres...
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/...
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)
***自前の軸 [#he47dc08]
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)
***プロットの中心に軸を持つ三次関数を描く [#o7036612]
# Draw the basic curve (limits were established by tri...
# The axes are turned off blank axis labels used.
curve(x * (2 * x - 3) * (2 * x + 3), from = -1.85, to =...
xlim = c(-2, 2), ylim = c(-10, 10), axes = FALSE, xla...
# Draw the axes at the specified positions (crossing at...
# The ticks positions are specified (those at (0,0) are...
# Note the use of las=1 to produce horizontal tick labe...
# 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...
title(main = expression(italic(y==x * (2 * x - 3) * (2 ...
(from Ross Ihaka)
***3次元曲面を等高線プロットの様に塗分ける、例えば大きな ...
# 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 fa...
hgt <- 0.25 * (z[-nx,-ny] + z[-1,-ny] + z[-nx,-1] + z[-1...
hgt <- (hgt - min(hgt))/ (max(hgt) - min(hgt))
*** 各面の色を指定した曲面を描く [#u6306b4c]
persp(x, y, z, col = gray(1 - hgt), theta = 35)
persp(x, y, z, col = cm.colors(10)[floor(9 * hgt + 1)], ...
(from Ross Ihaka)
*** polygon 関数使用例 [#lcf2ddcf]
y <- c(1.84147098480790, 1.90929742682568, 1.14112000805...
0.243197504692072, -0.958924274663133,0.7205845...
1.65698659871879, 1.98935824662338, 1.4121184852...
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 による並置出力 [#oe4ef49c]
(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]] へ移動
終了行:
COLOR(red){SIZE(25){Paul Johnson 氏の R Tips 集訳 1/2}}
後半 [[P_Johnson_tips_2]] へ移動~
2003/08/09 現在~
注：これは P. Johnson さんが r-help の記事から役にたちそ...
注意：あまりに長くて表示に時間がかかるので二つに分割しま...
----
(私はこれを "StatsRus" と名付けるのが気がきいていると思っ...
もしあなたが R の Tips 集が必要なら、これがそうです。これ...
Brian D. Ripley の言葉を注意しましょう、
"One enquiry came to me yesterday which suggested that so...
distributions do not know that R comes with two Guides in...
directory plus an FAQ and the help pages in book form. I...
distributed with all the binary distributions (they are n...
installed by default). Windows-specific versions are ava...
rw1001d1.zip and rw1001d2.zip"
FaqManager ソフトに対しその著者の Stas Bekman に謝辞を捧...
#contents
~
* 1 データの入出力 [#k0b83b2c]
**1.1 生の数値を R に取り込む [#t84e9115]
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
これは本当に簡単。ファイル "myData"で、空白区切りの数値が...
//This is truly easy. Suppose you've got numbers in a spa...
myDataFrame<-read.table("myData",header=TRUE)
このコマンドは素晴らしく動く。
//command works great.
"?read.table" とタイプすれば、他の区切り記号を使っている...
//If you type "?read.table" it tells about importing file...
タブ区切りのデータを持っていて、空白を欠損値を示すのに使...
//Suppose you have tab delimited data with blank spaces t...
> myDataFrame<-read.table("myData",sep="\t",na.strings="...
(R1.3 の新しいフィーチャー) データが gzipで圧縮されたフ...
//(new R1.3 feature) Suppose your data is in a gzipped fi...
myDataFrame <- read.table(gzfile("myData.gz"), header=T)
別々のファイルの列を読み込む際には、以下のようにして、こ...
//If you read columns in from separate files, combine int...
variable1 <- scan("file1")
variable2 <- scan("file2")
mydata <- cbind(variable1, variable2)
もしくは、以下を使えば、同様なことができる。
//or use the equivalent:
mydata <- data.frame(variable1, variable2)
R オブジェクトのデータフレームを、以下のコマンドで保存す...
//Optionally save dataframe in R object file with:
write.table(mydata, file="filename3")
**1.2 データアクセスに関する基本的注意 [#y50a7128]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
データフレーム "x" の列を列番号でアクセスするには、例えば...
//To access the columns of a data frame "x" with the colu...
注：原文にも抜けているのだが，x["pjVar1"] は x$pjVar1 と...
注の注：「第1列の名前が"pjVar1"であることを知っていれば」...
もし，そのデータフレームに対して attach(x) を行うと，単に...
//If you attach that data frame with attach(x) then you c...
回帰分析において，data=x オプションを含めることにより，あ...
//In a regression, include the option data=x and then you...
リストの要素を得るためには x[[1]] と...
//もとの和訳ではそこを考慮していなかったので，意味不明の...
//Grab an element in a list as x[[1]]
例えば，もしデータフレームが "nebdata" ならば，
//For instance, if a data frame is "nebdata" then
nebdata[1, 2]
あるいは
//or
nebdata[nebdata$Volt=="ABal", 2]
あるいは，データフレームにアタッチすることで，
//or by attaching the data frame
attach(nebdata)
nebdata[Volt=="ABal", 2]
detach(nebdata)
(from Diego Kuonen)
**1.3 新しい Data Import/Export マニュアルのチェック [#r9...
(13/08/2001 Paul Johnson <pauljohn@ukans.edu>)
R-1.2 において，R チームは R にデータを読み込んだり R か...
//With R-1.2, the R team released a manual about how to g...
**1.4 データを R と他のプログラム (Excelなど)の間で交換す...
(04/05/2001 Paul Johnson <pauljohn@ukans.edu>)
私はこの覚え書きを新しい R のデータインポート・データエク...
//I made these notes before the new R Data Import/Export ...
Excel -> R
エクセルで，例えば table.txt などという名前のテキストファ...
//Use 'save as space delimited text (.prn)' in Excel to c...
R -> Excel
テキストファイルに書き出すときには，sink("table") とする...
//Use sink("table"), where 'table' is the name of the R o...
(from Griffith Feeney)
R-Excel インターフェースは試したかい？
http://lib.stat.cmu.edu/R/CRAN/contrib/extra/index.html ...
(from Charles Raux)
//Did you try R-Excel interface ? see http://lib.stat.cmu...
私の学生達に教えようと思っていることを書いておく（インポ...
//Here's what I'm planning to teach my students (Import o...
今のところ，最も簡単なやり方は，えとね，データをテキスト...
//Currently the easiest way is to to ... export data as a...
現在の最も普遍的な三つの PC プログラムにおけるテクニック...
//最小公分母というのは，日本ではあまり使わないな
//The current techniques for the three most common PC pro...
SPSS：データエディタの File, Export メニューにより，タブ...
//SPSS: Use |File, Export| menu in the data editor and ex...
このやり方はカンマが小数点として使われているかどうか loca...
//This refers to locales where the comma is used as decim...
Excel: File, Save as... を使い，タブ区切り（.txt）粗選択...
//Excel: Use |File,Save as...| and choose tab-sep (|.txt|...
SAS: カンマ区切りかタブ区切りでエクスポートできる。SAS は...
//SAS: It is possible to export in a comma-separated or t...
(from Peter Dalgaard)
**1.5 データフレームのマージ [#y72d434c]
(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)
experiment<-data.frame(times=c(0,0,10,10,20,20,30,30),ex...
simul<-data.frame(times=c(0,10,20,30),simul=c(3,4,5,6))
以下のように、マージされたデータセットが欲しい。
//I want a merged datatset like:
> 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
答え
//Answers
merge(experiment, simul)
が，貴方の望む全てだろう。
//does all the work for you.
(from Brian D. Ripley)
あるいは，以下を検討する
//Or consider:
exp.sim<-data.frame(experiment,simul=
simul$simul[match(experiment$times,simul$...
(from Jim Lemon)
このようなデータフレームを持っている。
//I have dataframes like this:
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
以下が欲しい
//And I want
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 のレスは以下のとおり
//In response, Ben Bolker said
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を利用した結果は，以下のようにな...
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
（訳注その２）次のようにしたらうまくいくけど，どうでしょ...
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...
merge(data1,data2,by=1)
**1.6 一時に一つの行を加える [#q885c512]
(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)
Question: I would like to create an (empty) data frame wi...
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 ...
If you know the number of results in advance, say, N, do ...
df <- data.frame(time=numeric(N), temp=numeric(N), press...
df[1, ] <- c(0, 100, 80)
df[2, ] <- c(10, 110, 87)
...
または
//or
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 on...
If you don't know the final size you can use rbind to add...
[top]
**1.7 データフレームに対するもう一つのマージ法 [#v32711f4]
(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 オブジェクトがあるかどうか検査する [#p66901e7]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
R がロードする際に、R はロードするものがないとき警告は出...
//If you load things, R does not warn you when they are n...
is.null(list$component)
上のコマンドは list という名のものにある component と名づ...
//to check a thing named component in a thing named list.
存在しないデータフレームのカラムを "[" を使ってアクセスす...
//Accessing non-existent dataframe columns with "[" does ...
>data(trees)
>trees$aardvark
NULL
>trees[,"aardvark"]
Error in [.data.frame(trees, , "aardvark") : subscript o...
(from Thomas Lumley)
**1.9 乱数を発生する [#c21fab27]
(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 ...
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 ...
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 平均/分散を固定して乱数を発生 [#bbb2a147]
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
任意のツールで乱数を生成する場合、指定したちょうどの平均...
//If you generate random numbers with a given tool, you d...
そこで、平均０のサンプルをどのようにして強制的に作成する...
//So, how to force the sample to have a mean of 0? Take a...
> 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
そして平均５、標準偏差２のサンプルができる。
//and now create your sample with mean 5 and sd 2:
# 以下は同じ
> x <- x*2 + 5
> mean(x)
[1] 5
> var(x)
[1] 4
(from Torsten.Hothorn)
**1.11 重みつきの複製 [#z9bcb072]
(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...
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 freque...
If you have frequency weights you may have to expand your...
expanded.data<-as.data.frame(lapply(compressed.data,
function(x) rep(x,compressed.data$weig...
should do it.
(from Thomas Lumley)
**1.13 分割表をデータフレームに変換する [#s9c0f947]
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
もし８次元のクロス表があるなら、表にセルの８つのファクタ...
//Given a 8 dimensional crosstab, you want a data frame w...
R1.2 はこれを扱うのに as.data.frame.table() 関数を導入し...
//R1.2 introduces a function as.data.frame.table() to han...
この処理はマニュアルでもできる。ここに関数がある（expand....
//This can also be done manually. Here's a function(it's ...
dfify <- function(arr, value.name = "value", dn.names =
names(dimnames(arr))) {
Version <- "$Id: dfify.sfun,v 1.1 1995/10/09 16:06:12 d3...
dn <- dimnames(arr <- as.array(arr))
if(is.null(dn)) stop("Can't data-frame-ify an array with...
names(dn) <- dn.names
ans <- cbind(expand.grid(dn), as.vector(arr))
names(ans)[ncol(ans)] <- value.name ans }
関数名は "data-frame-i-fy" を縮めたものである。
//The name is short for "data-frame-i-fy".
質問の例では、多重配列が適切な dimnames を持っていると推...
//For your example, assuming your multi-way array has pro...
my.data.frame <- dfify(my.array, value.name="frequency")
(from Todd Taylor)
**1.14 書き込み: アスキーファイルのデータ [#o0598ccd]
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
たとえば 28 x 28 データマトリックスを生成するコマンドを持...
//Say I have a command that produced a 28 x 28 data matri...
write.table(mat,file="filename")
Pesl Thomas は「"tab"という名前のデータフレームの内容をフ...
//Pesl Thomas said, "I want to output a file to my disc w...
> tab
or 1 2 3 4 5 6 1 57 80 25 23 46 23 2 106 35 59 8 51 4
これが使える
//You can use
write.table(unclass(tab))
これも使える
//You could also do
write(t(tab),ncol=NCOL(tab))
なぜなら、 write() は行と列を転置する。(from Thomas Lumley)
//since write() transposes rows and columns. (from Thomas...
注目してもらいたいのは、MASS ライブラリには、write.matrix...
//Note MASS library has a function write.matrix which is ...
*2 データフレームを使った作業: 記録、合体 [#d7e9a175]
**2.1 変数をデータフレーム(又はリスト)に加える [#l5ca6160]
(02/06/2003 Paul Johnson <pauljohn@ku.edu>)
In r-help, I asked " I keep finding myself in a situation...
> 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 ...
eval(substitute(myList$newColumn<-aColumn),
list(newColumn=as.name(varName)))
as an alternative to parse().
-(Thomas Lumley)
**2.2 変数名を「その場で」作る [#oe2e9c8b]
(10/04/2001 Paul Johnson <pauljohn@ukans.edu>)
変数名とそれにセットする値を一緒に貼り付ける。assign を使...
//Paste together a variable name, set it to a value. Use ...
> assign(paste("file", 1, "max", sep=""), 1)
> ls()
[1] "file1max"
(Brian Ripley, June 18, 2001)
**2.3 一つの列を記録する、値を別の列に出力する [#ofcc96a7]
??Recode one column, output values into another column
(20/06/2001 Paul Johnson <pauljohn@ukans.edu>)
Aleksey Naumov がこの質問をした。例えば、ここに簡単なデー...
//Aleksey Naumov asked this question. For example, here i...
var1
1 1
2 2
3 3
どのようにして、上の "var1" をグループ化した値 'var2' を...
//How do I add a new column to data - 'var2' in which I g...
var1 var2
1 1 4
2 2 4
3 3 5
data$var2 <- c(4,4,5)[data$var1]
または
//or
data <- transform(data, var2=c(4,4,5)[var1])
(from Peter Dalgaard)
I admit the previous seems complicated, and I'm inclined ...
x <- rnorm(10000)
x[x > 1.5] <- NA
And if you don't want to replace the original variable, c...
You can put other variables inside the brackets, so if yo...
x[y == 1] <- 234
So, to do the v1 v2 example above, I think it is clearest...
> 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 bel...
For example, suppose I get a bunch of variables coded on ...
1= no
6 = yes
8 = tied
9 = missing
10 = not applicable.
Recode that into a new variable name with 0=no, 1=yes, an...
It seems like the replace() function would do it for sing...
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...
test <- replace(test, test == 1, 0)
test <- replace(test, test == 2, 1)
test
I suppose a better function would take from and to lists ...
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 wa...
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...
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...
**2.4 指示(ダミー)変数をつくり出す [#lb43f1a5]
//Create indicator (dummy) variables
(20/06/2001 Paul Johnson <pauljohn@ukans.edu>)
以下に2例ある。
c is a column, you want dummy variable, one for each vali...
> 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 )
質問；５つのカテゴリーのアル変数があり、各カテゴリに対す...
//Question: I have a variable with 5 categories and I wan...
答え；行の添え字または model.matrix を使う。
//Answer: Use row indexing or model.matrix.
ff <- factor(sample(letters[1:5], 25, replace=TRUE))
diag(nlevels(ff))[ff,]
または
//or
model.matrix(~ ff - 1)
(from Brian D. Ripley)
**2.5 時系列回帰のための変数の遅延値をつくり出す [#vf0da2...
//Create lagged values of variables for time series regre...
(02/06/2003 Paul Johnson <pauljohn@ku.edu>)
Peter Dalgaard explained, "the simple way is to create a ...
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 関連する観測値を持たない因子水準をデータセットから...
//How to drop factor levels for datasets that don't have ...
(08/01/2002 Paul Johnson <pauljohn@ukans.edu>)
もっともいいのが、水準を落とすことである。たとえば、以下...
//The best way to drop levels, BTW, is
problem.factor <- problem.factor[,drop=TRUE] (from Brian...
**2.7 データフレームから観測値を選択する、部分集合を取り...
//Select/subset observations out of a dataframe
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
ある変数が欠損値である行を省きたいなら、John Fox の言うよ...
//Want to drop out observations that are mising on a vari...
x2 <- na.omit(x)
これにより、欠損値を含む行すべてが除かれたデータフレーム ...
//produces a copy of the data frame x with all rows that ...
欠損値の除外については、Peter Dalgaard は以下のコマンドを...
//For exclusion of missing, Peter Dalgaard likes
subset(x,complete.cases(x)) or x[complete.cases(x),]
//"is.na(x)”を追加することよりも、 x != "NA" の方が望ま...
x != "NA" を追加するよりも、 "is.na(x)” の方が望ましい。
//adding "is.na(x) is preferable to x != "NA"
Suppose you want observations with c=1 in df1. This makes...
df2 <- df1[df1$c == 1,]
and note that indexing is pretty central to using S (the ...
Or use "match" select values from the column "d" by takin...
> 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 lik...
> 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 ov...
maleOver40 <- subset(d, sex %in% c("m","M") & age > 40)
maleOver40 <- d[(d$sex == "m" | d$sex == "M") & d$age > ...
**2.8 最初の観測値をそれぞれから抹消する [#k5d732bf]
//Delete first observation for each
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
Given data like:
1 ABK 19910711 11.1867461 0.0000000 ...
2 ABK 19910712 11.5298979 11.1867461 111
6 CSCO 19910102 0.1553819 0.0000000 ...
7 CSCO 19910103 0.1527778 0.1458333 ...
remove the first observation for each value of the "sym" ...
newhilodata <- hilodata[-c(1,6,13),]
To solve the more general problem of omitting the first i...
newhilodata <- subset(hilodata, diff(c(0,as.integer(sym)...
(actually, the as.integer is unnecessary because the c() ...
(from Peter Dalgaard)
Alternatively, you could use the match function because i...
> 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 ...
subset(hilodata, duplicated(sym))
has got to be the simplest variant.
(from Peter Dalgaard)
**2.9 データのランダム標本を選ぶ [#r4bf29d7]
//Select a random sample of data
(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 関数を忘れるな [#j1cad700]
//Selecting Variables for Models: Don't forget the subset...
(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)
直接行その他を削除することでデータを管理できるが、subset(...
//You can manage data directly by deleting lines or so fo...
Peter Dalgaard が「組み込み」データセット airquality を使...
//Peter Dalgaard gave this example, using the "builtin" d...
data(airquality)
names(airquality)
lm(Ozone~.,data=subset(airquality,select=Ozone:Month))
lm(Ozone~.,data=subset(airquality,select=c(Ozone:Wind,Mo...
lm(Ozone~.-Temp,data=subset(airquality,select=Ozone:Mont...
lm の RHS の "." は「すべての変数」を意味しており、rhs の...
//The "." on the RHS of lm means "all variables" and the ...
**2.11 すべての数値変数を処理し、文字変数を無視するには？...
??Process all numeric variables, ignore character variabl...
(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),...
mean<-sapply(data[,indices],mean) sd<-sapply(data[,indic...
**2.12 二つ以上の変数に関しソート [#x55a44eb]
//Sorting by more than one variable
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
Can someone tell me how can I sort a list that contains d...
name M
1234 8
1234 8.3
4321 9
4321 8.1
I also would like to set a cut-off, so that anything belo...
I take it that the cutoff is on the value of M. OK, suppo...
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 im...
If your two variables are in a data frame you can combine...
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 ...
dat <- dat[order(dat$name, pmax(coff, dat$M)), ]
(from Bill Venables)
I want the ability to sort a data frame lexicographically...
Here's how:
spsheet[order(name,age,zip),]
(from Peter Dalgaard)
**2.13 ある因子により定義されるサブグループ内で順位を取る...
//Rank within subgroups defined by a factor
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
by() のヘルプを読みなさい
//Read the help for 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 ...
//Use na.omit to screen missings from a data frame:
(22/08/2000 Paul Johnson <pauljohn@ukans.edu>)
欠損値が "insure" という名前のデータセットより除去される...
//To make sure missings are omitted from a dataset called...
thisdf <- na.omit(insure[,c(1, 19:39)])
body.m <- lm(BI.PPrem ~ ., data=thisdf, na.action=na.fail)
**2.15 各行から一つの集約値を計算 [#d9eb624c]
//Aggregate values, one for each line
(16/08/2000 Paul Johnson <pauljohn@ukans.edu>)
質問:200行で，各行に32個の実数データを含むテキストファイ...
//Question: I want to read values from a text file - 200 ...
Peter Dalgaard は以下を提案した。
//Peter Dalgaard says do thi:
y <- apply(as.matrix(read.table("myfile")), 1, mean)
x <- seq(along=y)
（ファイル形式によっては，read.table に，いくつかのオプシ...
//(possibly adding a couple of options to read.table, dep...
**2.16 各因子に対する集約値を保持する新しいデータフレーム...
//Create new data frame to hold aggregate values for each...
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
これでどう
//How about this:
Z<-aggregate(X, f, sum)
(X の変数はすべて合計できるものとする)
//(assumes all X variables can be summed)
または、（もし）Xが因子も含んでいる場合、要約統計の計算に...
//Or: [If] X contains also factors. I have to select vari...
Z <- data.frame(f=levels(f),x1=as.vector(tapply(x1,f,sum...
(from Wolfgang Koller)
**2.17 データフレームを選択的に集約 [#d43804cb]
//Selectively aggregate a data frame
(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),...
//note the first two columns are excluded by [,(1:2)] and...
最初の2列は [,-(1:2)] で除外され，by オプションで出力のた...
**2.18 実数から各桁数値を一時に一つずつ取り出す [#yc94b1e9]
//Rip digits out of real numbers one at a time
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
//I want to "take out" the first decimal place of each ou...
小数点以下1桁目の数値を取り出し出現度数をプロットする。次...
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 リスト中の各副行列から一つの項目を取り出す [#t2307...
//Grab an item from each of several matrices in a List
(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)
Let Z denote the list of matrices. All matrices have the ...
lapply(Z, function(x) x[1,2])
should do this, giving a list. Use sapply if you want a v...
**2.20 データセット中の値を示すベクトルを得る [#y017ddba]
//Get vector showing values in a dataset
(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の命令を表す文字列の値を計算する [#t38d404f]
//Calculate the value of a string representing an R comma...
(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 m...
Or
eval(parse(text="ls()"))
Or
eval(parse(text = "x[3] <- 5"))
**2.22 "which" はあるテストをパスする観測値の添字値を抜き...
//"Which" can grab the index values of cases satisfying a...
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
大きなベクトルのデータを箱ひげ図で分析して、外れ値を見つ...
//To analyse large vectors of data using boxplot to find ...
which(x==boxplot(x,range=1)$out)
**2.23 行列/データフレーム中の一意的な行を見出す [#k57449...
//Find unique lines in a matrix/data frame
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
Jason Liao writes:
"I have 10000 integer triplets stored in A[1:10000, 1:3]....
Peter Dalgaard answers, "As of 1.4.0 (!):
unique(as.data.frame(A))"
**2.24 数値変数だけを選ぶ [#h40cdbff]
//Select only numerical variables
(12/11/2002 Paul Johnson <pauljohn@ku.edu>)
sapply(dataframe, is.numeric)
または
//or
which(sapply(data.frame, is.numeric))
Thomas Lumley 11/11/2002
* 3 行列とベクトルの演算 [#yb7940fe]
**3.1 単位行列をつくり出すには？ [#pa2417c8]
//How to create an identity matrix?
(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" を一つの長いベクトルに変換 [#z97db96f]
//Convert matrix "m" to one long vector
(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) をつくり出す [#e...
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
なぜこれが有用なのかは分かりませんが、いくつかの列と行列...
//I don't know why this is useful, but it shows some row ...
結局は Brian Ripley の方法を使うことになった。この方法は...
//I ended up using Brian Ripley's method, as I got it fir...
A <- matrix(1, n-1, n-1)
rA <- row(A)
rA[rA + col(A) <= n]
しかしながら、 その後、Andy Royle の御蔭で、もっと簡単で...
//However, thanks to Andy Royle I have since discovered t...
// sequence(n-1:1) もとからこうなっていたのだろうが，良...
sequence((n-1):1)
[1] 1 2 3 4 1 2 3 1 2 1
Karen Kotschy
**3.4 すべての n 個めの項目を選ぶ [#o8322b1e]
//Select every n'th item
(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)
非常に長いベクトルよりｎ個目の要素をすべて取り出す、これ...
//extract every nth element from a very long vector, vec?
vec[seq(n, length(vec), n)] #(from Brian Ripley)
seq(1,11628,length=1000) はインデックスとして使える 1:116...
//seq(1,11628,length=1000) will give 1000 evenly spaced n...
My example:
vec <- rnorm(1999)
newvec <- vec[1, length(vec), 200]
> newvec
[1] 0.2685562 1.8336569 0.1371113 0.2204333 -1.2798...
[7] -0.2366385 0.5060078 0.9680530 1.2725744
This shows the items in vec at indexes (1, 201, 401, ...,...
**3.5 ベクトル中の 1.5 に最も近い値の添字を見出す [#b9a08...
//Find index of a value nearest to 1.5 in a vector
(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 ベクトル中の非負値項目の添字を見つける [#ha85c609]
//Find index of nonnegative items in vector
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
which(x!=0) (from Uwe Ligges)
rfind <- function(x) seq(along=x)[x != 0] #(from Brian ...
スピードとメモリーの効率性に関しては、以下がある。
//Concerning speed and memory efficiency I find
as.logical(x)
これは以下よりもよく
//is better than
x!=0
そして
//and
seq(along=x)[as.logical(x)]
これは以下よりもよく
//is better than
which(as.logical(x))
このため
//thus
which(x!=0)
は最も短いが
//is shortest and
rfind <- function(x)seq(along=x)[as.logical(x)]
以下はコンピュータではもっとも効率的であるようである。
//seems to be computationally most efficient
(from Jens Oehlschl?gel-Akiyoshi)
**3.7 欠損値の添字を見出す [#j42f9f57]
//Find index of missing values
(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)
//Suppose the vector "Pes" has 600 observations. Don't do...
Pes というベクトルが600の観測値を持つとする。欠損値の添え...
(1:600)[is.na(Pes)]
//The `approved' method is
正しいやり方は，次の通り。
seq(along=Pes)[is.na(Pes)]
//In this case it does not matter as the subscript is of ...
この場合にはどちらでも良いことではあるが，添え字の長さが0...
(from Brian Ripley)
//However, the solution I gave
しかし，私は以下の解を与える
which(is.na(Pes))
is the one I still really recommend; it does deal with 0-...
**3.8 ベクトル中の最大値を持つ項の添字を見出す [#hdd3d6af]
//Find index of largest item in vector
(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 行列中の値を変更する [#ad45fc77]
//Replace values in a matrix
(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 がデータフレームなら、データを強制的にマトリックス形式...
//If la is a data frame, you have to coerce the data into...
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
以下を試してみよう
//Try this:
la <- ifelse(la == 0, 1, la)
(from Martyn Plummer)
日本語訳を見てのコメント~
la がデータフレームのとき，何もせずいきなり la[la==0] <- ...
**3.10 行列から特定の行を取り除く [#rdc52cd6]
//Delete particular rows from matrix
(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 ある基準を満たす項目の総数を数える [#a8a81810]
//Count number of items meeting a criterion
(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)
以前の質問で述べた which() の結果に"length()" を以下のよ...
//Apply "length()" to results of which() described in pre...
length(which(v<7))
または
//or
sum( v<7,na.rm=TRUE)
**3.12 相関行列から偏相関係数を計算する [#pcb94fcf]
//Compute partial correlation coefficients from correlati...
(08/12/2000 Paul Johnson <pauljohn@ukans.edu>)
多変量変数間の相関係数から，二変数以外の変数を制御して二...
//I need to compute partial correlation coefficients betw...
実際問題として，これはきわめて簡明直裁。R が相関係数行列...
//Actually, this is quite straightforward. Suppose that R...
Rinv<-solve(R)
D<-diag(1/sqrt(diag(Rinv)))
P<- -D %*% Rinv %*% D
P の非対角要素は他の変数の影響を除いた（コントロールした...
//The off-diagonal elements of P are then the partial cor...
(from John Fox )
一般的に，分散・共分散行列の逆行列を求めて，対角成分が1に...
//In general you invert the variance-covariance matrix an...
pcor2 <- function(x){
conc <- solve(var(x))
resid.sd <- 1/sqrt(diag(conc))
pcc <- - sweep(sweep(conc, 1, resid.sd, "*"), 2, resid...
return(pcc)
}
pcor2(cbind(x1,x2,x3))
see J. Whittaker's book "Graphical models in applied mult...
以下は，私が今使っているもので，それぞれの偏相関係数の有...
//This is the version I'm using now, together with a test...
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...
"*")
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 -...
df = t.df))
}
return(pcc)
}
(from Kaspar Pflugshaupt)
**3.13 正準変数を計算 [#y693ff4f]
//Compute Canonical variates
(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)
同じ観察対象について，データ行列 A と B があるとする。
//Suppose you have data matrices A and B on the same obse...
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=Bp...
}
この関数には，トリックがある。正準変量の平均は0で，分散は...
//should do the trick. The canonical variates are zero-me...
**3.14 多次元行列(R 配列)をつくり出す [#nf518cfb]
//Create a multidimensional matrix (R array)
(20/06/2001 ? <?>)
//Brian Ripley said:
Brian Ripley の話では
my.array<-array(0,dim=c(10,5,6,8))
//will give you a 4-dimensional 10 x 5 x 6 x 8 array.
上のようにすれば、 10 x 5 x 6 x 8 の 4次元の配列が与えら...
//Or
または以下のような方法もある。
array.test <- array(1:64,c(4,4,4))
array.test[1,1,1]
1
array.test[4,4,4]
64
**3.15 多数の行列を結合する [#vf71c9d3]
//Combine a lot of matrices
(20/06/2001 ? <?>)
100個の n×2 行列のリストがあるとする。どうやったら，これ...
//Given a list with 100s of (nx2) matrices, how do you co...
list は行列のリストである。
//list is the list of matrices.
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 特定の論理に基づき "近傍" 行列をつくり出す [#kd19f...
//Create "neighbor" matrices according to specific logics...
(20/06/2001 ? <?>)
あるセルがある位置に近傍があるかどうかの有無を０と１で示...
//Want a matrix of 0s and 1s indicating whether a cell ha...
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(...
bishop.case <- function(y,z,mat){coldiff(y,z,mat)==1 & r...
queen.case <- function(y,z,mat){rook.case(y,z,mat) | bis...
matrix(as.numeric(
sapply(x,function(y)sapply(x,rook.case,y,mat...
matrix(as.numeric(
sapply(x,function(y)sapply(x,bishop.case,y,m...
matrix(as.numeric(
sapply(x,function(y)sapply(x,queen.case,y,ma...
(from Ben Bolker)
**3.17 ある "key" により二つの数値列をマッチング [#s3f21f...
//"Matching" two columns of numbers by a "key"
(20/06/2001 ? <?>)
The question was:
I have a matrix of predictions from an proportional odds ...
小さな例として，以下の x, y
//As a toy example, if I have
> 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] (t...
Is
x[cbind(seq(along=y), y)]
what you had in mind? The key is definitely matrix indexi...
**3.18 上・下三角行列をつくり出す [#idb327b7]
//Create Upper or Lower Triangular matrix
(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...
> 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 の逆行列を計算 [#j6878b82]
//Calculate inverse of X(20/06/2001 ? <?>)
solve(A)は A の逆行列を返す
**3.20 行列添字の面白い使い方 [#r4d2915e]
//Interesting use of Matrix Indices
(20/06/2001 ? <?>)
Here's a similar problem.
Mehdi Ghafariyan said "I have two vectors A=c(5,2,2,3,3,2...
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...
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 固有値の例 [#p3eb8b69]
//Eigenvalues example
(20/06/2001 ? <?>)
Tapio Nummi asked about double-checking results of spectr...
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 の適用等 [#wedb7571]
**4.1 一つの関数から複数の値を返す [#dfe54df3]
//Return multiple values from a function
(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)
Below is a summary on returning variable from subfunction...
To access the variables within a function, return a list ...
"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" 値を掴み出す [#v0bd0df1]
//values out of a list of significance tests
(22/08/2000 Paul Johnson <pauljohn@ukans.edu>)
//Suppose chisq1M11 is a list of htest objects, and you w...
chisq1M11 が htest オブジェクトのリストだとする。そして，...
> apply(cbind(1:1000), 1, function(i) chisq1M11[[i]]$p.v...
//And Peter Dalgaard showed a more elegant approach:
Peter Dalgaard はもっとエレガントなやり方を示した。
sapply(chisq1M11,function(x)x$p.value)
このどちらも，"function" と呼ばれる関数を作り，それぞれの...
（訳注：本当は無名の関数を作っているのだが）
//In each of these, a simple R function called "function"...
**4.3 ifelse の用法 [#zca930f7]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
//If you want to calculate something, but only if y is *b...
y が 0 から 1 の範囲にあるときにのみある計算をしたいとき...
ifelse(y==0,0,y*log(y))
訳注：上のは，y が 0 であるときとそれ以外で別の計算をする...
(from Ben Bolker)
**4.4 apply を用いて各セル毎の確率の行列を作る [#je988350]
//Apply to create matrix of probabilities, one for each "...
(14/08/2000 Paul Johnson <pauljohn@ukans.edu>)
ガンマ関数のスケールパラメータと形状パラメータの様々な値...
//Suppose you want to calculate the gamma density for var...
gridS <- expand.grid(scales, shapes)
survLilel <- function(ss) sum(dgamma(Survival,ss[1],ss[2...
Likel <- apply(gridS,1,survLilel)
実際には私は以下のように使う
//Actually, I would use
sc <- 8:12; sh <- 7:12
args <- expand.grid(scale=sc, shape=sh)
matrix(apply(args, 1, function(x) sum(dgamma(Survival, s...
log=T))), length(sc), dimnames=list(scale=sc, shape=s...
(from Brian Ripley)
**4.5 外積 outer [#j2426a97]
(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)
outer(x,y,f) does just one call to f with arguments creat...
**4.6 あるものが公式か関数か検査する [#vd03b825]
//Check if something is a formula/function
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
公式（モデル式）は "formula" クラスを持つので，inherits(o...
//Formulae have class "formula", so inherits(obj, "formul...
もし，Y~X+Z+W ではなくて ~X+Z+W であることが分かっている...
//And if you want to ensure that it is ~X+Z+W rather than...
inherits(obj,"formula") && (length(obj)==2)
(from Thomas Lumley)
**4.7 変数のベクトルに関する最適化 [#h99431a3]
//Optimize with a vector of variables
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
The function being optimized has to be a function of a si...
opt.func <- function(arg)
t(Y-(X[,1] * arg[1] + X[,2] * arg[2])^delta) %*% covar...
(Y-(X[,1] * arg[1] + X[,2] * arg[2])^delta)
(from Douglas Bates)
**4.8 S+ のような slice.index [#t3cf7951]
(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)
slice.index<-function(x,i){k<-dim(x)[i];sweep(x,i,1:k,fu...
(from Peter Dalgaard)
*5 グラフ [#l8916f99]
**5.1 グラフ描きの前に par で特性を変更 [#m4ae1348]
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
par() はグラフの品質に影響を与える多目的コマンド。例を以...
//par() is a multipurpose command to affect qualities of ...
***1. 以降のプロットを一つの絵にまとめる: [#s88a8fc8]
par(mfrow=c(3,2))
プロットの 3x2 のマトリックスを作成する。
//creates a 3x2 matrix of plots.
***2. 余白を次の例のように変更 [#ka7a90a6]
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 グラフ出力を保存 [#c7347e49]
(15/08/2001 Paul Johnson <pauljohn@ukans.edu>)
グラフ保存には２つの方法がある。グラフの保存前に、出力形...
//There are two ways to save and/or print graphs. Before ...
There are other devices representing kinds of file output...
>?device
Rがシステムで見つけるデバイスのリストが見れる。
//to see a list of devices R finds on your system.
各デバイスについて、より多くを見つけられる。
//For each device, you can find out more about it.
?x11
または
//or
?png
また MS Windowsでは
//or, in MS Windows
?windows
(等々)
//(and so forth)
Before you start using a device, you can configure it. Th...
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 wh...
> dev.set(3)
Now here is the big thing. The devices are configured dif...
Suppose we want png output. We create an instance of a pn...
> png(filename="mypicture.png", width=480, height=640, p...
Note, if no filename is here, it will print to the defaul...
> dev.off()
Some print devices can accept multiple plots, and dev.off...
The postscript device has options to output multiple grap...
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 us...
So what is the alternative?
If you already have a graph on the screen, and did not pr...
dev.copy() opens the device you specify. It is as if (but...
dev.copy(device=png, file="foo", width=500, height=300)
dev.off()
The dev.copy() function takes a device type as the first ...
If you forget the last line, it won't work. It's like a w...
If you use dev.copy or dev.print, you may run into the pr...
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 ...
x11(height=6, width=6)
There are some specialized functions that can do this in ...
The dev.print() function, as far as I can see, is basical...
dev.print(device=postscript,file="yourFileName.ps")
Before you run the graphs, you can do use the postscript(...
If you want pdf output, as of R 1.3 there is a pdf device...
If you just need to test your computer setup, Bill Simpso...
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, poin...
The format of png use is the same.
As of R1.1, the dev.print() and dev.copy2eps() will work ...
> ps <- function(file="Rplot.eps", width=7, height=7, .....
dev.copy2eps(file=file, width=width, height=heigh...
}
> 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 user...
The mismatch of "size" between devices even comes up when...
> dev.print(height=6, width=6, horizontal=FALSE)
You might want to include pointsize=20 or whatever so the...
One user observed, "Unfortunately this will also make the...
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-pla...
The rescaling of a windows() device works differently, an...
**5.3 プロット出力を別個のファイルに保存するには [#z501e9...
(10/04/2001 Paul Johnson <pauljohn@ukans.edu>)
postscript(file="test3.%d.ps",onefile=FALSE)
ファイル名を test3.1.ps, test3.2.ps などと付ける。(from T...
//gives files called test3.1.ps, test3.2.ps and so on. (f...
**5.4 用紙サイズをコントロールする [#o2fe077f]
(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 [#p68d2627]
(20/06/2001 Paul Johnson <pauljohn@ukans.edu>)
Advice seems to be to output R graphs in a scalable vecto...
Keith Reckdahl's Using EPS in LaTeX documents is freely a...
For a description of EPSF and various other vector and bi...
I personally swear by the LaTeX Companion and the LaTeX G...
Thanks to Prof. Ripley for the hint, maybe someone other ...
-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...
?psfrag{wrho}{$W(?varrho)$} <-- this replaces...
?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 c...
**5.6 "Snapshot" グラフとそれらの閲覧 [#ma9b3933]
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
Ordinarily, the graphics device throws away the old graph...
Have a look to ?recordPlot to see how to keep snapshots. ...
On Windows you can choose "history"-"recording" in the me...
**5.7 scatterplot の行列の上三角部分を選択する [#s7b37d6b]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
必要なのはプロットする点を指定するだけ:
//You only need to specify the points to be plotted:
## 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 i...
x <- col(mmtop94.2)[temp] # col indices
scatterplot3d(x, y, z, type="h")
(from Uwe Ligges)
**5.8 密度関数をプロットする (例えば正規分布) [#bff44199]
(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 エラーバー付きでプロット [#ufb394cd]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
以下にやり方を示す。
//Here is how to do it.
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)
//or
または
segments(x,ucl,x,lcl)
(from Bill Simpson)
**5.10 推定密度関数付きのヒストグラム [#n8f63224]
(14/08/2000 Paul E Johnson <pauljohn@ukans.edu>)
同じスケールで密度推定とヒストグラムが欲しければ、以下の...
//If you want a density estimate and a histogram on the s...
> IQR <- diff(summary(data)[c(5,2)])
> dest <- density(data, width = 2*IQR) # or some smalle...
> hist(data, xlim = range(dest$x), xlab = "x", ylab = "d...
probability = TRUE) # <<<--- this is the vit...
> 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 いくつもの線プロットを重ね描きするには？ [#rc63ca52]
(15/08/2000 Paul Johnson <pauljohn@ukans.edu>)
//Use par("new"=TRUE) to stop the previous output from be...
前の出力がかき消されないようにするためには，以下のように ...
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="")
//Note the par() is used to overlay "high" level plotting...
par() は，高水準のプロット命令を重ねるものであるが，ablin...
Here's another idea: form a new data.frame and pass it th...
plot(cbind(ts1.ts,ts2.ts),xlab="time", ylab="", plot.typ...
//or better something like:
以下のようにすれば，もう少しまし。
plot(cbind(ts1.ts,ts2.ts),plot.type="single",
col=c("yellow","red"),lty=c("solid","dotted")) #co...
//It can also be helpful to contrast 'c(ts1.ts,ts2.ts)' w...
'c(ts1.ts,ts2.ts)' と 'cbind(ts1.ts,ts2.ts)' を比較すると...
(from Guido Masarotto)
//For time series data, there is a special function for t...
時系列データについては特別な関数がある。
ts.plot(ts1.ts, ts2.ts) # same as above
//See help on plot.ts and ts.plot (from Brian D. Ripley)
plot.ts と ts.plot のヘルプを見よ。 (from Brian D. Ripley)
//Here's an example for histograms:
以下は，ヒストグラムの例である。
> hist(rnorm(100, mean = 20, sd =12), xlim=range(0,100),...
> par(new = TRUE)
> hist(rnorm(100, mean = 88, sd = 2), xlim=range(0,100),...
//The label at the bottom is all messed up. You can put t...
横軸ラベルがめちゃくちゃになる。これを避けるには，xlab=""...
（訳注：タイトルもめちゃくちゃになるので，main="" も加え...
//Here's another example that works great for scatterplots
散布図を書くときにうまい方法がある。
（蛇足の訳注：xlim, ylim で xlim=c(min(x1,x2),max(x1,x2))...
> 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 グラフを行列状に配置する [#f3268d7d]
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
//The standard way to array several images across a page ...
ページに複数のグラフを配置する標準的な方法は，以下のよう...
par(mfrow=c(2,3))
//And then the next 6 plot/image commands will be laid ou...
そうすれば，次の 6 つのグラフは2行3列に配置される
The layout option can give a powerful set of controls for...
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(...
x2=as.factor(round(runif(100,1,4))) ) attach(x...
by(y,list(x1,x2),function(x)print(x)) xn <- as.integer(x...
## 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(x...
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 r...
And to label the whole page, use the "mtext" function
**5.13 線と棒グラフを結合する [#d1da03b8]
(07/12/2000 Paul Johnson <pauljohn@ukans.edu>)
David James was kind enough to help me out and enlighten ...
xbarplot <- function(y, col=par("col"), border = par("fg...
{
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 回帰散布図: グラフに当てはめ回帰直線を加える [#v5c...
(17/08/2001 Paul Johnson <pauljohn@ukans.edu>)
suppose you have
reg <-lm( y ~ x)
plot(x,y)
then do
abline(reg)
**5.15 散布図中のプロット文字を制御する [#q6e48599]
(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...
If you specify pch, only the first character is taken as ...
Matt Wiener suggested creating a color vector and then us...
> 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 は，プロット記号を選択するために整数を用いることを覚...
//Remember pch uses integers to select plotting characters.
x <- 1:10
y <- 1:10
res <- -4:5
plot(x, y, pch = ifelse(res < 0, 20, 1))
もし論理式が真であれば，記号20が，さもなくば記号1が使われ...
//If true, use character 20. Else use character 1. Experi...
Then note you can create a vector of integers to control ...
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 ...
**5.17 サイズと色を修正した散布図 [#d297ed49]
(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="Siz...
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 ...
(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 散布図：三番目の変数に関してサイズを調整する [#h03...
(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 散布図: 点を結ぶ線を滑らかにする [#ua9136e8]
(02/06/2003 Paul Johnson <pauljohn@ku.edu>)
//Rado Bonk asked how to smooth a line connecting some po...
Rado Bonk がいくつかの点をなめらかに結ぶ曲線を求めている...
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 回帰散布図: 推定値をプロットに加える [#a01ad6ac]
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
I don't know how to get the estimate of the line printed,...
I used this to indicate the R-squared value on a regressi...
text(5, 0.6, as.expression(substitute(R^2 == r, list(r=r...
where R2.the was computed a few lines above. (from Emmanu...
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>)
チックマークは表示するが、チックマークの数値ラベルは表示...
//Ticks, but no numbers:
plot(1:10, xaxt="n")
axis(1, labels=FALSE)
X軸およびY軸のタイトルを省くが、各チックの数値ラベルは残...
//This leaves off the x and y titles. You still have the ...
plot(x,y,ann=FALSE)
もっと多くのグラフの制御をしたいなら、以下のようなことが...
//If you want a lot of control you can do
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)
//do ?axis
? axis と入力してみよ (from Bill Simpson)
Control range of axes through the graph command itself. T...
> 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 軸: ラベルの回転 [#d9e30912]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
一度で、par(las=) は見かけ上これをする。入力は以下の通り。
//At one time, apparently par(las=) would do this, as in
par(las=2) or par(las=3)
上では軸のラベルを90度回転する。
//to make the 90 degrees to axes labels.
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 ...
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?nb...
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 軸: 軸に整形した日付を加える [#mf83c8fb]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
CRANで入手できる"date" パッケージを試してみる。
//Try package "date", available at CRAN.
作図例は以下のとおり。
//Example for plotting:
library(date)
plot.date(.....)
## if x-values are dates, you can also use plot(...) after
library(date)
(from Uwe Ligges)
**5.24 軸: プロット中の軸を逆転する [#w0080bf7]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
R1.1.1 にはバグがあるが、後のバージョンでは、Ross Ihaka ...
//In R1.1.1, there is a bug, but in later versions Ross I...
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 so...
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 ... (althoug...
Another person said try (if you don't want axes):
plot(x, -y, axes=FALSE)
axis(1)
axis(2)
**5.25 軸: 日付によるラベル付き軸 [#e3bab1b9]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
以下のようにデータセットをプロットする必要がある。
//I have to plot the data set like the following.
> x <- c("05/27/00","06/03/00","08/22/00",...)
> y <- c(10,20,15,...)
"date"パッケージを試すべきだ。このパッケージは CRAN で入...
//You should try the package "date", available at CRAN. E...
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 を利用したものだ。
//Here is one using 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日、す...
//which will work if you were in the USA and your first s...
**5.26 軸: 軸ラベルの上付き添字 [#c283699d]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
plot(c(1:10), ylab= expression(x^"++")) (from Peter Dalg...
**5.27 軸: 位置の補正 [#b93bfb1b]
(31/12/2001 Paul Johnson <pauljohn@ukans.edu>)
質問：プロットするとx軸はグラフの下に，y軸はグラフの左に...
//Question: When I plot them default x axis takes place o...
以下のようにすればよいでしょう
//Maybe this helps:
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 散布図に "エラー矢印" を加える [#x4c9351f]
(30/01/2001 Paul Johnson <pauljohn@ukans.edu>)
arrows(..., code=3, angle=90, ...) を使えば、まったく簡単...
//I think using arrows(..., code=3, angle=90, ...) is qui...
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, le...
arrows(x-se.x, y, x+se.x, y, code=3, angle=90, le...
最初の arrows() は y のエラーバーで、2番目のは x のエラ...
//The first arrows() draws the error bars for y, and the ...
Emmanuel Paradis
**5.29 時系列: 一つのグラフに複数の線をプロットする [#d80...
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
同じグラフで２つを描画するのは以下のようにする。
//To do the two in the same graph:
plot(x)
lines(y)
または
//or
points(y)
または
//or
matplot(cbind(x,y),type="l")
グラフを上下に分割して描画するには以下のようにすればよい。
//To do separate graphs one above the other:
par(mfrow=c(2,1))
plot(x)
plot(y)
色の設定、範囲、軸のラベルなどもできる。
//You can do other things like set colors, ranges, axis l...
(from Jon Baron 12/29/2001)
**5.30 時系列: 当てはめ値と実際のデータをプロット [#hc154...
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
> library(ts)
> data(LakeHuron)
> md <- arima0(LakeHuron, order=c(2,0,0), xreg=1:98) pl...
> lines(LakeHuron-md$resid,col="red")
(from Adrian Trapletti)
**5.31 プロットにテキストを挿入 [#x4c67d7b]
(22/11/2000 Paul Johnson <pauljohn@ukans.edu>)
次の文字列をプロットに入れたい。"The average temperature ...
//I want this in a plot: "The average temperature is 23.1...
text(x,y,paste("The average temperature is ", variable, ...
(from Thomas Lumley)
**5.32 はみ出る変数をプロット [#r765e81f]
(07/12/2000 Paul Johnson <pauljohn@ukans.edu>)
Graphs blow up if a variable is unbounded, so if you are ...
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 d...
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)...
dotplot(1:2, labels = my.names)
**5.34 プロットのタイトルに数式等の複雑なものを使う [#l2f...
(11/11/2002 Paul Johnson <pauljohn@ku.edu>)
substitute() を使う。なにかこのようになる。
//Use substitute(). Something like
e<-substitute(expression(paste("Power plot of ", x ^ alp...
alpha == ch.a)), list(ch.a=formatC(alpha,wid=1)))
plot(x, x^alpha, main=e)
(from Peter Dalgaard)
以下のようにしているとすると
//Suppose you have
> that <- 1
> plot(1:10)
そして、タイトルに "that" を指定する必要がある。
//and you need to specify "that" in the title.
expression(paste(hat(theta),'= ',that)))
Generally, to get an expression either use parse on your ...
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 45...
lat<-45
plot(1,main=substitute("Data at "*lat*degree*" North",li...
(from Thomas Lumley)
You can use different way: First export plot for R as an ...
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 P...
plot(x[!is.na(z)],y[!is.na(z)],xlim=range(x),ylim=range(...
points(x[is.na(z)],y[is.na(z)],col=2) (from Ross Darnell)
**5.36 lattice: その他の例 [#bc8fdfd6]
//misc examples
(12/11/2002 Paul Johnson <pauljohn@ku.edu>)
The lattice library came along recently, I've not explore...
library(lattice)
# 10 rows
x <- matrix(rnorm(100), ncol=10)
lp1 <- levelplot(x, colorkey=list(space="bottom"), aspe...
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...
print(lp2)
This one from Deepayan Sarkar (deepayan@stat.wisc.edu) is...
library(lattice)
xyplot(rnorm(100) ~ rnorm(100) | gl(2, 50),
strip = function(factor.levels, ...) {
strip.default(factor.levels = expression(1 * degree, 2 *...
...)
})
**5.37 三次元の散布図 [#x867c42a]
(11/08/2000 Paul Johnson <pauljohn@ukans.edu>)
以下をタイプすると
//Type
?persp
?contour
ベースパッケージにあるこれらの特徴を参照できる。
//to see about features in the base package.
CRAN には、scatterplot3d パッケージがある。
//There is a scatterplot3d package on CRAN.
以下も参照のこと：
//See also:
http://www.statistik.uni-dortmund.de/leute/ligges.htm (f...
以下もまた参照のこと：
//See also:
Get Xgobi (an X windows graphics library) and the R contr...
**5.38 値を反映した線種を用いた三次元等高線図 [#gd556439]
(06/04/2001 Paul Johnson <pauljohn@ukans.edu>)
正の値が実線で、負の値が点線のコンターを描画したい場合,
//Suppose you want a contour plot solid contour lines for...
種明かしをしてみると levels で間隔を指定して、 add= 引数...
//The trick is to specify the levels and to use the 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), lt...
contour(z, levels = 0, lty = "dashed", add = TRUE)
contour(z, levels = seq(0.1, 1, by = 0.1), add = ...
(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...
このようにしておけば，例えば色も変えられる
contour(z, levels = k, lty = c("dotted", "dashed"...
**5.39 グラフのアニメーション [#d937e3ff]
(13/08/2000 Paul Johnson <pauljohn@ukans.edu>)
ImageMagick (www.imagemagick.org)を使えば本当に容易にでき...
//It is quite easy to do with ImageMagick (www.imagemagic...
Rで、１５個のヒストグラムを作成する。
//In R, create 15 histograms:
> for(i in 1:5) {
+ jpeg(paste("fig", i, ".jpg", sep = ""))
+ hist(rnorm(100))
+ dev.off()
+ }
それから、コマンドライン（わたしはLinuxで試したが、他のプ...
//Then from the command line (I tried it using Linux, tho...
% convert -delay 50 -loop 50 fig*.jpg animated.gif
これは animated.gif を作成した。これは私の一連のファイル...
//This created animated.gif, a nice animation of my seque...
**5.40 グラフのユーザー領域の背景色を余白と異なったものに...
(06/09/2000 Paul Johnson <pauljohn@ukans.edu>)
"par" を使って、 boundaries/axes の描画位置を見つけ出し、...
//Use "par" to find out where the boundaries/axes are, an...
//For example:
例えば，
## 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(co...
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 実際に動くグラフコード幾つか (その他) [#x24d401b]
(18/06/2001 Paul Johnson <pauljohn@ukans.edu>)
*** 点プロット、直線上書き、凡例、PS 出力 [#ldad56be]
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 関数の使用例 [#h...
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 =...
xlab = "x", ylab = "sin(x)", main = "Y = sin(x)...
axis( 2, at = seq( -1, 1, by=1 ),las = 2 ) box(lty="dott...
plot( cos(x), type = "l", xaxs = "i", yaxs = "i", axes =...
xlab = "x", ylab = "cos(x)", main = "Y = cos(x)...
axis( 2, at = seq( -1, 1, by=1 ),las = 2 )
box(lty="dotted")
abline( h = 0, lwd = 1 )
(from Ko-Kang Wang)
***正多角形のプロット [#hb2462cc]
Below is a function that generates regular polygons, fill...
borders, of n sides (n>8 => circle), with "diameter" pres...
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/...
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)
***自前の軸 [#he47dc08]
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)
***プロットの中心に軸を持つ三次関数を描く [#o7036612]
# Draw the basic curve (limits were established by tri...
# The axes are turned off blank axis labels used.
curve(x * (2 * x - 3) * (2 * x + 3), from = -1.85, to =...
xlim = c(-2, 2), ylim = c(-10, 10), axes = FALSE, xla...
# Draw the axes at the specified positions (crossing at...
# The ticks positions are specified (those at (0,0) are...
# Note the use of las=1 to produce horizontal tick labe...
# 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...
title(main = expression(italic(y==x * (2 * x - 3) * (2 ...
(from Ross Ihaka)
***3次元曲面を等高線プロットの様に塗分ける、例えば大きな ...
# 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 fa...
hgt <- 0.25 * (z[-nx,-ny] + z[-1,-ny] + z[-nx,-1] + z[-1...
hgt <- (hgt - min(hgt))/ (max(hgt) - min(hgt))
*** 各面の色を指定した曲面を描く [#u6306b4c]
persp(x, y, z, col = gray(1 - hgt), theta = 35)
persp(x, y, z, col = cm.colors(10)[floor(9 * hgt + 1)], ...
(from Ross Ihaka)
*** polygon 関数使用例 [#lcf2ddcf]
y <- c(1.84147098480790, 1.90929742682568, 1.14112000805...
0.243197504692072, -0.958924274663133,0.7205845...
1.65698659871879, 1.98935824662338, 1.4121184852...
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 による並置出力 [#oe4ef49c]
(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]] へ移動
ページ名:
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
