初心者のための R および RjpWiki に関する質問コーナー
注意:このコーナーは重量オーバーで表示に時間がかかり過ぎるようになりました。今後は過去記事へのコメントを除き、新設の Q&A (初級者コース) コーナーを利用下さい。
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)色々試してみましたが、、やり方がわかりません。
ご教授の程よろしくお願いします。
伊藤 (2005-04-17 (日) 15:59:59)
お願いします。
私は生物学の研究をしており、2種類の寄主植物上が発育期間とその後の産卵数との関係に及ぼす影響について解析しています。具体的には、横軸に発育日数をとり縦軸に産卵数をとります。この関係をそれぞれの寄主植物別で解析すると、両方とも二次関数へのあてはまりがよさそうでした。そこで次に、寄主植物の種類が二次関数のパラメータに影響しているかを評価するために、以下のモデルをたてました。
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も試みてみましたが、ランダム因子がないので不可なようです。代案としてどういうコマンドが考えられるか、ご教唆いただければと思います。よろしくお願いします。
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ここしばらくは,細部にわたった解答をして,質問者をスポイルすることもいとわずと言う解答スタイルを取っています。 -- 青木繁伸 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というダミー変数を作る必要があると理解しました。ありがとうございました。 -- 伊藤 2005-04-17 (日) 21:27:00
thorman (2005-04-16 (土) 21:51:53)
最近Rで統計学を始めたのですが(以前はExcelを使っていました)、Rでデータのソートってできるのでしょうか。
例えば、x<-c(1,3,4,4,6,5,2)とxにデータを入れた後、このデータを降順(あるいは昇順)にならべたりするのはできるのでしょうか。
まだ、Rについてよく分からないのですが、C言語でデータをソートするようなプログラムを作るように、じぶんでfor構文でソートするためのプログラムを作らなければダメなのでしょうか?
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このソートアルゴリズムは,優れたものではないですが,わかりやすいもののひとつでしょうね。 -- 青木繁伸 2005-04-17 (日) 20:39:56
ごろう (2005-04-14 (木) 13:06:21)
クラスター解析で樹形図を作成しているのですが、項目の日本語を入れたまま作図しようとするとエラーになるので、項目を削除して作図をしているのですが、そうすると図の下には数字しか表示されません。R上で図の下に表示される数字と項目を対応させる方法はあるのでしょうか?
自分の作図方法を下に示します。 よろしくお願いします。
使用OS Microsoft windows XP
R Version 2.0.0library(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)
beginner for R (2005-04-13 (水) 16:31:53)
いくつのデータをバッチ処理によってプロット出力したいのですが、そのときに入力ファイル名をプロットのmain titleに指定したいのですが、どのようにやればいいでしょうか?
どうぞよろしくお願いします。
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 に渡してやることもできますよね。 -- 青木繁伸 2005-04-13 (水) 16:56:53
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"が入っているのですが。。。
> hoge.test <- function(l="hoge") { + list("fuga"=c("PV.IV","IV"),l=l) + } > hoge.test("fugafuga") $fuga [1] "PV.IV" "IV" $l [1] "fugafuga"ちゃんとしたプログラムの一部としては意味のあるパーツなんでしょうかね? -- 青木繁伸 2005-04-07 (木) 15:04:45
> 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℃ 2005-04-07 (木) 16:04:46
> 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
初心者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 日本語版を使っています。
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"あまり,有用な関数とも思えません。 -- 青木繁伸 2005-04-05 (火) 19:10:27
fmt <- sprintf("%%se%%+0%ii", as.integer(exponent))としたら、希望通りになりました。スペルミス、失礼いたしました。青木先生のサイトは よく 参考にさせていただいています。 -- 初心者11242号 2005-04-05 (火) 23:22:06
初心者 (2005-03-31 (木) 01:48:22)
シュミレーションをしていたら以下のようなエラーが出たのですが、Error: subscript out of boundsってなんでしょうか?
> y <- matrix(1:4, 2,2) > y [,1] [,2] [1,] 1 3 [2,] 2 4 > y[3,2] Error: subscript out of bounds
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")
> x <- 1:100/50*pi > y <- cbind(sin(x), cos(x)) # 同じx座標に対する二つの関数値を合わせた行列 > matplot(x,y, pch=".") # 二つの関数グラフを同時にかく
plot(x,y) abline(lm(y~x))
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は7つの回帰式を持っており、散布図にはそれぞれが重ねて描画されると思っていたのですが、うまくいかないのです。 -- yuta 2005-04-01 (金) 15:04:01
> 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ということになったんですが,どうですか? -- 投稿ルール 2005-04-01 (金) 16:18:57
> 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号 2005-04-01 (金) 16:56:28
Rをインタプリタでなくバッチで実行すると何百倍も早くなったりすんでしょうか? -- 豊田 2005-03-29 (火) 08:21:06
超初心者 (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 をしようとすると上記エラーがでます。
IE上に直接入力するとうまくいくのに、理由がわかりません。
かなりの初心者ですが、だれかお教え願いますでしょうか????
環境:win2k
初心者 (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ユーザーです。
> 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
future user of Mac mini (2005-03-16 (水) 19:28:26)
R Mac 版はMac mini (メモリー256MB/HDD 40GB)でも、問題なく動作するのでしょうか?
さち (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"))) }困った事に、得られた線分が、単色(”赤”)になってしまいました。。。
どうしたら良いのでしょうか。
ご指導をお願いします。
> 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))ここまでは、前回と同じで、非線形回帰分析をさせた結果を表示させている部分です。
> 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")))}色分けさせるように記述したつもりでしたが、結果は、線分が単色(赤)で表示されてしまいました。
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])}同時刻に修正しているので,編集の衝突が起きていますね(^_^;)。 -- 青木繁伸 2005-03-16 (水) 17:31:36
> 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 の結果を使えば良かったのですね。
> newpar <- result$m$getPar(); a=newpar[1]; b=newpar[2]; c=newpar[3]をします。 頑張って勉強します。-- さち 2005-03-17 (木) 10:31:08
初心者 (2005-03-16 (水) 15:17:13)
Excelの表をread.csvで読み込もうとしたら以下のメッセージが出たのですが、どうしたらよいのでしょう?なお、Rは2.01、Mac OS X環境です。
incomplete final line found by readTableHeader on `Data.csv'
% 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 から,書き出しを行う。 -- 青木繁伸 2005-03-16 (水) 15:58:45
というものです。 -- 青木繁伸 2005-03-19 (土) 10:05:03NeoOffice/JはMac OS X用の機能豊富なオフィス・ソフト(ワープロ、表計算、プレゼンテーション、作図)です。OpenOffice.orgというオフィス・スイートに基づいたもので、数多くのMac OS Xネイティブの機能を提供し、MicrosoftTM Officeなど他の人気のオフィス・ソフトの書類の読み込み、編集、交換ができます。
GNU一般公衆利用許諾契約書(GPL)の下で無料のオープンソース・ソフトとして公開しています。NeoOffice/Jは充実した機能と日常的な用途に充分の安定性を備えています。
nama (2005-03-16 (水) 12:21:56)
すみません。教えていただきたい事があります。
散布図X軸を打点があるところだけにしたいのですがどうすればいい
のでしょうか? xaxt="n" で消去することは出来たのですが、その後
3点のみ表示させたいのです。よろしくお願いいたします。* * * ---+----+----+----+----+----+---- 40 60 90
x <- c(40, 60, 90) y <- c(10, 4, 19) plot(x, y, xaxt="n") axis(side=1, at=x)
おさる (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を使ってある程度のことはできるのですが,いまひとつ爽快感が得られません. -- おさる 2005-03-14 (月) 03:10:36
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-- 青木繁伸 2005-03-14 (月) 10:35:45
ピンバッジ (2005-03-09 (水) 18:02:31)
全く何も分からないので何を書いていいのかも分からないのですが、異なる餌を与えた昆虫の生存率の比較を、RのTukeyで統計しました。しかし、比率を比べる場合はカイ二乗検定後、角変換(アークサイン変換)をした後にTukeyをしなければならないことを知りましたが、この一連の操作をRでどのようにやるかわかりません。また、どういった形のデーターで比較するかもわかりません。
本当にわからないことばかりですので、お手数ですができる限り詳しく教えていただけると有 難く思います。足りない情報がありましたら教えてください。
当初昆虫匹数 生き残った昆虫匹数 生存率 餌A a1 a2 a2/a1 餌B b1 b2 b2/b1 餌C c1 c2 c2/c1みたいなものなんですか?
ぶー (2005-03-08 (火) 16:28:53)
こちらに書き込むのはまずいでしょうか?
Rを使うのにESSから使用しているんですが、ESS-5.2.5だと、source や read.table の補完機能が効きません。ESS-5.2.3だとOK。
elisp を見ていないのでなんなんですが、、、同じ現象で回避されている方がいらっしゃいましたらご指導願えないでしょうか?
あ、LInux 2.4.7 で使用しています。
石川誠 (2005-03-06 (日) 20:43:09)
お世話になります。
重回帰式を作成するとき、有効な変数を選択する方法を教えて下さい。
総当り法はあったのですが、変数が多いと時間がかかってしまいます。
summary(lm1 <- lm(Fertility ~ ., data = swiss)) slm1 <- step(lm1) summary(slm1) slm1$anova
さち (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のようなもので色づけができたら、とてもうれしいです。
どうかよろしくお願いします。
points(predict.c,col=ifelse(predict.c>=15.1,"blue",ifelse(predict.c>=10.1,"green","red")))ということなんでしょうか??(黄色というのは見にくいので緑にしましたが)
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")))これは、2次元のグラフの記されている各点を色分けしていますね。3次元版になったら、きっと必要になると思うので、参考にさせて頂きます。
filled.contour(1:length(predict.c),1:length(predict.c), matrix(predict.c,nrow=length(predict.c),ncol=length(predict.c)))とか. -- takahashi 2005-03-03 (木) 12:37:08
highvalley (2005-03-01 (火) 19:05:38)
Rを勉強中のhighvalleyと申します。
現在多項式による回帰分析(最小2乗法)を行いたいのですが、皆さんはどのように行っていますでしょうか?
私はlm()を使って以下のように行っています。例えば目的変数をy、説明変数をxとして2次の多項式で回帰分析をする場合は、
result <- lm( y ~ 1+x+I(x^2) )
これで目的を達する事は出来ているのですが、この方法だと次数が多くなると大変です。もう少しスマートな方法があるような気がして、質問を致しました。良い方法をご存知でしたら教えて下さい。よろしくお願いします。
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) }ヒデー,ラッパーだ -- 青木繁伸 2005-03-01 (火) 21:19:19
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.
しま (2005-02-20 (日) 15:29:58)
新しい関数を追加したいのですが,どのようにすればいいのでしょうか?
初歩的な質問で申し訳ないです..
ひょうひょう (2005-02-10 (木) 18:48:00)
グラフィックス参考実例集:ラティスグラフィックスで以下のような説明がありますが、具体的に教えてください。
「lattice 関数は lightgrey の背景色等幾つかの固有の既定値を採用しているので、必要に応じ par 関数で変更する。」
par(bg="white")
としても変更できませんでした。よろしくお願いします。
library(lattice) trellis.par.get() trellis.par.set(background = list(col="white")) show.settings() ?trellis.par.get質問の基本的な姿勢が,以下のサイトに書かれてあります.
ちょろちゃん(;_;) (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のスクリプトの書き方が
悪い以外のなにものでもないような気がいたしますが、
何かお気づきの点等ございましたら、アドバイスいただければ幸いです。
どうぞよろしくお願いいたします。
--max-vsize=N Set vector heap max to N bytes; --max-nsize=N Set max number of cons cells to Nどうでしょうね。vsize か nsize か,それともどっちも関係ないか。 -- 2005-02-04 (金) 15:05:07
bob3 (2005-01-31 (月) 00:03:57)
bob3と申します。
Rを使った対数線形モデル(log-linear models)による分割表(クロス集計表)の分析をしようとしています。
ところが2点ほど不明な部分があり、ご相談にあがりました。
ご相談したいのは以下の2点です。
1)ANOVAコーディングよる対数線形モデルのやり方
2)関数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コーディングで掲載されています。
まず、Rによる対数線形モデルの手順を調べ、以下のようにしてみました。(長くなるので出力は省略しています。)
# 既存の分割表を入力 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と、なっております。
よろしくお願いいたします。
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 2005-02-01 (火) 00:33:31
loglin(rocket,margin=list(c(1,2),c(1),c(2)),param=TRUE) loglm(freq ~ kyori * zure, data = rocket.df, family = poisson)$param-- なかの 2005-02-04 (金) 22:11:04
実は、私がANOVAコーディングで得たいと思っているのは「標準化係数(標準効果)」なのです。
これは重回帰分析における標準化偏回帰係数に相当するものですが、loglin や loglm で得られるのは「推定値(効果)」だけのようです。
標準化係数(標準効果) = 推定値(効果) ÷ 標準偏差(標準誤差) なので、標準偏差(標準誤差)が得られれば何とかなりそうなのですが……
なかなか難しいものですね。
……ん、待てよ。AICでモデルの選択をするんだったら、そのあとで変数やカテゴ リごとに検定をするのはおかしいのかな?
重回帰分析でAICを使って変数選択するときは、有意でない変数も含めたモデルを 採用することもありますね。
ということで、ちょっと分からなくなってきました。すみません。
参考にしている文献ですが、これも実は私のほうが教えていただきたいぐらい で、「質的情報の多変量解析」http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/index.htmlぐらいしかありません。
Rでの対数線型モデルについてということでは、主に以下の文書を参考にしています。
http://www.ci.tuwien.ac.at/~zeileis/teaching/Biostatistics03/examples4.pdf
http://www.stat.ohio-state.edu/~tjs/865/handout-3.pdf -- bob3 2005-02-07 (月) 00:14:15
竹内 (2005-01-30 (日) 09:06:30)
Rでメタアナリシスを行いたいのですが、どなたか参考になるもの(教科書、webなど)をご存知でしょうか?google、amazonなどで一応探してはみましたが、「これ」というものがみつかりませんでした。英語のものやS言語のものでも結構ですのでどなたかご存知のかた、よろしくお願いいたします。
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 の引数で与えるだけでよいことになるでしょう。
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 20forなどを使えば出来そうかとは思うのですが、プログラム経験も乏しいので良く分かりませんでした。
どなたか御教授よろしくお願いします。
また、参考URL or 文献なども教えていただけるとうれしいです。
> 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 のオブジェクトを作る } }
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をうまいこと分解していけばよいのかな?と考えていました。
listoftables[[12, 34]]等と添字操作する方が簡便なのでは。リスト listoftables は
dim(listoftables) <- c(51,51)とすれば行列として操作できます。 -- 2005-01-25 (火) 20:43:26
初心者 (2005-01-17 (月) 20:18:02)
どうも以前にgamma関数を上書き定義してしまったらしく、使えなくなってしまいました。一度Rを閉じても元に戻りません。初期定義に戻すにはどうすればいいでしょうか。
Jack (2005-01-16 (日) 17:00:47)
次のようなデータフレームがあります。
V1, V2は処理の種類(カテゴリー)を表し、1,2,3,4,5は計測値1〜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を使わない方法があるのでしょうか。 よろしくおねがいします。
> 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,] # 実際の行数に切り詰め
# テストデータ作成 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))これがどういうことかは,
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))いきなり,これが出てくるとわけわからんですね。 -- 青木繁伸 2005-01-16 (日) 21:33:44
おやじっち (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
の部分を参照したのですが、なにぶん素人なものでどこが間違っているのかよく分かりません。
どうか教えて下さい。
scan("/Users/foo/Desktop/bar/baz.dat")とすればいいのだが,毎回そんなの書くのいやだということなら,command+Dでワーキングディレクトリを bar ディレクトリにすれば
scan("baz.dat")だけで読める。 -- 青木繁伸 2005-01-13 (木) 21:37:45
青木繁伸 (2005-01-11 (火) 14:32:44)
Windows ユーザから以前,
http://aoki2.si.gunma-u.ac.jp/R/excel.html
の動きがおかしいという問い合わせがありました。
複数列をコピーしてもなぜか一行中の全部の数字を連結した文字列になる(タブで区切られない)ので,書かれているような動きをしないということでした。その人は1.9でやっているからでしょうかと言っておりましたが,今日たまたまWindows版のRをさわっていて2.0.1でもちゃんと期待通り動かないことが分かりました。
なんででしょうか。
エディタなどにペーストすると列の間にはちゃんとタブコードが入っているのですが。
/* * 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では上手くいくかも. -- なかま 2005-01-11 (火) 16:08:17
ファイルの中身 $ 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なんとオバカな仕様で,書式が「標準」のときには「タブ」のみ,「数値」のときには「空白+タブ」が入っているようです。というか,数値の後に空白を挿入しているみたいで,行末でも空白一個が入ってます(みっともない)。
はいじま (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-----------------------------
## on Unix with enscript available ps <- pipe("enscript -o tempout.ps","w") capture.output(example(glm), file=ps) close(ps)
winga (2005-01-05 (水) 20:57:16)
windowsXPでR2.0.1を使用しているものです。パッケージfSeriesでaparch分析をしようとしています。aparchSimというコマンドでシミュレーションを行っているのですが、出力された結果には負の値が含まれてしまうことが良く分かりません。aparchSimは自分で設定したパラメータをもとに将来の条件付標準偏差を求めていると思うのですが、この値が負になることはないと思うのですが・・・。help(aparchFit)を読んでも書いていないので、どなたか分かる方がいらっしゃったら教えていただければと思います。
> 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. >こうなります。 -- 青木繁伸 2005-01-09 (日) 01:48:07
QDU (2005-01-04 (火) 08:13:30)
R とは直接関係ありませんが、RjpWiki を見るブラウザーについて前から困っていることがあります。Mozilla で見ると半角英数字が空白になることがある。例えば、このページのトップ見出し中の 「R」が消えてしまいます。一方「RjpWiki」は問題無く表示されています。しかたなしにいつも Konqueror で見ていますが、このページ(だけ)が左右の欄が一部重なりあって見にくくなります(一方 Mozilla では問題無し)。何か解決のヒントをご存知の方いらっしゃいませんか。使用環境は Linux (Knoppix 3.6) です。
ゆき (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 を使用して色をつけていますが、これを、50Mや100M間隔で色を自分で割り振って、色付けを行いたいのですが、どうすれば良いのでしょうか。
## 色分け範囲を自前で指定 -> 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)
## 色分け範囲を自前で指定 -> 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)
## なおこの例の引数 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"
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] }でベクトルを作っています。
chara.list <- substr(chara.a, start=1, stop=1)としました。 -- Akira 2004-12-27 (月) 12:01:39
chara.list <- substr(chara.a, start=nchar(chara.a), stop=nchar(chara.a))あっていますか? -- Akira 2004-12-27 (月) 14:56:28
> test XorY A 1000X X TRUE 1001X X FALSE 1000Y Y FALSE 1001Y Y TRUEというdata.frameなのでXorYのところにrownamesの末尾をもってきました。-- Akira 2004-12-28 (火) 19:04:34
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値だけを得る方法はあるのでしょうか?
# おそらく質問者がやったこと > 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. 氏があっと驚く一行コードを紹介してくれるかも、--- 彼は一体どういう素性の人なのかしら?)
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しか使えないのでしょうか?
> 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"
col = c("gray","blue","black")[ match( Type, c(1,-20000,0) ) ]
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"))
> 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 (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
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.
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
cricket (2004-12-18 (土) 00:14:58)
データフレームで、あるラベルに対応する値だけを変更するにはどうすればいいのでしょうか。たとえば、data$treatmentが"control"の場合だけに、data$observationに10を加えたいのです。初歩的なことだと思いますが、よろしくお願いします。
> 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など -- 青木繁伸 2004-12-18 (土) 00:47:45
> 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
> 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 = 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
nakanaka (2004-12-17 (金) 19:31:57)
他の計量経済学の掲示板にも質問させていただきましたが,なかなかレスがないようなのでここでも質問させてください.現在,単位根検定のADF検定をパッケージ”fSeries”で行っているのですが,トレンド項,ドリフト項の取捨に困ってます.
トレンド項,ドリフト項のt(τ)値を見れば,有意かどうか判断できると思うのですが,”R”では,それらの推計式が示されない(つまり,トレンド項,ドリフト項のt値が見れない)ので推計モデルにトレンド項やドリフト項を含んでいいのかどうかがわかりません.このような場合どういう風に推計式を選択すればよいのでしょうか?fSeriesのマニュアルを見てもこの推計式の選択についての検定は書かれてなかったように思われます.推計式のトレンド項,ドリフト項のt値を見る方法,もしくは,他の解決策があれば教えていただけると幸いです.
よろしくお願いします.
Mari (2004-12-17 (金) 16:48:30)
テーブルを読み込んだときに、そのデータのリストの一覧を入手するような関数ってあるのでしょうか?
例えば、Test というデータのカラムにColumn ------- type1 type3 type2 type1 type3
というような項目があり、このデータが何万件にも及ぶ場合、
何種類のtype があるか、という情報を入手したいのですが。。
show(Test$Column)
ですと、ただ一覧がでてきてしまいます。
調べ方も甘いのかもしれませんが、ご教授いただけますと
幸いです。
よろしくお願いいたします。
fuji (2004-12-14 (火) 13:51:37)
正準相関係数に関して質問です.
cancor()関数では正準相関係数は出力されますが,重み係数は出力されないのでしょうか?
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)+..
のパラメータがどれにあたるのかわかりません.
とても初歩的な質問で申し訳ありませんが,教えていただけないでしょうか?
よろしくお願いいたします.
shiGe (2004-12-13 (月) 07:28:32)
アメリカ50州のデータに関して,それぞれの州の2変量OLSフィットの図(つまり全州において同じ2変量間の散布図及び回帰直線)を1つの図にまとめて描きたい(50個の図が並ぶ)のですが,plot()に関する情報を検索しても,分かりませんでした(というよりどのように検索すべきかもいまひとつ分からないのですが)。
仕方ないので当てずっぽうでplot(Y ~ X, data=hoge, for=STATE)や,for=STATEの部分をby=STATEなどに代えてやっても当然のようにエラーになります。ご存知の方がいらっしゃいましたらご教示いただけますでしょうか。