COLOR(green){SIZE(20){初心者のための R および RjpWiki に関する質問コーナー}}~
~
COLOR(red){SIZE(18){注意：このコーナーは重量オーバーで表示に時間がかかり過ぎるようになりました。今後は過去記事へのコメントを除き、新設の [[Ｑ＆Ａ (初級者コース)]] コーナーを利用下さい。}}
~
#contents
**Callの中身の値を変数に取得したいです。 [#w2d0b4d3]
>[[sally]] (2005-04-26 (火) 11:57:46)~
~
こんにちは。国際化おめでとうございます。~
今、グラフプロットに直線回帰をしています。~
そこででてくるパラメータの中身をグラフに書く為に、値を取得したいと思っています。
plot(logCy5, logCy3)
reg <-lm( logCy3 ~ logCy5)
abline(reg)
ここで、Rのコンソールで、reg リターンとすると、以下の値がコンソール中に出てきます。
Call:
lm(formula = logCy3 ~ logCy5)
Coefficients:
(Intercept) logCy5
-0.01426 0.22954
これを、どうにかして変数に取得し、グラフ上に書きたいのです。
a <- XXX(reg) -変数に取得？？
text(1000000, 2, labels = "a", pos=2, cex=1)
色々試してみましたが、、やり方がわかりません。~
ご教授の程よろしくお願いします。~
//
-いろいろやり方はありますが，例えばtext(max(logCy5),min(logCy3),sprintf("logCy3=%6.5f+%6.5f logCy5",coef(reg)[1],coef(reg)[2]),pos=2)とすれば書けます。 -- &new{2005-04-26 (火) 12:33:54};
-ありがとうございました！かけました。 -- [[sally]] &new{2005-04-26 (火) 13:04:09};
**カテゴリー変数の入った非線形回帰はどうすればいいでしょうか。 [#o634ff56]
>[[伊藤]] (2005-04-17 (日) 15:59:59)~
~
お願いします。~
~
私は生物学の研究をしており、２種類の寄主植物上が発育期間とその後の産卵数との関係に及ぼす影響について解析しています。具体的には、横軸に発育日数をとり縦軸に産卵数をとります。この関係をそれぞれの寄主植物別で解析すると、両方とも二次関数へのあてはまりがよさそうでした。そこで次に、寄主植物の種類が二次関数のパラメータに影響しているかを評価するために、以下のモデルをたてました。~
egg=a*(development)^2+b*(development)+c*host+d~
hostの効果を入れたモデルとないモデルの間でAICを比較しようと思いました。~
~
前置きが長くなりましたが本題です。非線形のあてはめを行うときにはnlsを使うことを知っていたので、上の式を~
nls(egg~a*(development)^2+b*(development)+c*host+d,start=list(a=2,b=1,c=1,d=1),data)~
として解析しようとしたら、
Error in numericDeriv(form[[3]], names(ind), env) :
Missing value or an Infinity produced when evaluating the model
というエラーが出て止まってしまいました。"host"（寄主の種類；factor)を抜くとうまくいくので、nlsではカテゴリー変数を評価できないのかと思います。~
~
nlmeも試みてみましたが、ランダム因子がないので不可なようです。代案としてどういうコマンドが考えられるか、ご教唆いただければと思います。よろしくお願いします。~
//
-詳しいことは文面だけからはわかりませんが、エラーメッセージは目的関数（非線形最小自乗誤差）がパラメータに関して微分できないと言っているように見えます。ということは hosts が因子ということとは無関係なのでは？確認のために、optim 等で直接自乗誤差を最小化してみるともう少し様子がはっきりするかも知れません。 -- [[MKR]] &new{2005-04-17 (日) 18:08:01};
-hosts 変数は本当に因子（内部表現で整数）になっていますか？ 文字のままということはないですね？ -- [[MKR]] &new{2005-04-17 (日) 18:15:23};
-データを示すのは差し障りがあるでしょうが，一部分とか，本当のデータを操作したものとかを提示することはできるのではないかと思います。そのようなことも不可能ということでしたら，そもそもあなたの革新的アイデアをこのような場に露出すること自体不適切でしょう。相談するからには，要するに，ほかの質問の場合も同じですが，「データも含めた実行可能なプログラムの提示」が必須です。ただ，「現象を再現できる仮想データとプログラムでもよい」ということはあるのです。要するに，他の人が実際に問題を再現できて対処法がこれだというのを確定できるだけの実際のデータが必要なんです。その基準に満たない質問は，解答を切実に要求していないと見なされても仕方ないでしょう。最近，この掲示板では，質問しっぱなし，解答もらいっぱなしという傾向が見られますが，そのようなことは潜在的回答者のやる気をそぐことになると思います。ま，やる気をそがれない回答者もいるでしょうし，あなたは一回限りの質問者かもしれないのだが。 -- [[青木繁伸]] &new{2005-04-17 (日) 18:48:37};
-追加：「寄生植物の種類が二次関数のパラメータに依存」ということは，そもそも変だと思いますが？たとえて言えば，「4種の血液型の二乗が何とかを予測する」というようなことを言うのと同じようなことではないのですか。そうでないとしたら，説明が不十分だと思います。「寄生植物の種類が二次関数のパラメータに影響しているか」ということですが，寄生植物の種類というのがどういう風に表されていて，その二乗というのがどういうことなのか。まさか，寄生植物A,B,Cを1,2,3で表していて，その二乗，1,4,9がモデルに組み込まれると言うことですか。そうではないのでしょう（と，希望的観測）。ついでに言えば，記事の投稿ルールも把握されていませんね。直しておきます。どのように直されたか，検討されると幸せでしょう。 -- [[青木繁伸]] &new{2005-04-17 (日) 19:03:20};
-読んでも読んでもよく分かりません。「寄主植物の種類が二次関数のパラメータに影響」って，host は一次しか入ってないし，factor って書いてあるし。？？初期値が悪いだけなんじゃない？ -- [[青木繁伸]] &new{2005-04-17 (日) 19:21:20};
-エラーメッセージから見れば，host に NA とか Inf とか，モデル式を評価するときに NA や Inf を発生したと言うことなんでしょ？元のデータに NA や Inf は，ないのでしょうね。あ〜，まだるっこしいですね。-- [[青木繁伸]] &new{2005-04-17 (日) 19:23:13};
-こういう風にやればできるんじゃないのという，見本を作りました。host は factor で，0,1,2という値をとると仮定しています。as.factor をそのままnls では使えません。
set.seed(777)
development <- rnorm(10)
host <- sample(0:2, 10, replace=TRUE)
host1 <- as.integer(host==1)
host2 <- as.integer(host==2)
egg <- development^2+2*development+3*(host==1)+4*(host==2)+rnorm(10,sd=0.01)
data <- data.frame(development=development, host=host, egg=egg)
result <- nls(egg~a*(development)^2+b*(development)+
c1*host1+c2*host2+d,data,start=list(a=1,b=1,c1=1,c2=1,d=1))
plot(egg, result$m$fitted())
summary(result)
# summary の結果
Formula: egg ~ a * (development)^2 + b * (development) + c1 * host1 +
c2 * host2 + d
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 1.016296 0.010428 97.462 2.16e-09 ***
b 1.975588 0.010862 181.876 9.53e-11 ***
c1 3.006908 0.011833 254.116 1.79e-11 ***
c2 4.001636 0.011679 342.622 4.02e-12 ***
d -0.002181 0.011612 -0.188 0.858
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 0.009724 on 5 degrees of freedom
Correlation of Parameter Estimates:
a b c1 c2
b + 1
c1 , . 1
c2 . . , 1
d , . * +
attr(,"legend")
[1] 0 ` ' 0.3 `.' 0.6 `,' 0.8 `+' 0.9 `*' 0.95 `B' 1
#ref(nls.png)
ここしばらくは，細部にわたった解答をして，質問者をスポイルすることもいとわずと言う解答スタイルを取っています。 -- [[青木繁伸]] &new{2005-04-17 (日) 20:01:19};
- 不快な思いをさせて大変申し訳ありませんでした。
お聞きしたかったことは、線形回帰における共分散分析のように、寄主植物間で切片が異なるか、といった分析をしたい、という意味でした。表現が正確でなかったことをおわびします。
~データフレームは以下のような構造です。developmentは発育期間、eggは産卵数、hostは寄主の種類です。
development egg host
17 73 A
18 33 A
19 41 A
20 28 A
17 23 A
17 11 A
18 51 A
19 56 A
19 10 A
18 25 B
20 25 B
19 18 B
18 9 B
19 25 B
21 1 B
19 21 B
行ったのは次のようなことです。
y<-nls(egg~a*development^2+b*development+c*host+d,
start=list(a=2,b=1,c=1,d=1),data)
皆さんがご指摘のとおり、hostのA, Bをそのままnlsに入れたのがいけなかったようです。青木先生のコードを見て、hostA, hostBというダミー変数を作る必要があると理解しました。ありがとうございました。 -- [[伊藤]] &new{2005-04-17 (日) 21:27:00};
-fec というのは egg?? カテゴリー変数をダミー変数に変換するとき，必要なダミー変数の個数はカテゴリー数マイナス1です。従って，host が A,B の二つならダミー変数は1つで良いです。ダミー変数一つで，A なら 1, Bなら0 （逆でも良いけど）とすればよいわけです。これ，重要。試験に出ます(^_^) -- [[青木繁伸]] &new{2005-04-17 (日) 22:05:01};
-すみません。直しました。 -- [[伊藤]] &new{2005-04-17 (日) 22:12:32};
-変数をくくるかっこも不要。初期値 e もformulaにない。直しておきます。 -- [[青木繁伸]] &new{2005-04-17 (日) 22:24:40};
-nls 関数は "A" を実数倍することを迫られて、困ったと見えますね。数値かどうかのチェックがあっても良いような気もしますが、さすがにそれは想定外だった？ -- &new{2005-04-18 (月) 14:11:26};
**データのソートについて （なんでも掲示板から移動） [#n1ab57b9]
>[[thorman]] (2005-04-16 (土) 21:51:53)~
~
最近Rで統計学を始めたのですが（以前はExcelを使っていました）、Rでデータのソートってできるのでしょうか。~
例えば、x<-c(1,3,4,4,6,5,2)とxにデータを入れた後、このデータを降順（あるいは昇順）にならべたりするのはできるのでしょうか。~
まだ、Rについてよく分からないのですが、C言語でデータをソートするようなプログラムを作るように、じぶんでfor構文でソートするためのプログラムを作らなければダメなのでしょうか？~
//
-投稿場所が違っていますよ。必要な関数はすべて揃っています。？sort, ?rev, ?rank, ?order 等で関連関数の使い方を調べてみて下さい。 -- [[QDU]] &new{2005-04-16 (土) 22:00:37};
-探すのが面倒なときは，プログラムを書けばいいのです。あなたの技量によりどっちがより早く解答にだどりつけるかちがうでしょうけど。以下は，プログラムの一例。Excel もソートは作りつけの機能が有るのでわざわざVBAでプログラムすることもないわけですが，場合のよっては作りつけのソート機能が性能の悪いものだったり，標準的なソート機能では要求を満たせないとかの場合には，遅くても仕方ないから自分でプログラムしなきゃならんことも有るわけです。
x<-c(1,3,4,4,6,5,2)
sort <- function(x)
{
n <- length(x)
for (i in 1:(n-1)) {
min <- x[i]
pos <- i
for (j in (i+1):n) {
if (x[j] < min) {
min <- x[j]
pos <- j
}
}
if (i != j) {
x[pos] <- x[i]
x[i] <- min
}
}
return(x)
}
sort(x)
# 結果
[1] 1 2 3 4 4 5 6
このソートアルゴリズムは，優れたものではないですが，わかりやすいもののひとつでしょうね。 -- [[青木繁伸]] &new{2005-04-17 (日) 20:39:56};
**クラスター解析の樹形図の作成について [#rd964b0a]
>[[ごろう]] (2005-04-14 (木) 13:06:21)~
~
クラスター解析で樹形図を作成しているのですが、項目の日本語を入れたまま作図しようとするとエラーになるので、項目を削除して作図をしているのですが、そうすると図の下には数字しか表示されません。R上で図の下に表示される数字と項目を対応させる方法はあるのでしょうか？~
自分の作図方法を下に示します。 よろしくお願いします。~
使用OS Microsoft windows XP~
R Version 2.0.0
library(cluster)
A<-read.table("gorou.txt")
A1<-A[,-1]
A.euclid<-dist(A1, method="euclid")
Euclid.complete<-hclust(A.euclid, method="complete")
plclust(Euclid.complete, hang=-1)
//
-意味が不明です。まず日本語化Rを使って問題が解決しませんか。日本語に拘らないなら（オリジナル版なら拘りようがない）、項目を英語にしておけば済まないのですか。 それと、試行条件を丁寧に書いてくれたのは良いのですが、gorou.txt がなければ、結局私には「あなたが困っている状況」が再現できません。この場で簡単に紹介できるようなミニチュアデータを使うか、R の組み込みデータを使って紹介して頂くと、回答者は幸せになれます。-- &new{2005-04-14 (木) 15:10:23};
-まだはじめたばかりで、質問の仕方がよく分かっていませんでした。申し訳ありません。 アドバイスを参考にして、いくつか自分で試してみたいと思います。 ありがとうございました。 -- [[ごろう]] &new{2005-04-14 (木) 19:03:18};
**plotでのmain titleの指定 [#x3eafa57]
>[[beginner for R]] (2005-04-13 (水) 16:31:53)~
~
いくつのデータをバッチ処理によってプロット出力したいのですが、そのときに入力ファイル名をプロットのmain titleに指定したいのですが、どのようにやればいいでしょうか？ ~
どうぞよろしくお願いします。~
//
-今の状態（つまり main タイトルが書けない点は不満だが，ちゃんと動いているプログラム）を，見せて頂くといいと思いますよ。そうしないと，手間がかかるから，回答してやろうと思う人が少なくなる。~
ファイル名をどのように得ているのかとか，，~
でも，以下のプログラムを見れば検討はつきますか？見当が付かないようでしたら，あきらめてください。
file.name <- c("data1.dat", "data2.dat") # ファイル名のベクトル
par(ask=T)
for (i in 1:2) {
fn <- file.name[i] # ファイル名
dt <- read.delim(fn) # データの入力
plot(dt, main=fn) # main タイトル付きの描画
}
ようするに，ファイル名を文字変数に代入し，read.delim の file パラメータ や plot の main パラメータにはその変数を渡してやれば良いのです。複雑な加工をしてから main に渡してやることもできますよね。 -- [[青木繁伸]] &new{2005-04-13 (水) 16:56:53};
**関数とlistコマンドに関する質問 [#wdd5d9f8]
>[[21℃]] (2005-04-07 (木) 14:56:40)~
~
hoge.test <- function(l="hoge") {
list("fuga"=c("PV.IV","IV"),l=c("Yes","No"))
}
hoge.test("fugafuga")
とすると、lが評価されずに~
$fuga
[1] "PV.IV" "IV"
$l
[1] "Yes" "No"
となります。$lが$fugafugaとなるようにするにはどうすればよいのでしょうか。print(l)を入れるとlに期待通り"fugafuga"が入っているのですが。。。~
//
-やろうとしてらっしゃることの意義が分かりませんが，以下のようにすればあなたの思うようになるというのでしょうか？~
結論は，以下の l=l の部分で，= の左右の l は意味が違うと言うことで。ちゃんちゃん！
> hoge.test <- function(l="hoge") {
+ list("fuga"=c("PV.IV","IV"),l=l)
+ }
> hoge.test("fugafuga")
$fuga
[1] "PV.IV" "IV"
$l
[1] "fugafuga"
ちゃんとしたプログラムの一部としては意味のあるパーツなんでしょうかね？ -- [[青木繁伸]] &new{2005-04-07 (木) 15:04:45};
-もし質問の趣旨が、リスト返り値の成分名を引数で与えたいということなら、例えば -- &new{2005-04-07 (木) 15:32:08};
> hoge2 <- function(name1="a", name2="b")
{ tmp <- list(rnorm(2), runif(2))
names(tmp) <- c(name1, name2)
tmp
}
> hoge2()
$a
[1] -1.627812797 -0.005876096
$b
[1] 0.8406733 0.6035658
> hoge2("A","B")
$A
[1] 1.111073 1.008163
$B
[1] 0.3292604 0.1637242
-ちゃんとしたプログラムのパーツのつもりです。解決しました。ありがとうございました。
> hoge.test <- function(l="hoge") {
+ x <- list("fuga"=c("PV.IV","IV"),l=c("Yes","No"))
+ names(x) <- c("fuga",l)
+ print(x)
+ }
> hoge.test("fugafuga")
$fuga
[1] "PV.IV" "IV"
$fugafuga
[1] "Yes" "No"
期待通りに、$lではなく$fugafugaになりました。c()の中ならlが評価されて、list()の中ならlはlのままというのが釈然としませんが、とにかく解決しました。ありがとうございました。 -- [[21℃]] &new{2005-04-07 (木) 16:04:46};
-$fugafuga になるようにということだったんですね。前のコメントにも書きましたが，l=l というときの左辺の l はリスト要素の名前ですから，names で与えないといけないのですね。 -- [[青木繁伸]] &new{2005-04-07 (木) 16:20:49};
- 成分名の一部をかえるだけで良いなら次のようにします。-- &new{2005-04-07 (木) 17:36:59};
> x <- list(a=runif(2), b=rnorm(2))
> x
$a
[1] 0.4666695 0.3686140
$b
[1] 0.1382112 -0.3268364
> names(x)[2] <- "c"
> x
$a
[1] 0.4666695 0.3686140
$c
[1] 0.1382112 -0.3268364
- list("fuga"=c("PV.IV","IV"),l=c("Yes","No")) というように、左辺を "fuga" のようにダブルクオートでくくって書いてしまうような所に、間違いの根元があったと思いますが。 -- &new{2005-04-07 (木) 22:48:18};
-試してみればわかるように、二重引用符の有無の問題ではないですね。 -- &new{2005-04-08 (金) 06:53:58};
-付いていてもいなくても同じ結果になるが、二重引用符を付けなければならないと思ったこと自体が問題。=の左に書くのは、変数名や、引数名や、文字定数ではないでしょう？すくなくとも、引数名ではない。~
文字定数として"付きで書くと思っていたから、変数名も書けるのではないか、更に、引数名も書けるのではないかと思ったのだろう。 -- &new{2005-04-08 (金) 07:24:06};
-好きなように文字列を生成してパースして実行も可能. eval(parse(text=paste("list(fuga=c(\"PV.IV\",\"IV\"),",l,"=c(\"Yes\",\"No\"))"))) 二重引用符の問題は,"hoge foo"<-1:10 などとしても, `hoge foo`で参照可能(後退引用符?), listの場合は現状単純にstring=stringで左辺には名前を期待しているに過ぎない. もちろん "foo boo"など,間に空白などが入ってもちゃんと処理が行える. 私はパズルは苦手なので, もし l=c("Yes","no") の lの展開によって, パズルが隠匿(かんがえなくてもよい)されるなら,むしろこういった表現の追求は実は良いのではないかと. さらに毎年恒例の4/1発表作品に向けパズル度を高めるのも良し. どんな難しいパズルでもほんの少しのコメントを追加すれば, 後で自分がわからない事は無いでしょう. -- [[なかま]] &new{2005-04-08 (金) 10:54:17};
-みなさんありがとうございました。よく理解でき、すっきりしました。私の愚考の筋道はご指摘の通りです。eval(parse(text=paste()))も実は試みたのですが、書き方が悪かったみたいです。 -- [[21℃]] &new{2005-04-08 (金) 13:54:38};
**sprintf の指数表示について [#g21d28cd]
>[[初心者11242号]] (2005-04-05 (火) 14:48:32)~
~
sprint のフォーマットで e を指定したときに、R2.0.1 では e の後ろが3桁になって表記されますが、これを2桁や4桁にするには どうしたらよいのでしょうか？~
現在は
> sprintf("%8.5e", 0.0001)
[1] "1.00000e-004"
となっているのを "1.00000e-04" や "1.00000e-0004" としたいのです。~
2桁の場合は formatC を使い, 4桁の場合は
sprintf("%4.2e0", 0.0001) # おかしいですね。これでは期待される結果は得られないと思いますが
などと無理やりフォーマットを変えていますが、もっと良い方法があれば教えてください。~
OS は Windows XPで、R2.0.1 日本語版を使っています。~
// ●● 校正しておきますね。。。全く手が掛かるっちゃありゃしない。
//
-ないものは作る。ただ，本気で作るほどの関数でもない。~
もう少しオプションを増やしてもいいかもしれないが，参考になれば。~
ま，言ってみれば，「こそくな関数」です。~
exponent は指数部の桁数（符号も含む）
extend.sprintf <- function(x, exponent=4)
{
a <- unlist(strsplit(sprintf("%10.5e", x),"e"))
fmt <- sprintf("%%se%%0%ii", as.integer(exponent))
sprintf(fmt, a[1], as.integer(a[2]))
}
> x <- 0.00012345
> extend.sprintf(x, exponent=1)
[1] "1.23450e-4"
> extend.sprintf(x, exponent=2)
[1] "1.23450e-4"
> extend.sprintf(x, exponent=3)
[1] "1.23450e-04"
> extend.sprintf(x, exponent=4)
[1] "1.23450e-004"
> extend.sprintf(x, exponent=5)
[1] "1.23450e-0004"
> extend.sprintf(x, exponent=15)
[1] "1.23450e-00000000000004"
> x <- 123456789.0123
> extend.sprintf(x, exponent=4)
[1] "1.23457e0008"
あまり，有用な関数とも思えません。 -- [[青木繁伸]] &new{2005-04-05 (火) 19:10:27};
-ありがとうございます。私にとっては有用です。
fmt <- sprintf("%%se%%+0%ii", as.integer(exponent))
としたら、希望通りになりました。スペルミス、失礼いたしました。青木先生のサイトは よく 参考にさせていただいています。 -- [[初心者11242号]] &new{2005-04-05 (火) 23:22:06};
-おっしゃるとおり，%+04i で正の場合にも + の符号が付きますね（補足） -- [[青木繁伸]] &new{2005-04-05 (火) 23:41:17};
-おお！きれいに編集されて見やすくなっている... お手数をおかけしてしまい、申し訳ありませんでした。 -- [[初心者11242号]] &new{2005-04-07 (木)};
**しつもん [#q8d604eb]
>[[初心者]] (2005-03-31 (木) 01:48:22)~
~
シュミレーションをしていたら以下のようなエラーが出たのですが、
Error: subscript out of bounds
ってなんでしょうか？~
//★★ 投稿法くらい調べてから投稿したら？タイトルが「しつもん」というレベルでは，所詮無理なんだろうけど。
//
-どういうシュミレーションしたのか，きちんと追試可能なように書いた方が良いですよ．あと，「しつもん」というタイトルは相応しくありません．もう少し，生きるうえで必要なマナーを守りましょう． -- [[がくせー]] &new{2005-03-31 (木) 03:05:21};
- 次の例を参考に考えてみて下さい。-- &new{2005-03-31 (木) 08:00:08};
> y <- matrix(1:4, 2,2)
> y
[,1] [,2]
[1,] 1 3
[2,] 2 4
> y[3,2]
Error: subscript out of bounds
-回答は質問のレベルによるという見本ですね。~
"Error: subscript out of bounds" は，「添え字が範囲を超えました」ということです。~
2005-03-31 (木) 08:00:08 の回答例を解説すると，「行列は2×2の大きさなので一番目の添え字に3を指定してしまうと，3行2列目の要素ってのはありませんよ」ということです。~
このような初心者だと，英語のエラーメッセージが日本語化によって「添え字が範囲を超えました」なんていう日本語メッセージになっても，理解できるのかどうかわからないですね。 -- &new{2005-03-31 (木) 09:10:37};
-もしかして,始めて学ぶコンピュータ言語がRと言う幸せ(不幸?)な人が出てきている?(^o^)/? -- [[なかま]] &new{2005-03-31 (木) 13:20:10};
**関数群を重ねて描画する方法について [#y0182bba]
>[[yuta]] (2005-03-30 (水) 14:56:18)~
~
以下のように複数の回帰直線を得る関数を定義したあとで、
taus <- c(0.99,0.95,0.9,0.85,0.8,0.75)
MyRq <-c(rq(y~x, tau = taus,MyData))
散布図上にその複数の関数を描画しようと思い、以下のようなスクリプトを単純に組みましたが、回帰直線が一つだけしか描画されません。恐らく上書きされてしまっているのだと思うのですが、単純に重ね合わせをするような指定の方法はあるでしょうか。お知恵をお借りしたく思います。宜しくお願いします。
plot(x,y)
new = t
abline(c(LineRq),col="blue")
//★★ 投稿の仕方くらい調べてから投稿しなさいな！
//
-自己レスで申し訳ありません。c(rq ・・・はおかしいですね。これを修正しても状況には変化はありませんでしたが。 -- [[yuta]] &new{2005-03-30 (水) 15:10:35};
-一寸準備が面倒ですが、その分確実で世話がない方法は matplot 関数を使うことです。関数毎の座標範囲の差異を手動で調整するのは結構面倒です。?matplot で微妙な調整法を調べてみて下さい。 -- &new{2005-03-30 (水) 15:44:00};
> x <- 1:100/50*pi
> y <- cbind(sin(x), cos(x)) # 同じx座標に対する二つの関数値を合わせた行列
> matplot(x,y, pch=".") # 二つの関数グラフを同時にかく
-見出しにつられてお手つきコメントをしましたが、散布図に回帰直線を重ね画きしたいという趣旨でしょうか。命令 rq というのが何のことかわかりませんが、例えば lm 関数であれば単に次のようにすれば良いのでは。 -- &new{2005-03-31 (木) 17:06:08};
plot(x,y)
abline(lm(y~x))
-回答ありがとうございます。ｒｑは分位点回帰の関数で、複数の分位点をtausとして代入しているので,MyRqは複数の回帰直線を含むオブジェクトになっているはずなのですが、このスクリプトでは１つしか描画されないのです。うーむ。 -- [[yuta]] &new{2005-04-01 (金) 11:32:21};
-「複数の分位点をtausとして代入しているので,MyRq?は複数の回帰直線を含むオブジェクトになっているはずなので」~
そうはおっしゃいますが，追試もできないのです。あなたが実際にやったことをちゃんと書かないと。~
「c(rq ・・・はおかしいですね。これを修正しても」~
というのも，どこをどのようになおしたのかさえわかりません。~
回答が欲しくないなら，質問する必要もないでしょう。 -- [[質問のマナー]] &new{2005-04-01 (金) 13:35:11};
-申し訳ありません。背景説明を省略しすぎました。問題のスクリプトはquantregパッケージを読み込んだ後で
taus <- c(0.99,0.95,0.9,0.85,0.8,0.75,0.50) #複数の分位点を指定
Rq <-rq(stack.loss~Air.Flow, tau = taus,stackloss)
plot(Air.Flow,stack.loss)
NEW = TRUE
lines(Rq,col ="blue")
とするとできるかと思います。この状態で、Rqは７つの回帰式を持っており、散布図にはそれぞれが重ねて描画されると思っていたのですが、うまくいかないのです。 -- [[yuta]] &new{2005-04-01 (金) 15:04:01};
-「複数の回帰直線を含むオブジェクトになっているはずなのですが」=> 本当ですか？-- &new{2005-04-01 (金) 15:18:46};
-プログラムや出力結果を含む質問を投稿するときは，実際にRに入力したもの・Rが出力したものをコピー・ペーストすべし。~
上のものは全角空白で区切られて文章の中にあったものを，投稿規程にそうように編集したのだけど。~
投稿規程をよく読んでくださいね。~
NEW = TRUE というのは，なんでしょうね。~
で，ま，書かれたようにやってみたんですが，
> library(quantreg)
[1] "quantreg package loaded"
[1] "quantreg now includes the nlrq and nprq packages"
> Rq <-rq(stack.loss~Air.Flow, tau = taus,stackloss)
> plot(Air.Flow,stack.loss)
Error in plot(Air.Flow, stack.loss) : Object "Air.Flow" not found
> NEW = TRUE
> lines(Rq,col ="blue")
Error in plot.xy(xy.coords(x, y), type = type, col = col, lty = lty, ...) :
plot.new has not been called yet
ということになったんですが，どうですか？ -- [[投稿ルール]] &new{2005-04-01 (金) 16:18:57};
-そもそも，Rq の前に attach(stackloss) が必要なのか。それと，最後の lines は abline と書きたかったんとちゃうか？だからちゃんと動くプログラムをペーストすべしということ。ヒマだから付き合っているんだけどね。世の中，暇人ばかりじゃないのよ。 -- &new{2005-04-01 (金) 16:31:39};
-ま，一番最後の abline のつもりの lines を以下のようにすれば，とりあえずは，あなたの目的は達せられるのかもね。
> library(quantreg)
[1] "quantreg package loaded"
[1] "quantreg now includes the nlrq and nprq packages"
> attach(stackloss)
> taus <- c(0.99,0.95,0.9,0.85,0.8,0.75,0.50)
> Rq <-rq(stack.loss~Air.Flow, tau = taus, stackloss)
> plot(Air.Flow,stack.loss)
> for(i in 1:ncol(Rq$coefficients)) {
+ abline(Rq$coefficients[,i],col=i)
+ }
上の解は，エープリル・フールじゃないよ -- [[暇人32号]] &new{2005-04-01 (金) 16:56:28};
-なぜ，yutaさんは，「不完全な質問」をしたのでしょうか（決して責めているわけではないです）．その原因を，正直に書いていただくと，今後の参考になると思います． -- [[aa]] &new{2005-04-01 (金) 17:12:33};
-まず、皆さまのお手を煩わせてしまったことをお詫びします。最初はただ単純に、引数指定などで足りない点があるのかと考えてコードの記述を簡単にしてしまった次第です。いずれにせよ第３者が検証できないのではどうにもならないのですよね。以後気をつけます。なお、暇人32号さんのスクリプトで問題は解決しました。どうもありがとうございます。 -- [[yuta]] &new{2005-04-01 (金) 17:48:05};
~
**トップページから記事一件移動しました [#t4eecb91]
Ｒをインタプリタでなくバッチで実行すると何百倍も早くなったりすんでしょうか？ -- [[豊田]] &new{2005-03-29 (火) 08:21:06};
-質問する場所が違いますよ．バッチとインタプリタの比較は多分行われていないのでは？参照：[[R Site Search:http://finzi.psych.upenn.edu/R/Rhelp02a/archive/28051.html]] ．どういう目的が存じませんが，batch mode じゃなくても，かなり速いですよ．参照：[[Speed comparison of various number crunching packages:http://www.sciviews.org/benchmark/]] -- [[me]] &new{2005-03-29 (火) 13:40:28};
-バッチモードとコンパイルモード（これはまだ無い）を混同していませんか。バッチモードも結局はインタプリターを動かしているだけで、画面表示がない分わずかは早くなる可能性もある？ -- &new{2005-03-29 (火) 14:48:45};
----
**パッケージのＵＰＬＯＡＤでエラーがでます。 [#j06538b4]
>[[超初心者]] (2005-03-28 (月) 13:27:38)~
~
trying URL `http://cran.md.tsukuba.ac.jp/bin/windows/contrib/2.0/PACKAGES'
Error in download.file(url = paste(contriburl, "PACKAGES", sep = "/"), :
cannot open URL `http://cran.md.tsukuba.ac.jp/bin/windows/contrib/2.0/PACKAGES'
In addition: Warning message:
unable to resolve 'cran.md.tsukuba.ac.jp'.
package のupload をしようとすると上記エラーがでます。~
ＩＥ上に直接入力するとうまくいくのに、理由がわかりません。~
~
かなりの初心者ですが、だれかお教え願いますでしょうか？？？？~
環境：win2k~
//
-エラーメッセージも重要ですが、ヒントを得るには、どういう命令を実際に実行したのかも大事ですよ。 -- &new{2005-03-28 (月) 15:22:39};
-upload ではなくて，download ですね。。。 -- &new{2005-03-28 (月) 18:51:49};
-proxy が必要になる環境からなら，[[CRAN国内ミラーの使い方]]を参照してみてください -- [[青木繁伸]] &new{2005-03-28 (月) 18:58:03};
**エラーが取れません。。。 [#u8ec8994]
>[[初心者]] (2005-03-27 (日) 16:16:20)~
~
ファイナンスのシュミレーションをするのにRを使ってやってみたのですが、
Error in if ((800/K[i, j]) * Dil * (S[i, j] - K[i, j]) > C[i, j])
missing value where TRUE/FALSE needed
と言うエラーが出てしまって困っています。~
該当箇所は
if( (800/K[i,j]) * Dil *(S[i,j] - K[i,j]) > C[i,j] ){
C[i,j] <- (800/K[i,j]) * Dil *(S[i,j] - K[i,j])
}
と言うもので、K, Dil, S,C などのオブジェクトは全て前に定義しています。~
この部分のどこかにまずいところがあるでしょうか?~
教えてください、よろしくお願いします。~
~
windowsユーザーです。~
//
-エラーが起きた箇所での (800/K[i,j]) * Dil *(S[i,j] - K[i,j]) と C[i,j] の値を確認するのが先決でしょうね。次の例が参考になるでしょうか。 -- &new{2005-03-27 (日) 16:34:13};
> if(NA > 3) cat("A")
Error in if (NA > 3) cat("A") : missing value where TRUE/FALSE needed
> if(NaN > 3) cat("A")
Error in if (NaN > 3) cat("A") : missing value where TRUE/FALSE needed
-早速の返事ありがとうございます。 -- [[初心者]] &new{2005-03-27 (日) 17:10:56};
-上のご指摘の件ですが、これはつまり if( ) の中身が変な値になったと言うことでしょうか? -- [[初心者]] &new{2005-03-27 (日) 17:12:33};
-数値として解釈できない値がどこかで生じたということでしょうね。 -- &new{2005-03-27 (日) 19:19:43};
**R for Mac mini [#l7a86be1]
>future user of Mac mini (2005-03-16 (水) 19:28:26)~
~
R Mac 版はMac mini (メモリー256MB/HDD 40GB)でも、問題なく動作するのでしょうか？
//
-試して，ぜひその結果を報告してください。というか，まだMac mini を持ってませんか。。。。明日までお待ち頂ければ，お返事できるかもしれませんが。一般的に言えば，「動かないはずないと思いますけど。。。」ということでしょうか。ああ，それと，「問題なく」というのは限りなく主観的な問題ですので，貴方が問題ないとする状況の定義を是非お知らせ下さい。「CUI でしか使えないなんて思ってなかった」なんて言われても困るので(^_^;) -- [[青木繁伸]] &new{2005-03-16 (水) 21:47:51};
-私も知りたいです。どうなったのかなど、是非、教えてください。 -- &new{2005-03-17 (木) 10:40:28};
-Mac mini を持っている人に，R をインストールさせてもらいました。そのマシンでは問題なく動作しました。 -- [[青木繁伸]] &new{2005-03-17 (木) 11:59:31};
-さっそく、ためしてみようっと -- &new{2005-03-17 (木) 11:59:31};
-動作情報ありがとうございます。もしもよければ動作Ｍａｃの環境（メモリーとＨＤＤ）教えていただけませんか？増設しなくても動かせるのか知りたいのですが。 -- [[future_user_of_Mac_mini ]] &new{2005-03-17 (木) 17:11:29};
-特に大量のメモリは不要でしょう。多い方が，いろいろなソフトを同時に立ち上げておけるとかあるでしょうが。そもそも，メモリ増設しないと動かないようなソフトは少ないでしょう。 -- [[青木繁伸]] &new{2005-03-17 (木) 17:29:48};
**線分の色分け [#y6e86c5a]
>[[さち]] (2005-03-16 (水) 16:14:04)~
~
またお世話になります。以前”射影した回帰曲線の色の付け方”でお世話になった、さちです。またご指導をよろしくお願いします。~
前回の質問から時間が経ってしまいましたので、新しく質問をたてることにしました。~
~
青木さんのご指導をもとに、フォーミュラ y= f(x) の逆関数 x=g(y) によって，x1=g(a), x2=g(b) を計算して lines(c(x1, x2), c(0,0), col="color you want")で線分の色分けを行おうとしました。~
フォーミュラ y= f(x) の逆関数 ：
a=50;b=20;c=0.1
g<-function(d) (-(1/c)*log( (a-y)/(b*y) ))
x2<-c(g(y))
このようにして、逆関数 x=g(y) を求めて、下記のようにして、線分を描かせてみましたところ、、、
X<-matrix(0,20:3)
X[,1]<-x<-1:20
X[,2]<-x2
<-matrix(0,20:1)
for(i in 2:20) {
lines(c(X[(i-1),1],X[i,2]), c(Y[i,1],Y[i,1]),
col = ifelse( predict.c>= 15.1, "blue",
ifelse( predict.c>= 10.1, "green", "red")))
}
困った事に、得られた線分が、単色（”赤”）になってしまいました。。。~
どうしたら良いのでしょうか。~
ご指導をお願いします。~
//
-初めから書き直します。~
非線形回帰分析を行って得られた、回帰曲線 y=f(x) の各点からx軸へ垂線を下ろし、その交点を結んで得られる線分を、y の値に従って、色分けすることを考えています。~
非線形回帰分析：
> f <-function(x) 20/(1+50*exp(-0.8*x))
> x<-1:20
> y<-f(jitter(x))
> (result <-nls(y~a/(1+b*exp(-c*x)), start=c(a=50,b=20,c=0.1)) )
Nonlinear regression model
model: y ~ a/(1 + b * exp(-c * x))
data: parent.frame()
a b c
20.039274 51.908061 0.788543
residual sum-of-squares: 0.3331449
> predict.c<-predict(result)
> plot(x,y,ann=F,xlim=c(0,20),ylim=c(0,20))
> par(new=T)
> plot(predict.c,type="l",xlim=c(0,20),ylim=c(0,20))
ここまでは、前回と同じで、非線形回帰分析をさせた結果を表示させている部分です。~
逆関数 x=g(y) を求めるために、パラメータ、a,b,c を f(x) と同じように設定しました。
> a=50;b=20;c=0.1
逆関数 g(y)を関数を以下のように求めました。
> g<-function(y) (-(1/c)*log( (a-y)/(b*y) ))
グラフを表示させる際、もっと複雑な場合に備えて、 for を使った繰り返し文を使いたかったので、下記のようにしました。
> x2<-c(g(y))
> X<-matrix(0,20:3)
> X[,1]<-x
> X[,2]<-x2
> Y<-matrix(0,20:1)
先に求めた predict.c に従って線分に色をつけようと考えて、下記のようにして、描かせました。
> for(i in 2:20) {lines(c(X[(i-1),1],X[i,2]), c(Y[i,1],Y[i,1]),
col = ifelse( predict.c>= 15.1, "blue",
ifelse( predict.c>= 10.1, "green", "red")))}
色分けさせるように記述したつもりでしたが、結果は、線分が単色（赤）で表示されてしまいました。~
どのようにすれば、問題が解決できるようになるでしょうか。ご教授をお願いします。 -- [[さち]] &new{2005-03-16 (水) 16:50:15};
-X <-matrix(0,20:3)は X<-matrix(0,20,3)だろうとか。それはおいておきましょう。 -- [[青木繁伸]] &new{2005-03-16 (水) 17:29:15};
-どうにかよめるようになりました。。。お騒がせを致しました。<(__)> -- [[さち]] &new{2005-03-16 (水) 17:30:24};
-最後の for ループの lines 文の col パラメータがちゃんと機能していないのです。たとえば以下のようにすれば，色分けはできます。ただ，変な色分けになっているのかもしれませんね。
col <- ifelse( predict.c>= 15.1, "blue", ifelse( predict.c>= 10.1, "green", "red"))
for(i in 2:20) {lines(c(X[(i-1),1],X[i,2]), c(Y[i,1],Y[i,1]),col = col[i])}
同時刻に修正しているので，編集の衝突が起きていますね(^_^;)。 -- [[青木繁伸]] &new{2005-03-16 (水) 17:31:36};
-青木さん！！有り難うございます！！やっと、それらしくなってきました(^-^) また頑張ってみます！！有り難うございました。<(__)> -- [[さち]] &new{2005-03-16 (水) 17:45:00};
-a=20.039274 b=51.908061 c=0.788543 でないと変？ -- &new{2005-03-16 (水) 17:47:49};
-重複して見苦しいところを削除しました。<(_ _)> -- [[さち]] &new{2005-03-16 (水) 17:48:38};
-もしや。。。と思って、ご指摘のパラメータで確かめてみました。
> a=20.039274; b=51.908061 ; c=0.788543
> g<-function(y) (-(1/c)*log( (a-y)/(b*y) ))
> x2<-c(g(y))
> x2<-c(g(y))
> X<-matrix(0,20,3)
> X[,1]<-x
> X[,2]<-x2
> Y<-matrix(0,20,1)
> col <- ifelse( predict.c>= 15.1, "blue", ifelse( predict.c>= 10.1, "green", "red"))
> for(i in 2:20) {lines(c(X[(i-1),1],X[i,2]), c(Y[i,1],Y[i,1]),col = col[i])}
確かにその通りでした。 result の結果を使えば良かったのですね。~
有り難うございました。-- [[さち]] &new{2005-03-16 (水) 17:53:28};
-念のため。newpar <- result$m$getPar(); a=newpar[1]; b=newpar[2]; c=newpar[3] としないと（少しだけ）損しますよ。 それと、あまりに長い式を一行に書くと、ブラウザーによっては画面が極めて見にくくなりますから、適宜改行して下さい。その方がみやすくもあるし（私が改行しておいたのを、また元に戻されたようなので）。--- &new{2005-03-16 (水) 19:17:11};
-さからったのではなく，そのほかの所を直しただけでは(^_^;) -- &new{2005-03-16 (水) 23:40:26};
-そうだったんですか？？てっきり、自分だと思って。そのままにしておけばよかった。。。 -- &new{2005-03-17 (木) 10:26:23};
-使い方と newpar の意味が今ひとつよく分かりませんが、損は、少しであっても、嫌なので、
> newpar <- result$m$getPar(); a=newpar[1]; b=newpar[2]; c=newpar[3]
をします。 頑張って勉強します。-- [[さち]] &new{2005-03-17 (木) 10:31:08};
-Q＆Aにヒントをみつけました。 -- [[さち]] &new{2005-03-17 (木) 12:34:48};
-str(result) で nls オブジェクトの構造を見て下さい。リストである result の成分 result$m には非線形回帰結果から必要な情報を後から取り出すためのさまざまな専用関数が（それ自身リストの形で）まとめて放りこんであります。その中で result$m$getPar という関数が当てはめパラメータを取り出すための関数です（名前から大体想像がつきますね、実は私もそうした関数があるはずと見当をつけて見つけました）。少し損すると書いたのは、パラメータは本来15,16桁で求められているはずだからです。実質的な違いが出るとはおもえませんが。 -- &new{2005-03-17 (木) 12:49:57};
-なるほど。。。そういうことだったんですね。確かに、str(result) をしてみると、沢山出てきました。もう少し、よく勉強してみます。 -- [[さち]] &new{2005-03-17 (木) 12:56:21};
**Warning message [#s8b18f27]
>[[初心者]] (2005-03-16 (水) 15:17:13)~
~
Excelの表をread.csvで読み込もうとしたら以下のメッセージが出たのですが、どうしたらよいのでしょう？なお、Rは2.01、Mac OS X環境です。~
incomplete final line found by readTableHeader on `Data.csv'~
//
-Excel で書き出したファイルを，エディタか何かでいじりませんでしたか（特に，ファイルの最後）。~
例えば，ファイルの最後の行の行末が 0x0a でない場合に，そのようなエラーメッセージになるでしょう。
% hexdump -C test.csv # このファイルは，数字6で終わっている。エラーが出る
00000000 31 2c 32 2c 33 0a 34 2c 35 2c 36 |1,2,3.4,5,6|
0000000b
% hexdump -C test2.csv # 普通のファイルは，改行コードで終わっている。エラーはない
00000000 31 2c 32 2c 33 0a 34 2c 35 2c 36 0a |1,2,3.4,5,6.|
0000000c
対処法は簡単。エディタで読み込んで，ファイルの最後の行末で，リターンキーを押して，上書き保存する。か，もう一度 Excel から，書き出しを行う。 -- [[青木繁伸]] &new{2005-03-16 (水) 15:58:45};
-行ったのはExcelに入力してCSVファイルで保存しただけなのですが、いかがでしょう？ -- [[初心者]] &new{2005-03-16 (水) 18:45:18};
-ですから，そう言う場合は，Excel のワークシートファイルとそれを csv ファイルにセーブしたファイルの両方を提示しないと，だれも現象を再現できないので，自分で再現しようとしてみたらうまくいくのだけどとか，こういう場合にはそうなりますねとしか，いえないじゃないですか。もっとも，大きなファイルを提示されてもこまるので，現象を再現できる最小限のファイルを提示してみてください。というかそれ以前に，何回やっても同じ結果（エラーが出る）ということでいいのですね。あるいは，hexdump してみたら確かに0x0a で終わっていないファイルだった（なんで 0x0a で終わらないんだろう。。。Excelが悪い可能性），0x0a で終わっているんだけど？（R が悪い可能性）。何ででしょう，何ででしょうといっているのではことが進まない。エラーの可能性を提示したなら，それについて，確認して，それではないですくらいいったらいかが？ああ，それと，R のバージョンは分かりましたが，Excel のバージョンは書いていませんね。 -- [[青木繁伸]] &new{2005-03-16 (水) 21:57:05};
-青木先生、説明が足りずにご迷惑をおかけしました。実は他の表計算ソフト（NeoOffice）で同じことをやってみたら無事読み込めました。Excel（2004,）の問題なのでしょうか？ -- [[初心者]] &new{2005-03-17 (木) 01:10:49};
-その，エラーを生じた，Excel 2004が生成した CSV ファイルは，確かに 0x0a では終わっていないんですね？こちらでも，明日調べてみます。 -- [[青木繁伸]] &new{2005-03-17 (木) 01:13:56};
-Excel 2004 を持っている人に頼んで送ってもらった CSV ファイルは，確かに 0x0a で終わっては＊い＊ま＊せ＊ん＊ね＊。またしても，黒い魔ソフトのちょんぼでしょう。と思いましたが，2004 でなくても，同じですね。ただ，そのような 0x0a で終わらないファイルであっても，問題なく read.csv で読める場合もあるようですね。
read.csv の引数との関係なんかがあって，一筋縄ではいかないのかややこしそうです。 -- [[青木繁伸]] &new{2005-03-17 (木) 12:02:34};
-[[NeoOffice/J:http://www.planamesa.com/neojava/ja/index.php]] だと，ちゃんと 0x0a で終わるファイルを作ってくれるので，これで完全にマイクロソフトとおさらばできそうです。~
[[NeoOffice/J:http://www.planamesa.com/neojava/ja/index.php]] は，~
> [[NeoOffice/J:http://www.planamesa.com/neojava/ja/index.php]]はMac OS X用の機能豊富なオフィス・ソフト（ワープロ、表計算、プレゼンテーション、作図）です。[[OpenOffice.org:http://www.openoffice.org/]]というオフィス・スイートに基づいたもので、数多くのMac OS Xネイティブの機能を提供し、[[MicrosoftTM Officeなど:http://aoki2.si.gunma-u.ac.jp/Hanasi/excel/index.html]]他の人気のオフィス・ソフトの書類の読み込み、編集、交換ができます。~
GNU一般公衆利用許諾契約書（GPL）の下で無料のオープンソース・ソフトとして公開しています。[[NeoOffice/J:http://www.planamesa.com/neojava/ja/index.php]]は充実した機能と日常的な用途に充分の安定性を備えています。
<というものです。 -- [[青木繁伸]] &new{2005-03-19 (土) 10:05:03};
**軸を指定したい [#pc9533ec]
>[[nama]] (2005-03-16 (水) 12:21:56)~
~
すみません。教えていただきたい事があります。~
散布図Ｘ軸を打点があるところだけにしたいのですがどうすればいい~
のでしょうか？ xaxt="n" で消去することは出来たのですが、その後~
３点のみ表示させたいのです。よろしくお願いいたします。
*
*
*
---+----+----+----+----+----+----
40 60 90
//
-すみません、変な「＊」が入ってしまいました。要は３つの打点があって、40,60,90の -- [[nama]] &new{2005-03-16 (水) 12:24:29};
-本当にもうしわけありません。axisで出来るって教えてもらえました。掲示板汚してすみませんでした -- [[nama]] &new{2005-03-16 (水) 12:30:00};
-ここでは「編集」ボタンを押せば，汚してしまったものを自分で片付けることができますよ． -- [[学生]] &new{2005-03-16 (水) 13:06:38};
-あと、どのように解決したかも書いていただけたら幸いです。namaさんはこんな感じで処理されたのでしょうか？ -- [[舟尾]] &new{2005-03-16 (水) 15:06:38};
x <- c(40, 60, 90)
y <- c(10, 4, 19)
plot(x, y, xaxt="n")
axis(side=1, at=x)
-「散布図Ｘ軸を打点があるところだけにしたい」というのは、「散布図のＸ軸の目盛り数値を、データの打点があるところだけに表示したい」ということだったんですね。あまり省略するとよく分からなくなりますね。 -- &new{2005-03-17 (木) 01:38:40};
-みなさんすみませんでした。掲示板の使い方もよくわからないで投稿してしまって。教えてもらったのは正に舟尾さんのおっしゃるようでした。axisの中でatを使って指定しました。また、わからないことがあったらよろしくお願いしたいと思います。 -- [[nama]] &new{2005-03-18 (金) 18:54:54};
**変化率の計算 [#cf6ca79e]
>[[おさる]] (2005-03-14 (月) 02:25:43)~
~
お世話になります．~
~
変化率の計算方法で，うまいものは無いでしょうか．~
~
変化率2= (x[2]-x[1])/x[1]~
変化率3= (x[3]-x[2])/x[2]~
~
といった感じの計算を，データフレームの変数xの各行に対して（簡単スマートに）行いたいのですが．．．~
例えば，次のような株価のデータがあるとき，
x y z
905 728 589
950 773 597
930 737 580
960 770 606
1080 845 688
1030 808 669
965 733 687
992 759 691
956 745 702
これら企業の株価それぞれについて，上記のような変化率を求めたいのです．これらのデータは，データフレームに入っておりまして，できれば新たなデータフレームを作成できれば良いなぁと思っています．forを使ってある程度のことはできるのですが，いまいち爽快感が得られません． -- [[おさる]] &new{2005-03-14 (月) 03:10:36};~
これら企業の株価それぞれについて，上記のような変化率を求めたいのです．これらのデータは，データフレームに入っておりまして，できれば新たなデータフレームを作成できれば良いなぁと思っています．forを使ってある程度のことはできるのですが，いまひとつ爽快感が得られません． -- [[おさる]] &new{2005-03-14 (月) 03:10:36};~
//
- m<-data.frame(x=runif(10),y=runif(10),z=runif(10));attach(m);(x[2:length(x)]-x[1:length(x)-1])/x[1:length(x)-1] -- &new{2005-03-14 (月) 10:35:27};
-爽快感が得られるかどうか，以下のようなやり方もあるでしょう。~
d <- as.matrix(read.table("temp2.data", header=TRUE))
d2 <- data.frame(apply(d, 2, function(x) {n <- length(x); (x[-1]-x[-n])/x[-n]}))
d2
x y z
2 0.04972376 0.06181319 0.013582343
3 -0.02105263 -0.04657180 -0.028475712
4 0.03225806 0.04477612 0.044827586
5 0.12500000 0.09740260 0.135313531
6 -0.04629630 -0.04378698 -0.027616279
7 -0.06310680 -0.09282178 0.026905830
8 0.02797927 0.03547067 0.005822416
9 -0.03629032 -0.01844532 0.015918958
-- [[青木繁伸]] &new{2005-03-14 (月) 10:35:45};
-素晴らしい！！！こういう一行野郎が欲しかったのです．すっきり爽快です．ありがとうございました． -- [[おさる]] &new{2005-03-14 (月) 12:44:13};
-みっつ上の x[1:length(x)-1] は、x[1:(length(x)-1)] 良くある誤り。動かしてみればすぐ分かることだけど。 -- &new{2005-03-17 (木) 00:10:07};
**Tukeyで比率の統計をしたいのです [#s6caaeb4]
>[[ピンバッジ]] (2005-03-09 (水) 18:02:31)~
~
全く何も分からないので何を書いていいのかも分からないのですが、異なる餌を与えた昆虫の生存率の比較を、RのTukeyで統計しました。しかし、比率を比べる場合はカイ二乗検定後、角変換（アークサイン変換）をした後にTukeyをしなければならないことを知りましたが、この一連の操作をRでどのようにやるかわかりません。また、どういった形のデーターで比較するかもわかりません。~
本当にわからないことばかりですので、お手数ですができる限り詳しく教えていただけると有 難く思います。足りない情報がありましたら教えてください。~
//
-1. R の Tukey とは，具体的にはどの関数ですか~
2. 異なる餌とは，何種類ですか~
3. 昆虫は，それぞれの餌で何匹ずつ飼育したのですか~
4. 一定期間後に，生き残っている昆虫の割合を生存率としたのですか~
5. 生存率（割合）のデータの平均値に差があるか検定したのですか~
6. 「比率を比べる場合はカイ二乗検定後」というからには，平均値の差の検定ではないのですね~
7. 「一連の操作をRでどのようにやるかわかりません」ということは，角変換の公式は分かるのですね。どこが分からないのか，もう少し詳しく。~
8. ひょっとして，あなたのデータは，
当初昆虫匹数 生き残った昆虫匹数 生存率
餌A a1 a2 a2/a1
餌B b1 b2 b2/b1
餌C c1 c2 c2/c1
みたいなものなんですか？~
分からないことだらけで，答えることもできません(^_^;) -- [[青木繁伸]] &new{2005-03-09 (水) 19:07:50};
-御丁寧な質問をありがとうございます。~
質問に関してですが（答えがかみ合ってなかったら申し訳ありません）~
1. どの関数？という質問に対して当てはまらないかもしれませんが、RにTukeyHSDと入力してさせるとうけいです。~
2. 餌の種類は7種です~
3. 6種類の餌は30匹飼育し、1種類の餌は10匹を飼育しました。~
4. はい。何日後に何匹生き残っているかをしらべ、それを生存率としました。~
5. はい。餌の違いによる生存率の違いに差があるかどうかをしらべたいです。~
6. はい。たぶん平均値の差の検定ではありません~
7. 角変換の公式も分かりません。データーをどんな形（比率にしたもの 等）にして読み込ませるかも分かりません。わかることの方が少ないので、ほとんどわかってないと考えていただけたほうがいいと思います。~
8. はい。ほぼそれと同じです。~
御迷惑お掛けしますがどうぞよろしくお願いします。 -- [[ピンバッジ]] &new{2005-03-11 (金) 15:26:47};
-8. の質疑応答からいくと，chisq.test( ) では。それに引き続いての多重比較と言うことなら，pairwise.prop.test( ) かな。pairwise.prop.test の example も，ちょっとは，参考になるかな。~
それにしても，TukeyHSD 関数を，どのように使ったんですか。。。 -- [[青木繁伸]] &new{2005-03-11 (金) 20:05:57};
-ありがとうございます。 よく分かっていなかったので使い方が違った様ですね。大変お恥ずかしいことに、TukeyHSDは、生存しているものを1、死亡したものを0として、日ごとに餌の種類によっての生存の違いを見ているつもりでした。
ありがとうございました。やり直します。-- [[ピンバッジ]] &new{2005-03-14 (月) 13:42:05};
**ESS-5.2.5で補完機能が効かない [#z1d172da]
>[[ぶー]] (2005-03-08 (火) 16:28:53)~
~
こちらに書き込むのはまずいでしょうか？~
Ｒを使うのにESSから使用しているんですが、ESS-5.2.5だと、source や read.table の補完機能が効きません。ESS-5.2.3だとOK。~
elisp を見ていないのでなんなんですが、、、同じ現象で回避されている方がいらっしゃいましたらご指導願えないでしょうか？~
あ、LInux 2.4.7 で使用しています。~
//
-すみません。Ｒは2.0.1でLinuxは2.4.27でした。 -- [[ぶー]] &new{2005-03-08 (火) 16:32:49};
-わたしの所は大丈夫でした.Vine3.1,Rは自家製,カーネルはディストリビューションの最新です. -- [[なかま]] &new{2005-03-08 (火) 18:55:22};
-そうですか、、、私の設定がまずかったかな？時間をみて再チャレンジしてみます。 -- [[ぶー]] &new{2005-03-10 (木) 11:58:18};
**重回帰式の変数選択の方法 [#p16c59f7]
>[[石川誠]] (2005-03-06 (日) 20:43:09)~
~
お世話になります。~
重回帰式を作成するとき、有効な変数を選択する方法を教えて下さい。~
総当り法はあったのですが、変数が多いと時間がかかってしまいます。~
- いろいろあると思いますが，関数 step() などはどうでしょう？関数の意味・使い方等はヘルプをご覧下さい．[[舟尾]] &new{2005-03-06 (日) 20:52:08};
summary(lm1 <- lm(Fertility ~ ., data = swiss))
slm1 <- step(lm1)
summary(slm1)
slm1$anova
//
-そう当たりで法で時間がかかっても結果が出るなら step 法よりも良いのでは？ step 法は必ずしも最良な組合せを選ぶとはかぎりません。データが多すぎる、変数が多すぎる等の理由で、時間が途方もなくかかって実行できないというなら、むしろ問題そのものを見直すべきかも知れません。つまり、すべてのデータ・変数を使うべきか。 -- &new{2005-03-07 (月) 21:16:09};
-総当たり法にしろ，ステップワイズ法にしろ，それらを使って得られた答えが，その学問分野で好意をもって迎えられるかどうか一度先行研究をサーベイなさることをお勧めします。モデルに含めたくない変数が含まれている，あるいはその逆の場合，貴方はその結果をどのように評価するでしょうか。貴方の論文を読む人（査読者，一般読者）はどのように評価するでしょうか。これは，学会の分野にもよるでしょうが。 -- [[青木繁伸]] &new{2005-03-07 (月) 21:24:36};
-非常に勉強になります！コンピュータに自動で変数選択をさせるのは昔から議論があるようで，私も「変数の意味も考えずに選択させるべきでない！」と習った記憶があります．青木先生のおっしゃるとおり，分野にもよるでしょうね．例えばデータマイニングでは変数選択がよくやられているようです（データマイニング自体が良い！悪い！という議論もあるでしょうが^^; ） -- [[舟尾]] &new{2005-03-07 (月) 22:02:20};
-ありがとうございます。薬物設計（ドラッグデザインの探索）では、うまく相関がでれば良いといった風潮です。確かに、GAやニュラールネットといったデータマイニングでよく使われる手法にて変数選択することが多いです。しかし、開発段階になると説明のつかない変数をもってきて説明するのはだめといった状況です。分野によるのは確かです。 -- [[石川誠]] &new{2005-03-08 (火) 23:07:28};
**射影した回帰曲線の色の付け方 [#jb54a256]
>[[さち]] (2005-03-02 (水) 15:37:36)~
~
Rの初心者です。~
方々探したのですが、こんなことはあまりされていないのか、それとも、探し方が悪いのか、ヒントになりそうなものも見つかりませんでしたので、こちらに投稿させて頂きました。分かりにくい文章ですみませんが、どなたか、教えてください。
> f <-function(x) 20/(1+50*exp(-0.8*x))
> x<-1:20
> y<-f(jitter(x))
> (result <-nls(y~a/(1+b*exp(-c*x)), start=c(a=50,b=20,c=0.1)) )
Nonlinear regression model
model: y ~ a/(1 + b * exp(-c * x))
data: parent.frame()
a b c
20.039274 51.908061 0.788543
residual sum-of-squares: 0.3331449
> predict.c<-predict(result)
> plot(x,y,ann=F,xlim=c(0,20),ylim=c(0,20))
> par(new=T)
> plot(predict.c,type="l",xlim=c(0,20),ylim=c(0,20))
この回帰曲線をx軸に射影（投影）したものに、yの値に従って、色を付けることを考えています。~
例えば、
y=0.0 ~10.0 赤
y=10.1~15.0 黄色
y=15.1~20.0 青
のように。collor.paretsにあるheat.colorsのようなもので色づけができたら、とてもうれしいです。~
どうかよろしくお願いします。~
//★★投稿書式くらい，確認してから投稿しなさいね
//
-「x軸に射影（投影）したものに、yの値に従って、色を付ける」というのが，ややこしくて，認知症老人には解釈できなかったのですが，単に，予測値 predict.c にある y の値によって色を変えると言うことでしたら，
points(predict.c,col=ifelse(predict.c>=15.1,"blue",ifelse(predict.c>=10.1,"green","red")))
ということなんでしょうか？？（黄色というのは見にくいので緑にしましたが）~
本来，lines も points も，引数 col は「ベクトルであっても良いが最初の値しか使わないよ」と書いてあるのだが，points はちゃんと色分けできるみたい。lines はだめだ。~
ま，考えてみれば，n番目の点が 14 で n+1 番目の点が 16 なら，それを結ぶ線の色を何にしたらよいか，認知症老人でなくても R も悩むのではないかと思うなぁ。~
どうしても色分けをやりたければ，各点を結ぶ線分と上記色分けの水平線との交点を計算して，その交点を境に，指定した色分けで二つの線分を描くということをすれば良いわけですね。~
色鉛筆で，グラフ用紙にグラフを描くときのことを考えればよいわけだ。それよりは，ある意味めんどうだが。~
全て同じだが，コンピュータに何かやらせようとしたとき，もしそれを誰かに頼んで（命じて）やらせようとしたときに，もし相手が小学生くらいだったらどういう風に説明すれば良いだろうかと考えればよい。コンピュータなんて，言われたことはちゃんとやるけど，ちゃんと教えられないと何にもできない。ちゃんと教えるというのはかなり難しい。日本語で小学生に指示するのと，R の文法にそってコンピュータに指示するのと，基本的には同じ。 -- [[青木繁伸]] &new{2005-03-02 (水) 23:42:44};
-レスが遅れてすみませんでした。丁寧なご指導をありがとうございます。質問者である私が、問題の趣旨を理解できていないのですから、青木さんが混乱されても、仕方ありません。私の方こそ、若年性認知症ですよ。説明が足りない上に、混乱させてしまうような最悪な文章で、本当にすみません。 ~
各点からx軸に下ろした垂線とx軸との交点を繋いだ曲線（この場合は、直線になりますよね？）そのものに色を付けることを考えていました。
色付けの条件としては、前にも書きました通り、yの値に対応して、色を変化させます。
y=0.0 ~10.0 赤
y=10.1~15.0 黄色
y=15.1~20.0 青
こうすると、色分けされたバーコードのようになるかも、、、と思ったのです。~
各点を結ぶ線分と上記色分けの水平線との交点を計算して，その交点を境に，指定した色分けで二つの線分を描く方法につては、これからやってみます。
points(predict.c, col = ifelse(predict.c >= 15.1, "blue",
ifelse(predict.c >= 10.1, "green", "red")))
これは、２次元のグラフの記されている各点を色分けしていますね。３次元版になったら、きっと必要になると思うので、参考にさせて頂きます。~
ありがとうございました。 -- [[さち]] &new{2005-03-03 (木) 11:18:34};
-adhocですがfilled.contourで代用されたらいかがですか?
filled.contour(1:length(predict.c),1:length(predict.c),
matrix(predict.c,nrow=length(predict.c),ncol=length(predict.c)))
とか． -- [[takahashi]] &new{2005-03-03 (木) 12:37:08};
-takahashi さんありがとうございます。早速ためしてみました。とても奇麗な色ですね。~
数値軸上に、少し太めの線として描かせる事しか頭になかったので、うれしくなりました。これを使って色々試してみます。-- [[さち]] &new{2005-03-03 (木) 12:50:32};
-フォーミュラ y= f(x) の逆関数 x=g(y) によって，x1=g(a), x2=g(b) を計算して lines(c(x1, x2), c(0,0), col="color you want") で良いのでは？ -- [[青木繁伸]] &new{2005-03-03 (木) 22:51:36};
-レスが遅れてすみません。青木さん、有難うございます。目から鱗が落ちました。~
早速やってみます!! -- [[さち]] &new{2005-03-04 (金) 16:06:56};
**多項式による回帰分析 [#z83b850b]
>[[highvalley]] (2005-03-01 (火) 19:05:38)~
~
Rを勉強中のhighvalleyと申します。~
現在多項式による回帰分析（最小２乗法）を行いたいのですが、皆さんはどのように行っていますでしょうか？~
私はlm()を使って以下のように行っています。例えば目的変数をy、説明変数をxとして２次の多項式で回帰分析をする場合は、~
~
result <- lm( y ~ 1+x+I(x^2) )~
~
これで目的を達する事は出来ているのですが、この方法だと次数が多くなると大変です。もう少しスマートな方法があるような気がして、質問を致しました。良い方法をご存知でしたら教えて下さい。よろしくお願いします。~
//
-nls()が使えます。[[このpdfファイル:http://phi.med.gunma-u.ac.jp/medstat/it09.pdf]]に簡単な使用例を書いてあります。ご参考まで。 -- [[中澤]] &new{2005-03-01 (火) 19:13:13};
-早速のお返事どうもありがとうございます。資料も読ませて頂きました、大変参考になります。nlsで非線形の回帰ができることがわかったのですが、この関数だとやはり次数が大きくなった場合（例えば６次とか）は、７項並べて書く事になるのですよね？ これは仕方の無い事なのでしょうか。多項式の次数をどんどん上げていった場合、やはりlm()と同じようにスマートではなくなってしまいます。 -- [[highvalley]] &new{2005-03-01 (火) 19:48:50};
-「次数が高くなると大変」ということですが，formula を書く手間はどうってことないでしょう。問題は，各項の相関ではないでしょうか。poly を使うのが良いと思います。~
どうしても lm を使う（書く）手間を省くことだけが目的と言うことなら，簡単なラッパーを書けば良いだけでしょう。しかし，繰り返しますが，次数が高くなるとこの方法はお勧めできません。
takousiki <- function(x, y, k) # パラメータは，独立変数，従属変数，多項式の次数
{
n <- length(x)
z <- matrix(x, n, k)
for (i in 2:k) {
z[,i] <- x^i
}
lm(y ~ z)
}
ヒデー，ラッパーだ -- [[青木繁伸]] &new{2005-03-01 (火) 21:19:19};
-どうもありがとうございます。おっしゃる通り、業務全体から見るとformulaを書く手間は最初のほんの少しなので大した事は無いのですが、それもラッパーを書けば解決できてしまうのですね。だんだん、関数の調べ方や、見つからなかった時の対処の仕方等Rの流儀が見えてきました。どうもありがとうございました！ -- [[highvalley]] &new{2005-03-02 (水) 10:06:23};
//★★重複投稿したのに気づいたら，自分で消してね（気づいてもいないのだろうけど）
//-どうもありがとうございます。おっしゃる通り、業務全体から見るとformulaを書く手間は最初のほんの少しなので大した事は無いのですが、それもラッパーを書けば解決できてしまうのですね。だんだん、関数の調べ方や、見つからなかった時の対処の仕方等Rの流儀が見えてきました。どうもありがとうございました！ -- [[highvalley]] &new{2005-03-02 (水) 10:28:59};
-まず R site search (トップページ参照)を使いましょう（日の下に新しいものはない...旧約聖書）。例えば次のようなパッケージ、関数が見つかります。 -- [[QDU]] &new{2005-03-05 (土) 11:56:32};
polyreg {mda} R Documentation
Polynomial Regression
Description
Simple minded polynomial regression.
Usage polyreg(x, y, w, degree = 1, monomial = FALSE, ...)
Arguments
x predictor matrix.
y response matrix.
w optional (positive) weights.
degree total degree of polynomial basis (default is 1).
monomial
If TRUE a monomial basis is used (no cross terms). Default is FALSE.
... currently not used.
Value
A polynomial regression fit, containing the essential ingredients
for its predict method.
-どうもありがとうございます！ R site searchからpolynomialで検索して、polyregにたどり着く事が出来ました！「多項式」や「回帰」で検索していたのですが、キーワードを英語にすることも必要なのですね。どうもありがとうございました。 -- [[highvalley]] &new{2005-03-07 (月) 13:33:31};
**関数の追加 [#b3e15d07]
>[[しま]] (2005-02-20 (日) 15:29:58)~
~
新しい関数を追加したいのですが，どのようにすればいいのでしょうか？~
初歩的な質問で申し訳ないです．．~
//
-初等以前の問題のような気がします。関連文献を少しは読んでから、再質問して下さい。 -- &new{2005-02-20 (日) 16:27:58};
-簡単なことほど難しいし，何通りも実現方法があったり，少しずつ違っていたりする。まず，一番簡単には，関数の定義はセーブ付きで終了すれば、.Rdata ファイルに定義が残るので、実行のたびに読み込まれ，以後いつでも使えるようになる。~
[初歩的な質問で申し訳ない]と思うなら，質問しない方がよい。 -- &new{2005-02-20 (日) 20:13:38};
-話題の新刊 The R Tips を買って読めば大丈夫. わたしも初心者です. :-) -- [[なかま]] &new{2005-02-21 (月) 13:23:05};
-初心者の答えです。単純には、関数名<-function(引数){関数の中身}で良いのではないでしょうか？こんな簡単なことではないのかもしれませんが。 -- [[MS]] &new{2005-02-21 (月) 20:57:18};
-うむ。それは，関数の定義法では？~
そもそも，「しま」さんが書いた「新しい関数を追加したい」というのはどういうことを指していたのか？？~
私は，「関数を追加して，いつでも使えるようにしたい」と読んだのだが，深読みですか。わはは~
日本語の初心者じゃないのだろうから，もう少しちゃんと説明して質問すべきでしょうねぇ~
というか，私の読解力が低いということで，私の日本語能力に問題があるという，根本的な問題が明白になったとか。わたし，にほんごのしょしんしゃです，やさしくおしえてください -- &new{2005-02-21 (月) 21:17:12};
-参考になるかわかりませんが, RFC1855を読んでみましょう. にほんごのしょしんしゃでも, たいていの言語に訳されていると思います. :-) -- [[なかま]] &new{2005-02-21 (月) 22:08:20};
-蛇足ながら，google などで RFC1855 を検索すると色々出てきます。たとえば，[[http://www.cgh.ed.jp/netiquette/rfc1855j.html:http://www.cgh.ed.jp/netiquette/rfc1855j.html]] など -- &new{2005-02-21 (月) 22:19:57};
**levelplotの背景色の変更 [#s2505fa8]
>[[ひょうひょう]] (2005-02-10 (木) 18:48:00)~
~
グラフィックス参考実例集：ラティスグラフィックスで以下のような説明がありますが、具体的に教えてください。~
「lattice 関数は lightgrey の背景色等幾つかの固有の既定値を採用しているので、必要に応じ par 関数で変更する。」~
par(bg="white")~
としても変更できませんでした。よろしくお願いします。~
//
- [[学生]] (2005-02-10 (木) 20:32:34)
library(lattice)
trellis.par.get()
trellis.par.set(background = list(col="white"))
show.settings()
?trellis.par.get
質問の基本的な姿勢が，以下のサイトに書かれてあります．~
参考:http://myu.daa.jp/osiete/index.html
// ★ 書かれているないようについては申し分ありません
// ★ どうせなら，記事の書き方については ヘルプ を参照のこと
-Q&A, 初心者用Q&A コーナーへの質問テンプレートがあるといいかもしれませんね。最低明らかにしておく要件を入力しないと投稿できないような。いいえ、岡田さんに申し上げているのではありません。 -- [[MKR]] &new{2005-02-10 (木) 22:21:49};
// ★ 何を用意しようとも，それを利用しないのが無法者（初心者）というものです
// ★ その証拠となるのが，善意のちゃんとした意志を持った投稿者さえ，フォーマットに従えないという事実が証明しているのです(^_^;)
// ★ このようなコメントに気づく人も少ないというのが，そのような見解を補強していと思いますが(^_^;)(^_^;)
-ご回答ありがとうございました。~
ぶっきらぼうな質問に答えていただいて、恐縮です。解決しました。質問の背景をもう少し説明すべきというご指摘かと思います。以後気をつけます。過去に同様の質問がないか？などはかなりの時間をかけて調べたつもりだったのですが、もしかして既にどこかに回答があったのでしょうか？ -- [[ひょうひょう]] &new{2005-02-11 (金) 21:12:18};
-トップページにある R site search でキーワード lattice background color で検索してみて下さい。RjpWiki に投稿する手間の何倍も早くヒントがえられますよ(おまけに余計なこともいわれずにすむ :-)) -- &new{2005-02-12 (土) 00:19:22};
-グラフィックス参考実例集の解説は間違いですね(もしくは仕様が変わったか)。「学生」さんがいいたかったことはおそらく、RjpWiki の記事等だけを参考にせず、lattice パッケージのヘルプドキュメントを一応チェックして下さいということかとおもいます。library(help=lattice) とすれば含まれる関数の一覧がえられますから、個別にヘルプドキュメントを通読すれば、作者が必要とおもったこと（そして関連事項へのポインターが）は書かれているはずです。研究の上、グラフィックス参考実例集にまとめをお願いします。皆さんがひょとして想像されているよりもはるかに実質的投稿・回答者は小数です。 -- &new{2005-02-12 (土) 07:32:02};
-おっしゃる通り、R site searchでまったく同じ質問を発見しました。最近Rを使い始めたばかりですが、RjpWikiのおかげで敷居が低く初心者にとっては心強い存在です。グラフィックス参考実例集に投稿できるくらいのレベルになりたいです。 -- [[ひょうひょう]] &new{2005-02-12 (土) 17:15:35};
**Error cannot allocate ..というエラーがいつも出てしまいます [#z147e2d4]
>[[ちょろちゃん(；＿；)]] (2005-02-04 (金) 09:34:32)~
~
みなさま始めまして、R超初心者のちょろと申します。~
初心者のため、ご相談内容にも説明不足の点もあるかと思います。~
その際はご指摘ください。~
~
□相談内容~
~
Windowsで動くRスクリプトがLinux上で稼動すると~
エラーが出てしまう。~
~
□Rのエラー内容~
Error: cannot allocate vector of size 16478 Kb~
Execution halted~
~
vector of size ○○○。。。というサイズは~
32MBの場合もありますし、上記のように16MBの場合もあります~
~
□実行環境~
(スクリプトが動く環境)~
OS:Windows 2000~
Memory:512MB~
R:1.9.1~
(スクリプトが動かない環境)~
OS:Linux Redhat 7~
Memory:2GB~
R:1.9.0~
~
WindowsとLinuxでは、Rを動かす際のメモリーのとり方など~
違いがあるのでしょうか。~
このようなエラーが出てしまうのはRのスクリプトの書き方が~
悪い以外のなにものでもないような気がいたしますが、~
何かお気づきの点等ございましたら、アドバイスいただければ幸いです。~
どうぞよろしくお願いいたします。~
//
-肝心の作業内容 script が無いと判断困難ですね。普通は超巨大なオブジェクトを作ろうとした時にでるエラーですが、お話ではそうでもないような。一つの解決策候補は R のバージョンを最新に変えてみることでしょうか。 -- &new{2005-02-04 (金) 09:50:28};
-アドバイス有難うございます。確かにScriptを提示せず質問するとは私かなりのフトドキものでした。すみません。ただ色々ありましてScriptのどの部分をどのように提示すればよいか迷っております。あと個人的に思っていることとしては、Windows版だとMemory.limit(2000000)などとメモリーを一定量確保しその中でRを動かすという事ができると思いますが、Linuxの場合memory.limitコマンドは見当たらないのですが、、、色々分からないことばかりで大変です。自分でも調査を続けます。お手数をかけますが引き続きアドバイスいただければ幸いです。よろしくお願いいたします。 -- [[ちょろちゃん]] &new{2005-02-04 (金) 14:39:16};
-R をコマンドラインから起動するんでしょうね。R --help とやってみると，以下のようなオプションがあるかもしれません。この数値をいじってみるとか。
--max-vsize=N Set vector heap max to N bytes;
--max-nsize=N Set max number of cons cells to N
どうでしょうね。vsize か nsize か，それともどっちも関係ないか。 -- &new{2005-02-04 (金) 15:05:07};
-現在の R ではよほどクリティカルな作業でない限りメモリー管理を考える必要はないはずです。Linux 版には mem.limit() 関数がありますが、これは MSW 版の Memory.limit() 関数と同じものかな？ 現在の R のメモリー使用状態は gc() 関数の副作用として、表示されます。もしかして、作業途中にとんでもなく大きなオブジェクトが知らぬ間に使われているということはないですか。オブジェクトのサイズは object.size() 関数でわかりますが一時ファイルだと困りますね。大部のファイル読み込みを伴う際には、一時ファイルのサイズが大きくなりすぎることがあるようです。？Memory でいくつかの注意が得られます。 -- &new{2005-02-04 (金) 18:31:51};
-.RDataに巨大なオブジェクトが残っているとか(んでもって毎回読み込んでるとか). ls(all=T)の結果は? -- [[なかま]] &new{2005-02-05 (土) 01:23:28};
- rm(list=ls(all=TRUE)) で隠しオブジェクトを含むすべてのオブジェクトを抹消してから再実行したらどうなりますか？、ということです。 -- &new{2005-02-05 (土) 09:14:40};
-みなさま、色々とあたたかいアドバイスありがとうございます。m(_ _)m。頂いたアドバイスを元にちょっと頑張ってみます。また不明な点がありましたらご相談させていただくかと思いますが、その際はどうぞよろしくお願いいたします。 -- [[ちょろちゃん]] &new{2005-02-09 (水) 13:42:43};
-本件進展がありましたので、ご連絡致します。後々気が付いたのですが、error allocate...というエラーが出ると申し上げておりました実行条件ですが、PerlからRCMDにてRのプログラムを呼び出し実行した時に生じておりました。そこで、該当LINUX上でRによってRのプロンプトを表示し、同じプログラムを走らせて見たところ、WindowsのRと同じように無事処理が終了いたしました。よってLINUX上でプログラムを実行する場合は、「R cmd BATCH ***.R」といったバッチ処理にてプログラムを実行することと致しました。私良く分からないのですが、Perl＋RCMDという組み合わせはメモリーを大量に消費してしまうなどと言ったことがあるのでしょうか？なお、該当LINUXにインストールしていたPERLは5.6.0でした。 -- [[ちょろちゃん]] &new{2005-02-09 (水) 15:54:18};
-「PerlからRCMDにてRのプログラムを呼び出し実行した時に生じておりました」-> 最初に注意すべきでしたね。 -- &new{2005-02-10 (木) 22:18:24};
**Ｒで分割表の対数線形モデル [#j516ba78]
>[[bob3]] (2005-01-31 (月) 00:03:57)~
~
bob3と申します。~
~
Ｒを使った対数線形モデル（log-linear models）による分割表（クロス集計表）の分析をしようとしています。~
ところが２点ほど不明な部分があり、ご相談にあがりました。~
ご相談したいのは以下の２点です。~
~
１）ANOVAコーディングよる対数線形モデルのやり方~
２）関数loglinによる「特定の単一の変数の主効果と誤差のみ」を認めるモデル、および「誤差のみ」を認めるモデルの指定の仕方。~
~
例題として http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/pdf/ch03-01.pdf の表3.2「ロケットの発射試験」を使っています。この例題の飽和モデルによる推定値は http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/pdf/ch03-04.pdf の表3.11にダミーコーディング、表3.12にANOVAコーディングで掲載されています。~
~
まず、Ｒによる対数線形モデルの手順を調べ、以下のようにしてみました。（長くなるので出力は省略しています。）~
~
# 既存の分割表を入力
rocket <- array(data=c(5,7,8,9,3,21,7,9,6), dim=c(3,3))
dimnames(rocket) <- list(c("A1","A2","A3"),c("B1","B2","B3"))
# モデル選択
model0 <- loglin(rocket, list(c(1, 2)), param=TRUE)
model1 <- loglin(rocket, list(1, 2), param=TRUE)
model0
model1
p0 <- 1-pchisq(model0$lrt, model0$df)
p1 <- 1-pchisq(model1$lrt, model1$df)
p0
p1
AIC0 <- model0$pearson-2*model0$df
AIC1 <- model1$pearson-2*model1$df
AIC0
AIC1
# ここでは飽和モデルを採用し分析してみる
rocket.df <- as.data.frame.table(rocket)
colnames(rocket.df) <- c("航続距離","横方向のずれ","度数")
glm0 <- glm(度数 ~ 航続距離 * 横方向のずれ, data = rocket.df, family = poisson)
anova(glm0, test = "F")
summary(glm0)
~
ところが、これで出力されるのはダミーコーディングによる推定のみです。~
~
ANOVAコーディングによる推定値を得るにはどのようにすればよいのでしょうか。~
~
また、関数loglinでモデルを指定する際、[AB]と[A][B]というモデルについては問題ないのですが、[A]のみ、[B]のみ、[ ]のみ（誤差のみ）のモデルについてはどのように指定すればよいのでしょうか。~
~
なお、環境は~
> version
_
platform i386-pc-mingw32
arch i386
os mingw32
system i386, mingw32
status
major 2
minor 0.1
year 2004
month 11
day 15
language R
と、なっております。~
~
よろしくお願いいたします。~
//
-模範的な質問ですね。当てはめ結果の anova 関数によるまとめが R 開発者の考えるところの ANOVA コーディングなのでしょうが、ある文献に出ている例と全く同じ形式での出力を期待するのは虫がよすぎます。引用例のどういう項目が特に必要なのでしょうか？また、lm, glm のモデル式の指定の例は以下のようで大丈夫とおもいますが？ -- &new{2005-01-31 (月) 07:34:46};
glm0 <- glm(度数 ~ 航続距離, data = rocket.df, family = poisson)
glm0 <- glm(度数 ~ 横方向のずれ, data = rocket.df, family = poisson)
glm0 <- glm(度数 ~ 1, data = rocket.df, family = poisson)
- ご返答、ありがとうございます。
~
~ lmやglmでのモデル式の指定方法を教えていただき、ありがとうございます。助かります。これで他のモデルとAICを使った比較が出来ます。
~ （やはりloglinではこれらの式は指定できないのだろうか……）
~
~ 私がANOVAコーディングで知りたいのは、各変数の全カテゴリと全交互作用の標準効果（標準化係数）です。この例題の飽和モデルでいえば、A1、A2、A3、B1、B2、B3、A1*B1、A1*B2、A1*B3、A2*B1、A2*B2、A2*B3、A3*B1、A3*B2、A3*B3、それぞれの標準効果です。それによって、どのような要因がどの程度影響を与えているのか、また有意であるのかどうかを知りたいと思っています。
~
~ よろしくお願いします -- [[bob3]] &new{2005-02-01 (火) 00:33:31};
-トップページの R site search を使って、例えばキーワード loglin anova で検索してみて下さい。MASS パッケージの loglm 関数は lm, glm 風のインタフェイスを持つ loglin 関数へのインタフェイスです。formula 引数が可能で、anova メソッドも存在するようです。 -- &new{2005-02-01 (火) 09:12:51};
-前の方も指摘されていますが、loglin または loglm でパラメータを出力すると、anova コーディングでの結果が得られます。
loglin(rocket,margin=list(c(1,2),c(1),c(2)),param=TRUE)
loglm(freq ~ kyori * zure, data = rocket.df, family = poisson)$param
// ★★ こうやるのよ！
-- [[なかの]] &new{2005-02-04 (金) 22:11:04};
//
-返答を下さった皆様、ありがとうございます。勉強になります。
~実は、私がANOVAコーディングで得たいと思っているのは「標準化係数（標準効果）」なのです。
~これは重回帰分析における標準化偏回帰係数に相当するものですが、loglin や loglm で得られるのは「推定値（効果）」だけのようです。
~標準化係数（標準効果） ＝ 推定値（効果） ÷ 標準偏差（標準誤差） なので、標準偏差（標準誤差）が得られれば何とかなりそうなのですが……
~なかなか難しいものですね。
-- [[bob3]] &new{2005-02-05 (土) 02:40:35};
-純粋に質問です。bob3 さん、loglinear における「標準化係数」というのはどのような意味があるのでしょうか？もしくは、それが解説されている本などを教えて頂けるとさいわいです。 -- [[774]] &new{2005-02-05 (土) 10:09:26};
- 標準化係数（標準効果）それぞれのカテゴリとその交互作用が有意かどうかの検
定をしたいと思っています。
~ ……ん、待てよ。AICでモデルの選択をするんだったら、そのあとで変数やカテゴ
リごとに検定をするのはおかしいのかな？
~ 重回帰分析でAICを使って変数選択するときは、有意でない変数も含めたモデルを
採用することもありますね。
~ ということで、ちょっと分からなくなってきました。すみません。
~ 参考にしている文献ですが、これも実は私のほうが教えていただきたいぐらい
で、「質的情報の多変量解析」http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/index.htmlぐらいしかありません。
~ Ｒでの対数線型モデルについてということでは、主に以下の文書を参考にしています。
~ http://www.ci.tuwien.ac.at/~zeileis/teaching/Biostatistics03/examples4.pdf
~ http://www.stat.ohio-state.edu/~tjs/865/handout-3.pdf -- [[bob3]] &new{2005-02-07 (月) 00:14:15};
-御指摘の「標準化係数」というのはパラメータ推定値のz得点のことなのですね。カテゴリ毎のパラメータを検定してモデルを構築するという方法をとったことがなかったので、なぜ「標準化係数」というのが必要なのか、一瞬理解できませんでした。普段は、モデルの適合度と残差のチェックからモデルを構築しています。対数線形モデルでは個別のパラメータが有意であるかどうかよりも、関係（オッズ比）に注目することの方が個人的には多いです。 -- [[774]] &new{2005-02-07 (月) 06:09:08};
**Rでメタアナリシス [#y120ea9f]
>[[竹内]] (2005-01-30 (日) 09:06:30)~
~
Rでメタアナリシスを行いたいのですが、どなたか参考になるもの（教科書、webなど）をご存知でしょうか？google、amazonなどで一応探してはみましたが、「これ」というものがみつかりませんでした。英語のものやS言語のものでも結構ですのでどなたかご存知のかた、よろしくお願いいたします。~
//
-探して出てきたけど，これというものではなかったもののリストを書いておくというのはいかが。そうでないと，自分で勧めようと思おうものが，これというものでないと判断されたものだったらいやですものね。 -- &new{2005-01-30 (日) 11:07:40};
-丹後俊郎先生著の，メタアナリシス入門（朝倉書店）はいかがでしょうか？S-Plusのコードがついてます．また，packageは，rmetaが利用できます． -- &new{2005-01-30 (日) 13:01:32};
-[[R(フリーの統計ソフト)によるメタアナリシス手順:http://www.genome.med.kyoto-u.ac.jp/ra/rmeta/menu.htm]]~
[[Microsoft Excel spread sheet for Meta-analysis:http://www.kdcnet.ac.jp/hepatology/technique/statistics/meta.htm]]~
[[ここからはじめるメタ・アナリシス?Ｅｘｃｅｌを使って簡単に:http://bookweb.kinokuniya.co.jp/htm/4880037001.html]] -- &new{2005-01-30 (日) 14:30:28};
-"R(フリーの統計ソフト）によるメタアナリシス手順”をホームページに掲載した者ですが、Rのバージョンアップに伴ってか、記載中のmeta.DSL関数はエラーが出るようになっている模様です。遺伝子関連解析（SNPアレル頻度比較）のメタアナリシスを想定した、簡易関数を http://d.hatena.ne.jp/ryamada22/20051126/1132980110 に掲載しました。ローカルファイルを読み込んでMantel-Haenszel法の結果オブジェクトとプロット図を返します。 -- [[ryamada]] &new{2005-11-26 (土) 15:49:27};
**n×nマトリックス状にplotを配列する方法 [#q5c5cba6]
>[[kejuyan]] (2005-01-26 (水) 17:28:06)~
// ★★ 元発言者から何の反応もないですが，もう十分だと思いますので，編集しておきますワ
昨日は、大変有益なアドバイスを頂きました。ありがとうございました。~
今日は、例えば、2列のデータX00Y00からX19Y19までのデータをプロットし，20×20のマトリックス状に配置したいと考えております。~
2×2のマトリックスであれば，以下のようにも書くものです。
par(mfrow=c(2,2))
par(mar=c(0,0,0,0))
plot(y~x,X00Y00)
plot(y~x,X00Y01)
plot(y~x,X01Y00)
plot(y~x,X01Y01)
昨日教えていただいたプログラムをアレンジして，以下のプログラムを作成しました。
par(mfrow=c(20,20))
par(mar = c(0, 0, 0, 0))
for (i in 0:19) {
if (i < 10) I <- paste("0",i,sep="") else I <- as.character(i)
for (j in 0:19) {
if (j < 10) J <- paste("0",j,sep="") else J <- as.character(j)
oname <- paste(“X”,I,”Y”,J, sep="")
y <- paste(“ Y_”,oname)
x <- paste(“X_”,oname)
plot(y~x, oname, col="red",lty=1, type="l",axes=FALSE, frame = TRUE)~
}
}
しかし，
Error in eval(expr, envir, enclos) : invalid second argument
というエラーが出てきました。~
そこで、xの中身を調べてみたら，以下のように ”（二重引用符）が付いていました。
> x
[1] "I_ X00Y00"
なので、この " を無くしてやれば動作するのではと考えているのですが、" の取りかたが分かりませんでした。~
" の取り方、またはより良い方法をご存知の方、御教授よろしくお願いします。~
-一つの解答案を示しておきましょう。でもそれを示すだけでは，やさしい対応でもなんでもないのですけどね。やむをえない。~
par(mfrow=c(2,2), mar=c(0,0,0,0))
for (i in 0:1) {
for (j in 0:1) {
plot( read.table(sprintf("X%02iY%02i.txt", i, j)))
}
}
読み込んだデータフレームの列名を変えたりしていたのは，plot の引数で与えるだけでよいことになるでしょう。~
結局，これだけならば，データフレームを全部読み込んで保持しておく必要はさらさらないということですね。~
たいていの場合，読み込んだデータフレームを全部保存しておく必要はないでしょうね。統計処理にしても，処理結果を保存しておくだけでたいていの場合は大丈夫なはずです。というか，そういうようにプログラムを書く方が幸せになれるでしょう。 -- &new{2005-01-26 (水) 18:47:19};
-（参考まで）paste もしくは sprintf 関数は「文字列」を作るだけです。文字列 oname を「名前」に持つ変数を作り、それに値 x を代入するのが assign(oname, x) です。 -- &new{2005-01-26 (水) 21:45:35};
**複数ファイルの一括読み込み [#z3529c14]
>[[kejuyan]] (2005-01-25 (火) 15:09:18)~
~
// ★★ どっちもどっちがだ，書き込みのルールや書き込み方くらい，読んでから書き込もうよね。
座標(x,y)で得られたn行2列データファイルX00Y00.csv〜X50Y50.csv【X00Y00は(x,y)=(0,0)でのデータ】があります。~
~
1つのファイルであれば以下のようにX_<ファイル名>,Y_<ファイル名>の形式で読み込むことが出来ました。~
これを多数のファイルで一括読み込みしたい場合はどのようにしたらよろしいでしょうか？
> X00Y00 <- read.csv("C:/X00Y00.csv")
> X00Y00
X Y
1 1 2
2 2 4
3 3 6
4 4 8
5 5 10
6 6 12
7 7 14
8 8 16
9 9 18
10 10 20
> colnames(X00Y00) <- c("X_X00Y00","I_X00Y00")
> X00Y00
X_X00Y00 I_X00Y00
1 1 2
2 2 4
3 3 6
4 4 8
5 5 10
6 6 12
7 7 14
8 8 16
9 9 18
10 10 20
forなどを使えば出来そうかとは思うのですが、プログラム経験も乏しいので良く分かりませんでした。~
どなたか御教授よろしくお願いします。~
また、参考URL or 文献なども教えていただけるとうれしいです。~
//
//★★ リンク間違ってましたので修正しておきました
-まずチップス集は最低みましょうね。R に[[ファイルを読み込む tips 集(暫定版)]] あたりは参考になりませんか。 -- &new{2005-01-25 (火) 15:53:42};
-そんなにたくさんの（csv）ファイルを用意できませんので実地には確かめていませんが、次のようなやりかたはどうでしょうか。 ファイル X00Y00.txt, X01Y00.txt, X00Y01.txt, X01Y01.txt (実体は同じ)があると仮定しています。-- &new{2005-01-25 (火) 18:01:21};
// ★★ ここにはもともと空行があったのだけど，空行があると次の引用部分のインデントが変になる
// ★★ 書き込みの結果を見ればおかしいのが分かるのだから，直せばいいのに
> for (i in 0:1) {
if (i < 10) I <- paste("0",i,sep="") else I <- as.character(i)
for (j in 0:1) {
if (j < 10) J <- paste("0",j,sep="") else J <- as.character(j)
oname <- paste("X",I,"Y",J, sep="") # 文字列 "X00Y00" 等を作る
fname <- paste(oname, ".txt", sep="") # ファイル名文字列 "X00Y00.txt" 等を作る
x <- read.table(fname) # ファイルを読み込み
colnames(x) <- c(paste("X_",oname,sep=""), paste("I_",oname,sep=""))
assign(oname, x) # 文字列 oname を名前に持ち、中身が x のオブジェクトを作る
}
}
> X00Y00
X_X00Y00 I_X00Y00
1 1 2
2 2 4
3 3 6
4 4 8
5 5 10
6 6 12
7 7 14
8 8 16
9 9 18
10 10 20
> X01Y01
X_X01Y01 I_X01Y01
1 1 2
2 2 4
3 3 6
4 4 8
5 5 10
6 6 12
7 7 14
8 8 16
9 9 18
10 10 20
===== 別の記述方法
for (i in 0:1) {
for (j in 0:1) {
oname <- sprintf("X%02iY%02i", i, j) # 文字列 "X00Y01" 等を作る
x <- read.table(sprintf("%s.txt", oname)) # ファイルを読み込み
colnames(x) <- c(sprintf("X_%s",oname), sprintf("I_%s",oname))
assign(oname, x) # 文字列 oname を名前に持ち、中身が x のオブジェクトを作る
}
}
//-返信を下さった方々！ありがとうございました。 -- [[kejuyan]] &new{2005-01-25 (火) 18:22:51};
//★★ 発言の重複もあるし，改行できてないし。
//★★ 直しておくけど，気づくこともないのだろうな。
-返信を下さった方々！ありがとうございました。~
まず最初にアドバイス通り、ファイル読み込みTipを再読し、RプログラミングTipと組み合わせることで
filenames <- list.files(path = ".", pattern = NULL, all.files = FALSE,
full.names = FALSE, recursive = FALSE)
listoftables <- lapply(filenames, read.csv)
names(listoftables) <- filenames
listoftables
という記述をすれば、listが作製できるので、後はlistをうまいこと分解していけばよいのかな？と考えていました。~
ここでしばらく悩んでいたところ、プログラムを書き込んで頂いたので、それを試してみました。~
そしたら、、、私の欲しかったプログラムはこれだ！！というくらい完全に動作しました。~
返信を下さった方々本当にありがとうございました。 -- [[kejuyan]] &new{2005-01-25 (火) 18:35:59};
-(老婆心ながら) 51x51 個のオブジェクトを作っても後が大変なのでは。むしろリストを行列にして
listoftables[[12, 34]]
等と添字操作する方が簡便なのでは。リスト listoftables は
dim(listoftables) <- c(51,51)
とすれば行列として操作できます。 -- &new{2005-01-25 (火) 20:43:26};
-本来の解析内容によりますが、もし X..Y.. がすべて同じサイズならばむしろ4次元配列 XY として読み込む方がアクセス速度は早くなりそうですね。つまり XY[1,1,.,.] が X00Y00 になるようにする。添字が 0-50 ではなく 1-51 に なるのがすこし面倒ですが。 -- &new{2005-01-26 (水) 07:20:30};
-そもそも，全部読み込んで保持しておく必要があるのかどうなのか。 -- &new{2005-01-26 (水) 09:29:50};
-保持しておく必要は，とりあえずはなさそうですね。 -- [[青木繁伸]] &new{2005-01-26 (水) 17:49:11};
**関数を上書きしてしまった？ [#z8e67a89]
>[[初心者]] (2005-01-17 (月) 20:18:02)~
~
どうも以前にgamma関数を上書き定義してしまったらしく、使えなくなってしまいました。一度Rを閉じても元に戻りません。初期定義に戻すにはどうすればいいでしょうか。~
//
-関数の再定義をしてしまって，R を終わるときに q() の問い合わせに yes を答えたのでしょうね。~
R を起動した後，rm(gamma) で再定義を消去して，q() でもう一度 yes を答えましょう。もう一度立ち上げると，本来のgamma関数が復活します。たぶん。~
この手順で駄目ならば，上書きしてしまった経緯・状況をもっとちゃんと説明してください。 -- [[青木繁伸]] &new{2005-01-17 (月) 20:29:23};
-ありがとうございました！回復しました。 -- [[初心者]] &new{2005-01-18 (火) 16:44:34};
-てっ取り早い方法としてワーキングディレクトリのファイル .RData を消去すれば私的定義のオブジェクトはすべて消えますから、本来のていぎが復活します。 -- &new{2005-01-19 (水) 08:25:40};
-.RData を消去すれば，消去されたくない必要なものがあっても，もろともですから，お勧めはしませんでした。 -- [[青木繁伸]] &new{2005-01-19 (水) 09:03:33};
**分割表形式のデータフレームを通常のデータ形式に戻したい [#t000dd1f]
>[[Jack]] (2005-01-16 (日) 17:00:47)~
~
次のようなデータフレームがあります。~
V1, V2は処理の種類（カテゴリー）を表し、1,2,3,4,5は計測値１〜５（間隔尺度）の頻度（個数）を表しています。
V1 V2 1 2 3 4 5
1 1 1 3 1 0 2
1 2 3 1 4 5 5
2 1 2 1 9 8 7
2 2 1 0 8 8 1
解析に際し、これを生データに戻したいのです。つまり
V1 V2 V3
1 1 1
1 1 2
1 1 2
1 1 2
1 1 3
1 1 5
1 1 5
1 2 1
1 2 1
1 2 1
1 2 2
.
.
というようにです。~
forループを使えばできそうな気がしたのですが、肝心の計測値を抽出するところや、抽出した値をどう反映させればいいのかわかりませんでした。~
あるいは、forを使わない方法があるのでしょうか。
よろしくおねがいします。~
//
-例えば -- &new{2005-01-16 (日) 18:04:01};
> x
V1 V2 X1 X2 X3 X4 X5
1 1 1 1 3 1 0 2
2 1 2 3 1 4 5 5
3 2 1 2 1 9 8 7
4 2 2 1 0 8 8 1
> xx <- as.matrix(x) # 一旦行列に直す方が操作しやすい
> xx
V1 V2 X1 X2 X3 X4 X5
1 1 1 1 3 1 0 2
2 1 2 3 1 4 5 5
3 2 1 2 1 9 8 7
4 2 2 1 0 8 8 1
> y <- numeric(0)
> for (i in 1:dim(xx)[1]) {
for (j in 3:dim(xx)[2]) {
if (xx[i,j] != 0) {
for (k in 1:xx[i,j]) {
y <- c(y,xx[i,1:2],j-2)
}
}
}
}
> y <- matrix(y, nc=3, byrow=TRUE) # 行列に直す
> colnames(y) <- c("V1","V2","V3") # 列名を付けたければ
> y <- as.data.frame(y) # データフレームにしたければ
> y
V1 V2 V3
1 1 1 1
2 1 1 2
3 1 1 2
4 1 1 2
5 1 1 3
6 1 1 5
7 1 1 5
8 1 2 1
9 1 2 1
10 1 2 1
(途中省略)
66 2 2 4
67 2 2 4
68 2 2 4
69 2 2 4
70 2 2 5
==========================
# 別法
> y <- matrix(0, nr=1000, nc=3) # 十分大きな行数の行列を用意
> n <- 0
> for (i in 1:dim(xx)[1]) {
z <- rep(1:5, xx[i,3:7]) # 例えば i=1 ならこれはベクトル 1 2 2 2 3 5 5
N <- n + length(z)
y[(n+1):N, 1:2] <- xx[i, 1:2]
y[(n+1):N, 3] <- z
n <- N
}
> y <- y[1:n,] # 実際の行数に切り詰め
-早速ありがとうございました。行列にいったん直すのですね。検討してみます。 -- [[Jack]] &new{2005-01-16 (日) 20:46:31};
-すんごく，短い解があると思います。for ループも不要です。~
# テストデータ作成
x <- data.frame(V1=c(1,1,2,2), V2=c(1,2,1,2), X1=c(1,3,2,1), X2=c(3,1,1,0),
X3=c(1,4,9,8), X4=c(0,5,8,8), X5=c(2,5,7,1))
# 以下の4行が解 =============上とは出現順が違うが=================
f <- as.matrix(x[,3:7])
v1 <- x[,1]
v2 <- x[,2]
cbind(rep(v1[row(f)], f), rep(v2[row(f)], f), rep(col(f), f))
これがどういうことかは，~
[[http://aoki2.si.gunma-u.ac.jp/R/tenkai.html:http://aoki2.si.gunma-u.ac.jp/R/tenkai.html]]~
を見た方がわかりやすいかも。その応用です。~
行の名前が変になるけど，3行解として
f <- as.matrix(x[,3:7])
v12 <- as.matrix(x[,1:2])
cbind(v12[rep(row(f), f),], rep(col(f), f))
というのでもよいかも。~
単に行数を減らすだけなら，
f <- as.matrix(x[,3:7])
cbind(as.matrix(x[,1:2])[rep(row(f), f),], rep(col(f), f))
いきなり，これが出てくるとわけわからんですね。 -- [[青木繁伸]] &new{2005-01-16 (日) 21:33:44};
-ありがとうございました･･･しかしまだ理解できていません。あとでじっくり追ってみたいと思います。青木先生のサイトは大変勉強になるので、よく利用させてもらっています。 -- [[Jack]] &new{2005-01-17 (月) 08:30:01};
**shapefileが読めません。 [#x356afb4]
>[[おやじっち]] (2005-01-12 (水) 00:35:36)~
~
spdepでmoran's Iを計算しようと思い、vectorworksで作成したファイルをdxf形式になおし、ArcMAPで読み込んだものをshp形式で取り出しました。その際属性でデータを入れました。~
Rでmaptools,spdep,tripackをロードし
// ★★ 記述法をよく理解してくださいよ
x <- read.shape(system.file("shapes/ファイル名.shp", package="maptools")[1])
と入力しましたが
Error in read.shape(system.file("shapes/ファイル名.shp", package = "maptools")[1]) : ~
unable to open file
と表示され読めませんでした。~
http://web.sfc.keio.ac.jp/~maunz/~
のサイトの~
R language/空間重み付け行列とMoran’s I~
の部分を参照したのですが、なにぶん素人なものでどこが間違っているのかよく分かりません。~
どうか教えて下さい。~
//
-"system.file("shapes/ファイル名.shp", package="maptools")[1]"の示すパスにシェープファイルはあるのでしょうか？お使いのＯＳは？Ｒのバージョンは？maptoolsのバージョンは？ -- &new{2005-01-12 (水) 05:08:02};
-OSはmac OS10.3.7でRは2.01でmaptoolsは0.4-8です。パスを指定するやり方がよく分かりませんでした。ちなみにこの"ファイル名.shp"というファイルはデスクトップの"ファイル名"のフォルダの中に入れてあるのですが。本当に初歩的ですみません。 -- [[おやじっち]] &new{2005-01-12 (水) 16:36:47};
-パスの指定法~
Mac OS X の場合，ユーザ名 foo のデスクトップに bar フォルダを作って，その中に baz.dat というファイルを作ったとすると，そのファイルを読むには，
scan("/Users/foo/Desktop/bar/baz.dat")
とすればいいのだが，毎回そんなの書くのいやだということなら，command+Dでワーキングディレクトリを bar ディレクトリにすれば
scan("baz.dat")
だけで読める。 -- [[青木繁伸]] &new{2005-01-13 (木) 21:37:45};
-多分、exampleを見られてsystem.fileとされたのだと思いますが、ご自身のshapeファイルの場合は、そのまま read.shape("パス名+ファイル名")で、OKです。頑張って下さい。 -- [[谷村]] &new{2005-01-17 (月) 21:03:23};
-read.shape("パス名+ファイル名") で読み込めました。有り難うございました。 -- [[おやじっち]] &new{2005-01-19 (水) 05:46:59};
**Windows版Excelからのセル範囲コピー [#vacd4c6c]
>[[青木繁伸]] (2005-01-11 (火) 14:32:44)~
~
Windows ユーザから以前，~
http://aoki2.si.gunma-u.ac.jp/R/excel.html~
の動きがおかしいという問い合わせがありました。~
複数列をコピーしてもなぜか一行中の全部の数字を連結した文字列になる（タブで区切られない）ので，書かれているような動きをしないということでした。その人は1.9でやっているからでしょうかと言っておりましたが，今日たまたまWindows版のRをさわっていて2.0.1でもちゃんと期待通り動かないことが分かりました。~
なんででしょうか。~
エディタなどにペーストすると列の間にはちゃんとタブコードが入っているのですが。~
//
-実験してみました。おそらく範囲指定をした EXCEL データの「セル書式設定」が「数値」になっていないのが原因かと思われます。
-- 標準：一行中の全部の数字を連結した文字列になる~
-- 数値：行列としてちゃんと読み込まれる~
「標準」だけでなく「数値」以外の形式（「日付」や「文字列」など）はちゃんと読み込んでくれないようです。~
それにしても、非常に便利な汎用関数ですねぇ〜！ -- [[舟尾]] &new{2005-01-11 (火) 15:40:58};
-多分1.9.0以降はRguiではtabの入力がスキップされる
/*
* Filter R commands out of a string that contains
* prompts, commands, and output.
* Uses a simple algorithm that just looks for '>'
* prompts and '+' continuation -- won't work when
* other prompts are used (e.g., as in a debugging
* session.)
* Always return the length of the string required
* to hold the filtered commands.
* If cmds is a non-null pointer, write the commands
* to cmds & terminate with null.
*/
からかと.なのでRtermでは上手くいくかも. -- [[なかま]] &new{2005-01-11 (火) 16:08:17};
-じゃなかった,Rguiでは元々tabの入力を考慮してなかったりしました.Rtermもたまたまのような気がします. -- [[なかま]] &new{2005-01-11 (火) 16:37:36};
-ありがとうございました。Macintosh版（V.X版って，へんな表記だなぁ）ではなんの問題もなかったので疑っても見なかったです。~
それで，以下のような実験をしました。驚くべき結果でした(^_^)~
WIndows の Excel で以下のようなファイル test.txt を作りました。4行2列で，1，2行は「標準」，3，4行は「数値」の書式です。
ファイルの中身
$ cat test.txt
1 1
2 1.414213562
3.00000 1.73205
4.00000 2.00000
ダンプ
$ hexdump -c test.txt
0000000 1 ?t 1 ?r ?n 2 ?t 1 . 4 1 4 2 1 3 5
0000010 6 2 ?r ?n 3 . 0 0 0 0 0 ?t 1 . 7
0000020 3 2 0 5 ?r ?n 4 . 0 0 0 0 0 ?t
0000030 2 . 0 0 0 0 0 ?r ?n
なんとオバカな仕様で，書式が「標準」のときには「タブ」のみ，「数値」のときには「空白+タブ」が入っているようです。というか，数値の後に空白を挿入しているみたいで，行末でも空白一個が入ってます（みっともない）。~
これをコピーして，R に持ってくるとき，Windows版 R(Rgui) では，確かにタブだけだと無視されるようです。「数値にしたら読めるというのは，Excelが余計に付加した空白の「おかげ」であるようです。~
Macintosh版のExcelでも，「数値」のときには「空白+タブ」という仕様は同じですが，Macintosh版のRでは，「タブ」は有効です。つまり，「標準」でも，「数値」でも同じように読めます。 -- [[青木繁伸]] &new{2005-01-11 (火) 17:27:45};
-なるほど！勉強になりました！Windows 版では「数値」の他に「会計＋記号なし」にすると問題なく読み込むことが出来るので『なんで？』と思っていたのですが、この場合は「空白＋タブ＋空白」が入っていました（行末は「空白＋改行文字」）。 -- [[舟尾]] &new{2005-01-11 (火) 18:40:56};
-matrix(scan("clipboard"),byrow=T,nc=2)にすればWindowsでも大丈夫なので，Windowsではコンソールデバイスからのリダイレクトの扱いに問題があるんでしょうねぇ。 -- [[中澤]] &new{2005-01-11 (火) 18:58:42};
-数値の前の空白はもしや符号だったりして... -- [[なかま]] &new{2005-01-11 (火) 19:15:23};
-なるほど,Winのクリップボードはバイナリーの場合もあるので,"clipboard"指定でも,textかどうか判別して読み込んでいますね.stdin代りの部分だとtabを抑止しているのは単にカーソル制御が面倒なだけかもしれませんが... -- [[なかま]] &new{2005-01-11 (火) 19:32:05};
-ちなみに，Windowsで青木先生の関数の""の部分を"clipboard"にすれば，手順は1,2だけで終わりますから，その方が楽かもしれません。 -- [[中澤]] &new{2005-01-12 (水) 11:19:00};
-なるほど。Mac でも clipboard という装置名を試してみたけどだめだったんです。Windows なら，その方がいいですね。 -- [[青木繁伸]] &new{2005-01-12 (水) 13:56:45};
**Xでない環境でpng()やjpeg()を使いたい [#y2108482]
>[[はいじま]] (2005-01-10 (月) 23:33:31)~
~
// ＊＊＊＊＊ ~ を使って不必要な物理的な行替えをしないように
LinuxでXの無い環境で、ver2.0.0を使っています。~
R Bookを読みつついろいろ試していますが、
Windows上ではpng()が利用できるのですが、
Linuxで、コマンドラインからRを実行するとエラーが出ます。~
Xはインストールされていないサーバーなのですが、
この場合、画像ファイルの作成はどのように行うのでしょうか。~
それとも行えないのでしょうか。
// ＊＊＊＊＊ ヘルプを見たり，テキスト整形のルールを表示させたりして，投稿法を学んでから投稿しましょうね。
// ＊＊＊＊＊ 今回はなおしておきます。
src-----------------------------
x<-1:100
y<-sin(x)
png('/tmp/sin.png')
plot(x,y)
dev.off()
src-----------------------------
err-----------------------------
Error in X11(paste("png::", filename, sep = ""), width, height, pointsize, :
unable to start device PNG
In addition: Warning message:
unable to open connection to X11 display`'
err-----------------------------
//
-RjpWiki にある R-FAQ の訳の 7.21 を見て下さい。 pictex ドライバを使うのも一つの手かも知れません。-- &new{2005-01-11 (火) 08:15:19};
-pdf ファイルだとXはなくてもいいようですが，pdf を使う（必要なら png に変換する）ことで切り抜けられませんか？ -- [[青木繁伸]] &new{2005-01-11 (火) 15:12:57};
- X が使えない環境を利用できないので確認不可能ですが、?capture.output にでている次のようなやりかたはうまくいきませんか？一度 ps ファイルができれば Imagemagick 等で他フォーマットに変換するのは簡単なはずですが。 -- [[MKR]] &new{2005-01-11 (火) 18:20:49};
## on Unix with enscript available
ps <- pipe("enscript -o tempout.ps","w")
capture.output(example(glm), file=ps)
close(ps)
-ありがとうございます。 PDFについては設定無しに実行できました。とてもきれいです。前もってPDFLibがインストールしていた環境だからできたのでしょうか?? psについては残念ながらcapture()がNULLと返しました。webで使いたいと考えていますので、次はxvfdを試してみようと思います。 -- [[はいじま]] &new{2005-01-11 (火) 22:43:12};
-無事Xvfbでpngファイルを作成できました。Rを起動する前にexport DISPLAY=":1.0"; /usr/X11R6/bin/Xvfb :1 -screen 0 1024x768x24& でXを起動しました。 -- [[はいじま]] &new{2005-01-12 (水) 01:23:27};
**パッケージfSeriesのaparchSimについて [#e666acf8]
>[[winga]] (2005-01-05 (水) 20:57:16)~
~
//あのね，文章を物理的に ~ で区切って一行の文字数を決めるんじゃないのよ！！なおしておくから，これからはちゃんとしなさいね。by 青木繁伸
windowsXPでR2.0.1を使用しているものです。パッケージfSeriesでaparch分析をしようとしています。aparchSimというコマンドでシミュレーションを行っているのですが、出力された結果には負の値が含まれてしまうことが良く分かりません。aparchSimは自分で設定したパラメータをもとに将来の条件付標準偏差を求めていると思うのですが、この値が負になることはないと思うのですが・・・。help(aparchFit)を読んでも書いていないので、どなたか分かる方がいらっしゃったら教えていただければと思います。~
//
-example(aparchSim)でやってみると，たくさんの数値が出力されるようですが，どこの数値が負の値になるといっておられるのでしょうか。負になりそうなところはResiduals:でしょうが，それは負がでてもあたりまえというか。ま，私にはよくわからない関数（分析手法）なんで，よけいな口出しは無用なんですが。 -- [[青木繁伸]] &new{2005-01-05 (水) 22:28:18};
-どんなパラメータを設定したのかをaparchSimのコマンドで書いてみてください。時間がある時に見てみます。 -- [[矢野]] &new{2005-01-05 (水) 23:55:21};
-青木先生、矢野様、ありがとうございます。もう少し自分でも考えてみます。また、tseriesについてなのですが、R1.9.1まではlibrary(tseries)が実行できたのですが、R2.0.1にアップデートしてからは「Error: package 'quadprog' could not be loaded」となってtseriesを使うことができません。どなたか解決法を教えて頂けないでしょうか。。宜しくお願いいたします。 -- [[winga]] &new{2005-01-09 (日) 01:03:46};
-別の質問は，別のスレッドを起こす方が良いと思いますが。~
プラットホームが書いてありませんが，Windows なんでしょうね。~
Mac だとちゃんと動くようですが（というか，そんなエラーはでませんが）。なんででしょうね。~
バージョンダウンするしかないかも。
> library(tseries)
Loading required package: quadprog
'tseries' version: 0.9-24
'tseries' is a package for time series analysis and
computational finance.
See 'library(help="tseries")' for details.
>
こうなります。 -- [[青木繁伸]] &new{2005-01-09 (日) 01:48:07};
-quadprogというパッケージがインストールされていないか、インストールされていても（何らかの理由で）壊れてしまっているのではないかと思います。CRANからquadprogをインストールしてみてください。 -- [[矢野]] &new{2005-01-09 (日) 14:46:55};
-ありがとうございます。quadprogをインストールした結果、うまくいきました。以降違う質問は別スレッドを起こすように致します。 -- [[winga]] &new{2005-01-10 (月) 11:36:18};
**ブラウザーに関する質問二件 [#c1bbb8a4]
>[[QDU]] (2005-01-04 (火) 08:13:30)~
~
R とは直接関係ありませんが、RjpWiki を見るブラウザーについて前から困っていることがあります。Mozilla で見ると半角英数字が空白になることがある。例えば、このページのトップ見出し中の 「R」が消えてしまいます。一方「RjpWiki」は問題無く表示されています。しかたなしにいつも Konqueror で見ていますが、このページ(だけ)が左右の欄が一部重なりあって見にくくなります(一方 Mozilla では問題無し)。何か解決のヒントをご存知の方いらっしゃいませんか。使用環境は Linux (Knoppix 3.6) です。~
//
**自前のカラーバーを作る方法 [#a2e9bbdd]
>[[ゆき]] (2005-01-03 (月) 18:17:22)~
~
自前でカラーバーを作りたいのですが、どうすればよいの困り果てています。~
~
例えば、グラフィックス参考事例集にある火山の地形図のイメージでは、~
image2 <- function () {
data(volcano)
x <- 10*(1:nrow(volcano))
y <- 10*(1:ncol(volcano))
png("image2.png") # png デバイスを開く
# 地形図色調で色分けしてイメージ表示
image(x, y, volcano, col = terrain.colors(100), axes = FALSE)
# 等高線を重ねる
contour(x, y, volcano, levels = seq(90, 200, by=5), add = TRUE, col = "peru")
# 下部に軸、チックマークを描く
axis(1, at = seq(100, 800, by = 100))
# 左部に軸、チックマークを描く
axis(2, at = seq(100, 600, by = 100))
# 全体を囲む枠を描く
box()
# タイトル
title(main = "Maunga Whau Volcano", font.main = 4)
dev.off() # デバイスを閉じる
}
として、 terrain.colors を使用して色をつけていますが、これを、５０Ｍや１００Ｍ間隔で色を自分で割り振って、色付けを行いたいのですが、どうすれば良いのでしょうか。~
//
-御質問の趣旨「自前のカラーバー」が今一わかりませんが、おそらく (1) 同じ色を割り当てる範囲を自分で指定したい、(2) 範囲を表す色を自分で指定したい、のどちらか（もしくは双方）かと思います。このためには image 関数より filled.contour 関数の方が便利だと思います。以下の二つの例は example(filled.contour) 中の例をいじくったものです -- [[QDU]] &new{2005-01-03 (月) 21:11:51};
//引用（先頭一字半角開け）の前には空白を置くな
//注釈はおいても良いけど，空行は置くな
//インデントが変になる
//ちなみに言えば，インデントの前の文末に ~ おいても変なことになる
//経験則でわかるだろうか？
## 色分け範囲を自前で指定 -> level 引数を使用
## image 関数なら breaks 引数で範囲の分割点のベクトルを与える
> range(volcano)
[1] 94 195
> x <- 10*1:nrow(volcano)
> y <- 10*1:ncol(volcano)
> filled.contour(x, y, volcano, color = terrain.colors,
level=c(90,110,130,150,170,190,210),
plot.title = title(main = "The Topography of Maunga Whau",
xlab = "Meters North", ylab = "Meters West"),
plot.axes = { axis(1, seq(100, 800, by = 100)),
axis(2, seq(100, 600, by = 100)) },
key.title = title(main="Height?n(meters)"),
key.axes = axis(4, seq(90, 190, by = 10)))
> mtext(paste("filled.contour(.) from", R.version.string),
side = 1, line = 4, adj = 1, cex = .66)
#ref(Volcano1.png)
//ここにも置くな
## 色分け範囲を自前で指定 -> level 引数を使用
## 範囲色を自前で指定 -> col 引数で色名文字列ベクトルを指定
## 但し色の自前の指定は一人よがりになり勝ちですから、お勧めできません。
## すでに用意されている視覚的に慎重にデザインされたものを使うことをお勧めします。
> x <- 10*1:nrow(volcano)
> y <- 10*1:ncol(volcano)
> filled.contour(x, y, volcano,
level=c(90,110,130,150,170,190,210),
col = c("red", "blue", "yellow", "black", "cyan", "green"),
plot.title = title(main = "The Topography of Maunga Whau",
xlab = "Meters North", ylab = "Meters West"),
plot.axes = { axis(1, seq(100, 800, by = 100)),
axis(2, seq(100, 600, by = 100)) },
key.title = title(main="Height?n(meters)"),
key.axes = axis(4, seq(90, 190, by = 10)))
> mtext(paste("filled.contour(.) from", R.version.string),
side = 1, line = 4, adj = 1, cex = .66)
#ref(Volcano2.png)
//ここも同じ
## なおこの例の引数 color = terrain.colors の意味は、col 引数に与える色名ベクトルを
## カラーパレット関数 terrain.colors を用いて col = terrain.colors(10) のように
## しろという意味です。結果は以下のように色名を RGB 表記で与えたものになります。
> terrain.colors(10)
[1] "#00A600" "#2DB600" "#63C600" "#A0D600" "#E6E600" "#E8C32E" "#EBB25E"
[8] "#EDB48E" "#F0C9C0" "#F2F2F2"
-カラーバーに関するご指摘、すみませんでした。作成したかったのは、ご回答につけて下さった図(特に右横のもの)、まさにそのものです。自前の色付けの件は、確かにその通りだと思いましたので、そのようにします。有り難うございました。 -- [[ゆき]] &new{2005-01-04 (火) 10:32:48};
**ExcelのRight、Leftみたいな機能 [#n2644845]
>[[Akira]] (2004-12-27 (月) 10:38:26)~
~
Excelには文字列の右端、左端を得る関数があります。~
"abcd"なら、right("abcd", 1) -> "d" となります。~
Rには
strsplit("abcd", split="")[[1]][4]
とか~
文字数が不明な場合は
lapply(strsplit("abcd", split=""), rev)[[1]][1]
で文字を得られますが、~
文字ベクトルの全ての要素に適用して、listでなく右端の文字だけのベクトルを得る方法はあるのでしょうか？
chara.a <- c("abcd", "bcda", "cdab")
chara.list <- lapply(strsplit(chara.a, split=""), rev)
chara.b <- numeric(0)
for(i in 1:length(chara.list)){
chara.b[i] <- chara.list[[i]][1]
}
でベクトルを作っています。~
//
-substr という関数を調べるとよいでしょう。 -- [[青木繁伸]] &new{2004-12-27 (月) 11:03:01};
-ありがとうございます。
chara.list <- substr(chara.a, start=1, stop=1)
としました。 -- [[Akira]] &new{2004-12-27 (月) 12:01:39};
-ncharを使えばrevしなくても文字列末尾の参照はできますよ -- &new{2004-12-27 (月) 14:29:13};
-こういうことでしょうか?
chara.list <- substr(chara.a, start=nchar(chara.a), stop=nchar(chara.a))
あっていますか？ -- [[Akira]] &new{2004-12-27 (月) 14:56:28};
-テストデータで期待通りの結果が得られれば，あっているのでしょう。 -- [[青木繁伸]] &new{2004-12-27 (月) 21:07:46};
-回答以前に末端の文字ベクトルを求めてそれで何をしたいのか？が気になるケースですね。最終目的次第では更に気のきいたやりかたがありそうな雰囲気が。 -- &new{2004-12-28 (火) 18:30:51};
-全然凄いことではありません。~rownameに群分け情報があるので、取り出したかったのです。Excelで抽出していたのですが、バッチ処理にすれば手間が省けると思った次第です。~
> test
XorY A
1000X X TRUE
1001X X FALSE
1000Y Y FALSE
1001Y Y TRUE
というdata.frameなのでXorYのところにrownamesの末尾をもってきました。-- [[Akira]] &new{2004-12-28 (火) 19:04:34};
-そのようなデータが，機械的に作られてそれを解析する必要があるならやむを得ない操作になるのでしょうが，そうでないならデータ設計の失敗でしょうね。~
そうでないとしても，統計解析の前に適切なフィルターを噛ましてやれば良いだけでしょう。~
フィルターは R で構成する必要はなくて，より適切なスクリプト（perl とかなんとか）でも言い訳です。~
そのフィルター部分を Excel に担当させようとしたのが一つの敗因では？~
R の得意な分野とそうでない分野があるわけですから，R でない部分を何でやるかによって R の負担も変わるでしょうが。~
極端なことをいえば， Excel で全部やっても良かったんですから。~
処理の全体がわかりませんので，あれこれ言うのは控えておきましょう（十分よけいなことを言ってしまったか） -- [[青木繁伸]] &new{2004-12-28 (火) 21:19:29};
**R的なapply, lapply, sapplyの活用 [#xfd2395f]
>[[Akira]] (2004-12-27 (月) 10:24:04)~
~
「意味がある時は(ない時も)常にベクトル化せよ」を実践すべく頑張ってますが、困っています。次の処理はforを使わずにできますでしょうか?~
~
a.listとb.listがあります。a.listはnumericデータのdata.frameをlistでまとめ、b.listはa.listのdata.frameを2つに分類するためのnumericベクトル組をlistでまとめています。~
つまり、
length(a.list)=3
dim(a.list[[1]])
[1]20 200（これは[[2]]、[[3]]も同じ）
length(b.list)=2
b.list[[1]]
$g1
[1]1 2 3 4 5 6 7 8 9 10 11 12
$g2
[1]13 14 15 16 17 18 19 20
b.list[[2]]
$g1
[1]1 2 3 4 5 6 7 8 9 10
$g2
[2] 11 12 13 14 15 16 17 18 19 20
今このデータについて、b.listにある2つのカテゴリを使って、a.listの3つのdata.frameをそれぞれ2群検定したいと思っています。~
今は、
for(i in 1:length(b.list)){
for(j in 1:length(a.list){
g1 <- a.list[[j]][b.list[i], ]
g2 <- a.list[[j]][b.list[i], ]
ttest.p <- numeric(0)
for(k in 1:dim(a.list[[1]])[2]){
ttest.p[k] <- t.test(g1, g2, var.equal=T, alternative="two.sided",
na.action=na.omit)$p.value
}}}
の様なことをしてしまっています。~
せめて、apply(g1, 2, t.test, g2)みたいな方法でg1とg2のcol毎にt.testをしてP値だけを得る方法はあるのでしょうか？~
//
-複雑なデータ構造のようで，よく把握できませんが，期待通り動いているならいいのではないでしょうか。 -- [[青木繁伸]] &new{2004-12-27 (月) 15:49:00};
-こうした複雑なケースはミニチュアでも良いですから実際に自分の意図通りに動く例を添えないと回答する気になりませんよ。もし私の解釈が正しければ次のようになるでしょうか。-- [[QDU]] &new{2004-12-28 (火) 11:55:56};
# おそらく質問者がやったこと
> set.seed(31415)
> a1 <- as.data.frame(matrix(rnorm(4000),20,200))
> a2 <- as.data.frame(matrix(rnorm(4000),20,200))
> a3 <- as.data.frame(matrix(rnorm(4000),20,200))
> a.list <- list(a1,a2,a3)
> b.list <- list(list(g1=1:12,g2=13:20),list(g1=1:10,g2=11:20))
> an <- length(a.list); bn <- length(b.list)
> ttest.p <- matrix(vector("list", an*bn), bn, an) # 結果を入れるリストの行列
> for(i in 1:bn){
for(j in 1:an){
G1 <- a.list[[j]][b.list[[i]]$g1, ]
G2 <- a.list[[j]][b.list[[i]]$g2, ]
temp <- numeric(0)
for(k in 1:200){
temp[k] <- t.test(G1[,k], G2[,k], var.equal=TRUE,
alternative="two.sided",
na.action=na.omit)$p.value
}
ttest.p[i,j] <- list(temp)
}}
> str(ttest.p)
List of 6
$ : num [1:200] 0.869 0.155 0.458 0.768 0.362 ...
$ : num [1:200] 0.907 0.176 0.799 0.296 0.861 ...
$ : num [1:200] 0.226 0.691 0.465 0.704 0.363 ...
$ : num [1:200] 0.166 0.104 0.918 0.419 0.126 ...
$ : num [1:200] 0.7886 0.3279 0.0029 0.5534 0.4714 ...
$ : num [1:200] 0.5945 0.4196 0.0560 0.0134 0.4812 ...
- attr(*, "dim")= int [1:2] 2 3
# 内側のループを sapply で処理
> for(i in 1:bn){
for(j in 1:an){
G1 <- a.list[[j]][b.list[[i]]$g1, ]
G2 <- a.list[[j]][b.list[[i]]$g2, ]
ttest.p[i,j] <- list(sapply(1:200,
FUN=function(k) t.test(G1[,k], G2[,k],
var.equal=TRUE,
alternative="two.sided",
na.action=na.omit)$p.value))
}}
> str(ttest.p)
List of 6
$ : num [1:200] 0.869 0.155 0.458 0.768 0.362 ...
$ : num [1:200] 0.907 0.176 0.799 0.296 0.861 ...
$ : num [1:200] 0.226 0.691 0.465 0.704 0.363 ...
$ : num [1:200] 0.166 0.104 0.918 0.419 0.126 ...
$ : num [1:200] 0.7886 0.3279 0.0029 0.5534 0.4714 ...
$ : num [1:200] 0.5945 0.4196 0.0560 0.0134 0.4812 ...
- attr(*, "dim")= int [1:2] 2 3
外側の二重ループを消すこともできるのでしょうが、やりかたを考えるよりは、もっと大事なことが他にあるはず。(どうしても知りたければ r-help に質問するときっと Gabor G. 氏があっと驚く一行コードを紹介してくれるかも、--- 彼は一体どういう素性の人なのかしら？)
-ありがとうございます。sapplyの使い方の勉強になりました。ご指摘の通りもっと大事なことがあります。反省します。~
forが続くとメモリを消費するのかな?と思ったので、なるべく避けたいと思っていました。forの使用とメモリの消費は関係ないのでしょうか? -- [[Akira]] &new{2004-12-28 (火) 19:13:29};
-実は apply 関数は for loop を隠すだけで実際の速度向上は期待できません。メモリー消費も同じでは（未確認）。コードが簡潔になるというメリットは大きいですが。 -- [[QDU]] &new{2004-12-28 (火) 20:52:53};
-コードを簡潔にするために，あれこれ時間を掛けるのが良いことかどうかは議論の余地があるのではないでしょうか。~
簡潔なコードをいつでも書けるように修行するのは良いのでしょうが。~
取りあえずの目的を達するために書くスクリプトは，技巧を凝らすのではなく，確実を期せば良いと思います。何が確実かと言うことは，各人の技量に依存することなので，一概に論じることはできないでしょう。~
うまい書き方や，効率的なプログラムにこだわれば，結果としてはむだな時間を費やすことにもなりかねません。~
プログラミングの経験は十分あるが，R の経験はそれほどでもないという場合には，for ループでもなんでも，それまでに身につけているプログラミングテクニックを最大限利用すれば良いと思いますが。~
要するに，トータルのコストの問題ですね。R に習熟すれば，R の特性を生かして，より一層の時間短縮ができるとは思いますが。また，そうなるよう，修行すべきでしょうが。~
for だろうがなんだろうが，R 的に汚いやり方であっても，その時点でユーザが最小の努力で最大の効果が得られるような方法をとるべきだろうというのが，私の考えです。繰り返しますが，R 的な技量が高まれば，取るべき方法は徐々に変わってくるでしょう。~
メモリーを浪費するといっても，所詮個人のマイクロコンピュータでシングルユーザでの話ですし，メモリーを浪費して計算時間が数十倍になるわけでもないし（なるかもしれんが），計算時間が数十倍になっても数時間の問題になるのでなければ，問題ないでしょう？コンピュータに計算させておいて，その間に貴方の好きなことをやっていればいいのです。読書でも，スポーツでも，食事でも。。。~
それでもなお，もっと早く分析したいなら，それはきっと R のような言語に頼らない方がいいのです。C ででもプログラムを書くべきものなんでしょう。~
このような議論を受け入れてくれる人は少ないようなんですが，なんででしょうね。不思議なことです。不思議だと思う私の考え方の方が不思議なんでしょうね。 -- [[青木繁伸]] &new{2004-12-28 (火) 21:36:48};
- いちばん理想的なのは「他人が読んでも読みやすく、しかも計算速度が速いプログラムソース」なのでしょうが、実際にはこの 2 つを兼ね備えるのは難しいですよねぇ。~
「比較的簡単にプログラムが出来、それでいてそこそこ計算が速い R を使っているのだから、速度よりも読みやすさに力を入れた方が良いのでは？」というのが私のスタンスです。読みやすい・分かりやすいプログラムは誤解を招きにくいので、青木先生の「確実を期す」ことに繋がると思います（いくら速いからといっても、人様が書いた C 言語のプログラムは、プログラミングが苦手な人（例えば私）は読むのにメチャメチャ苦労しますし…）。「計算速度のことは完全に無視しろ！」などとは申しませんが、まずは間違いのない確実な計算をするようなプログラムを作るのが目標でしょう。そして、出来れば作っている際に読みやすい（誤解の少ない）プログラムを書くように心がけることが出来たら理想的ですよね。で、COLOR(red){もし}人様にプログラムを見せるのであれば、出来上がった（確実な計算をする）プログラムを見て、「この部分はこうすれば読みやすくなるんじゃないかな？」と編集すればよいのではないでしょうか？兎にも角にも、まずは「確実を期す」ことが一番です。計算速度を考えるのは、プログラムを読みやすいものに編集し終わった後にどうしても必要なのであれば…ってのが私のスタンスです。せっかく「最小の努力で最大の効果が得られる」R を使っているのですから、その恩恵を蒙りましょうよ。~
ここまで書いて思ったのですが、私のようなど素人が意見を述べるのもトンチンカンな気がしますねぇ＾＾；。まぁど素人の戯言だと聞き流して頂ければ幸いです。 -- [[舟尾]] &new{2004-12-29 (水) 12:52:09};
**プロットの3色の色分け [#r39ed7ed]
>[[Mari]] (2004-12-21 (火) 18:09:28)~
~
何度も本当にすみません。。。~
プロット上のスポットを値ごとに3色に分けたいのですが、以下のスクリプトではだめなようでした。
> plot(A2,M2,col = (if(Type==1){"gray"}else if(Type==-20000) {"blue"} else "black"))
Warning message:
the condition has length > 1 and only the first element will be used in: if (Type == 1) {
col をif文の中に入れてみたりしたのですが、plotではifelseしか使えないのでしょうか？~
//
-一歩ずつ前進しているようですね。以下が参考になるでしょうか。途中で y <- x$C として省エネしています。 -- &new{2004-12-21 (火) 19:00:58};
> x <- data.frame(A=rnorm(5), B=runif(5), C=c(1,2,2,3,1))
> plot(x$A, x$B, col= ifelse((y <- x$C)==1,"red", ifelse(y==2, "green", "blue")))
なお、試みられたコードが失敗する理由も、初心者の典型的な「躓きの石」のひとつです。この場合 Type はベクトルですが、構文 if(Type==1) はスカラ変数 Type しか受け付けません。ベクトルを与えると if(Type[1]==1) とされ、注意がでます。一方 ifelse (Type==1 , , ) はベクトル変数 Type を受け入れます。
> Type =c(1, -20000, 1, 0, -20000, 1)
> y = ifelse(Type==1, "gray", ifelse(Type==-20000, "blue", "black"))
> y
[1] "gray" "blue" "gray" "black" "blue" "gray"
-本当にありがとうございました!!! 今日の一連の流れをfunction関数で自動でできるようにするところまでできました!!! 今日は大変な進歩の一日でした。ありがとうございます。 -- [[Mari]] &new{2004-12-21 (火) 19:36:23};
-発展問題。関連して w <- foo(x,y,z) とすると x[a]==y[b] なら w[a]=z[b] となる関数は？もしかして既にそんざい？こうした関数があれば上の問題は col=foo(Type, c(1,-20000,0), c("gray","blue","black")) とすれば良いことになりますが。例えばこんなの ==> [[多重比較マッチング関数]] -- &new{2004-12-21 (火) 19:53:32};
-よくよく考えたらもっと簡単な解がありました。-- &new{2004-12-23 (木) 07:52:44};
col = c("gray","blue","black")[ match( Type, c(1,-20000,0) ) ]
**プロットの色分け [#a83c1f51]
>[[Mari]] (2004-12-21 (火) 16:54:14)~
~
いつもお世話になっております。今回はプロットを作成して、色分けを行いたいと思っております。~
~
インポートしたデータフレーム（dat2）のA2,M2 というカラムの値を使って、プロットを書くのですが、それぞれにControlType というものが数値として入っており、これごとに色を変えたいと思っております。
> unique(dat2$ControlType)
[1] 1 -20000 0
とりあえず、2色でifelseをプロット関数に入れて行ってみましたが、エラーがでました。
> plot(A2,M2, col = ifelse(dat2$ControlType=0,"gray","blue"))
Error: syntax error
別の関数に変換するほうがよいのかなと思い、下記のようにしましたが、~
やはりエラーがでました。
> Type <- dat2$ControlType
> unique(Type)
[1] 1 -20000 0
> plot(A2,M2, col = ifelse(Type=0,"gray","blue"))
Error in ifelse(Type = 0, "gray", "blue") :
unused argument(s) (Type ...)
某HPの以下のスクリプトを参考にしたのですが、何が悪かったのでしょうか？
> plot(x, y, col = ifelse(y>0.5, "red", "blue"))
//
-つぎの二つの例を見比べて下さい。初心者が良くひっかかる落し穴です。教訓：「一番最後までのこるバグは、一番つまらないバグである」 -- &new{2004-12-21 (火) 17:04:05};
> x=1; Y <- ifelse(x==1, 2,3)
> Y
[1] 2
> x=1; Y <- ifelse(x=1, 2,3)
Error in ifelse(x = 1, 2, 3) : unused argument(s) (x ...)
-ありがとうございました!!! たしかに・・・ -- [[Mari]] &new{2004-12-21 (火) 17:35:02};
**read.table で最初の9行をスキップしたい [#t28ecfd0]
>[[Mari]] (2004-12-21 (火) 11:09:47)~
~
エクセルで作成したデータのインポートを行いたいのですが、1行目から9行目までがコメントになっており、10行目からのインポートを行うために、~
read.table のオプションのskip を使ったのですが、下記のようなエラーが返ってきました。~
~
具体的には、下記のように行いました。どうやら、1行目のデータ数と、10行目以降のデータ数の違いによるもののようですが、どのようにすればよいでしょうか？~
~
データは1行目から9行目まではランダムなコメントが書いてあり、10行目がヘッダー、11行目からが実データになります。 11行目のデータは87列ほどのデータですが、1行目には50列ほどしかコメントで使用されていません。~
~
アドバイスいただけますと助かります。よろしくお願いいたします。
> dat<- read.table("Test.txt",skip=9)
Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
line 1 did not have 87 elements
//
-おそらく ?read.table の次の箇所が参考になるのでは。 -- &new{2004-12-21 (火) 11:37:23};
The number of data columns is determined by looking at the first
five lines of input (or the whole file if it has less than five
lines), or from the length of 'col.names' if it is specified and
is longer. This could conceivably be wrong if 'fill' or
'blank.lines.skip' are true, so specify 'col.names' if necessary.
-次のようなテストファイルを読み込んでみると問題は起きませんが、1 行目（ヘッダ部分が）完全でしょうか？ -- &new{2004-12-21 (火) 12:02:51};
comment 1 a b c d e f g h
comment 2
comment 3
comment 4
comment 5
A B C D E
1 2 3 4 5
6 7 8 9 10
11 12 13 14 15
> x <- read.table("test.txt", skip=5, header=T)
> x
A B C D E
1 1 2 3 4 5
2 6 7 8 9 10
3 11 12 13 14 15
-ありがとうございます。 コメント部分を削除しても、read.table で読むことはできず、read.csv で読んでやるとうまく行くよう -- [[Mari]] &new{2004-12-21 (火) 12:44:55};
-すみません、切れました。コメントをはずして、read.csv で読むとうまくいきました。どうやら、本体部分になにか欠損値のようなものがあるのかと思います。分かりましたら、こちらにUpdateしたいと思います。 -- [[Mari]] &new{2004-12-21 (火) 13:02:12};
-以下のように行うと、うまく行きました。本体部分のタブ区切りに何か問題があったのかもしれません。Warning が出ているのが気になりますが、names でカラムの内容を確認してみるとちゃんと読めていました。> tab5 <- read.table("test5.txt", skip = 9, header = TRUE , sep = "?t") -- [[Mar]] &new{2004-12-21 (火) 13:52:14};
-コメント部分を#でコメントアウトしてもexcel上では#が"#（"が追加されている）となります。（以前これをやってしまいました） -- [[Akira]] &new{2004-12-27 (月) 10:57:02};
**条件にあうものだけに、値を変更したい。 [#lfbf8f67]
>[[cricket]] (2004-12-18 (土) 00:14:58)~
~
データフレームで、あるラベルに対応する値だけを変更するにはどうすればいいのでしょうか。たとえば、data$treatmentが"control"の場合だけに、data$observationに10を加えたいのです。初歩的なことだと思いますが、よろしくお願いします。~
//
-色々やりかたはあるでしょうが例えば素朴に -- &new{2004-12-18 (土) 00:32:57};
> x <- data.frame(A=1:4, B=c("a","b","b","c"))
> x
A B
1 1 a
2 2 b
3 3 b
4 4 c
> x$A[x$B=="b"] <- x$A[x$B=="b"] + 10
> x
A B
1 1 a
2 12 b
3 13 b
4 4 c
> x$A[x$B=="b"] <- 0
> x
A B
1 1 a
2 0 b
3 0 b
4 4 c
-上の方と同じく，いろいろやり方はあるでしょうが
> data <- data.frame(
treatment=c("control", "treatment", "treatment", "control", "control"),
observation=c(1,3,2,1,4))
> data
treatment observation
1 control 1
2 treatment 3
3 treatment 2
4 control 1
5 control 4
> data$observation <- data$observation+ifelse(data$treatment=="control", 10, 0)
> data
treatment observation
1 control 11
2 treatment 3
3 treatment 2
4 control 11
5 control 14
など -- [[青木繁伸]] &new{2004-12-18 (土) 00:47:45};
-なるほど。確かにいわれてみればたしかに。難しく考えすぎていたようです。ありがとうございました。 -- [[cricket]] &new{2004-12-18 (土) 10:59:20};
-こんなのもあり -- &new{2004-12-18 (土) 11:34:33};
> x=data.frame(A=1:4, B=c("a", "b","b","c"))
> x$B == "b"
[1] FALSE TRUE TRUE FALSE
> x$A <- x$A + 10 * (x$B == "b")
> x
A B
1 1 a
2 12 b
3 13 b
4 4 c
-発展問題。逆に文字列部分を変更したい時は？データフレームでは文字列は因子とされてしまうので、単純に文字列として置き換えるとエラーになるようだが？ (こんなこと普通はしない？) -- &new{2004-12-18 (土) 11:42:08};
> x=data.frame(A=1:4, B=c("a", "b","b","c"))
> x$B[x$A==2] <- "c" # すでにある文字列には問題無く置き換えられるようだ
> x
A B
1 1 a
2 2 c
3 3 b
4 4 c
> str(x$B) # "a","b","c" は内部的には数値 1,2,3 と表現されている(?)
Factor w/ 3 levels "a","b","c": 1 2 2 3
> x$B[x$A==2] <- "d" # 最初に無い文字列に置き換えるとエラーで NA とされる
Warning message:
invalid factor level, NAs generated in: "[<-.factor"(`*tmp*`, x$A == 2, value = "d")
> x
A B
1 1 a
2 2 <NA>
3 3 b
4 4 c
-x$A <- x$A + 10 * (x$B == "b") のような使い方もできるんですね！目からウロコです。 -- [[cricket]] &new{2004-12-18 (土) 14:24:24};
-存在しないファクターを加える（強引な)方法。元の質問の趣旨から恐らく脱線気味でごめんなさい。
> x = data.frame(A=1:4, B=c("a","b","b","c"))
> x = rbind(x[1,], data.frame(A=x$A[2], B=c("d")), x[3:4,])
> x
A B
1 1 a
11 2 d # <- 確かに変わったが、行ラベルが何故かおかしくなった
3 3 b
4 4 c
> x[2,] # しかし添字操作では問題無し
A B
11 2 d
**単位根検定についてｐａｒｔ２ [#mf19ce3a]
>[[nakanaka]] (2004-12-17 (金) 19:31:57)~
~
他の計量経済学の掲示板にも質問させていただきましたが,なかなかレスがないようなのでここでも質問させてください.現在,単位根検定のADF検定をパッケージ”fSeries”で行っているのですが,トレンド項,ドリフト項の取捨に困ってます.~
トレンド項,ドリフト項のｔ（τ）値を見れば,有意かどうか判断できると思うのですが,”R”では,それらの推計式が示されない（つまり,トレンド項,ドリフト項のｔ値が見れない）ので推計モデルにトレンド項やドリフト項を含んでいいのかどうかがわかりません.このような場合どういう風に推計式を選択すればよいのでしょうか？fSeriesのマニュアルを見てもこの推計式の選択についての検定は書かれてなかったように思われます.推計式のトレンド項,ドリフト項のｔ値を見る方法,もしくは,他の解決策があれば教えていただけると幸いです.~
よろしくお願いします.~
//
-いろんなところに匿名でマルチポストするのは感心しない。すでにそっちの匿名掲示板で正しい回答をもらっている。[[マルチポストされた質問と回答:http://www.ichigobbs.net/cgi/15bbs/economy/0703/303-307]] -- &new{2004-12-18 (土) 04:46:45};
-[[R-sig-finance:https://stat.ethz.ch/pipermail/r-sig-finance/2004q4/000229.html]]のadfTestのコード例を書き換えて、lmの結果を取り出すようにすればOKです。よい勉強になると思うので、チャレンジしてみてください。 -- [[矢野]] &new{2004-12-20 (月) 08:50:38};
**テーブル内のリスト一覧 [#bf96ebd4]
>[[Mari]] (2004-12-17 (金) 16:48:30)~
~
テーブルを読み込んだときに、そのデータのリストの一覧を入手するような関数ってあるのでしょうか？~
例えば、Test というデータのカラムに~
Column
-------
type1
type3
type2
type1
type3
~
というような項目があり、このデータが何万件にも及ぶ場合、~
何種類のtype があるか、という情報を入手したいのですが。。~
~
show(Test$Column)~
ですと、ただ一覧がでてきてしまいます。~
~
調べ方も甘いのかもしれませんが、ご教授いただけますと~
幸いです。~
~
よろしくお願いいたします。~
//
-因子なら,levels(Test$Column),数値ならunique(Test$Column) -- &new{2004-12-17 (金) 17:20:08};
-できました！ありがとうございます -- [[Mari]] &new{2004-12-17 (金) 17:33:58};
-参考。unique は文字列ベクトルでも使えます。 -- &new{2004-12-17 (金) 22:56:35};
**正準相関係数に関して [#g73d5f1b]
>[[fuji]] (2004-12-14 (火) 13:51:37)~
~
正準相関係数に関して質問です．~
cancor()関数では正準相関係数は出力されますが，重み係数は出力されないのでしょうか？~
//
-[[青木先生作 cancor 関数:http://aoki2.si.gunma-u.ac.jp/R/cancor.html]] とその解説をじっくり読めばいかがでしょう。 -- &new{2004-12-14 (火) 23:16:26};
-参考図書の解答が間違っているのかもしれません．他の例題ではうまくいきました． -- [[fuji]] &new{2004-12-15 (水) 07:56:00};
**VECMのパラメータについて [#bd8fcf13]
>[[Yo]] (2004-12-13 (月) 17:37:37)~
~
VECMのパラメータ推計を，パッケージ"urca"の"ca.jo"で行っているのですが，出力された値のうちどれがどの式（長期均衡式・EC式・階差式）に対応したパラメータなのかわからずに困っています．~
"ca.jo-class"の説明を見て"PI""GAMMA"が~
~
EC式：EC=Y(t)-(a0+a1*X(t-1)+a2*Z(t-1)+..)~
~
のパラメータだろうと思っているのですが，~
~
VECM:?Y(t)=b0+b1*EC(t-1)+b2*Y(t-1)+b3*X(t-2)+..~
~
のパラメータがどれにあたるのかわかりません．~
とても初歩的な質問で申し訳ありませんが，教えていただけないでしょうか？~
よろしくお願いいたします．~
//
-VECMってなんですか？いきなり独りよがりに式を並べられても誰も分からないのでは？ -- &new{2004-12-14 (火) 18:55:43};
-VECMとは多変量誤差モデルのことで、VARに共和分分析が組みこまれたものです。式は一般的なものだと思い書いたのですが・・説明不足ですみません。計算はできるのですが、出力される値の意味がわからず困っています。ご存知の方がいらっしゃれば教えていただけないでしょうか？ -- [[Yo]] &new{2004-12-15 (水) 13:11:27};
-一般的なお話として。~
使ったプログラムが合っているか（R なので，たいていは大丈夫でしょうが），使い方を間違えていないか（これはあり得る），をチェックするために，答えのわかっているデータを食わせてみるということはよくやることでしょうね。そうすれば，出力されるどの数値がどれに対応するかはすぐにわかるのではないでしょうか。~
もっとも，この上の（正準相関分析 cancorについての）質問のように，教科書に書いてある答えが間違っているということもあるし，違っているように見えても仮定の違いだとかはありますが。だとしても，少なくとも違う答えが出たら，何でだろうと調べますよね。複数のソースからの例題を食わせてみるというのが最も慎重な振る舞いでしょう。 -- [[青木繁伸]] &new{2004-12-15 (水) 13:37:09};
-ありがとうございます。さっそく試してみます。 -- [[Yo]] &new{2004-12-15 (水) 16:05:41};
-VECM (Vector Error Correction Model)には詳しくありませんが、実際にVECMでの結果を出すのはca.joではなく、cajoolsではないでしょうか？使い方はhelp(cajools)で見てください。 -- [[矢野]] &new{2004-12-20 (月) 08:47:58};
**カテゴリーごとのプロット [#x299ddcd]
>[[shiGe]] (2004-12-13 (月) 07:28:32)~
~
アメリカ50州のデータに関して，それぞれの州の２変量OLSフィットの図（つまり全州において同じ２変量間の散布図及び回帰直線）を１つの図にまとめて描きたい（50個の図が並ぶ）のですが，plot()に関する情報を検索しても，分かりませんでした（というよりどのように検索すべきかもいまいち分からないのですが）。~
アメリカ50州のデータに関して，それぞれの州の２変量OLSフィットの図（つまり全州において同じ２変量間の散布図及び回帰直線）を１つの図にまとめて描きたい（50個の図が並ぶ）のですが，plot()に関する情報を検索しても，分かりませんでした（というよりどのように検索すべきかもいまひとつ分からないのですが）。~
~
仕方ないので当てずっぽうで
plot(Y ~ X, data=hoge, for=STATE)
や，for=STATEの部分をby=STATEなどに代えてやっても当然のようにエラーになります。ご存知の方がいらっしゃいましたらご教示いただけますでしょうか。~
//
-pairs(USJudgeRatings) なんてどうでしょう？関数 pairs() については help(pairs) でヘルプを見ることが出来ます． -- [[舟尾]] &new{2004-12-13 (月) 08:05:35};
-Tips 集の関数の項 [[pairs用 回帰直線つきpanel]] がずばりでは？しかし50州はひとつに並べるにはあまりに多過ぎるような。 -- &new{2004-12-13 (月) 08:27:42};
-あ！そのような記事がありましたか！ -- &new{2004-12-13 (月) 18:00:33};
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
