![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
COLOR(green){SIZE(20){初心者のための R および RjpWiki に関する質問コーナー}}~ COLOR(blue){注意}：新規記事用の入力欄は以下の目次の直後にあります。各記事への追加コメントの入力欄は、各記事の後にあります。質問が長くなる場合は、まず出だしだけ投稿し、次に編集ボタンでそれを修正、追加するのが便利です。>>[[簡単な投稿法]] もし新規に質問をされる場合は、効率的にコメントを得るためにも、善意のコメント者の時間を無駄にしないためにも、以下に特にご注意下さい。[[R のメイリングリストへの投稿記事のガイド]] も参考になります。 - 広い意味で R に関係する話題に限定 - 学校の宿題は自分で考えるべきです - 具体的な状況がわかる様に背景説明をけちらない - [[google などの検索エンジン:http://www.google.co.jp]] や [[CRAN の検索エンジン:http://finzi.psych.upenn.edu/search.html]]、 [[マニュアル:http://buran.u-gakugei.ac.jp/~mori/LEARN/R/]] を調べた上での質問かどうかを書き添えれば、空振りコメントが減るでしょう。せっかくのコメントに「実はそれはすでに調べたことなのですが」と後から言われたときの回答者の不快感を想像すべきです。 - 自分には馴染みでも、他の人には馴染みの無い手法・概念もあります。日本語は得てしてわかりにくい(英語名をそえると案外わかるものです)ものです。 - もし問題が解決した(結局解決しない場合を含め)場合はその旨報告していただくとコメント者や他の人に参考になります。また後で関連情報をこのWikiに要約記事としてご投稿いただくと今後の参考になります。 - 後で他の人が参照する際の便宜のため、簡潔で内容が推察できるようなタイトルをつけて下さい~ - 回答者は神通力を持ってはいません。COLOR(magenta){「これこれしたが、うまく動きません」といった質問}には、あなたの試みたことが回答者に再現できるような、データ、コード、エラー出力(そして、R のバージョン、使用 OS、使用非標準パッケージ名)を添えることをお勧めします。あまりに長くなる場合は適宜編集して下さい。~ - もし,初心者のあまりに稚拙な質問に苛々してしまう方の為には,[[Ｑ＆Ａ (上級者コース)]]があります. ここなら, 稚拙な質問は来ないはずですので精神衛生上にもよろしいかと思います. ~ ~ ---- COLOR(red){過去の記事のアーカイブ} 。アーカイブ中の記事への追加コメントは自由です。ただしあまり古い記事の場合は、このページに新規に項目を立てられた方が良いかも知れません。~ -[[初級Ｑ＆Ａ アーカイブ(1)]] (元記事が 2004-08-03 より 2004-11-30 まで)~ -[[初級Ｑ＆Ａ アーカイブ(2)]] (元記事が 2004-12-13 より 2005-04-26 まで)~ -[[初級Ｑ＆Ａ アーカイブ(3)]] (元記事が 2005-05-02 より 2005-11-07 まで)~ ---- #contents ---- COLOR(red){注意}：投稿失敗した記事をページ先頭の「新規｜編集｜差分｜添付」の「編集」から消すときは、#article と書いた行まで消さないように注意しましょう。新規の記事投稿ができなくなります。もし、編集方法がわからなければ [[簡単な投稿法]] を見てください。 #article **R（UNIX）でグラフを作成後、illustratorで見るとき >[[時雨]] (2006-02-08 (水) 19:31:46)~ ~ こんにちは。時雨と申します。~ UNIX版のRで、plotのグラフを作り、そのグラフをadobe illustrator10で加工しようと思っていたのですが、なぜか、生成されたファイルをillustrator10で開こうとすると、メモリとグラフ本体が90度ずれてしまうというバグが発生します。開いたときに警告文で「所在不明のフォントは初期設定のフォントで代･･･」と出てきて、「書類「〜」には、システムにないか、書式中での指定とは形式が異なるフォントが使用されています。」~ というメッセージが出てきます。~ グラフを作成するときは、 postscript("filename.eps", horizontal = FALSE, onefile = FALSE, height = 9, width = 14) plot(1:10) dev.off() で作成しました。ちなみに、win版でもこの方法でグラフを作ると、「このファイルは、MPS解析で解析できません」と出てきてしまいます。~ いろいろ調べたのですが、さっぱりわからないので、どなたかお教えいただければ幸いです。~ // -family のデフォルト指定に問題があるのかな，postscript("filename.eps", height = 9, width = 14, family="URWHelvetica") とすれば一応問題ないのかも。。。 -- &new{2006-02-08 (水) 19:54:23}; -PDFで編集した方がいいかも. -- [[なかま]] &new{2006-02-09 (木) 02:40:38}; -返信ありがとうございます。一部間違いがあり、horizontal = FALSE は付け加えていませんでした。教えていただいた -- [[時雨]] &new{2006-02-09 (木) 16:08:39}; -フォントを使用したところ、なんとかエラーなしで使うことができました。しかしながら、horizontal=FALSEをつけると、グラフが正方形になってしまいます。heightとwidthを設定しているのに、大きさを変更できません。何か解決法はありますでしょうか？ -- [[時雨]] &new{2006-02-09 (木) 16:09:58}; -horizontal = FALSE, onefile = FALSE, paper = "special", と三つの呪文を唱えて下さい. ?postscript すると書いてあります. -- [[なかま]] &new{2006-02-09 (木) 16:35:50}; -時雨さん！！発言は，編集できますよ。！！ -- &new{2006-02-09 (木) 18:07:46}; #comment **平滑化 >[[ワース]] (2006-02-08 (水) 16:55:39)~ ~ filter等を使ったデータの平滑化をしたいと思っていますが、Webで検索してもいまいち用例が見つかりません。~ 「Rの基本パッケージ中の時系列オブジェクト一覧」にも、filterの説明は3行くらいしかないし…~ Butterworthのようなローパスフィルタを実現する用例を教えていただけないでしょうか。~ // -こうした情報を探すには [[必要な関数を探す]] が参考（関連して、トップページを少し模様替えしました）になるでしょう。しかし、R site search でキーワード検索しても質問があるだけで解答はなし。どうもなさそうです。信号処理特有の手法でしょうか。おそらく同じような効果が得られる別の方法が R にはあると思うのですが、はて？ ちなみに Google で "filetype:r Butterworth" で検索するとそれらしいコードが見付かりますが、これは Ratfor のサブルーティンです。自作の参考にはなるかもしれません。 -- [[MKR]] &new{2006-02-08 (水) 18:34:07}; -site:r-project.org lowpass filter で検索するともう少しましなのが得られたみたいですが。。。 http://cran.r-project.org/doc/packages/rimage.pdf http://cran.r-project.org/doc/packages/rwt.pdf え？見たけど，役立たずだったから聞いた？用例は載ってないか。。。 -- &new{2006-02-08 (水) 18:49:27}; -ありがとうございます。そうですか、無いんですか…。うわー、そうするとデジタルフィルタの設計理論から勉強しなくちゃいけないのか…すごい時間かかりそうだし、無理だ…。 -- [[ワース]] &new{2006-02-09 (木) 09:08:05}; #comment **plotの点の大きさ >[[かず]] (2006-02-06 (月) 17:51:51)~ ~ plot関数で点プロットをしたときの、点の大きさを変えるにはどうしたらいいのでしょうか？~ // -? par でヘルプを見てみましょう。たくさんありますが，そのうちの cex かな？~ cex A numerical value giving the amount by which plotting text and symbols should be scaled relative to the default. Note that some graphics functions such as plot.default have an argument of this name which multiplies this graphical parameter. って，書いてあります。~ [[R Tips:http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html]]なら易しく読めますね。 -- &new{2006-02-06 (月) 18:01:45}; -できました。ありがとうございます。 -- [[かず]] &new{2006-02-07 (火) 09:28:29}; #comment **検定関数の中身を見る関数 >[[souta]] (2006-02-02 (木) 10:38:48)~ ~ 正規分布に従う対応のある2標本平均値の差の信頼区間を求めようとしています。「工学のためのデータサイエンス入門 」P89にある検定関数「t.test(x,y,alternative="two.sided",var.equal=TRUE)」の中身を見たいのですがどなたか教えていただけないでしょうか？このソースを基本に自作の検定関数を作成したいと思いまして・・・。~ ~ 過去の質問に「関数の中身を見る関数」について解説がありましたのでやってみましたが > t.test function (x, ...) UseMethod("t.test") <environment: namespace:stats> と表示されるだけでした。~ // -何回も出てくる FAQ だし，まとめのページもあるし。[[簡単な投稿法]]も読んで欲しいしねぇ。~ 中級者コースへの投稿を繰り返しておく。~ まず，関数名だけをタイプしてみる > deriv function (expr, ...) UseMethod("deriv") <environment: namespace:stats> なんだぁ？これだけか？と思ったら， > methods(deriv) [1] deriv.default deriv.formula と，やってみる。ここに示されるのが，実際に呼び出される候補の関数名。~ 単に関数名をコンソールにタイプするのではなくて，以前に出力されて namespase: の次にある名前をまず書いて，その後にコロン（：）を三つ書いてその後に，その関数名を書く > stats:::deriv.default function (expr, namevec, function.arg = NULL, tag = ".expr", hessian = FALSE, ...) .Internal(deriv.default(expr, namevec, function.arg, tag, hessian)) <environment: namespace:stats> もう一つ > stats:::deriv.formula function (expr, namevec, function.arg = NULL, tag = ".expr", hessian = FALSE, ...) { if ((le <- length(expr)) > 1) .Internal(deriv.default(expr[[le]], namevec, function.arg, tag, hessian)) else stop("invalid formula in deriv") } <environment: namespace:stats> これらの関数のどれがどのように使われるかは，help で検索のこと。~ ここで示されるのは，R の範囲で書かれるプログラムであって，複雑な計算処理プログラムは更にここから呼ばれる。上の例では .Internal(deriv.defarult.... などとある部分かな。~ そのような部分についてのプログラムの詳細は，（えい面倒だなぁ。。），普通のダウンロード・インストールでは実装できないことに依存するので，，，というか，そこまで望む人には，ソースを探索する必要がありますということだけを言えばいいだろうし，そういわれても分からない人は，ソースをみても特に情報は得られないのではないか。 -- &new{2006-02-02 (木) 10:47:59}; -ぜんぜん開発が進んでませんが,[[ごんずい:http://rgonzui.nakama.ne.jp/R/]]もよろしくと言いますか...(ほとんど誰も使ってない...) -- [[なかま]] &new{2006-02-02 (木) 10:59:00}; ->何回も出てくる FAQ だし，まとめのページもあるし。簡単な投稿法も読んで欲しいしねぇ。 -- [[souta]] &new{2006-02-02 (木) 12:27:02}; -上のコメントとの次に「ドウモです。気をつけます。」と入力したかっただけです。ありがとうございました。 -- [[souta]] &new{2006-02-02 (木) 12:28:51}; -編集もできるんですがね -- &new{2006-02-02 (木) 12:41:02}; #comment **read.fortran での入力行数 >[[アール]] (2006-01-30 (月) 22:14:14)~ ~ Mac の R 2.2.1 ですが，read.fortran でテキストファイルから入力するとき，1レコードが複数行の時には合計して2000行しか読んでいないようなのです。1レコード1行の時には，ファイルに存在する行の全部を読み込めます。~~ read.fortran のソースは見ていない（見てもわからないと思うので）のですが，なにか，勘違いしているのでしょうか。~ 5000 行のファイルを作ります > sink("temp") > for (i in 1:5000) print(i) > sink() そのファイルから4行単位に4変数を読み込むと > fn <- file(description="temp", open="r", encoding="SHIFT-JIS") > x <- read.fortran(fn, list(c("a5"),c("a2"),c("a5"),c("a2"))) > print(nrow(x)) [1] 500 # x は 500 行しかありません（元のファイルから2000行しか読んでいません） > close(fn) 2行単位に2変数を読むと > fn <- file(description="temp", open="r", encoding="SHIFT-JIS") > x <- read.fortran(fn, list(c("a5"),c("a2"))) > print(nrow(x)) [1] 1000 # x は 1000 行しかありません（やはり元のファイルから2000行しか読んでいません） > close(fn) > fn <- file(description="temp", open="r", encoding="SHIFT-JIS") 1行毎に1変数を読むと > x <- read.fortran(fn, c("a5")) > print(nrow(x)) [1] 5000 # 元のファイルの全ての行を読んでいます > close(fn) 取りあえず是を回避して，全てのデータを読む方法はありますでしょうか。~ // -Linux版でも同じ結果ですね... -- &new{2006-01-30 (月) 23:40:46}; -自己レスですが。~ ソースを見たところ，read.fortran は read.fwf を呼び，read.fwf のパラメータに buffersize というのがあり，一度に読みとる行数を設定することができるようで，それがデフォルトでは 2000 となっているようですね（すくねーー）。~ read.fortran のヘルプでは 「... Other arguments for read.table」とは書いてあるが，read.fwf を呼ぶよとは書いてないし，ましてや buffersize などという引数についても何にも書いてない。~ 困るじゃないのよねぇ。 > sink("temp") > for (i in 1:5000) cat(sprintf("%4i\n", i)) > sink() > fn <- file(description="temp", open="r", encoding="SHIFT-JIS") buffersize で十分な行数を指定する > x <- read.fortran(fn, list(c("i5"),c("i5"),c("i5"),c("i5")), buffersize=6000) > close(fn) > print(nrow(x)) [1] 1250 # 全部読めたようだ ということで，一件落着？~ それにしても，実際にデータ解析で使うと，いろいろ思いも掛けない落とし穴が見つかるモンですねぇ。 -- [[アール]] &new{2006-01-30 (月) 23:52:56}; -[[BTS:http://bugs.r-project.org/]]によろしくお願いします -- &new{2006-01-31 (火) 00:01:54}; -ま，バグとはいえない（マニュアルの不備）ので -- [[アール]] &new{2006-01-31 (火) 00:08:32}; -Documentationカテゴリのバグに相当するのではないでしょうか？国際的に共有する価値のある報告だと思いますので、ぜひ報告して下さい -- [[岡田]] &new{2006-01-31 (火) 11:21:28}; #comment **ファイル上の関数を使いたい。できればLinuxで >[[かわぐち]] (2006-01-30 (月) 15:08:29)~ ~ The R Bookを入手したので早速13章のロジスティック解析を試みようとしました。ところが、取り込みかたがよくわかりません。sjisのようなのでnkf -eで変換し、attachやsource、CMD INSTALLを試したのですが読み込んでくれません。OSはVinelinux3.2でRは2.2.1です。よろしくおねがいします。~ // -自己レスです。泥くさいですがコピーペーストを繰り返して使えるようになりました。このような使用法でよろしいのでしょうか。 もう一つ、この吉田氏の作成された関数を用いて出した結果を発表や論文に使わせていただきたいのですが、その場合何か使用制限や注釈は必要でしょうか。 -- &new{2006-01-30 (月) 18:39:23}; -[[Rを引用する:http://www.okada.jp.org/RWiki/index.php?R%A4%F2%B0%FA%CD%D1%A4%B9%A4%EB]]や[[なんでも掲示板:http://www.okada.jp.org/RWiki/index.php?%A4%CA%A4%F3%A4%C7%A4%E2%B7%C7%BC%A8%C8%C4#content_1_0]]の「RjpWiki への記事投稿へのご提案」を参照ください -- &new{2006-01-30 (月) 18:45:13}; -公表されているので、関数自体に著作権はありますが、それを利用した結果を発表することは（きちんと引用元を書けば）認められるのが一般的慣例では。自分でちょっと変更を加えたりしている場合など、一応作者に断ったほうがよいかなと思う場合は出版社に問い合わせれば著者に取りついでもらえますよ。 -- &new{2006-01-30 (月) 22:46:05}; -Vineな人ですね(^o^)/ 外部ファイルに書いた関数などを読み込むときは通常 source("ファイル名")を使います。読んでくれないという判断はどのような根拠で行ったのですか？ -- [[谷村]] &new{2006-01-30 (月) 23:29:29}; -source にも，encoding 引数があるので，shift-jis のままでも読めますよ。encoding="SHIFT-JIS" とする。読み込めないと言うが，そのときのエラーメッセージはどうなってましたか？読み込み対象のディレクトリ指定が合っていなかったのかな？ -- [[アール]] &new{2006-01-31 (火) 00:24:44}; -みなさんありがとうございます。わたしのミスだと思いますのでゆっくりやり直してみます。引用に付いては参考にさせていただきます。 -- [[かわぐち]] &new{2006-01-31 (火) 19:41:35}; #comment **グラフを速く描画する方法 >[[げのまー]] (2006-01-30 (月) 11:57:30)~ ~ たくさんのデータポイントを持つグラフを描くのは時間がかかります。~ 描画の更新を抑制するなどして，速くグラフを描く方法はありますか。~ // -あなたの使用しているコンピュータ（OS）の種類，R のバージョンは？~ また，どれくらいのデータポイントをどのように描いたら，どれくらいの時間がかかるのですか？~ 例えばそれが5秒程度であるとして，それでも遅いと思う理由はあるのですか？~ 具体的に質問して頂かないと，具体的に答えられません。 -- &new{2006-01-30 (月) 12:57:59}; -OSはWindowsXP, Rは2.1.0。各々約10000点のデータポイントで8枚のグラフを描画するのに2分ぐらいかかります。我慢できない程度ではありませんが，グラフを速く描画するための一般的な方法があれば知りたいなと思いました。RjpWiki内を探しましたがそういうTipsは見つけられませんでした。 -- [[げのまー]] &new{2006-01-30 (月) 13:21:16}; -描画の更新を抑制するのであれば、Rcmd BATCHを使うのはいかがでしょうか。少しは速くなるかも。他には、CPUの知識が必要ですが、ATLAS,libgotoを組み込むとか・・・。 -- [[okinawa]] &new{2006-01-30 (月) 13:56:43}; -私も無意味に細かい地図の描画でイライラすることがあります。そのようなときは画面に出力せずにファイルに出力します。劇的にはやくなります。png(), postscript(), pdf()などを試してみてください。 -- [[谷村]] &new{2006-01-30 (月) 23:39:47}; -Rcmd BATCHやファイル出力を試してみます。ありがとうございました。 -- [[げのまー]] &new{2006-01-31 (火) 10:02:45}; -（追伸）バッチを使う場合、当然ですが谷村先生が提示された画像のファイル出力が前提です。バッチコマンドが難しいようであれば、谷村先生の方法が簡単です。(^o^) -- [[okinawa]] &new{2006-01-31 (火) 10:41:00}; #comment **行列の抽出における drop オプション >[[アール]] (2006-01-28 (土) 00:00:24)~ ~ 行列の行・列の抽出時に，抽出された結果が1行または1列になる場合に行列ではなくてベクトルになるのを防ぐのが drop=FALSE ではなかったでしょうか？ データ解析のプログラムを書いていて，おかしなことに遭遇したのですが，私の誤解だったのでしょうか。 列を取り出すときにはうまく作動しているのですが，行を取り出すときには変な風になっているように思えます。 > x <- matrix(1:12, 3, 4) > x [,1] [,2] [,3] [,4] [1,] 1 4 7 10 [2,] 2 5 8 11 [3,] 3 6 9 12 行を取り出すときには drop を指定しなければ，その行が「ちゃんとベクトルとして」取り出されています。 > x[1,] [1] 1 4 7 10 しかし，drop を指定すると，drop の如何に関わらず，その行の第一要素しか取り出されていません。 > x[1, drop=FALSE] [1] 1 > x[1, drop=TRUE] [1] 1 列を指定するときにはちゃんと機能しているようなのですが。 > x[,2, drop=FALSE] [,1] [1,] 4 [2,] 5 [3,] 6 > x[,2, drop=TRUE] [1] 4 5 6 なんか変だと思って，こうやったらちゃんと行の場合も行列として取り出せました。 > x[1,, drop=FALSE] [,1] [,2] [,3] [,4] [1,] 1 4 7 10 あれ？？マニュアルのどこかに書いてあるんだろうか。~ 書いてありましたね。~ ? "[" の example の 5 行目。しかし，上の方の説明では，この example を見ろと書いてあるだけのような気がします。~ [[行列Tips大全]]，[[配列Tips大全]] にも [,] を次数の記述が完全になるように複数個おいた後 dorp 引数を書いている例がありますが，複数の [,] が必要ですよと明示的に書いていないように見えます。~ // #comment **polr.summary(MASS)をfunction()中で使う際のエラー >[[mass]] (2006-01-26 (木) 22:31:41)~ ~ library(MASS)にある、polrという関数を使う際に、疑問があります。~ ~ 【エラー内容】~ function()中で、summary.polrを使うと、エラーが出る。~ ~ 【やりたい事】~ 引数をデータセット名にして、polr.summary()で返される結果を利用した計算をしたい。~ ~ 【知りたいこと】~ 「やりたい事」で書いた要求を満たすためには、どうすれば良いのか知りたい。~ ~ 【エラー事例】 library(MASS) #polr data(housing) #dataset func1 <- function(temp){ polr1 <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = temp) return(polr1) } func2 <- function(temp){ polr1 <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = temp) polr2 <- summary(polr1) return(polr2) ###summary.polr()を返す } func1(housing) #エラーなし func2(housing) #エラーあり Re-fitting to get Hessian 以下にエラーeval(expr, envir, enclos) : オブジェクト "temp" は存在しません // -func2 でエラーが起きるのは、polr1 は内部情報としてデータセットが "temp" であると記録しており、summary 関数は "temp" というデータセットを元に必要な量を計算しようとするからです。polr1 に記録されたデータセット名 "temp" は引数ではありませんから関数引数 housing に置き換えられません。次のようにすれば一応ご希望の結果は得られますが、そもそも formula がデータセット housing 特有の変数名 Sat 等を使っていますから、何をされたいのか意味不明。 -- [[MKR]] &new{2006-01-27 (金) 00:12:25}; > func2 <- function(temp){ + attach(temp) + polr1 <- polr(Sat ~ Infl + Type + Cont, weights = Freq) + summary(polr1) + } > func2(housing) Re-fitting to get Hessian Call: polr(formula = Sat ~ Infl + Type + Cont, weights = Freq) Coefficients: Value Std. Error t value InflMedium 0.5663924 0.1046528 5.412110 InflHigh 1.2888218 0.1271561 10.135741 TypeApartment -0.5723552 0.1192380 -4.800107 TypeAtrium -0.3661912 0.1551733 -2.359885 TypeTerrace -1.0910195 0.1514860 -7.202113 ContHigh 0.3602834 0.0955358 3.771187 Intercepts: Value Std. Error t value Low|Medium -0.4961 0.1248 -3.9739 Medium|High 0.6907 0.1255 5.5049 Residual Deviance: 3479.149 AIC: 3495.149 #comment **outer 関数で引用する関数の定義法について >[[アール]] (2006-01-25 (水) 14:17:44)~ ~ こんなことをする人はいないので，不具合でも何でもないが，ちょっと焦ったので。~ outer 関数で引用する関数が返す値は，引数が絡んでいないとエラーになる。具体的には，定数を返そうとしたらエラーになった。そのようなときには，だましてやるとよいようだ。たまたま0を返す関数を作っていたので，x*0 とかにすればオーケーだった。~ > x <- y <- c(-1,0,1) > fun <- function(x, y) return(0) > outer(x, y, fun) 以下にエラーouter(x, y, fun) : dim<- : dims [product 9] は object [1] の長さに整合しません > fun2 <- function(x, y) {z <- x+y; return(0)} > outer(x, y, fun2) 以下にエラーouter(x, y, fun2) : dim<- : dims [product 9] は object [1] の長さに整合しません > fun3 <- function(x, y) return(x*0) > outer(x, y, fun3) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0 // -これは引数が絡んでないというか関数が返す値の構造の問題。同様のことは他の多くの関数でも生じる（plotとか)。 -- [[takahashi]] &new{2006-01-30 (月) 14:30:22}; #comment **関数のプロット >[[ying]] (2006-01-24 (火) 15:50:30)~ ~ 多変数関数 f(x1,x2,x3,x4)=a*x1+b*x2+c*x3+d*x4 (a,b,c,d:定数)を二次元画面で，x1とx2（任意の二つ変数）を横と縦軸として，グラフを描きたいですけど，どうすればいいですか？Rで描けますか？~ // -x3, x4 はどうやって表現するのですか？~ これは，いわゆるその「超平面」の式ですよねぇ。平面ではあるけど，超が付くので，描けないでしょう？~ この式は変数が5つなので，5次元空間です。~ 5次元空間は，我々3次元に住んでいる人間には想像もできないのだから描けないのでは？ ~ f(x1,x2)=a*x1+b*x2+c というのは，3次元空間での平面を表す式で，これを図に描いても，つまらない絵にしかならない。というか，3次元空間の平面をリアルに描くのは難しい。 -- &new{2006-01-24 (火) 16:10:03}; -x3とx4を任意の値を代入する．ここで、例えばx3=3,x4=4とする．後，1つを訂正する：f(x1,x2,x3,x4)=0のグラフを描きたい． -- [[ying]] &new{2006-01-24 (火) 17:06:54}; -f(x1,x2,3,4)=0 のグラフ? y=f(x1,x2,3,4)のグラフということ？ -- &new{2006-01-24 (火) 18:45:39}; -f(x1,x2,3,4)=0のグラフです． -- [[ying]] &new{2006-01-25 (水) 12:17:09}; -あなたもね，f(x1)=a*x1+b=0 のグラフがどんなものなのかがわかっていないんじゃないでしょうか。f(x1,x2)_a*x*x1+b*x2+c=0 だって，同じようなモンなんですよ。 -- &new{2006-01-25 (水) 12:25:34}; -こんな図になります。 #ref(mochi-yaki-ami.png) 黒く表されているのが結果として描かれる平面。これは，x,y平面に平行で，z 軸と 0 で交わる平面。persp 関数はグリッドを描くが，縮小したので黒くつぶれた。~ プログラムは y <- x <- seq(-1,1,length=51) f <- function(x, y) {(x+y)*0} z <- outer(x, y, f) persp(x, y, z, zlim=c(-0.5, 0.5)) のようになります。f の定義で単に0を返すように定義すれば良いところを，おかしなことをやっている理由については別件で。 -- &new{2006-01-25 (水) 13:54:50}; -「f(x1,x2,3,4)=0のグラフです」を解釈すると a*x1+b*x2+c*3+4*d=0 を満たす点 (x1,x2) のグラフ（直線）ということになるような気がするのだけれど。本当にそんなことが知りたいのかしら。 -- &new{2006-01-25 (水) 14:51:04}; #comment **直交表実験における２つの水準間の母平均の差の推定方法 >[[ExperimentalDesign]] (2006-01-24 (火) 14:29:09)~ ~ 直交表実験の結果は次のとおりです。xが特性値でA〜Hまでが2水準の因子です。Eの因子はありません。交互作用はA:B,A:C,B:C,C:D,D:Fのみ考慮しました。~ DATA <- structure(list(x = c(16, 38, 41, 25, 38, 58, 43, 45, 19, 27, 22, 38, 5, 41, 32, 28), A = c(-1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1), B = c(-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1), C = c(-1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1), D = c(-1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1), G = c(-1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, -1), H = c(-1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1), F = c(-1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1)), .Names = c("x", "A", "B", "C", "D", "G", "H", "F"), row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16"), class = "data.frame") 実験の結果から分散分析を次のようにしたところ~ summary(aov(x~.+(A+B+C)^2+C:D+D:F,DATA)) G,H,A:C,B:C,D:Fの効果が小さいことがわかったのでこれを無視することにしました。~ この結果から、A1B1C1D1F1水準(L1)とA1B2C1D2F2水準(L2)での母平均の推定を行うために~ mod <- lm(x~.-G-H+A:B+C:D,DATA) L1<-predict(mod,data.frame(A=-1,B=-1,C=-1,D=-1,F=-1,G=0,H=0), level=.95,interval="confidence",se.fit=T) L2<-predict(mod,data.frame(A=-1,B=1,C=-1,D=1,F=1,G=0,H=0), level=.95,interval="confidence",se.fit=T) として、望む結果を得られました。~ この結果を受け、L1とL2の差がどれくらいあるのか推定したいのですがどのようにすれば推定できるのでしょうか？望む結果は fit = -46.0 , lwr = -60.7 , upr = -31.3 です。~ // -lwr = -60.7 , upr = -31.3 の信頼度は 99% ということでしょうか？~ また，望む結果というのは，L1とL2の結果からどうやったら示す結果が得られるかと言うことでしょうか？~ そういうことであると仮定して， > (L1$fit-L2$fit)[1]+c(-1,1)*qnorm(0.995)*L1$residual.scale [1] -60.68449 -31.31551 かなぁ？ -- &new{2006-01-24 (火) 14:57:34}; -ありがとうございます。できました。信頼度は手元の本では95%で,信頼区間は{mu(L1) - mu(L2)}±t(φε,0.05)*sqrt(Vε/nd)で計算しているのですが、有効繰り返し数ndが因子の割り付け法によって異なってくるので、その影響も考慮して推定してくれる関数はないのかな？と思って質問させていただきました。言葉足らず、かつ後だしで申し訳ありません。 -- [[ExperimentalDesign]] &new{2006-01-24 (火) 17:54:35}; -どういうふうに「できました。」なのか？ -- &new{2006-01-24 (火) 18:07:21}; -質問のL16直交表による２水準のみの因子の実験では > (L1$fit-L2$fit)[1]+c(-1,1)*qnorm(0.995)*L1$residual.scale で信頼区間を得ることができたので、そういう意味で「できました」と書きました。何かおかしかったでしょうか？ -- [[ExperimentalDesign]] &new{2006-01-24 (火) 18:47:14}; -プログラムと実行結果の両方を提示しているので，「できる」のは当たり前なのに，「できました」とことさら書かれ，有効繰り返し数がどうのと書かれていたので，それ以外のやり方であなたが独自に包括的な問題解決にいたったのだろうと思いましたね。信頼度が私は99%なのかと聞いたのに，あなたは95％と書いたし，よくわからなかったのです。違和感を感じたと言うことです。 -- &new{2006-01-24 (火) 18:52:05}; -私の文章がわかりにくく、誤解させてしまったことにお詫び申し上げます。申し訳ありません。ごめんなさい、完全に勘違いしておりました。この実験ではqnorm(0.995)*L1$residual.scaleがt(φε,0.05)*sqrt(Vε/nd)に等しくなると「勝手に」思い込んでおりました。ですので「できた」というのは間違いでした。本当に申し訳ございません。ある水準AiBjCk...と別の水準AlBmCn...の母平均の差の推定を、 {mu(AiBjCk...) - mu(AlBmCn...)}±t(φε,0.05)*sqrt(Vε/nd) という形で、有効繰り返し数ndを自分で計算しなくても求めてくれる関数やパッケージがありますか？という質問のほうが適切でした。重ね重ねお詫びいたします。 -- [[ExperimentalDesign]] &new{2006-01-24 (火) 20:35:48}; -φε，Vε はなんなんでしょう？それがわかるとどうにかできる？（わからないとどうにもならない？）いずれにしろ，数式に出てくる記号の定義は明確にしておかないと。 -- &new{2006-01-24 (火) 21:02:23}; -お手数おかけして申し訳ありません。φεはプーリング後の分散分析表 summary(aov(x~.-G-H+A:B+C:D,DATA)) のResidualsのd.f.です。Vεは同じくResidualsのm.s.です。よろしくお願いいたします。 -- [[ExperimentalDesign]] &new{2006-01-25 (水) 00:15:51}; -そこまでわかっているなら，自分でプログラムするとよいのでは？ -- &new{2006-01-25 (水) 09:22:00}; -自分ではプログラムする技術がないので関数やパッケージを探していたのですが、ないみたいですね。どうも長い間お付き合いいただき、ありがとうございました。 -- [[ExperimentalDesign]] &new{2006-01-25 (水) 12:13:57}; #comment **２つのグラフィックスの軸の一致 >[[mitsu5]] (2006-01-20 (金) 16:55:33)~ ~ きわめて初歩的なことですが、よろしくお願いいたします。~ boxplotと散布図を合成したいと思いました。 a1=rnorm(30) a2=rnorm(30) a3=rnorm(30) x=c(rep(1,30),rep(2,30),rep(3,30)) boxplot(a1,a2,a3) par(new=T) plot(x,c(a1,a2,a3),xlim=par("usr")[1:2],axes=F) こうすればかなり良い線は行くのですが、X軸の1と3で微妙にずれてしまいます。Y軸はこのままで一致しているようです。しかしこれでは駄目かもしれません。~ 軸の一致の調整はどうすればよいのでしょうか。よろしくご指導ください。~ // -par はいりません。最後の plot を points にします（余分なパラメータは要りません） a1=rnorm(30) a2=rnorm(30) a3=rnorm(30) x=c(rep(1,30),rep(2,30),rep(3,30)) boxplot(a1,a2,a3) #par(new=T) points(x,c(a1,a2,a3)) #ref(boxplot.png) 今回はこちらで直しておきますが，今後は[[ヘルプ]]を読んで，きれいに投稿しましょう。 -- &new{2006-01-20 (金) 17:28:28}; -素早いレス、どうもありがとうございます。こんな簡単に出来るのですか。素晴らしい。~ こんなに簡単とは吃驚ですが、google検索 OKADA.JP.ORG を検索でも見つからないし~ RjpWiki->Tips紹介->グラフィックス参考実例集 R-Tips (r-tips.pdf) Statistics with R Graficos Estadisticos con R.pdf R graph gallery.pdf~ これらを見ても複雑すぎて、判らないし、~ ＞ 今後はヘルプを読んで，きれいに投稿しましょう。~ 何処のヘルプを読めばよいでしょうか。ご指導のほどよろしくお願いいたします。簡単明瞭な、解説書があるのでしょうか？ それが一番良いのですが。教えてください。 -- [[mitsu5]] &new{2006-01-20 (金) 22:27:36}; -あれ？あなたのお使いのブラウザがなんなのかわかりませんが「ヘルプ」と言うところの色がほかと違っていませんか。そこをクリックすればいいのですよ。ついでながら，あなたの今回の投稿も，自分では改行したつもりが実際には改行になっていなくって変なところに空白があるだけの記事になっているのは，投稿後に確認すればわかっていると思いますね。今回ももう一回なおしておきましょうか？？？ -- &new{2006-01-20 (金) 22:44:47}; -直すのも重複してしまったようですね。最終調整をよろしくです。。それにしても，投稿方法を確認することなく投稿するというのも無謀というか。。。 -- &new{2006-01-20 (金) 22:50:45}; -何のヘルプ->Rjpwikiの[[ヘルプ]]です（返信のところにリンクがあるでしょう？）。「習うより慣れろ」とも申します。Rの使用とともにwikiの使用もがんばってください。楽をするにもそれ相応の努力が必要です。直そうとしたらかぶりましたね。またもや。調整しました。ここ最近、ひどく初歩的なミスが多いですね。Rの裾野が広がってきたといいうことでしょうか？ -- &new{2006-01-20 (金) 22:43:31}; -追加：同じ座標系でプロットを追加するのは，(1)最初の描画は plot 関数で（プロット自体を描かない（座標系を決めるだけ）とか，座標軸とかその名前を抑制する必要があることも）(2)追加の描画は，点を追加するときは points 関数，線を追加するときは lines 関数，(3)その他の低水準関数を使ってグラフを完成させる。 -- &new{2006-01-20 (金) 23:05:48}; -私の理解不十分で申し訳ありません。どうも皆さんはLinuxですか？それともMac? >ついでながら，あなたの今回の投稿も，自分では改行したつもりが実際には改行になっていなくって変なところに空白があるだけの記事になっているのは，投稿後に確認すればわかっていると思いますね。 私はWindowsです。ここに投稿したＲの a1= .... etc をcopy and pasteすれば、私のR上でそのまま直ぐグラフが描かれますが。これでは、Wikiではいけないのですか？ ここで問題になったのは、LinuxとWindowsの改行のせいでしょうか？ ?r?n と?nの仕様の違いですか？ いや、ヘルプでは、改行は全部 「行末に~を書くと行末改行になります。」これにするのでしょうか？このあたりはどうしたらよいのでしょうか。OSにかかわらず、行末は ~ これであると書いていただければありがたいです。 何か複雑で難しい。 要するに、Windowsから、この投稿部位へ気軽に書くわけにはいかないということですね。 -- [[mitsu5]] &new{2006-01-21 (土) 11:45:13}; -mitsu5さん。私もwindowsです。投稿するためにはwikiの文書の成型法を少しくらい学んでからにしてくださいと言われているのだと思いますよ。wikiではメールのように単にTEXTを貼り付けただけでは、きれいに表示されなかったりします。私も最初は戸惑いましたが、教えていただいたり他のページの内容を見て学んだりして、少しはうまくなりました。「気軽に」とはいかないかもしれませんが、それがwikiの「仕様」でもありますね。 -- [[okinawa]] &new{2006-01-21 (土) 14:16:47}; -あなたの先ほどの記事，かなり読みづらいです。お節介ですが，以下のように直してみましたが同でしょう。~ 私の理解不十分で申し訳ありません。~ どうも皆さんはLinuxですか？それともMac?~ >ついでながら，あなたの今回の投稿も，自分では改行したつもりが~ >実際には改行になっていなくって変なところに空白があるだけの~ >記事になっているのは，投稿後に確認すればわかっていると思いますね。~ 私はWindowsです。ここに投稿したＲの a1= .... etc をcopy and pasteすれば、私のR上でそのまま直ぐグラフが描かれますが。これでは、Wikiではいけないのですか？~ ここで問題になったのは、LinuxとWindowsの改行のせいでしょうか？?r?n と?nの仕様の違いですか？~ いや、ヘルプでは、改行は全部 「行末に~を書くと行末改行になります。」これにするのでしょうか？このあたりはどうしたらよいのでしょうか。~ OSにかかわらず、行末は ~ これであると書いていただければありがたいです。何か複雑で難しい。~ 要するに、Windowsから、この投稿部位へ気軽に書くわけにはいかないということですね。 -- [[mitsu5]] &new{2006-01-21 (土) 11:45:13}; -＞ Windowsから、この投稿部位へ気軽に書くわけにはいかないということですね~ だあれも，そんなこと意ってませんが。~ 初心者だから難しいこと言ってくれるなというのんは，甘えかもしれませんがね。 -- &new{2006-01-21 (土) 15:15:18}; -[[練習用ページ:http://www.okada.jp.org/RWiki/index.php?cmd=read&page=%CE%FD%BD%AC%CD%D1%A5%DA%A1%BC%A5%B8&word=%CE%FD%BD%AC]]を活用しましょう -- &new{2006-01-22 (日) 10:17:12}; -私が書くとまた袋叩きと思いますが、一言。R helpにはほとんど全部 examples が付いています。ここのヘルプを見ると、段落、インライン、ブロックetcと難解なことばかりです。例示が一切ありません。数学の授業ならば、理論だけ教えれば学生は理解したとする方法です。今時、例題の無い授業など出来ないでしょう。例題があるのが普通です。ヘルプを改革してください。改革などと言うとまた大変ですか？ -- [[mitsu5]] &new{2006-01-22 (日) 10:45:05}; -「不満のあるかたが改良する」がこうしたサイトの運営精神です。期待しています（笑） -- [[間瀬茂]] &new{2006-01-22 (日) 11:09:15}; -ヘルプを改革してください→自分でしてください。あなたを含む誰でもヘルプを改変できるようになっているはずです。Open Sourceの世界ではそれを「言い出しっぺの法則」と言います。期待していますよ。 -- &new{2006-01-22 (日) 13:50:09}; -例題や十分なサポートが必要でしたら、「有料」のＳ言語を使いましょう。ここは、あくまでも「無料」の「有志」により運営されているサイトです。お間違えのないように。文面上、やもすると中傷されているように感じるかもしれませんが、それは、「袋叩き」ではなく「礼儀」を重んじているための反応ですので誤解なさらないようにしてください。 -- [[okinawa]] &new{2006-01-22 (日) 13:52:43}; -なにか，勝手にひがんで？いるようですね。~ ちなみにここの[[ヘルプ]]は，ここの掲示板のヘルプではあるものの，実際はこの掲示板がよって立っている「PukiWIki のヘルプ」なんですけどね。ここで文句をいうよりは，[[PukiWiki:http://pukiwiki.sourceforge.jp/]] で文句を言った方がよい（もっとも，そこで文句を言っても，ここで付けられたコメントほどの親切なコメントすら付かないのではないかと思いますが） -- &new{2006-01-22 (日) 18:59:01}; -「[[簡単な投稿法]]」というページを作りましたので，もう大丈夫ですね。もっともっと書けばよいのだろうけど，最小限で収めておくことも必要だろうと思いますので，まあ，これくらいで勘弁しておこう。 -- &new{2006-01-22 (日) 22:12:37}; -だんだん，「最小限」ではなくなっていくところがお宅の所以です。すみませんねぇ。 -- &new{2006-01-23 (月) 19:11:03}; -それにしても，前から感じていたんですが，[[ヘルプ]]の記載は実際と違うところがあるんですねぇ。 -- &new{2006-01-23 (月) 19:57:06}; #comment **polyroot関数について >[[にゃあ]] (2006-01-19 (木) 14:32:36)~ ~ Rでx^2−2x−3＝0という方程式を解こうと思います。~ Rではpolyroot関数を用いて、 polyroot(c(1,-2,-3)) で解けると思うのですが、出力が > polyroot(c(1,-2,-3)) [1] 0.3333333+0i -1.0000000+0i となって、解がおかしいような気がします。~ どなたか検算をお願いします。~ // -ヘルプをよくよみましょう。-- &new{2006-01-19 (木) 15:20:23}; > polyroot(c(-3,-2,1)) [1] -1+0i 3-0i -あ、逆でしたね。失礼致しました。 -- [[にゃあ]] &new{2006-01-19 (木) 15:22:25}; #comment **UNIX:R-2.0.1のバグ？でしょうか write.table出力不良 >[[ホームチーム]] (2006-01-19 (木) 10:15:26)~ ~ 現在、WindowsとUnixでRのプログラミングを行っております。~ バージョンは双方共にR２．０．１です。~ Windowsでは、テキストエディターでプログラムを書いて、RのGUIをダブルクリックし、コピペで動作確認をしています。~ Unixでは、Windowsで確認したプログラムを、FTPで設置し、R CMD BATCHで実行しています。~ ~ プログラムの中に write.table(summary_matrix,file="hoge.txt",sep="?t",col.names=NA) という部分があります。 プログラム自身は、自作関数による処理、クラスター解析や主成分分析など色々な処理を行い、途中結果も多数ファイルとして出力し、その最後にサマリーとして該当hoge.txtを出力して終了するような流れになっています。~ 開発環境であるWindowsXPでは、該当プログラムは全て問題なく動作し、サマリーもテキストファイルとして正しく出力されますが、Unix(ソラリス８)ではなぜか、サマリーとして出したいデータフレームのrow.names情報が、きちんと名前を正しく定義しているにもかかわらず、デフォルトの"1","2","3",,,,という行番号で出力されてしまうバグ？に苛まれております。~ UNIXで動作させたときに本当に正しくrow.namesがセットできているか、write.table直前で、cat(row.names(summary.matrix))をしてみても、Win,Unix共に正しく情報はセット出来ていることを確認しています。~ 本現象について、お気づきの点等ありましたらご教授いただきたく考えております。~ // -最新版にしてみて，なおかつエラーが発生するかどうか見てください。旧版にバグがあろうとなかろうとだれも興味を持ちません。 -- &new{2006-01-19 (木) 11:10:28}; -確かに仰る通りですね、考えてるくらいなら、最新版で確かめてみるべきでした。すいません。。。(；_；) -- [[ホームチーム]] &new{2006-01-19 (木) 11:35:26}; -「旧版にバグがあろうとなかろうとだれも興味を持ちません。」-> 日本語版 R fortunes があれば即採用（笑）。 -- &new{2006-01-19 (木) 16:14:35}; #comment **Rバグフィックス、機能追加、改良の履歴を追うには？ >[[やたま]] (2006-01-19 (木) 10:07:08)~ ~ 現在Windows環境でR-2.0.1を使っております。~ ただ、現状最新のVersionは2.2.1かと思いますので、バージョンのアップを考えています。~ そこで、R2.0.1→2.2.1に至るまでは数多くのバージョンアップがあったかと思いますが、その際数々のバグフィックスが等々あったと思いますが、その履歴等を追うにはどこのページを見ればよろしいでしょうか？~ Rをインストールしたフォルダ以下にある「CHANGES」テキストファイルを見るしか方法はないでしょうか？~ 稚拙な質問で恐縮ですが、アドバイスいただければ幸いです。~ // -[[RソースのNEWSファイル:http://cran.md.tsukuba.ac.jp/src/base/NEWS]]を見るのが最も正確ではないでしょうか。これは2.0.0以降ですが、より古い情報はソースアーカイブ中のONEWSとかOONEWSというファイルに書かれています。 -- &new{2006-01-19 (木) 15:45:09}; -ありがとうございます。NEWSファイル早速見てみたいと思います。誠にありがとうございます。 -- [[やたま]] &new{2006-01-19 (木) 18:00:27}; #comment **svmのプログラムについて >[[ying]] (2006-01-17 (火) 17:57:37)~ ~ 今svmについて理解している．でも，分からない所があって，教えていただけませんか？~ help(svm)からsvmについてたくさん書いてるけど，実際svmのprogramを知りたいです．どうすればいいですか？~ 実際このようにやってみたら： > edit(svm) function (x, ...) UseMethod("svm") <environment: namespace:e1071> これしかなかった．もっと詳しいプログラムについて知りたいので，教えてください．~ // -プログラムを見たいということなら，以下のようだけど。それで用が足りるのかどうか。 > library(e1071) > e1071:::svm.default > e1071:::svm.formula ? svm をして，ヘルプおよび，そこからリンクされているウエッブページの内容をよく読むことの方が先決ではないでしょうか。-- &new{2006-01-17 (火) 18:08:29}; -help(svm) の最後に参考プログラムがありますが、これではたらないという意味ですか。たしか、前に google で R のコードを検索する方法が紹介されていましたが、どこだったか？ -- &new{2006-01-17 (火) 18:48:30}; -トップ頁の R site search を使いキーワード "svm" で検索するとか。 -- &new{2006-01-17 (火) 20:12:43}; -ありがとうございます！！とても助かりました！！ -- [[ying]] &new{2006-01-17 (火) 20:15:12}; #comment **モンテカルロ法による円周率 >[[しゃけ？]] (2006-01-16 (月) 15:46:54)~ ~ 下記のような方法を見つけたのですが、一行目から何をやっているのか全く意味不明です。申し訳ありませんが、誰か詳しく解説していただけないでしょうか？ ULONG_MAX = 4294967295 以下ばっさり，削除（それにしても，投稿法くらい守ってね）~ // -R と関係ありませんね。宿題だったらここはおかど違いですよ。 -- &new{2006-01-16 (月) 16:16:20}; -rubyのRだと思ったのかな？rubyはここhttp://www.ruby-lang.org/ja/ -- &new{2006-01-16 (月) 16:57:34}; -ちなみに，R だと， n<-100000;sum((colSums(matrix(runif(n*2)^2,2))<1))*4/n 位でできるしィ -- &new{2006-01-16 (月) 17:23:28}; -もっと短くなる n<-100000;sum(runif(n)^2+runif(n)^2<1)*4/n もっともっと？ -- &new{2006-01-17 (火) 07:57:38}; #comment **plot.svmについて >[[ying]] (2006-01-14 (土) 16:42:59)~ ~ データirisについてサポートして，グラフを描きたいけど、多次元のデータだから、二次元のグラフを描こうと思う。どうすればいいですか？それと、Ｒの中にplot.svmについての例：~ plot(m2, iris, Petal.Width ~ Petal.Length, slice = list(Sepal.Width = 3, Sepal.Length = 4)) がある。これはいったいどんな意味ですか？~ 教えてください！！~ // -plot.svm が何なのか，私にはわかりません！！！！！！！！！！！~ 示された例でも，m2 はなんなのでしょうか！！！！！！！！！！！！！！！~ 一般的に多次元のデータを二次元平面にどのように描いたらよいのかということなら，バイプロットとか，主座標分析とか主成分分析とか正準判別分析とかいろいろありますね！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！ ！！！！ -- &new{2006-01-14 (土) 18:38:57}; // -yingさんの疑問が?plot.svmで得られるexampleを指していると解釈しますと、 ## a simple example library(MASS) data(cats) m <- svm(Sex~., data = cats) plot(m, cats) ## more than two variables: fix 2 dimensions data(iris) m2 <- svm(Species~., data = iris) plot(m2, iris, Petal.Width ~ Petal.Length, slice = list(Sepal.Width = 3, Sepal.Length = 4)) ですから、plot(m2,...)は「m2についてfomulaで指定する次元について2Dplotする、そのときの色分けはsliceで指定する」といっていると思います（推測）。-- [[Akira]] -パッケージ e1071 中にある svm (Support Vector Machine) 用のプロットメソッドのことでしょうか。SVM とは比較的新しい判別手法のことですが、ここで手軽にその意味を説明できるようなものではないですね。 -- [[MKR]] &new{2006-01-15 (日) 11:01:49}; -e1071に実装されているSVMは[[LIBSVM:http://www.csie.ntu.edu.tw/%7Ecjlin/libsvm/]]ですので[[guide:http://www.csie.ntu.edu.tw/%7Ecjlin/papers/guide/guide.pdf]]をご覧になってはいかが？ -- [[Akira]] &new{2006-01-15 (日) 22:13:25}; -ご返事ありがとうございます．ここで，slice = list(Sepal.Width = 3, Sepal.Length = 4)の意味はよく分からないです．何で３と４を取るの？他の値に変えて見たけど，正しいグラフが出ると出ない時がある．教えてください！お願いします！ -- [[ying]] &new{2006-01-17 (火) 17:28:04}; -?plot.svmに説明があるようです -- [[Akira]] &new{2006-01-18 (水) 00:06:23}; #comment **Rを使った学術論文 >[[Akira]] (2006-01-14 (土) 13:33:59)~ ~ Rの引用（Materials&Methodsとして）を何種類か見たのですが統一されていないようにも思います。それとも、お作法などあるのでしょうか？~ bioconductorとExcelでもできる簡単な処理ですが、投稿したいと思ってます。~ // -投稿先の雑誌が指定する投稿規定では引用の記述方法も決まっていると思います。これは雑誌のよって異なるため、論文（雑誌）によって引用の記述が異なるのは自然なことと思います。私の場合は、citation()得られるBibTeX形式の書誌情報を使って、雑誌指定のスタイルファイル（または自作のスタイルファイル）で処理しています。本文中には必要に応じてRのバージョンや使ったパッケージ・バージョンを付記するとよいかも。 -- [[谷村]] &new{2006-01-14 (土) 18:21:52}; -ありがとうございます。参考にさせていただきます。ただ、生物系でTexを使う人は少ないようです。 -- [[Akira]] &new{2006-01-16 (月) 08:25:03}; -補足させてください。citation()ではBibTeX形式以外にも一般的な形式での書誌情報が出力されます。また、BibTeX形式のファイルをお使いの文献データベースソフトインポートして使うという方法もあります。文献ソフト自体が対応していなくても、誰かが作った変換ソフトはあると思います。ちなみに医学系もTeXを使う人は希少です。それから口やかましくて申し訳ありませんが、(誤)Tex(正)TeXです。 -- [[谷村]] &new{2006-01-16 (月) 16:02:17}; -[[Rを引用する]]を作成しておきました -- [[谷村]] &new{2006-01-16 (月) 16:49:09}; -TeX を使う人は LaTeX を使う人より少ないと思いますよ。 -- &new{2006-01-17 (火) 18:14:48}; #comment **vectorの計算結果をmatrixまたはvectorに保存したい >[[柳]] (2006-01-13 (金) 13:02:58)~ ~ x<-c(1:3) y<-c(4:6) z<-matrix(0,length(y),length(x)) myfunction<-function(x,y){ for(i in 1:length(x)){ z[,i]<-x[i]+y } } をするとparse errorが出ますがなぜかぜんぜんわかりません~ 教えて下さい~ // -parse error というのは，出ませんでしたが。質問するときには，コンソールに入力したものをそのままコピー・ペースとしないとこちらで再現することができなくなりますので，ご注意のほどを。~ また，この後に myfunction(x,y) というように関数を引用するのでしょうが，その結果は z に付値することを意図しているのでしょうが，付値されませんね。それは z[,i] の後が <- になっているからで，大域変数に付値するには <<- を使わなければならないと言うことはご存じでしょうか。もっとも，<<- を使うのはお行儀の悪い方法なので，myfunction 内で局所変数に付値し，それを return で返し，返ってきたものを z に付値するというプログラミングの方がよいと思います。~ 質問（プログラム）の投稿についても，再度[[ヘルプ]]をご覧になり，今後はきれいに投稿してくださいね。今回はこちらで直しておきます。 -- &new{2006-01-13 (金) 13:09:31}; #comment **関数のソースコード >[[まいける]] (2006-01-12 (木) 17:35:58)~ ~ はじめまして。~ Rにoptimという関数がありますが、そのソースコードを入手することは可能でしょうか？~ また、そのソースコードはC等で書かれているものなのでしょうか？~ 大変初歩的な質問で恐縮ですが、なにとぞご教示よろしくお願い致します。~ // -コンソールのプロンプトに，optim と入力してみましょう -- &new{2006-01-12 (木) 17:49:31}; -早速のアドバイスありがとうございました。 -- [[まいける]] &new{2006-01-12 (木) 18:10:04}; -もっと先を知りたいのなら，src/main/optim.c を見ればよいでしょう。 -- &new{2006-01-12 (木) 18:13:54}; -すみません。もうひとつお教えください。src/main/optim.cは各人のPCにおけるRの環境下にあるということでしょうか？ -- [[まいける]] &new{2006-01-12 (木) 20:30:49}; -R をインストールすれば必ずソースもインストールされるわけでは無い（特に MSW ではバイナリーだけ）ですから、CRAN の R Sources から最新ソースをダウンロードし、解凍します。 -- &new{2006-01-12 (木) 21:34:43}; -何度もご親切にお教え頂きありがとうございました。無事入手できました。 -- [[まいける]] &new{2006-01-13 (金) 10:16:08}; -http://rgonzui.nakama.ne.jp/search?q=package:R-2.2.0%20optim&fm=c というのもあります． -- [[なかお]] &new{2006-01-14 (土) 17:13:59}; #comment **factorのlevelを削る方法 >[[Moz]] (2006-01-10 (火) 20:25:23)~ ~ x <- factor(c(rep("A",4),rep("B",2),rep("C",3))) x [1] A A A A B B C C C Levels: A B C x <- x[x!="C"] x [1] A A A A B B Levels: A B C levelからCを消すにはどうしたらよいのでしょうか。table()で集計する~ ときにじゃまになります。無理矢理~ x <- as.factor(as.character(x)) x [1] A A A A B B Levels: A B としましたが、もっと正統なやり方があるのではないかと思って質問しました。~ // -正統かどうか知りませんが x <- factor(x[x!="C"]) -- [[takahashi]] &new{2006-01-10 (火) 20:54:44}; -えぇと、 x <- x[, drop=T] ですね。 -- [[びい]] &new{2006-01-11 (水) 06:11:20}; -takahashiさん、びいさんありがとうございました。今後は、x <- x[x!="C",drop=TRUE] を使うことにします。Tという名前のオブジェクトを昔に作成したのがRの起動時に読み込まれていたので、drop=Tではうまくいかず、TRUEと書く必要がありました。 -- [[Moz]] &new{2006-01-11 (水) 19:40:41}; #comment **データベースからのデータを使いたい >[[ひらい]] (2006-01-09 (月) 00:10:48)~ ~ windowsでR2.2.0を使っています。~ MySQLの中のデータベースに入っているデータを使って、ユークリッド距離(関数dist)を求めたいのですがどうしたらよいでしょうか？~ どなたかご教授お願いします。~ // -RMySQL パッケージをインストールし、研究すべし。 -- &new{2006-01-09 (月) 00:37:00}; -ちなみにRからデータベース内のデータはすでに見ることができます。 -- [[ひらい]] &new{2006-01-09 (月) 00:56:52}; -あなたも，もう少し具体的に質問しなくちゃいけないし，回答しようと思う側も，回答してやっても良いかなという気があるなら，もうちょっと誘導してやっても良いじゃないの？回答者が期待するような知識・経験がないからこそ質問しているんじゃないのかなぁ -- &new{2006-01-09 (月) 21:24:59}; -老婆心から質問のしかたについて書いてみました。（気を悪くしないでくださいね） OS:windows(XP)、 R:ver2.2.0(J)、 DB:MySQL、RMySQLにてRとダイレクトリンクしています。、 R内に取り込んだdataframeは XXX、 dataframeのフィールド名は、A,B,C、 このデータを使ってユークリッド距離(関数dist)をもとめたいのですがどうしたらよいでしょうか？ このような形で質問していただけると、正確な回答がつくと思いますよ。回答がつかないのは質問の仕方に問題があることが多いです。（質問者の意図が、回答者にわからないため、パスされてしまいます） -- [[okinawa]] &new{2006-01-10 (火) 12:59:41}; -既に実施して上手く行かなかった例があると他の人の参考になりますよ。 -- &new{2006-01-10 (火) 14:40:28}; #comment **Rの画面に表示される速度 >[[aa]] (2006-01-05 (木) 00:18:37)~ ~ Windows XP で，Rを使っています。~ 長いプログラム（1000行ほど）をRのコンソールに一気に流した際，画面に表示される速度がversionによって全く異なります。~ ~ 余程複雑な計算でなければ，新しいversion（2.20以降）の方が，画面に表示される時間がかなり遅くなります。~ 計算の実行速度ではなく，表示時間が遅くなります。~ この原因をご存知の方，アドバイスを頂けると幸いです。~ // -バージョンこれこれではしかじか秒，違うバージョンかくかくではこれこの通りの秒数のように，具体的に書いて頂かないとね〜。数秒が数倍になったって，いいじゃないの。 -- &new{2006-01-05 (木) 00:31:02}; -version は，1.81-2.11 vs 2.20-2.21 です。具体的な秒数は，実行速度でないので，system.time()では測れません。プログラムの長さ/内容に依存されますが，数秒が数倍ではなく，数分の違いが出ます。 -- [[aa]] &new{2006-01-05 (木) 00:39:47}; -重箱の隅つつきで申し訳ないが。R のバージョン表記は 2.2.1 のごとし。表示時間はストップウオッチででも計るか。。 -- &new{2006-01-05 (木) 00:52:50}; -たとえば，a <- sqrt(3) というのを1000行入力して，次のプロンプトが出るまでにどれくらいの時間数がかかるのでしょうかね。私のマシンで2.2.0だと20秒くらいかかりました。古いバージョンでどれくらいかかるか，インストールし直さないと測れないのでやめておきます。~ [長いプログラム（1000行ほど）をRのコンソールに一気に流す]というのは，プログラムをコピーして，Rのコンソールにペーストするみたいなことを意味しているのでしょうかね？前述の1000行のプログラムも，ファイルに書き出しておいて source("foo.R") みたいにすると，2，3秒で次のプロンプトが出ますが。~ //修正: sourece -> source 遅くなる理由が知りたいというならしようがないですが，我慢の限界内で快適に実行するにはどうしたらよいかと言うことなら，そのような対処法を採るのも一法かと。 -- &new{2006-01-05 (木) 00:57:03}; -ウインドウのサイズを小さくすると速くなりません？ -- &new{2006-01-05 (木) 08:51:03}; #comment **変数が行として表されているファイルからのインポート >[[本太]] (2006-01-01 (日) 23:36:29)~ ~ 変数が列として表されているCSVファイル(d1.csv)からデータをxに読み込む場合、 x<-read.table("C:/R/d1.csv",header=F,sep=",") と記述しています。~ 変数が行として表されているCSVファイル(d2.csv)を読み込み上記と等しいデータフレームxを作成するためには、どのようにすればよいのでしょうか？~ すみませんがよろしくお願いします。 ----- d1.csv ------ height weight age 172 72.6 25 159 59.8 27 164 78.7 42 179 85.4 35 164 65.7 30 ------------------- ------------d2.csv------------------- height 172 159 164 179 164 weight 72.6 59.8 78.7 85.4 65.7 age 25 27 42 35 30 -------------------------------------- ~ // -たぶん,ヘッダがついたファイルなのでは?などと思いつつ, x <- t(read.table("c:/r/d1.csv" ,header=T,sep=",")) でいいかと. -- [[なかま]] &new{2006-01-02 (月) 01:12:43}; -data.frameがほしいということですので、 x <- data.frame(t(read.table("c:/r/d1.csv" ,header=T,sep=","))) でしょうか？ただ、要素の属性は変わると思います。 -- [[Akira]] &new{2006-01-02 (月) 13:32:36}; -あなたの記述が若干不正確ではないかと思うのですが， 変数が列として表されているCSVファイル(d1.csv)からデータをxに読み込む場合、 x<-read.table("C:/R/d1.csv",header=F,sep=",") は x<-read.table("C:/R/d1.csv",header=T,sep=",") ではないかと思うのはさておき，~ 対象とするファイルが d2.csv で，その内容が，まさしく以下のように変数名プラス変数値...と言うことならば~ sep="," ということは，カンマで区切られていると言うことですよね~ 実行結果で示すよりは，入力ファイルの形式で示した方が誤解が少ないでしょうね。 height,172,159,164,179,164 weight,72.6,59.8,78.7,85.4,65.7 age,25,27,42,35,30 さて，そうであれば，あなたの本当に望むことは（回答例に示されるように）そんなに簡単には求められないので，以下のようになるのではないかと思いますね。~ 回答案も，実際にやってみればすぐにわかるとはいうものの~ もっとも，回答者をせめるのではなく，あなたが十分な情報を提供できなかったのではないかと反省すべきではありますが。 > x <- t(read.table("d2.csv" ,header=F,sep=",")) > name <- x[1,] > x <- x[-1,] > x <- matrix(as.numeric(x), nrow(x)) > colnames(x) <- name > x <- data.frame(x) > x height weight age 1 172 72.6 25 2 159 59.8 27 3 164 78.7 42 4 179 85.4 35 5 164 65.7 30 d2.csv を読み込んだときと同じ結果ということは，このようなことを意味するのではないかと，私は考えましたのですが，深読みですか？ -- [[あけましておめでたかったろう]] &new{2006-01-02 (月) 23:17:57}; #comment **Rcmdrテキストファイルのインポート >[[にゃあ]] (2006-01-01 (日) 04:41:36)~ ~ Rcmdrでデータのインポート⇒テキストファイルからでCSVファイルを読んだ後、データを表示を選ぶと、データエディタを動かそうとすると、画面がひきつったようになって沢山の残像が出るような画面になってしまいます。~ これってバグなんでしょうか?~ また、同様に出力した筈のCSVファイルでさえ、データを読み込める場合と読み込めない場合もあります。これもバグなんでしょうか?~ // -CSVファイルの容量はどれくらいなのでしょうか？重すぎたりしませんか？？ -- &new{2006-01-01 (日) 11:20:30}; -Windows系ならたぶんデータが巨大過ぎてどこかのヒープが溢れたんでしょう. わたしは32Mぐらいで良いと思うのですが現状は2M程(Win9x系の関係だと思うのですが,詳しく調べていません)です. Unix系なら涙ながら(殆ど拷問)に(Xt使って書き直した方が良かったかも...てっきりlibXに拘ってると思いきや最近のはlibXtも使い出した...涙)国際化しましたdeなら問題は無いと思ってるのですが... -- [[なかま]] &new{2006-01-01 (日) 11:58:22}; -読み込みが成功したファイルは1.73MBで、失敗したファイルは1.38MBです。両方とも同じように出力したファイルですが、失敗した方では『〜はデータフレームではないので、アタッチできません』と言ったエラーメッセージが表示されます。また、成功したファイルでも、例えば統計量のプルダウンメニューから、『モデルへの適合』を選んだらポップアップメニューが出てきますけれども、変数選択のスクロールバーがおかしかったり、また、読み込んだデータがアクティブデータとして表示されない場合もあります。 -- [[にゃあ]] &new{2006-01-02 (月) 18:09:37}; -一つお聞かせ(OSやRバージョン)願えれば,わたしは通常の3倍!の速度でデバック出来ますが...(WinかLinuxかMacかわかればですが...) -- [[なかま]] &new{2006-01-02 (月) 21:59:10}; -抽象的なバグ報告は，何の情報も提供しない。どうも，ファイルサイズの問題でもなさそうなので，そのような場合には（この場合には難しいのかもしれないが），問題を再現できる最小限のデータセットを探索するのがよろしいかと。バグ・バスターズとしては，あいまいな状況報告では手の施しようもないのが本当のところです。対策を望む場合にはまず，そのエラーを再現することが大切なわけで，そのエラーを再現できるデータを提供することがなにより大事なんです。。と，他人事ながら（他人というのは，私がバグ・バスタースではないということですが） -- &new{2006-01-02 (月) 22:12:56}; -WindowsはXPです。Rは闇2.2.0です。 -- [[にゃあ]] &new{2006-01-02 (月) 22:45:57}; -それだけの情報では，まだ，不足ではないかなと。。。 -- &new{2006-01-02 (月) 23:18:39}; -しかしながらデータをここに載せるわけにはいかないので。行数が3,500行くらいあるのです。他に何が必要なのでしょう? -- [[にゃあ]] &new{2006-01-02 (月) 23:26:47}; -とりあえず 行数 は伺えましたので列数も教えてもらえますか。合間見て考えてみます 。 -- [[なかま]] &new{2006-01-02 (月) 23:48:27}; -あと、学生さんか社会人の方か存じませんが、タイトルは後で参考になるようなマトモな物に付け直してくださいね。 -- [[なかま]] &new{2006-01-02 (月) 23:55:12}; -成功した方の1.73MBのファイルは行数は3,945行105列です。失敗した1.38MBのファイルの行数は3,152行105列です。 -- [[にゃあ]] &new{2006-01-03 (火) 01:19:27}; -同じ現象は再現できませんでしたが、ちょっと気になる動作を発見しました。たまたまオーバランした部分が影響の出るときと出ないときの違い＊かも＊しれません。 （ちょっと時間がかかるかも...） -- [[なかま]] &new{2006-01-06 (金) 16:58:18}; -＞なかまさん わざわざありがとうございました。現時点では闇2.2.1にアップグレードしたのですが、問題は依然解決してません。こちらでももうちょっと調べてみたいと思います。 -- [[にゃあ]] &new{2006-01-06 (金) 18:47:21}; -失敗するファイルの列数を100以下にすれば問題ないかも... -- [[なかま]] &new{2006-01-07 (土) 02:25:18}; -分かりました。やはり直接CSVを読み込まないで、一旦Excel上でEditして列数減らしてみます。 -- [[にゃあ]] &new{2006-01-07 (土) 14:11:09}; -う〜ん、やっぱり具合が宜しくないようです。データが読み込めてもアクティブデータセットとして認識はしてくれないようです。 -- [[にゃあ]] &new{2006-01-07 (土) 15:56:52}; -害が無ければ(無保証,無責任でよければ)、一度私宛(eiji.nakama@gmail.com)に送ってもらえますか?. -- [[なかま]] &new{2006-01-07 (土) 21:20:22}; #comment **日付変数をデータフレームの要素に格納したい >[[じゃんぽけ]] (2005-12-29 (木) 11:43:38)~ ~ Tips等調べて１日格闘してみましたが、わからないポイントがあり投稿しております。~ ~ 【質問内容】~ 日付に関する演算を繰り返し行うプログラムを書く予定でして、~ ~ １．文字型で格納しておいて取り出す度にDATES()で変換して使う~ ２．文字型で格納しておいて文字型に変換したものを別の変数に保存して使う~ ３．最初から日付型で格納しておいて取り出す~ ~ で３．の方法で実行する方法がないかなと思い模索している所です。~ ~ 日付型の変数をデータフレーム内の要素に代入しようとすると、シリアル値が代入されてしまいます。~ これを日付型のまま代入する方法はございますでしょうか？~ （列をベクトル演算で変換した場合には日付型で代入できるようです。）~ ~ また、まだまだ勉強不足のため、下記のような挙動になる要因が何なのかが理解できていない状況です。~ ~ 何度も質問致しまして恐縮ではございますがご教授頂ければ幸いです。~ よろしくお願いします。~ ~ 【プログラム例】~ > DATE <- dates("05/12/31",format="y/m/d") > DATE2<- "05/12/31" > > test1 <- data.frame(a=rep(NA,3)) > test1$a[1] <- DATE > test1 a 1 13148 2 NA 3 NA > mode(test1$a) [1] "numeric" > > test2 <- data.frame(a=rep(NA,3)) > test2$a[1] <- DATE2 > test2 a 1 05/12/31 2 <NA> 3 <NA> > mode(test2$a) [1] "character" > > test3 <- data.frame(a=rep(NA,3)) > test3$a[1] <- dates(test1$a[1],format="y/m/d") > test3 a 1 13148 2 NA 3 NA > mode(test3$a) [1] "numeric" > > test4 <- data.frame(a=rep(NA,3)) > test4$a <- dates(test2$a,format="y/m/d") > test4 a 1 05/12/31 2 <NA> 3 <NA> > mode(test4$a) [1] "numeric" // -chronライブラリのdatesクラスですね。型はnumericでクラスがあるだけなので，要素の代入ではクラスが代入されないことから，そういう結果になるのだと思います。最初のdata.frame宣言時にa=dates(rep(NA,3))としておけば大丈夫です。chronライブラリを使わないでDateクラスを使った場合でも，基本的に同様です。 -- [[中澤]] &new{2005-12-29 (木) 12:58:06}; -中澤さん、ありがとうございます。事前にデータフレームで定義することで解決できました。参考までに、定義の仕方で出力が異なりました。 > data.frame(a=rep(dates("12/31/05"),3),b=0) a b 1 13148 0 2 13148 0 3 13148 0 > data.frame(a=dates(rep("12/31/05",3)),b=0) a b 1 12/31/05 0 2 12/31/05 0 3 12/31/05 0 ところで、上記のようにデータフレームにNAなど値を入れる場合であれば列の型を定義できるのですが、要素が全くない場合に事前に型を宣言する方法はあるのでしょうか？ 自己レスの修正をしました。 データフレーム作成時に、変数を日付型に宣言しておくことで 空データになっても日付型を保持していますが、 更新時に個別の行を更新しなければ日付型を保持できないようです。 > #データの定義をした後に空データにする > test <- data.frame(a=dates(c(NA,NA),format="y/m/d"),b=c(NA,NA)) > test <- test[FALSE,] > test [1] a b <0 rows> (or 0-length row.names) > > #空のデータフレームに行を追加するとちゃんと日付を保持 > test <- rbind(test,c("2000/01/03",0)) > test a b 1 00/01/03 0 > > #データフレームが１行の場合、全行をまとめて更新すれば日付型で更新 > test$a <- dates("2000/01/22",format="y/m/d") > test a b 1 00/01/22 0 > > #更に1行追加。（実験1用と実験2用のデータを作成） > test1 <- rbind(test,c("2000/01/05",0)) > test2 <- test1 > test1;test2 a b 1 00/01/22 0 2 00/01/05 0 a b 1 00/01/22 0 2 00/01/05 0 > > #実験１：データフレームが２行以上の場合、 #全行をまとめて更新すれば日付型でなく なる > test1$a <- dates("2000/01/22",format="y/m/d") > test1 a b 1 10978 0 2 10978 0 > > #実験１：このデータフレームに日付を文字型で #いれようとすると変数aは文字型に変化 > test1 <- rbind(test1,c("2000/01/03",0)) > test1 a b 1 10978 0 2 10978 0 3 2000/01/03 0 > mode(test1$a) [1] "character" > > #実験２：この場合は、データフレームの行を #それぞれ更新すればOK > > test2$a <- dates(rep("2000/01/22",length(test2$a)),format="y/m/d") > test2 a b 1 00/01/22 0 2 00/01/22 0 -- [[じゃんぽけ]] &new{2005-12-29 (木) 22:00:28}; #comment **条件に合致する行を抽出した上で、新たな変数（列）を追加・代入したい >[[じゃんぽけ]] (2005-12-28 (水) 11:51:41)~ ~ 一通りTips や Q&A を見てプログラムを書いてみたのですが、もっと簡単にできる方法がないかと思い投稿しております。~ ~ 以下のように、「FOR文を利用して、条件に合致する行を抽出した上で、新たな変数（列）を追加・代入したい」と考えています。~ ~ 考えたやり方では、~ 「１．事前に追加する変数（列）は定義しておく」~ 「２．抽出条件に当てはまる行がない場合は、IF文でエラーを飛ばす」~ ということで解決はできているようなのですが、もっと簡単にできる方法があるのではないかと思っています。~ ~ ご教授頂ければ幸いです。宜しくお願いします。~ ~ > #その１：新たなdという変数（列）に直接的に値を追加・代入 #テストデータを作成 > test <- data.frame(a=c(1,2,3),b=c(4,5,6),c=c(7,8,9)) > test a b c 1 1 4 7 2 2 5 8 3 3 6 9 > > #合致する行がなくても代入するとエラー > for(i in 4:3){ + test[test$a == i,]$d <- 2 + } 以下にエラー"$<-.data.frame"(`*tmp*`, "d", value = 2) : 置き換えは 1 列ですが、データは 0 列です > > #IF文で合致しない行の場合にスキップしてもエラー > for(i in 4:3){ + if(nrow(test[test$a == i,]) != 0){ test[test$a == i,]$d <- 2 } #こうやれば合致しないものにはエラーはでない + } Warning message: 4 個の変数を与えて、3 個の変数を置き換えようとしています in: "[<-.data.frame"(`*tmp*`, test$a == i, , value = list(a = 3, > > #その２：事前にdという変数（列）を作っておいて値を追加・代入 #テストデータを 作成 > test <- data.frame(a=c(1,2,3),b=c(4,5,6),c=c(7,8,9)) > test$d <- NA > test a b c d 1 1 4 7 NA 2 2 5 8 NA 3 3 6 9 NA > > #合致する行がなくても代入するとエラー > for(i in 4:3){ + test[test$a == i,]$d <- 2 + } 以下にエラー"$<-.data.frame"(`*tmp*`, "d", value = 2) : 置き換えは 1 列ですが、データは 0 列です > > #IF文で合致しない行の場合にスキップすればＯＫ > #ただし、文章としてはやや冗長感がある？ > for(i in 4:3){ + if(nrow(test[test$a == i,]) != 0){ test[test$a == i,]$d <- 2 } #こう やれば合致しないものにはエラーはでない + } > > test a b c d 1 1 4 7 NA 2 2 5 8 NA 3 3 6 9 2 > // -cbind(test,d=rep(2,length(test$a))[!!match(test$a,4:3)]) とか,Trueは1なので,cbind(test,d=(!!match(test$a,4:3))*2) でもいいです. -- [[なかま]] &new{2005-12-28 (水) 15:05:34}; -なかまさん、ありがとうございます。例文を基にいろいろ試してみたいと思います。 -- [[じゃんぽけ]] &new{2005-12-29 (木) 09:39:05}; -「FOR文を利用して・・・」ということなので、以下のようなことをするのは反則ですよね…。すみません・・・。 -- [[舟尾]] &new{2005-12-29 (木) 10:58:54}; > ( tmp <- subset(test, a==2) ) # a==2 となる行を抽出 a b c 2 2 5 8 > ( transform(tmp, d=0) ) # データフレーム tmp に変数 d を追加 a b c d 2 2 5 8 0 > transform(subset(test, a==2), d=0) # 上の 2 行をひとまとめに -なかまさん、舟尾さん、ありがとうございました。最終的に、test <- transform(test,d=rep(2,length(test$a))[!!match(test$a,5:4)]) としました。同一のデータフレームで列dを変更したい場合、cbindだと列dがどんどん増えてしまうのでtransformを利用してみました。 -- [[じゃんぽけ]] &new{2005-12-29 (木) 13:17:31}; #comment **欠損値のあるデータに対する主成分分析 >[[Shigu]] (2005-12-27 (火) 18:21:52)~ ~ はじめまして、データ分析のためにRの勉強を始めたShiguと申します。~ googleと過去ログで調べたのですが、見つからなかったので質問させていただきます。~ ~ 欠損値がそこらじゅうにあるcsvデータ(400サンプル×130属性)に対して、主成分分析を行いたいのですが、欠損値をna.omitで省いてしまうと、4サンプル×130属性とほとんどデータがなくなってしまいます。~ このような欠損値のあるデータに対して、各属性の第一主成分の寄与の大きさを見るために主成分分析を行うにはどのようにすればいいのでしょうか？~ ~ どなたかご教授どうぞよろしくお願いします。~ // -すざまじい NA の嵐ですか。NA がどうしてそんなにたくさんあるのか（あるようになったのか）理由をよく考えないと，場合によっては，そのデータの信頼性自体がないということになるかもしれませんね。130属性の中に際だってNAの多いものはありませんか。そのようないくつかの変数を除くとサンプル数が確保できるようなら，そのようにするか。いよいよとなったら，ペアワイズ除去（cor または cov のパラメータの use で pairwise.complete.obs を指定）して求めた相関係数行列または分散・共分散行列からスタートする主成分分析を行う。 -- &new{2005-12-27 (火) 18:33:56}; -NAがゼロでも，130属性に対して400サンプルは少なくありませんか？ 項目（属性）の10倍位のサンプル数ってのがPCA（あるいは多変量解析の諸技法）のお約束では？ そもそも130属性ってのが... -- [[PCA]] &new{2005-12-27 (火) 20:07:33}; -理想的にはそういうことですが。安定性のある分析結果が得られるためにはということで。現実はなかなかそうもいかないことも。特に，卒論程度だとね。理論的には130変数なら，131サンプルあれば十分(^_^;)。なお，付け足しておくと，ペアワイズで求めた相関係数行列や分散共分散行列は非負定値行列にならないことがあるので，分析中エラーメッセージが出るかもね。 -- &new{2005-12-27 (火) 21:11:42}; -みなさまコメントありがとうございます。おっしゃるとおり卒論で解析を行っています。データは医療データで個々の患者さんに関する生活習慣・血液成分等がそろっているのですが、問診の未回答などが散在しており、NAが多くなってしまっている状況です。求まっているデータはより多く用いたほうがいいかと思っていたのですが、データの集まり具合を見て属性除去を行ったほうがいいのですね。ペアワイズ除去とそこからの主成分分析については初見ですが、教えていただいたコマンドのhelpとgoogleで調べてみます。一つ思ったのですが、他の値が埋まっているところとの関係から、NAのところを補間するようなことはRはできるのでしょうか？ また -- [[Shigu]] &new{2005-12-28 (水) 10:57:09}; -そのようなことをしてしまうと、主成分分析の意味は薄らいでしまうのでしょうか？ ぱっと思いついたことなので、的をはずしていましたらすみません。どうぞよろしくお願いします。 -- [[Shigu]] &new{2005-12-28 (水) 10:58:20}; -欠損値補完法について調べてみてはいかがでしょうか？（でも、本来ならば欠損値を埋める努力をするのが先決ですが） -- [[okinawa]] &new{2006-01-06 (金) 20:05:48}; -コメントありがとうございます。とりあえず、ある閾値以上の欠損がある属性は取り除き、それ以下のものについては他の部分の期待値を求めて全て埋めることにしました。データはもらい物ですので、そこを埋めるということは難しいので、とりあえずこれで解析してみます。 -- [[Shigu]] &new{2006-01-09 (月) 15:13:35}; -問診の欠損値ですと、尺度が名義尺度か順序尺度でしょうから、補完する場合の理由付けが難しいところですね。 -- [[okinawa]] &new{2006-01-09 (月) 15:23:40}; -そうなのです。もっと厳密な補間法もあるとは思うのですが、卒論の提出にそろそろ間に合わなそうなのでこれでがんばってみたいと思います。 -- [[Shigu]] &new{2006-01-20 (金) 12:22:06}; #comment **barplot関数でのErrorの処理 >[[鍋]] (2005-12-26 (月) 22:24:05)~ ~ > mydata [1] "15" "17" "17" "22" "11" "22" "13" "20" "19" "18" "11" "11" "18" "18" "24" [16] "13" "21" "17" "19" "11" "10" "19" "13" "6" "5" "5" "5" "3" "2" "2" [31] "0" "2" > barplot(mydata) Error in -0.01 * height : non-numeric argument to binary operator 何でこういうErrorが発生するのでしょうか？どう対応すればいいのでしょうか？わかる人教えてください！！~ // -どうやってデータを入力したのかわかりませんが、各数字がダブルクォーテーションで囲まれているのは、数字が「数」でなく「文字列」として扱われているということです。そのため、棒グラフの棒の高さを計算する際に 数値 * 文字列 という計算をしようとしたのでエラーが出たのです。 mydata <- c(15,17,17,22,11,22,13,20,19,18,11,11,18,18,24, 13,21,17,19,11,10,19,13,6,5,5,5,3,2,2,0,2) barplot(mydata) とすれば大丈夫。-- [[soysoy]] &new{2005-12-26 (月) 23:16:45}; -どうやって mydata に付値したのかを書かないと，教えようがないでしょう。取りあえずはデータを入力し直さなくても，barplot(as.numeric(mydata)) で結果を出すことはできますが。 -- &new{2005-12-26 (月) 23:38:25}; -errorメッセージはmydataがnumericになっていないと言っています。そもそもmydataはなぜcharacterになっていたのでしょうか？鍋さんはデータの型を勉強してはいかがでしょうか。また、後を考えると mydata <- as.numeric(mydata) した方が楽かも。 -- &new{2005-12-27 (火) 08:23:10}; -わかりました。ありがとうございます！！ -- [[鍋]] &new{2005-12-27 (火) 14:10:27}; -ひとりよがりだけではだめですね。「わかりました。ありがとうございます！！」だけじゃなく，回答の中に共通して逆質問があったでしょう？礼儀としても，それに答えましょうよ。あとからこのやりとりを読む人や，他の人の参考になるかもしれませんからね。 -- &new{2005-12-27 (火) 14:29:30}; #comment **"を含む文字列を作りたい >[[S. Sugawara]] (2005-12-25 (日) 14:37:16)~ ~ hoge1.sav-hoge50.savと連番になっているspssファイルをまとめてcsvにしたいのですが~ library(foreign)~ tmp="hoge"~ for(i in 1:50)~ eval(parse(text=paste(tmp,i,"<-read.spss(", """ , tmp, """ ,i,".sav, use.value.labels=TRUE, to.data.frame=TRUE)",sep="")))~ ~ のようにすると"""の部分で文法エラーになります。しかしこれをはずすとreadの中のhoge.savをファイル名と認識してくれないのでエラーとなります~ "を含む文字列を作成できれば解決すると思うのですが、その方法を御存じなら教えてください。または他の方法があればよろしくお願いします~ // -複雑な文字列を作るのには，sprintf 関数をお勧めします。以下のプログラムで出力される文字列中の引用符の前にバックスラッシュが付いて表示されるはずですが，実際の文字列中にはちゃんと引用符だけが入っているはずです。って，プログラム中と実行結果中のバックスラッシュが円記号で表示されちゃいます。。(^_^;)バックスラッシュですからね。 > tmp="hoge" > i=12 > text=sprintf("%s%i<-read.spss(?"%s%i.sav?", use.value.labels=TRUE, to.data.frame=TRUE)", tmp, i, tmp, i) > text [1] "hoge12<-read.spss(?"hoge12.sav?", use.value.labels=TRUE, to.data.frame=TRUE)" 二つ目の引用符の位置がおかしいと思いましたので直しましたが。 -- &new{2005-12-25 (日) 15:08:54}; -別法というか，このページ（目次）の中程にある「自動でファイル名を作成してファイルを読み込む場合に」を参照すると良いでしょうね -- &new{2005-12-25 (日) 15:22:49}; -解決しました。ありがとうございました。余談ですが、Rでは"もメタキャラクタになる、と解釈してよいのでしょうか? このWikiの正規表現のページには"がメタに含まれていなかったもので -- [[S.Sugawara]] &new{2005-12-25 (日) 22:07:37}; -ダブルクオートはメタキャラです. -- [[なかま]] &new{2005-12-26 (月) 11:01:10}; #comment **plot(density(x))のy軸をFrequencyにしたいです。 >[[Akira]] (2005-12-24 (土) 00:01:14)~ ~ 複数のデータに対してヒストグラムを重ねて、データの違いを視覚化しようとしています。ヒストグラムの柱の頂点を結ぶ曲線を描きたいので、plot(density(x))のy軸をFrequencyに変えたいです。~ #例えばmatというdata.frameがあるとして mat <- data.frame(rnorm(1000), rnorm(1000), rnorm(1000)) for(i in seq(ncol(mat))){ dat <- mat[, i] dat <- density(dat) dat$y <- dat$y * nrow(mat) plot(dat, xlim=c(-10, 10), ylim=c(0, 100), col=i) if(i != ncol(mat)) par(new=TRUE) } などとすればdensity(x)のy軸がFrequencyに変わると思ったのですが、うまくいきません。基本的なことかもしれませんがよろしくお願いします。~ // -5,6 行目は dat$y <- dat$y plot(dat, xlim=c(-10, 10), ylim=c(0, 0.5), col=i) ではないでしょうか？ -- &new{2005-12-24 (土) 09:51:25}; -ylim は適当に決めました。dat$y <- dat$y * nrow(mat) は density(mat) が返す mat$y がprobabilityだと思いましたので、データの要素数nrow(mat)をかければFrequencyになると思ったのですが、意味を勘違いしていますか？ -- [[Akira]] &new{2005-12-24 (土) 21:46:03}; -縦軸の名前が Density なんだけど，それを Frequency に変えたいと言うことだったんですか？それなら， plot(dat, xlim=c(-10, 10), ylim=c(0, 400), col=i, ylab="Frequency") でしょうけど。 -- &new{2005-12-24 (土) 22:03:24}; -ylab="Frequency"にすればdat$yの値も確率から頻度に変換されるのでしょうか？~ hist(x)やtruehist(x)にはそれぞれfreqとprobのオプションがあって変換できそうですが、plot(density(x))を調べてもそのような記述が見つかりません。~ やりたいこととしては、hist(x)かtruehist(x)で柱の頂点を結ぶ曲線がほしいのですが。 -- [[Akira]] &new{2005-12-25 (日) 10:37:55}; -曰く「ylab="Frequency"にすればdat$yの値も確率から頻度に変換されるのでしょうか？」~ そんな便利な機能はないですね。~ やってみればわかるように，描画される文字だけが変わるんですよ。 -- &new{2005-12-25 (日) 15:13:03}; -確かにそんな便利な機能は見たことがなかったです。そうなると先のご回答の意図がわかりません。~ で、やはりplot(density(x))のy軸は変わらないのでしょうか？bwとhが異なりますが、mat$y * nrow(mat)のplot(density(x))は、hist(x)やtruehist(x)で重ねた図と形状は似ているようです。でも、重なり部分があって視覚的にわかりにくいのが欠点です。 -- [[Akira]] &new{2005-12-27 (火) 08:17:10}; -内容をまだ理解できていませんが、[[ヒストグラムと密度の推定:http://www.okada.jp.org/RWiki/index.php?%A5%D2%A5%B9%A5%C8%A5%B0%A5%E9%A5%E0%A4%C8%CC%A9%C5%D9%A4%CE%BF%E4%C4%EA]]が関係しそうなのでリンクを張っておきます。ありがとうございました。 -- [[Akira]] &new{2006-01-16 (月) 08:29:57}; #comment **交互作用項を含む多重回帰の総当たりのAIC比較 >[[velvet]] (2005-12-23 (金) 13:13:52)~ ~ 現在以下の環境でRを勉強しているvelvetと申します~ OS:windows XP、R:Version2.2.0~ ~ 現在ある生物の出現が、どのような環境要因によって説明できるについて解析しています。~ 具体的には以下のようなデータです~ 100地点において、ある生物が存在したかどうか（A、0（出現しなかった）と1（出現した）のダミー変数です）を従属変数に、その100地点での3つの環境の変数（B、C、D、すべて連続変数です）を独立変数として、ロジスティック回帰モデルで、考えられる全ての独立変数の組み合わせのAICを比較しようと思っています。~ このような方法はweb上でいくつか存在しましたが~ 私の場合、3つの環境変数の間に相関が存在するため、これらの交互作用項を含むモデルも比較しなければなりません。~ つまり~ glm(A~B, family=binomial) から~ glm(A~B + C + D + B*C + B*D + C*D + B*C*D, family=binomial)というモデルまでのAICを計算したいわけです。~ どなたか具体的な方法を御存知でしょうか？~ 宜しく御願い致します~ // -[[総当たり法による重回帰分析:http://aoki2.si.gunma-u.ac.jp/R/All_possible_subset_selection.html]]をちょいと書き換えるか，書き換えるのが面倒なら，以下のプログラムで B, C, D, B*C, B*D, C*D, B*C*D の7つの項を独立に含む含まないと言うのを自動生成してやるとか。 name <- c("B","C","D","B*C","B*D", "C*D","B*C*D") nv <- length(name) str <- character(nv) n <- 2^nv-1 result3 <- character(n) for (i in 1:n) { k <- i m <- 0 for (j in 1:nv) { if (k%%2) { m <- m+1 str[m] <- name[j] } k <- k%/% 2 } result3[i] <- paste(c("glm(A ~", paste(str[1:m], collapse=" + "), ", family=binomial)"), collapse=" ") } print(result3,quote=FALSE) これを，sink でファイルに保存して，しかるべき時に source で読み込んで実行するとか。 [1] glm(A ~ B , family=binomial) [2] glm(A ~ C , family=binomial) : [127] glm(A ~ B + C + D + B*C + B*D + C*D + B*C*D , family=binomial) 127通りありますからねぇ。~ ああ，でも，glm(A ~ B*C*D , family=binomial) みたいなのは，ほかのと重複するんだけど，それはご愛敬と言うことで。-- &new{2005-12-23 (金) 17:12:56}; -早速のアドバイスありがとうございます。何とかなりそうです。 -- [[velvet]] &new{2005-12-23 (金) 19:39:49}; #comment **Unix上のRでのdevice出力(jpeg,png)に関して >[[サイハテ]] (2005-12-20 (火) 18:00:30)~ ~ 現在、以下の環境でRの勉強をしているサイハテと申します。~ ~ OS:Unix(Sun solaris8)~ R:Version２．０．１~ Rの使用形態：簡単なスクリプトを”R CMD BATCH”で実行する~ ~ 非常に稚拙な質問で申し訳ないのですが、Unix上のRでグラフ(plotコマンド)を出力する際、デバイスの制限があるようで、現状PDF形式でしか出力をすることができません。。。~ ~ できればR CMD BATCHで、WindowsではメジャーなJPEGやPNG形式で出力したいと考えているのですが、何か方法等ありますでしょうか？~ 教えていただけると非常に助かります。~ よろしくお願いいたしますm(_ _)m~ ~ 以下、?deviceの引用 Devices package:grDevices R Documentation List of Graphical Devices~ Description:~ The following graphics devices are currently available: * 'postscript' Writes PostScript graphics commands to a file * 'pdf' Write PDF graphics commands to a file * 'pictex' Writes LaTeX/PicTeX graphics commands to a file * 'xfig' Device for XFIG graphics file format * 'bitmap' bitmap pseudo-device via 'GhostScript' (if available). The following devices will be available if R was compiled to use them and started with the appropriate '--gui' argument: * 'X11' The graphics driver for the X11 Window system * 'png' PNG bitmap device * 'jpeg' JPEG bitmap device None of these are available under 'R CMD BATCH'. ~ // -直接の回答でなくてすみませんが，jpeg はグラフィックには不適でしょう。png はグラフィックには適していますが，拡大などの場合にジャギーが目立ちます。pdf は LaTeX や Word に取り込む場合でも，拡大してもスムーズできれいですよね。 -- &new{2005-12-20 (火) 18:31:57}; -Imagemagick をインストールし COLOR(red){convert foo.pdf foo.jpg} ではいかが。大抵の画像フォーマットは表示可能ですし、相互に変換可能です。 -- [[MKR]] &new{2005-12-20 (火) 18:35:26}; -左より，pdf, jpeg, png。プレゼンソフトに200%でペーストして，それをさらに200パーセント拡大表示して，スクリーンショットをとったもの。 #ref(hikaku.png) 実際のプレゼン画面上での pdf はもっときれい。-- &new{2005-12-20 (火) 19:44:09}; -アドバイスありがとうございました。大変参考になります<(_ _)> 丁寧なご教授本当に感謝しています。今後引き続き勉強、検討をしていきたいと思います。 -- [[サイハテ]] &new{2005-12-21 (水) 10:00:08}; #comment **ARモデルの当てはめ >[[itoken]] (2005-12-19 (月) 21:07:52)~ ~ ARモデルのあてはめをAICを使って、自分でプログラムし、次数選択したのですが、Rでもともと定義しているAICの次数選択と異なる結果になってしまいました。自分のモデルは~ (n-k)×log{Q/(n-k)}+2×k k:次数 n:データ数 Q:残差平方和~ というモデルで計算しました。~ RでのAICの式はどのようになっているのかを教えていただけると非常に助かります！よろしくお願い致します。~ // -どのライブラリを使って，どのようなプログラムで解析して，どののように異なる結果を得たというのでしょうか。質問の文面だけではさっぱりわけがわかりません。 -- &new{2005-12-19 (月) 21:45:15}; -何とか解決できました！！回答ありがとうございます！あと、もうひとつなんですが、あるデータがfに入っているとして、arモデルに当てはめるにはar(f)で当てはまるのですが、AICで次数決定して当てはめるのではなく、BICで次数決定をして当てはめたいのですが可能でしょうか？？よろしくお願いいたします。 -- [[itoken]] &new{2005-12-20 (火) 11:42:56}; -今後の参考のためにどういう問題でどう解決できたか記入していただけますと助かります。 -- [[Akira]] &new{2005-12-20 (火) 16:11:36}; #comment **Meadow＋ESSにおいて、”アンダースコア”のショートカットキーを無効にする方法は？ >[[三田]] (2005-12-19 (月) 18:18:31)~ ~ 現在、Rwikiのページを参考に~ Meadow+ESSでRのスクリプト作成を行っております。~ ~ Version:~ Meadow Emacs 21.4.1と表示されています~ ESS 2005-09-07版~ R 2.1.1~ ~ 上記のような環境において、Meadow上で変数名として＿(アンダースコア)を~ 入力したいのですが、ESSのショートカットキー~ ~ Rwikiのページ：~ http://www.okada.jp.org/RWiki/index.php?ESSUsage~ ~ において＿(アンダースコア)は自動的に”<-”に変更されてしまい、~ 何とも困っております。~ Googleや本Rwikiでも情報収集を試みたのですが、いきづまってしまいました。。ご教授いただければ幸いです(T_T)。~ // -どうしてもアンダースコアでなくてはならないのですか？こだわらないのなら，ドット（ピリオド）で代用すればよいと思います。アンダースコアは，使わない方が無難ではないでしょうか。 -- &new{2005-12-19 (月) 18:36:04}; -アンダースコアキーが付値記号にマッピングされているのは essl-s.el で定義されている仕様です。その中を詳しく見れば判りますが、アンダースコアキーを立て続けに2度入力するか、C-q _ （コントールキーと"q"を同時に押した後でアンダースコアキーを入力）すれば、アンダースコア "_" は入力できます。 -- &new{2005-12-19 (月) 19:14:22}; -諸先輩方、迅速、かつ真摯なるアドバイスを頂戴しましてまことにありがとうございました。頂きましたアドバイスを生かして今後のRスクリプト作成に挑戦していきたいと思います。この度は重ね重ねありがとうございましたm(_ _)m -- [[三田]] &new{2005-12-19 (月) 19:26:20}; #comment **検索用ベクトルの値を使ってデータフレームの条件抽出をしたい >[[小田]] (2005-12-15 (木) 10:04:39)~ ~ 一通りTips や Q&A を見てみたのですが、力不足で最適なプログラムを思いつかず投稿しています。~ どなたかご教授頂けるとうれしく思います。~ どうかよろしくお願いします。~ ~ 【実現したい事】~ 検索用ベクトルの値を使い、データフレーム内のある列で検索値に該当した全ての行を抽出したいと考えています。 >data<-data.frame(LABEL=c("A","B","C","D", "E","F","G","H","I","J"), VALUE=c(1:10)) >data LABEL VALUE 1 A 1 2 B 2 3 C 3 4 D 4 5 E 5 6 F 6 7 G 7 8 H 8 9 I 9 10 J 10 >#検索したい値 >search <- c("A","B","C") >search [1] "A" "B" "C" >#実現したいこと >data[data$LABEL =="A" | data$LABEL =="B" | data$LABEL =="C" ,] LABEL VALUE 1 A 1 2 B 2 3 C 3 いくつか方法を考えているのですが、最適なプログラムを思いつけず困っています。~ (1) 条件式をpasteで記述後、data[条件式,]を実行 > exec <- NA > for(i in 1:length(search)){ + exec <- paste(exec, ifelse(i > 1," | ",""), "data$LABEL =='",search[i],"'",sep="") + } > exec [1] "NAdata$LABEL == 'A' | data$LABEL == 'B' | data$LABEL == 'C'" > data[exec,] LABEL VALUE NA <NA> NA (2) IFELSEを使ってみた > data[ifelse(search == data$LABEL,TRUE,FALSE),] LABEL VALUE 1 A 1 2 B 2 3 C 3 Warning messages: 1: 長いオブジェクトの長さが短いオブジェクトの 長さの倍数になっていません in: is.na(e1) | is.na(e2) 2: 長いオブジェクトの長さが 短いオブジェクトの長さの倍数になっていません in: "==.default"(search, data$LABEL) どなたかご教授お願いしますm(_ _)m~ // -do.call("rbind",lapply(search,function(x){data[data$LABEL==x,]})) -- [[なかま]] &new{2005-12-15 (木) 10:17:59}; -クイズの答えの一つとして， data[apply(sapply(search, function(i) i == data$LABEL),1,sum)==1,] こちらは，表示順は元のデータファイルの出現順-- &new{2005-12-15 (木) 11:26:39}; -data[data$LABEL %in% search,] で充分か.(^_^; -- [[なかま]] &new{2005-12-15 (木) 11:37:25}; -data[is.element(data$LABEL, search),] も -- &new{2005-12-15 (木) 11:45:32}; -皆さんありがとうございます。無事解決できました！結局、search= list(label=c("A","B","C")) と置き、data[is.element(eval(parse(text=paste("data$",names(search[1]),sep=""))),search[[1]]),] と行うことで、列の指定と検索値の指定を行うようにしてみました。 -- [[小田]] &new{2005-12-15 (木) 15:41:16}; #comment **AIXへのRプログラムインストールについて >[[柴田]] (2005-12-14 (水) 17:44:56)~ ~ こんにちは。宜しくお願い致します。~ ~ OS:AIX5.2~ AIXにRをインストールしている段階なのですが、makeコマンド実行時に以下のメッセージが出力されてしまいます。~ 「make:make 1254-025 既存詳細ファイルが存在するか、またはターゲットを指定する必要があります」~ ~~ 作業内容としましては、CRANにて最新のソース「R-2.2.0.tar.gz」をダウンロードして、FTPにてAIXサーバの/home/operatorにコピーしました。それから、gunzipコマンド、tarコマンドにて上記ディレクトリにR_2.2.0というディレクトリが作成されました。そのディレクトに移動して./configureコマンドを実行した後に、makeコマンドを実行すると上記メッセージが出力されてしまします。~ この回避方法をご存知の方がいましたらご教授願いたいのですが、宜しいでしょうか。お願い致します。必要な情報などございましたらお送りいたします。~ どうぞ宜しくお願い致します。~ // -[[こちら:http://cran.r-project.org/doc/manuals/R-admin.html]]のC8セクションを御覧下さい。 少し癖(環境依存)がありますので, がんばってみて下さい. -- [[なかま]] &new{2005-12-14 (水) 17:59:16}; -早速のご回答有難うございます。がんばるしかないですね(＾＾;有難うございますm( )m -- [[柴田]] &new{2005-12-14 (水) 22:08:23}; #comment **stepAICによる変数増加法 >[[R初心者]] (2005-12-14 (水) 15:21:25)~ ~ 重回帰モデルの変数選択を行う上で、~ 変数増加法を用いたいのですが、~ stepAIC()の使い方が良くわかりません。~ helpによると、scope = list(upper=・・・~ としている辺りが増加法だと思うのですがよく理解できませんでした。~ ~ result<-glm（y~???)~ result2<-stepAIC(result,scope=list(upper=x1*x2*x3)）~ ~ とするのでは？とも思いましたが、???の部分に何を入れたら良いのか~ わかりませんでした。追加する1つ目の変数はこちらで、~ 設定しておいたほうが良いのでしょうか？~ ~ ~ どなたか、ご指導お願いいたします。~ // -最小のモデルにしておけば簡単です。 x1<-rnorm(100) x2<-rnorm(100) x3<-rnorm(100) y<-rnorm(100) g<-glm(y~1,gaussian) # 切片から始めたいなら。 library(MASS) stepAIC(g,direction="forward", scope=list(upper=~x1*x2*x3)) # scope=~x1*x2*x3でもよい。 upperの後の~を忘れないように。--[[KI]] &new{2005-12-14 (水) 17:14:48}; -例もつけて下さって、よく分かりました。どうもありがとうございました。切片からはじめれば良いのですね。早速、試してみたいと思います。 -- [[R初心者]] &new{2005-12-14 (水) 17:29:45}; #comment **数量化 I 類で総当り法による変数選択 >[[TOTO]] (2005-12-13 (火) 20:46:56)~ ~ 数量化 I 類で総当り法による変数選択（指標はAIC）をしたいのですが、どうすればよいのでしょうか？~ as.factorを使って説明変数を変換した後、mle.aic()を使ってみたのですが、アイテム変数単位ではなく、ダミー変数単位で選択されてしまいます。~ // -総当たり法ですから，全ての説明変数の組み合わせを発生させ，その都度目的の分析手法をコールして，結果を蓄積していって，最後にトップテンなどを表示させればよいのではないでしょうか。~ [[http://aoki2.si.gunma-u.ac.jp/R/All_possible_subset_selection.html:http://aoki2.si.gunma-u.ac.jp/R/All_possible_subset_selection.html]]に総当たり法の重回帰分析の実現法の一例があります。lm は説明変数が factor のときにもちゃんと手当てしてくれるので，殆どこのままでも実行できるのかな？うまくいったら教えて下さい。 -- &new{2005-12-13 (火) 21:52:18}; -紹介していただいたソースを改造し、うまくいきました。ありがとうございます！ -- [[TOTO]] &new{2005-12-14 (水) 17:23:14}; #comment **library(Hmisc); library(Biobase)だとBiobaseが読み込めません >[[Akira]] (2005-12-12 (月) 11:01:01)~ ~ bioconductorの質問はここで大丈夫でしょうか？~ 使用環境はXPpro+闇R2.2.0です。Hmiscを読み込んでからBiobaseを読み込むとBiobaseの読み込みに失敗します。逆は大丈夫でした。 # これだとエラーが出てbiobaseが読み込めません library(Hmisc); library(Biobase) # これだとOK。Biobase、Hmiscともに読み込めます。 library(Biobase); library(Hmisc) # エラーメッセージです > library(Hmisc) > library(Biobase) 要求されたパッケージ tools をロード中です 以下にエラーsetMethod("contents", "environment", function(object, all.names) { : no existing definition for function 'contents' エラー：.onLoad は 'Biobase' のための 'loadNamespace' に失敗しました エラー：'Biobase' に対するパッケージもしくは 名前空間のロードが失敗しました loadNamespaseが変更されているという意味だと思うのですが、読み込みに順番があるのでしょうか？あるとすれば、どこを調べると良いのでしょうか？~ ~ よろしくお願いします。~ // -各パッケージと R 本体とのバージョンはあっていますか（思いつきです）。 "Dependencies: R (>= 2.2.0), tools, methods" -- [[MKR]] &new{2005-12-12 (月) 11:24:11}; -library(tools)を単独で実行するとどうなりますか？ -- &new{2005-12-12 (月) 14:02:56}; -library(tools)は問題なくloadできます。library(tools); library(Biobase)とlibrary(methods); library(Biobase)はOKです。~ Hmiscのバージョンは最新（3.0-7）です。Biobaseも最新（1.7?）です。~ 他のパッケージでもこんな現象はあるのでしょうか？ -- [[Akira]] &new{2005-12-12 (月) 19:52:24}; #comment **Rをインストールしたのですが、CRANにアクセスできません。 >[[阿部勝延]] (2005-12-12 (月) 09:30:20)~ ~ R2.2.0（Windows版、Mac版）をインストールしたのですが、「R Book」P50にあるようにプロキシーサーバーの設定をしたのですが、アクセスできません。~ 状況は~ Windows Xp Professional＋R2.2.0~ プロキシーサーバーがあり、そこからファイアーウォールを抜けていくのでRのプロパティでリンク先"...Rgui.exe”--internet2と設定しています。しかし、Rを立ち上げ、ミラーサーバーの設定でJapan(Aizu)を選択して、パッケージの一覧の取得で落ちます。~ Macintosh（MacOSX10.3.9)+R2.2.0~ こちらは、ちょっと複雑で、先ほどのWindowsマシンをProxyとして、インターネットに出て行きます。~ せっていは、P50にあるとおり、~ Sys.putenv("http_proxy"="http://proxyhost:8080")~ を実行します。このコマンドは通りますが、パッケージの一覧の取得では、取得できません。~ もしいいアドバイスがありましたら、よろしくお願いいたします。~ // -プロキシサーバの名前が，本と同じく，proxyhost なんですか？ -- &new{2005-12-12 (月) 11:28:56}; -Japan(Aizu)は同時接続ユーザ数の制限がきついようですので、Japan(Tsukuba)ではどうなるでしょうか？ -- [[岡田]] &new{2005-12-12 (月) 12:49:37}; -お返事ありがとうございます。Japan(Tsukuba)はまだ試していませんでした。それと、proxyhostは、きちんと -- [[阿部勝延]] &new{2005-12-12 (月) 18:53:28}; -こちらのネットワークのホストネームにかえてあります。 -- [[阿部勝延]] &new{2005-12-12 (月) 18:54:31}; -会社でProxy越しにRを使っていますが、--internet2 を付けても外に出られませんでした。結局そのオプションは使わず、.Rprofileに options(CRAN="http://cran.md.tsukuba.ac.jp") と Sys.putenv("http_proxy"="http://プロキシーサーバ名:Port番号/") ←Port番号の後にもスラッシュを付けてます を書き加えると外に出て行けました。 -- [[船越]] &new{2005-12-12 (月) 19:22:37}; -"--internet2"オプションはproxyの設定としてWindows自体の設定を使うものなので，Internet Explorerのオプションで正しくproxy設定されていないと指定しても無駄になってしまいます．あと，proxyとして使いたいWindowsXPマシンはproxyサーバとしての動作を許可しているのでしょうか -- &new{2005-12-13 (火) 11:46:29}; -Xpの方は、Tsukubaを選択すると大丈夫でした。Macの方は、だめなので、サーバーを通さないでやることにしました。皆様ありがとうございました。 -- [[阿部勝延]] &new{2005-12-14 (水) 09:12:53}; #comment **subset= で複数の条件を指定したい >[[伊藤]] (2005-12-11 (日) 11:46:51)~ ~ ある条件に従うデータだけにおいて線形回帰をしたいと考えています。 x1<-gl(2,10) x2<-rep(gl(2,5),2) x3<-runif(20) y<-rnorm(20) lm(y~x3) このとき、{x1=1,x2=1}の場合だけを回帰したいときはどうすればよいでしょうか。~ 以下のようなことを試しましたが、うまくいきません。 lm(y~x3,subset=list(x1==1,x2==1)) #エラー lm(y~x3,subset=x1==1 && x2==1) # 全部のデータ？ lm(y~x3,subset=c(x1==1,x2==1)) # 最初の指定だけ有効？ お手数ですが、どうかよろしくおねがいします。~ // -& と && の違いに注意しましょう。それと，間違いがないことに確信が持てない場合には，途中結果を保持して，確認しましょう。 > data <- cbind(x1,x2,x3,y) > data x1 x2 x3 y [1,] 1 1 0.56633985 -0.86012923 [2,] 1 1 0.33230350 0.92239242 [3,] 1 1 0.99880629 0.39366260 [4,] 1 1 0.17977729 -0.36911252 [5,] 1 1 0.14474499 -0.74874814 [6,] 1 2 0.13511632 -0.74440885 [7,] 1 2 0.63733995 -0.28768662 [8,] 1 2 0.20289775 0.61580370 途中省略 [18,] 2 2 0.32580467 1.07307234 [19,] 2 2 0.01826360 1.16538160 [20,] 2 2 0.50006782 1.18701904 > data2 <- subset(data, x1==1 & x2==1) > data2 x1 x2 x3 y [1,] 1 1 0.5663398 -0.8601292 [2,] 1 1 0.3323035 0.9223924 [3,] 1 1 0.9988063 0.3936626 [4,] 1 1 0.1797773 -0.3691125 [5,] 1 1 0.1447450 -0.7487481 というところかなぁ。と。-- &new{2005-12-11 (日) 15:00:17}; -ありがとうございました。&と&&を完全に混同していました。Ｃ言語の癖をひきずってしまったようです。 しかし、MAC版で同じこと（subset(data, x1==1 & x2==1)）を試したところ、５×４の行列ではなく、20個のベクトルとして出力されてしまいました。これはどうしてでしょうか。-- [[伊藤]] &new{2005-12-11 (日) 18:48:44}; -MAC は使いませんが、行列属性が取り除かれている （drop=TRUE）ようなので、オプション drop=FALSE を付けてみたらどうでしょう（これが subset の既定動作のはずなので変ですが）。 -- &new{2005-12-12 (月) 09:54:43}; -最新バージョンにしましょうね。 -- &new{2005-12-12 (月) 10:09:47}; -ありがとうございます。後で試してみます。MACのほうのバージョンは2.0.0でした。 -- [[伊藤]] &new{2005-12-12 (月) 11:05:41}; -最新バージョンにしたら解決されました。ありがとうございました。 -- [[伊藤]] &new{2005-12-12 (月) 13:36:59}; #comment **一度に表示しきれないオブジェクトについて >[[Omar]] (2005-12-11 (日) 02:51:25)~ ~ 中身が多くて、スクロールバーを使っても最初の方の部分が表示されないようなオブジェクトを~ すべて見るためにはどうすればよいのでしょうか？~ 今見ようとしているのは線形モデルオブジェクトです。~ Win XP, R 2.2.0を使っています。~ // -help("page") とかしてみる。 -- [[（ただし）Linux ユーザー]] &new{2005-12-11 (日) 03:29:40}; #comment **R CMD BATCHの使い方 >[[g]] (2005-12-09 (金) 12:12:50)~ ~ 例えば，a.Rというファイルを作り，中にlibrary(Rcmdr)と記述して，コマンドプロンプト（WinXP SP2を使っています．）から， R CMD BATCH a.R と打つと，Rcommanderが起動するのですが，すぐにRも終了するので，実際は何もできません．何かオプションを加えて，Rを起動したままにすることができるのでしょうか．バッチ処理の性質上無理なのでしょうか．~ 実際は，R CMD BATCH a.Rをb.batというファイルに記述し，バッチファイルから起動しようと思ったのですが，うまくいきませんでした．~ ご意見いただけると幸いです．~ // -R Commanderを直接起動するバッチファイルを作成したい、という意図でしょうか？ -- &new{2005-12-12 (月) 14:06:11}; -実際は，Javaでパネルを表示させたいので，a.Rの中は，"自作関数名()"としたいのですが，未完成な点と，皆さんに再現して頂けるという点で，R Commanderを挙げました． -- [[g]] &new{2005-12-13 (火) 10:05:07}; -マイドキュメントの直下に.Rprofileをつくって、その中身を library(grDevices) library(Rcmdr) にすると、動作としては合っているものになるようです。-- [[岡田]] &new{2005-12-14 (水) 13:06:18}; -ありがとうございます．できました． -- [[g]] &new{2005-12-14 (水) 14:06:40}; #comment **フォントの設定を保存するには? >[[mura]] (2005-12-07 (水) 11:15:10)~ ~ OS: Mac OS X 10.4.3~ R: Version 2.2.0 (2005-10-06 r35749)~ ~ 「環境設定」で文字の色を変更すると、次回以降も有効なままでした。~ ところで、メニューバー→フォーマット→フォント→フォントパネルから文字の種類や大きさを変更したら、その場では有効なのですが、次回以降は、また初期設定に戻ってしまいます。~ フォントの設定を保存するにはどうしたらいいのでしょうか~ // -何回か設定したり再起動したりしているうちに、設定が保存されていました。お騒がせして大変申し訳ございませんでした。 -- [[mura]] &new{2005-12-07 (水) 12:08:52}; -解決済みですが、Mac版Rの環境設定ファイルは ~/Library/Preferences/org.R-project.R.plist です。Developer tools に入っているProperty List Editorで見ると中身を見たり編集できます。これが壊れていると設定が保存されなくなりますが、Property List Editorで開いてから保存しなおすと復活したりします。 -- &new{2005-12-12 (月) 14:08:57}; #comment **パッケージの更新 >[[にゃあ]] (2005-12-07 (水) 00:38:30)~ ~ パッケージの更新しようとしたら次のエラーメッセージが出てきます。~ ~ 警告：package 'lattice' is in use and will not be installed~ ~ “lattice”ってパッケージは更新できないのでしょうか?~ なお、Rは闇2.2.0で、Windows XPで使用しております。~ // -latticeのモジュールを既にロード済みなのでインストールはせんよと申しております(latticeをロードしてなくても,他にlatticeを要求するパッケージを使用していると思います)。どのような環境か(起動時に何かロードさせている?)存じませんが, 起動時のコマンドラインにオプション--vanillaを加えて起動後にインストールなりアップデートなりを行って頂けると、いいんではないかと。 -- [[なかま]] &new{2005-12-07 (水) 09:52:22}; -そうですね。MASS,graphics,grDevices,R2WinBUGS,BRugs,tcltk,car,Rcmdrが起動時にロードされる設定になってるんですが・・・。ちょっとアドヴァイス通りにやってみます。 -- [[にゃあ]] &new{2005-12-08 (木) 05:33:43}; -options(vanilla) : object "vanilla" not foundって出てくるんですが・・・・・・??? -- [[にゃあ]] &new{2005-12-08 (木) 05:46:28}; - > help(vanilla) No documentation for 'vanilla' in specified packages and libraries: you could try 'help.search("vanilla")' > help.search("vanilla") No help files found with alias or concept or title matching 'vanilla' って出てくるんですが･･････???-- [[にゃあ]] &new{2005-12-08 (木) 05:51:43}; -vanillaはパッケージの名前ではなくて起動時のoption指定のことです。R --vanillaしてください。 -- [[Akira]] &new{2005-12-08 (木) 08:35:39}; -WinXPでは，Rを起動するアイコンのプロパティを開いて，リンク先のところを，"C:?Program Files?R?R-2.2.0?bin?Rgui.exe"の後にスペースを空けて--vanillaと変更して「適用」してからOKして，そのアイコンでRを起動すればR --vanillaで起動したことになります。 -- &new{2005-12-08 (木) 10:41:47}; -無事解決いたしました。どうもありがとうございました。 -- [[にゃあ]] &new{2005-12-08 (木) 15:19:43}; #comment **MySQLとRが連携できない… >[[サイトウ]] (2005-12-05 (月) 18:02:07)~ ~ WinXPで闇R-2.2.0を使用しています。このバージョンだとRMySQLが対応していないということなのでMySQLにODBCドライバをインストールしてRと連携しようとしています。~ そこで、このサイトに載っている手順通りにDBIとRODBC双方のパッケージを読み込んでlibraryの定義をして、次の手順に進むと以下のようなエラーメッセージが出てしまいます。 > m <- dbDriver("MySQL") 以下にエラーdo.call(as.character(drvName), list(...)) : 関数 "MySQL" を見つけることが できませんでした どうしたらよいでしょうか？~ どなたかご教授お願いします。~ // -DBIとRODBCは別物のパッケージなのでは？RODBCからの接続だとデータソースを作らないといけないと思いますが・・・。 -- [[okinawa]] &new{2005-12-05 (月) 19:26:11}; -誠に申し訳ないのですが、データソースを作らないといけないというのはどういうことでしょうか？ -- [[サイトウ]] &new{2005-12-05 (月) 19:57:32}; -これはRとは関係ないのですが、MySQLやPostgreSQL,MSSQLserver,OracleなどをRDB(リレーショナルデータベース）と言います。ODBCとはMicrosoftが開発したRDBへ接続するための方法です。このODBCをRから使えるようにしたものがRODBCです。ODBCを使ってRDBへ接続する場合、RDBのドライバのインストール（例えばMySQLだとMySQLODBCドライバ)とRDBへ接続するためのデーターソースというファイルを作ります。（つなぐときのＩＤやパスワード、その他のパラメータを指定します）このデータソースが作れないとODBCでRDBへの接続はできません。 -- [[okinawa]] &new{2005-12-05 (月) 20:24:18}; -[コントロールパネル][管理ツール][データソース（ODBC)]を開いて、そこでデーターソースを作って保存します。詳しくはネットで調べてください。 -- [[okinawa]] &new{2005-12-05 (月) 20:30:04}; -ご丁寧にありがとうございます。早速、ODBCにデータソースを登録しましたが同じエラーが出てしまうのですが… -- [[サイトウ]] &new{2005-12-05 (月) 20:41:14}; -データソースは正常に登録できましたか？正常接続できたのであれば、RODBCでは odbcConnect("データソース名","ユーザー名", "パスワード","ケース") という指定になると思います。細かい指定内容はRjpWikiには無いようです。ＣＲＡＮ本家に行ってRODBCの英語のマニュアルを読むか、The R Bookには載ってます。 -- [[okinawa]] &new{2005-12-05 (月) 20:44:59}; -ケースはMySQLでは"mysql"を入れるようです。（RODBCパッケージマニュアルより） -- [[okinawa]] &new{2005-12-05 (月) 20:58:16}; -何度もすみません。okinawaさんの教えてくださっているところまでは無事に出来たのですが、やはり最初の質問と同じエラーが出てしまいます… -- [[サイトウ]] &new{2005-12-05 (月) 21:27:37}; -mysqlのクライアントで動作するか(低水準のネットワークやライブラリの問題),ODBC経由MS-Query,MS-Access等で動作するか(ODBCの問題),R上に来たら(dbDriverはODBCではないか?)等,順番に対処してください。データソースでつまずいているようですので最初に戻って順繰り詳細に考えてみましょう。(何が何に依存するか把握しないと決して幸せにはなれません) -- [[なかま]] &new{2005-12-05 (月) 22:05:44}; -まずはMySQLのテーブルをODBCでExcelから開くところあたりから練習してみたらいかがでしょうか。いきなりMySQLとRODBCではハードルが高すぎるようですね。 -- [[okinawa]] &new{2005-12-05 (月) 22:13:41}; #comment **jaccard係数によるクラスター分析 >[[ミズノ]] (2005-12-03 (土) 11:48:13)~ ~ jaccard係数を導いたのですが、その値を使ってクラスター分析のような樹形図が描けるのでしょうか？~ つまり、こちらで距離を指定して樹形図を描けるのでしょうか？~ どなたかご指導お願いします。~ // -hclust の第1引数が，"d : a dissimilarity structure as produced by dist."です。あなたが計算したものが，この定義に合致するなら，お望みのことができるでしょう。dist クラスにすることをお忘れなく。dist 関数の説明も見るといいです。 > x <- matrix(c(0,1,2,1,0,3,2,3,0),3) > x [,1] [,2] [,3] [1,] 0 1 2 [2,] 1 0 3 [3,] 2 3 0 > y <- as.dist(x) > y 1 2 2 1 3 2 3 > ans <- hclust(y) > plot(ans, hang=-1) こんな感じですね。 -- &new{2005-12-03 (土) 12:41:58}; -ご丁寧にありがとうございました。distクラスにすることで、hclust関数に引き渡せるのですね。 -- [[みずの]] &new{2005-12-03 (土) 15:30:31}; #comment **複数の生存曲線をプロットするには >[[かわぐち]] (2005-12-02 (金) 10:53:31)~ ~ 御世話になります。~ 生存曲線をプロットするには、[[Rによる統計処理:http://aoki2.si.gunma-u.ac.jp/R/index.html]]や中澤先生のRによる統計解析の本にlibrary(survival), Surv, survfitあるいはkm.survによる方法がかかれていますが、二つの群をpar(new=T)で重ね合わせようとするとx軸がずれます。~ そこでsurvfitを用いplot(res,xlim=c(0,20))とすると軸はあいますが95％信頼区間が邪魔です。~ 95%信頼区間の点線をのぞく方法はどうしたらよいでしょうか~ また、2群ではなく３群、4群となった場合にスマートにプロットさせる方法は~ ありますか。~ ~ データは”Rによる統計処理”から引用させていただきました。~ # 富永祐民「治療効果判定のための実用統計学 − 生命表法の解説 −」蟹書房(第3回改訂版）74 ページ，表 4.1 のデータ~ # 解析結果は 83-85 ページ~ # 1 は A 群，2 は B 群を表す~ > group <- c(1, 1, 2, 2, 1, 1, 2, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1)~ > event <- c(1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0)~ > time <- c(2, 20, 5, 1, 3, 17, 2, 3, 15, 14, 12, 13, 11, 11, 10, 8, 8, 3, 7, 3, 6, 2, 5, 4, 2, 3, 1, 3, 2, 1)~ // -survfit で，conf.type="none" を指定しておく。 -- &new{2005-12-02 (金) 11:59:17}; -conf.type="none" ありがとうございます。自己レスですがhelp(survfit)の例を見て下のようにしたら思ったような生存曲線が描けました。 res <- survfit(Surv(time,event)~group,conf.type="none") plot(res) ところで、グループ毎の直線を区別するためにイベント発生、打ち切り時のマークを変更するにはどうしたらよいでしょうか。 -- &new{2005-12-02 (金) 13:19:32}; -plot関数の指定をすればよいのでは？ plot(res, lty=c(1,2)) plot(res, col=c(1,2)) plot(res, lwd=c(1,2)) などでどうでしょうか？でも、pchはだめみたいですね。僕も教えてほしいです。~ それと、スクリプト部分は枠で囲んだほうが見やすいと思います。 -- [[Akira]] &new{2005-12-02 (金) 13:52:14}; -plot(fit,col=c("red","blue"),mark=c(1,2)) -- &new{2005-12-02 (金) 14:08:01}; -ありがとうございます -- [[Akira]] &new{2005-12-02 (金) 14:17:35}; -皆様ありがとうございました。 -- [[かわぐち]] &new{2005-12-02 (金) 15:35:43}; -先頭1桁に半角空白を付ければその行は薄い水色の枠内に入ります。そのような行の最後には「〜」をつけないこと。また，その前の文章の最後にも「〜」を付けないこと。前後に空白行を入れないこと。編集をクリックして出てくるウィンドウの内容を見れば，表現したいように表現するにはどのようにしたらよいかだいたいわかるでしょう。詳しくは，[[ヘルプ]]を見ましょう。このコメントは後で消去します。 -- &new{2005-12-02 (金) 18:20:16}; #comment **自動でファイル名を作成してファイルを読み込む場合に >[[蓮井]] (2005-11-30 (水) 17:50:27)~ ~ ファイル名が、1から40までの数字２つの組合せでできているファイル群があり（0101.txt、0102.txt・・0140.txt、0201.txt・・4040.txt）、~ これを順番に開いて計算結果をファイル出力する作業をしたいと考えています。~ ただ、このファイル群には、すべての組合せ分1600ファイルはなく、全体では800ファイル程度となっています。ファイル名のリストはありません。~ ~ そこで、 > for (i in 0:40) { > if (i < 10) I <- paste("0",i,sep="") else I <- as.character(i) > for (j in 1:40) { > if (j < 10) J <- paste("0",j,sep="") else J <- as.character(j) > fname0 <- paste(I,J,sep="") > fname1 <- paste(fname0, ".txt", sep="") として、ファイル名を自動的に作ってファイルを読み込もうとしたのですが、該当するファイルがない場合に「ファイル '0001.txt' を開くことができません 」となって止まってしまいます。~ ~ ファイル名がない場合に、その回をスキップするようにするには、どのようにすればよいのでしょうか？~ あるいは、別の方法を検討した方がよろしいでしょうか？~ ~ よろしくお願いいたします。~ // -邪道ならtry(読み込み処理) 真面目にやるなら dir(".","[0-4][0-9][0-4][0-9]?.txt") -- [[takahashi]] &new{2005-11-30 (水) 18:02:21}; -list.files(pattern="正規表現")とかtryとか -- [[なかま]] &new{2005-11-30 (水) 18:03:33}; -まじめかどうかわかりませんが，ちょっと変わったものを。 filenames <- outer(1:40, 1:40, function(i, j) sprintf("%02i%02i.txt", i, j)) for (str in filenames) { if (file.access(str) == 0) { ファイルが存在するので， 処理を行う } } 場合分けで paste というのは，避けたい。 -- &new{2005-11-30 (水) 19:24:09}; -無駄が無いという点では次の折衷案が良い？ 直前の例で for ( str in list.files("[:digits:]{4}.txt") ) -- [[MKR]] &new{2005-11-30 (水) 21:03:58}; -みなさま、ありがとうございます。まだ解決できないでおります、理解が遅くて申し訳ありません。~ 「ちょっと変わった」方法を試していたのですが、この場合その回のファイル名を取得するのにはどうすればよいでしょうか？~ ファイル名（の数字の部分）を変数として利用したいのです。~ あと、場合わけでpasteを避けた方が良い理由について、教えていただけますでしょうか。-- [[蓮井]] &new{2005-12-01 (木) 04:57:04}; -「ファイルが存在するので，処理を行う」の部分で，str がファイル名を表しています。paste を避けるというのは言葉足らずでしたが，sprintf だと一発で簡単に文字列を作成できますよ，ということです。~ 文字列にした後でまた数値を取り出すのは，as.integer(substr(str, 1, 2)) などですが，余計な作業なので，素直に二重ループにして，内側のループで sprintf でファイル名を作り，ファイルの存在を file.access で確認して，処理を行うときにはループの添え字もそのまま使うという具合で良いのではないでしょうかね。 -- &new{2005-12-01 (木) 08:46:11}; #comment **闇R-2.2.0でSweaveを使用したときのfont指定 >[[Akira]] (2005-11-30 (水) 13:53:30)~ ~ WinXPsp1で闇R-2.2.0を使用しています。Sweaveで画像を貼り付けたいのですが、pdf(..., cidfamily="")のような指定はどのようにすれば良いでしょうか？par(family="")でもうまくいきません。~ // -cidfamily引数はfamilyに統合されますので,Sweave中には書かないで,.Rprofile等に setHook(packageEvent("grDevices", "onLoad"), function(...) grDevices::ps.options(cidfamily="")) のように設定してください -- [[なかま]] &new{2005-11-30 (水) 15:21:41}; -なかまさま、ご返事ありがとうございました。早速、対応させていただきます。 -- [[Akira]] &new{2005-11-30 (水) 21:21:46}; #comment **R Commanderで非線形回帰したい >[[ikeike]] (2005-11-29 (火) 22:56:18)~ ~ R初心者＠勉強中です。~ R 2.2.0にパッケージRcmdr Version 1.1-2を入れていろいろといじっているのですが、~ 統計量→モデルへの適合、に「非線形回帰」を組み込む（？）~ ことは可能でしょうか？~ ~ Rでnlsを使えばよいのでしょうが、R Commanderの手軽さが気に入ってしまいました。~ ~ Fox教授のHPには、~ Fit models~ Linear regression~ Linear model~ Generalized linear model~ Multinomial logit model~ Proportional-odds logit model~ (http://socserv.mcmaster.ca/jfox/Misc/Rcmdr/)~ とあり、今のところ対応の予定はなさそうですが。~ ~ よろしくお願い致します。~ // -可能ですし、そんなに難しくはありません。適当に作ってみましたが恐ろしく使いにくい。やり方はlibrary/Rcmdr/RのlinearModel関数などをみれば解ると思いますよ。R以外にtcltkの知識が多少必要です。 -- [[takahashi]] &new{2005-11-30 (水) 19:56:32}; -takahashiさま、ありがとうございます。恐ろしく使いにくいですか。。nlsの扱いにも慣れてきたので、素直にR上で試行錯誤してみます。 -- [[ikeike]] &new{2005-11-30 (水) 22:26:06}; -言葉足らずですみません。恐ろしく使いにくいのは私のGUIのデザインセンスがまったく無いからなので(係数の初期値が扱いにくい)、うまくやれば使いやすく出来るかもしれません。 -- [[takahashi]] &new{2005-12-01 (木) 01:00:39}; -Rcmdr上のデータセットとかを使いたいという話なら以下のようなものでいけます Rcmdr.nls <- function(formula,data=NULL,start, control=NULL, algorithm=NULL, trace=NULL, subset=NULL, weights=NULL, na.action=NULL, model=NULL, lower=NULL, upper=NULL ){ UpdateModelNumber() modelValue <- paste("NonLinearRegressionModel.", getRcmdr("modelNumber"), sep="") if(is.null(data)) data <- ActiveDataSet() else data <- paste(substitute(data)) args <- paste(c(paste("start=", paste(c(substitute(start)))), paste("control=", paste(c(substitute(control)))), paste("algorithm=", paste(c(substitute(algorithm)))), paste("trace=", paste(c(substitute(trace)))), paste("subset=", paste(c(substitute(subset)))), paste("weights=", paste(c(substitute(weights)))), paste("na.action=", paste(c(substitute(na.action)))), paste("model=", paste(c(substitute(model)))), paste("lower=", paste(c(substitute(lower)))), paste("upper=", paste(c(substitute(model))))) [!sapply(list(start,control, algorithm,trace,subset,weights, na.action,model,lower,upper), is.null)],collapse=",") command <- paste("nls(", c(formula), ", data=", data, ",", args, ")", sep="") logger(paste(modelValue, " <- ", command, sep="")) assign(modelValue, justDoIt(command), envir=.GlobalEnv) doItAndPrint(paste("summary(", modelValue, ")", sep="")) activeModel(modelValue) } listNonLinearModels <- function(envir=.GlobalEnv, ...) { objects <- ls(envir=envir, ...) if (length(objects) == 0) NULL else objects[sapply(objects, function(.x) "nls" == (class(eval(parse(text=.x), envir=envir))[1]))] } Rcmdr.nls(density~a*conc,start=list(a=1)) #test code (active data set in #Rcmdr should be `DNase') -- [[takahashi]] &new{2005-12-01 (木) 11:19:33}; -takahashiさま。ありがとうございます。まだ？？な感じですが、努めて理解します。データは自分が解析しているものなんですが、確かに初期値は扱いにくいですね。 -- [[ikeike]] &new{2005-12-02 (金) 10:45:00}; #comment **分散分析 >[[oa]] (2005-11-29 (火) 22:20:26)~ ~ ~ 仮に以下のようにデータを読み込ませ、condition,sentence,primeを独立変数に、responseを従属変数にして分散分析したい場合には、どのように入力すればいいのでしょうか？ condition sentence prime response 1 1 2 1471 1 2 1 812 2 1 2 830 2 2 1 1498 3 1 2 792 3 2 1 1068 4 1 2 753 4 2 1 612 ~ // -help(aov) -- &new{2005-11-29 (火) 23:04:04}; -「単語検索」→「分散分析」 -- &new{2005-11-30 (水) 08:39:46}; -ありがとうございました。自分普段はSASを使っていまして、Rについてまだまだこれから精進するつもりです。 -- [[oa]] &new{2005-11-30 (水) 20:08:02}; #comment **ベクトルから０を取り出すこと >[[きたむら]] (2005-11-29 (火) 22:02:46)~ ~ ベクトルからある数値を取り出すには~ 例えば１を取り出すとき、~ > test<-c(1,2,3,4,5,0)~ > test[-1]~ [1] 2 3 4 5 0~ となりますが、~ ~ 0を取り出そうとすると、~ > test[-0]~ numeric(0)~ ~ となってしまいます。~ [-0]は取り出す以外の意味があるのでようか？~ ~ また、特定の値が複数ある場合、何個あるかを指定しないで~ それら全てを取り除くことは可能でしょうか？~ 例えば[1,2,0,0,0,3]→[1,2,3]~ ~ どなたか、ご指導よろしくおねがいします。~ // -test[-1] はベクトルの1番目を取り除くという意味。ベクトルの中の0を取り除きたい場合は test[test != 0] です。 -- &new{2005-11-29 (火) 22:42:28}; -ありがとうございます。全くの勘違いをしておりました。 -- [[きたむら]] &new{2005-11-29 (火) 23:23:14}; #comment **対数グラフに回帰直線と補助線を入れる方法は >[[かわぐち]] (2005-11-27 (日) 21:01:09)~ ~ 下のようなデータがあり、対数グラフを作成しました。 x <- c(3,8,11,15,16) y <- c(5,6,9.6,25.6,29.6) plot(x,y,log="y") 初めのデータy=5を基準値として横線を引き、2番目以降のデータ増加量を直線回帰して交点を計算するとx=7.38となりました。~ この2直線をグラフに追加するにはどうしたらよいでしょうか。 abline(5,0) abline(lm(log(y)~x)) ではうまく線が入りません。~ よろしくおねがいします。~ // -縦軸を実際の目盛りで描いたとき，回帰は指数曲線になることに注意。ということで， x <- c(3,8,11,15,16) y <- c(5,6,9.6,25.6,29.6) plot(x,y,log="y") x2 <- x[-1] y2 <- y[-1] ans <- lm(log(y2)~x2) x3 <- c(6, 18) y3 <- exp(ans$coefficients[1]+ans$coefficients[2]*x3) lines(x3, y3, col="red") abline(h=5, col="blue") #ref(exp.png) ってとこでどうでしょ。 -- &new{2005-11-27 (日) 23:31:41}; -望みどおりです。初めのデータの抜き方もスマートで完璧です。どうもありがとうございました。 -- [[かわぐち]] &new{2005-11-28 (月) 11:32:44}; -言葉の問題ですが、このグラフを表現するのに、データの増加量から得た回帰直線とすべきでしょうか。 -- [[かわぐち]] &new{2005-11-28 (月) 11:52:28}; -私も，そこのところで悩んだんですよ。縦軸は単に測定値なんでしょう？増加量じゃないでしょう。最初，データの差分を取ってやったりしたら，ずいぶん変な結果が出て，プロットしたら，なんだ単なる指数回帰じゃないかと思ったわけです。 -- &new{2005-11-28 (月) 12:22:49}; -では、yの値を指数回帰したところ、増加の開始時期はx=7.38であった。とすればいかがでしょうか。 -- [[かわぐち]] &new{2005-11-28 (月) 13:40:13}; -本当にデータ数は５つだけですか。それとデータは，水平な部分と指数的に増加する部分に区分できるというのは本当なのでしょうか。更にいえば、指数関数で回帰するなら、対数をとって直線回帰せずに、元のデータに直接指数関数を非線型回帰（例えば nls 関数を使い）するのが「普通はベター」なんですが。 -- [[QDU]] &new{2005-11-28 (月) 16:54:27}; -前に引き続いて， ans2 <- nls(y2~a*b^x2,start=list(a=1.1,b=1.2)) p <- ans2$m$getPars() y4 <- p[1]*p[2]^x3 lines(x3, y4) として，比較してみればよろしいかと。 -- &new{2005-11-28 (月) 21:37:33}; -御指摘ありがとうございます。もう少し具体的に申しますと治療データで、x軸は月数です。初めの治療が行われた後，再発が発見されるまででありこの間のデータ数は5個だけです。ちなみに2回目の治療後y=5.6となり今回の基準線と考えました．データの推定や解釈に大きな間違いがありましたら御指摘をお願いします。 -- [[かわぐち]] &new{2005-11-28 (月) 21:40:24}; #comment **学習セットと検証セットの範囲を切り分ける方法 >[[大軽貴典]] (2005-11-25 (金) 10:20:00)~ ~ はじめまして。~ Rでニューラルネットワークなどを使ってみようと思っています。~ いま、手元にデータセットがあります。このデータセットは事前にランダムに10個のグループ（各グループは同数）に分けております。グループの情報はカラムを追加させて1,2,..,10というふうに入力しています。~ これを用い、90%を学習セットに用い、残りの10%検証用に使おうと思っています。~ 一度だけなら分けて入力すればよいのですが、グループカラムの情報を使って、１回目はグループ1を学習セットに、２回目はグループ2と繰り返し最終的に１０回の作業をしようと思っています。~ ~ いきなりすべての工程をコード化できなくてもよいので、うまくあるデータセットからグループカラムの情報を使って学習セットと検証セットの切り出しをする方法はないでしょうか？~ ~ データセットは次のような構成になっています。~ 目的変数 グループ 説明変数１ 説明変数２ ３ ４ ５ ...~ A 1 1 -3.4~ B 3 0 3~ C 9 0 5.2~ B 10 1 7~ ~ どなたかアドバイスお願いできないでしょうか？~ よろしくお願いいたします。~ // -10％が検証セットですよね。以下のようにするのではダメでしょうか？ for (i in 1:10) { 学習セット <- データフレーム[データフレーム["グループ"] != i, ] 検証セット <- データフレーム[データフレーム["グループ"] == i, ] 何らかの分析 } どうでしょう。~ subset 関数を使うべし。以下のようにする方がモアベター。 for (i in 1:10) { 学習セット <- subset(データフレーム, データフレーム["グループ"] != i) 検証セット <- subset(データフレーム, データフレーム["グループ"] == i) 何らかの分析 } ですね。 -- &new{2005-11-25 (金) 10:22:04}; // -早速のアドバイスありがとうございます。教えていただいた方法でたしかに学習セット、検証セットはできるのですが、分析のときにNAがあってだめみたいなエラーが出ます。おそらくグループをまびいたのだけど行番号も減っておらず、分析時にあいだあいだにNAがあると認識しているのかなと思っています。とりあえずsvm(e1071)では、n-fold cross validationのオプションがあるので取り急ぎのデータ解析はオプションで対応しました。ですが、教えていただいたものを利用してsvmをかけることはできないでしょうか？教えていただいたsubsetを使う方が汎用性が高いと思います。 model<-svm(subset(training,select=-Object),training$Object) のような命令でNAがあるみたいなエラーが出ています。なお、trainingは、教えていただいた方法で作成した学習セット、Objectとは名義型の目的変数です。 -実際に NA があるのではないですか。というか，subset(training,select=-Object) と元のデータを比べてみて下さい。[[データフレームTips大全]]を参照しておいて下さい。 -- &new{2005-11-26 (土) 01:22:34}; -[[データフレームTips大全]]でわかりやすく教えていただきありがとうございます。無事に解決いたしました。 -補足すれば，どのようなことがあったために障害が生じていたのかということをお知らせ頂くと，今後同じことで悩む人がいないのではないかと思うわけです。 -- &new{2005-11-26 (土) 21:58:23}; #comment **PDFとTeX >[[aa]] (2005-11-23 (水) 17:04:35)~ ~ はじめまして。~ Rで作成したグラフィックスをPDFで書き出し，TeXでPDFを読み込む際に，エラーが発生します。~ 対応策をご存知の方，ご示唆を頂けると幸いです。~ ~ 【作成例 R】 setwd("C:/question20051123") #ディレクトリ plot(1:10, main="日本語の入ったPDF") #画像を作成 #画像 をPDFで保存 dev.copy(pdf, "fig1.pdf", cidfamily="Japan1Ryumin", width=12, height=5) dev.off() graphics.off() 【作成例 TeX】~ ?documentclass[12pt,a4paper]{jarticle} ?usepackage[dvips]{graphicx} ?begin{document} %図の挿入 ?begin{figure}[htbp] ?begin{center} ?includegraphics[keepaspectratio=true, height=100mm]{fig1.pdf} ?end{center} ?end{figure} ?end{document} コンパイル時のError Messageは， ! LaTeX Error: Cannot determine size of graphic in (ファイル名) (no BoundingBox) 当方の環境は以下の通りです。 OS: Windows XP R: 闇R2.20 TeX: LaTeX2e Ghostscript: Ghostscript 8.52 // -pdf ではなく，eps ファイルでは駄目なのでしょうか？？ -- &new{2005-11-24 (木) 01:00:16}; -大量の画像データを扱う場合，PDFでコンパイルした方が早いのです。ですので，可能であれば，R→PDF→TeXの連携を上手くしたいと思っています。 -- [[aa]] &new{2005-11-24 (木) 01:08:28}; -コマンドプロンプトから、 ebb graph_name.pdf というコマンドだったはず。-- [[蓮見]] &new{2005-11-24 (木) 01:23:04}; -蓮見さんのご指摘のように，bb ファイルがないのが原因でしょうね。ただ，Windows の場合，bb ファイルの作成コマンドは bmc -b graph_name.pdf だったような。ほかには mediabb.sty を利用するという方法も。いずれにせよ，R で作成した pdf が pdf として問題のないものになっているのであれば，R ではなく TeX の問題ですから，TeXWiki あたりを参照されては。 -- &new{2005-11-24 (木) 02:11:14}; -何がなくとも grep MediaBox hoge.pdf 等とすれば.BBoxの値は取得出来ると思います. -- [[なかま]] &new{2005-11-24 (木) 10:02:55}; -[[mediabb.sty:http://www.ns.musashi-tech.ac.jp/~inoue/Pages/TeX/mediabb.sty.html]]を使うことで，楽な解決ができました。ご紹介頂きありがとうございました。 -- [[aa]] &new{2005-11-25 (金) 17:48:17}; -解決済みのようですが、ご質問には[[Sweave:http://www.okada.jp.org/RWiki/index.php?cmd=read&page=Sweave%20%A5%DE%A5%CB%A5%E5%A5%A2%A5%EB%CC%F5&word=sweave]]を使っても良いと思います。以前、ここで質問させていただいて、[[関連ページ:http://www.okada.jp.org/RWiki/index.php?cmd=read&page=SweaveTips%BD%B8&word=sweave]]を作っていただきました。 -- [[Akira]] &new{2005-11-25 (金) 18:49:21}; #comment **イベントの少ない場合の生存分析 >[[かわぐち]] (2005-11-22 (火) 16:09:00)~ ~ はじめまして。最近Rを知って感激して統計とともに勉強中です。~ さて、生存分析にて型どおりKaplan-Meier,log-rank,Cox回帰分析をしたいと思います。例数は20ー30ほどですが、イベント発生数が群別で0のものや2ー3しかないものがあります。log-rankでは有意差があるのですが、Cox回帰はイベント発生が5個以上必要とのことです。このような場合の群別の比較はどうしたらよいのでしょうか。~ よろしくおねがいします。~ // -Rに限らないことですが、例数が少ない場合で無理にCoxを実施しても意味のある結果にはならないと思います。研究自体を例数を増やしたデザインに改変するか、log-rankの結果の範囲で言える結論のみにとどめて判断するべきでは。 -- &new{2005-11-26 (土) 11:27:12}; -ありがとうございました。イベントが少ないのは臨床データなので増やすと言うわけにもいきません。分割表分析に切替えて検討してみます。 -- &new{2005-11-27 (日) 11:44:43}; #comment **300行400列の行列をの各セル比較したい >[[ねこのて]] (2005-11-16 (水) 23:37:42)~ ~ 300行400列のテキスト形式の行列が２つあります。~ （細密数値情報と言うデータで、仕様は~ http://www.gsi.go.jp/MAP/CD-ROM/saimitu/htmls/format.html にあります。）~ ~ この２つのファイルの同じセルについて、数値（土地利用の分類コード）が何から何に変化したかを、場合毎に数えたいと思っています。~ （１から２がmケース、１から３がnケース、、、、、）~ ~ Ｒでのデータハンドリングの経験がなく、探せる範囲ではヘルプも探しましたがよくわかりませんでした。~ ~ 効率的な方法をお願いいたします。~ // -非定型的な問題を解決するには，自分で関数を書かねばならないでしょうね。~ 効率的であるかどうかは二の次のような気がします。~ 正確な答えを得るために，いかに素早くコーディングできるかでしょうね。~ コンピュータがいかに短時間に答えを出せるかは，通常は問題外でしょう。~ というのは，よっぽどまずいプログラムでなければ，それなりの時間内に答えを出せるでしょうから。~ というわけで，以下は回答例その1ということで。 > set.seed(1234) > x <- matrix(sample(1:5, 12, replace=TRUE), 4, 3) > y <- matrix(sample(1:5, 12, replace=TRUE), 4, 3) > # プログラムコンテストは，ここから開始 # 以下で 1000 は，適当な数。必要なら10000でも100000 でも > ans <- table(as.vector(x)*1000+as.vector(y)) > z <-as.integer(names(ans)) > ans <- cbind(floor(z/1000), z%%1000, ans) > colnames(ans) <- c("before", "after", "frequency") > ans before after frequency 1001 1 1 1 1002 1 2 1 2002 2 2 1 3001 3 1 1 3002 3 2 1 4001 4 1 1 4002 4 2 3 4005 4 5 2 5002 5 2 1 代替案がいくつか出た段階で，どれが一番計算時間が短いかというコンテストがあるかもしれませんが，そのコンテスト結果が出るまでには，第一案なり何なりが既に答えを出していることは間違いない。（答えが正しいかどうかは別問題。正しい答えを出さないプログラムは問題外）~ 代替プログラムを書こうと思う方は，このプログラムの4行目以降を書くとよろしいかも（1，2行はテストデータの生成） -- [[青木繁伸]] &new{2005-11-17 (木) 00:08:45}; -ありがとうございます~ 書いてあることは理解できたので、やってみようとしたのですが、~ 以下にエラーvector("integer", length) : 指定されたベクトルのサイズが長すぎます、となってしまい止まってしまいます。 -- &new{2005-11-17 (木) 01:57:31}; -おかしいですね。入力間違いではないですか？ -- [[青木繁伸]] &new{2005-11-17 (木) 09:18:49}; -各行は最初の7カラムがID、あとは2カラム×400なので > x <-read.fortran("S4_2617.TDU", c("i7", "400i2")) > dfx <-data.frame(x) > names(dfx) <-c("ID", paste("var", 1:400, sep="")) > dfx$ID <-NULL として読込み~ 教えていただいた > ans <- table(as.vector(x)*1000+as.vector(y)) を実施したところ、 > 以下にエラーvector("integer", length) : 指定されたベクトルのサイズが長すぎます、 となって実施できませんでした。 どこが間違っているのでしょうか？ お手数おかけしますがご指摘いただければ幸いです 。 -あなたのプログラムでは，x，y（dfx，dfyも）がデータフレームですね。データフレームに as.vector は期待する結果を返しません。ans <- table(as.vector(as.matrix(dfx))*1000+as.vector(as.matrix(dfy))) としてください(x,y を使うと ID も含まれますよ）。 -- [[青木繁伸]] &new{2005-11-17 (木) 11:03:25}; -どうもありがとうございました。試してみましたら無事カウントすることができました。もう少し気合を入れて勉強します。本当にどうもありがとうございました。 -- [[ねこのて]] &new{2005-11-17 (木) 11:36:40}; -ハンドル名でも「ねこのて」はまずいのではと。 -- &new{2005-11-17 (木) 17:35:25}; -何でまずいのか，私にはわかりませんねぇ(^_^) -- [[青木繁伸]] &new{2005-11-17 (木) 22:04:07}; #comment **共分散関数acf >[[Aki]] (2005-11-16 (水) 00:46:20)~ ~ 自己共分散関数を求めるacfで1次元の配列ではなくn×mの2次元配列を引数にもってくると~ グラフにはm×mのグラフが表示されます。~ これの非対角要素はccfで求められる相互共分散と同じなのでしょうか？~ // -ヘルプを読みましょう。~ Description~ ~ The function acf computes (and by default plots) estimates of the autocovariance or autocorrelation function. Function pacf is the function used for the partial autocorrelations. Function ccf computes the cross-correlation or cross-covariance of two univariate series. -- &new{2005-11-16 (水) 09:29:06}; -ヘルプを読んだつもりになっていました。~ acfは自己共分散関数と自己相関関数を計算するの理解できるのですが、自己であるのならグラフがm個表示されればよい気がします。~ しかしV1&V2といったグラフが出てくるのでこれは相互共分散関数or相互相関関数では？と思ったのですが、これはいったい何なのでしょう？~ 計算結果を見る感じではacfの非対角の要素がccfで計算したものと同じになっているように見えます。 -- [[AKi]] &new{2005-11-16 (水) 11:41:23}; -ソースが見られるときには，ソースを見るのがよいですね。~ > ccf function (x, y, lag.max = NULL, type = c("correlation", "covariance"), plot = TRUE, na.action = na.fail, ...) { type <- match.arg(type) if (is.matrix(x) || is.matrix(y)) stop("univariate time series only") X <- na.action(ts.union(as.ts(x), as.ts(y))) colnames(X) <- c(deparse(substitute(x)), deparse(substitute(y))) acf.out <- acf(X, lag.max = lag.max, plot = FALSE, type = type) ということで，引数をまとめて，acf を呼んでいるので acf と同じ結果になるわけでしょう。 -- &new{2005-11-16 (水) 13:02:30}; -なるほど、だいたい理解できました。~ いくつかの相互共分散をまとめて出したいときはccfでひとつずつ計算するより、acfでまとめてやったほうが楽ということですね。~ ありがとうございました。 -- [[AKi]] &new{2005-11-16 (水) 15:54:18}; #comment **ヒストグラムのy軸に中断線を入れたいです >[[Akira]] (2005-11-14 (月) 10:05:32)~ ~ ヒストグラムの作図においてy軸に中断線を入れることはできますでしょうか？"中断線"では検索できず、他に良いキーワードを知らないので、上手く見つけられませんでした。~ 例えば、y軸が1〜100の範囲で1〜20と80〜100で作図したいと思っています。~ よろしくお願いします。~ // -パッケージplotrixに、axis.breakという関数があります。これで上記の内容が実現できるかどうかは不明です。 -- &new{2005-11-14 (月) 10:52:17}; -example を実行させると，どうも期待する物ではないようですね。質問者がいうような図は，あまりお勧めではないので，そのような機能を持つ関数も提供しないというのが大人の態度なんでしょうね。 -- &new{2005-11-16 (水) 09:34:06}; -ご連絡ありがとうございます。よい作図でないのであれば仕方がないと思います。 -- [[Akira]] &new{2005-11-16 (水) 09:52:15}; -できないこともないようです。Google で "r-help axis break" で検索すると色々出てきます。 -- [[sei]] &new{2005-11-17 (木) 23:23:31}; -プログラムを書けば，どんなことでもできますね。できると言うことと，簡単にできるということと，そのようなことをしても良いと言うことはそれぞれ区別があるでしょうということです。 -- &new{2005-11-17 (木) 23:50:09}; #comment **両対数のヒストグラム >[[muramatsu]] (2005-11-11 (金) 15:48:14)~ ~ あるデータのヒストグラムを描きたいのですが、ばらつきが大きいので、両対数のヒストグラムを描きたいと思っています。~ "ヒストグラムと密度の推定"という記事をみましたが、両対数での描き方は載っていませんでした。~ ~ plot関数と同様に、 hist(data, log="xy") ともやってみましたが、 hist.default(data, log = "xy") : 'x' は数値でなければなりません というエラーが出てしまいました。しかし、ちゃんと数値データです。~ plot(table(data)) ではきちんとヒストグラムが描けています。~ ~ また、ヒストグラム作成時に区切り幅も対数で決める必要があると思いますが、それはとりあえずおいといて、 plot(table(data), log="xy") のようにやってみました。すると、 1: 軸の限界が有限ではありません [GScale(-1.#INF,2.18469,2, .); log=1] 2: Internal(pretty()) の範囲が大きすぎます.. 修正しました という警告がでて、plotに失敗しました。~ ~ どうしたらよいでしょうか？よろしくお願いします。~ // -「両対数のヒストグラム」って，どんな奴なんですか。~ 「両対数」というのは，二変数ともに対数を取るということで，両対数プロットというと，x軸，y軸を対数目盛りにした散布図のことですよね。~ data というのはベクトルデータですよね。plot(table(data)) でちゃんとヒストグラムが描けているというのがわからん。完全な連続値じゃなくて同じ値が多いデータかな。 x <- round(rlnorm(1000),0) plot(table(x)) みたいな感じかな。~ で，これを横軸だけを対数軸で（つまり片対数で）ヒストグラムを描くなら， x <- rlnorm(1000) hist(log(x)) # hist(x) と比較せよ みたいなのでよいのでは？横軸の目盛りも，対数目盛りになりますよ。~ しかし，対数目盛りも付け方に二通りのやり方があって，対数変換した後の数値を目盛る場合と，対数変換する前の数値を目盛るやり方。~ 示した方法では，前者になっている。後者のやり方は，[[どこか:http://aoki2.si.gunma-u.ac.jp/R/log-plot.html]]を参照する。 z <- x <- rlnorm(1000) hist(log(x), xlab="") z <- floor(log10(z)) z2 <- 1:10*10^(log.min <- min(z)-2) if ((n <- max(z)-log.min) > 0) { for (i in 1:n) { z2 <- c(z2, z2*10^i) } } log.z2 <- log10(z2) axis(1, at=log.z2, labels=z2, pos=-30) みたいな感じかな。描かれる横軸二種。上が変換後の値の数値に基づいて目盛った目盛り，下が変換前の数値に基づいて目盛った目盛り。（位置は決めうちしているので，汎用性はない） -- [[青木繁伸]] &new{2005-11-11 (金) 18:01:49}; #ref(pdf.png) -返信をしたいのですが、コメントの挿入のテキストボックスは、一行しかかけません。どうしたら複数行のコメントをかけるのでしょうか？ -- [[muramatsu]] &new{2005-11-24 (木) 22:27:09}; -このページの一番下に紙と鉛筆のアイコンがありますので、それをクリックしますと、このページ全体の文書データが表示されます。それを間違って消さないようにとても注意をしながら編集します。最後にページの更新ボタンを押すともとの画面に戻ります。 -- [[okinawa]] &new{2005-11-24 (木) 22:58:33}; -回答ありがとうございます。~ ~ ＞「両対数のヒストグラム」って，どんな奴なんですか。~ ~ 縦軸、つまり度数も対数軸で書きたいと思っています。~ 度数にすごくばらつきがあるのです。~ dataは [1] 83.553719 80.760000 100.750000 .... のようなベクトルデータです。~ 以下の方法で、確率分布を描くことはできたのですが、 dist <- density(data) plot(dist, log="xy") 下のようにヒストグラムを描こうとすると、 dist <- table(data) plot(dist, log="xy") やはり次のエラーが出ます。 1: 軸の限界が有限ではありません [GScale(-1.#INF,2.18469,2, .); log=1] 2: Internal(pretty()) の範囲が大きすぎます.. 修正しました 次のようにしても、 hist(data, log="xy") このようなエラーが出ます。 高水準 plot 関数で，パラメータ "log" を設定できません 両対数ヒストグラムを描くにはどうしたらいいのでしょうか？ -- [[muramatsu]] &new{2005-11-24 (木) 23:15:28}; -hist2 <- hist.default とし，hist2 の最後の方にある r <- structure( なんだらかんだらという行を r <- structure(list(breaks = breaks, counts = log10(counts), intensities = dens, density = dens, mids = mids, xname = xname, equidist = equidist), class = "histogram") に書き換える(counts=counts を counts = log10(counts) にする）。 次に上に示されてことを縦軸にも施す。さっき作った hist2 を呼ぶことに注意。 old <- par(mar=c(5,5,1,1), xpd=TRUE) z <- x <- rlnorm(1000) # hist ではなく hist2 を！ frq <- hist2(log(x), xlab="", ylab="") # 横軸を描く z <- floor(log10(z)) z2 <- 1:10*10^(log.min <- min(z)-2) if ((n <- max(z)-log.min) > 0) { for (i in 1:n) { z2 <- c(z2, z2*10^i) } } log.z2 <- log10(z2) axis(1, at=log.z2, labels=z2, pos=-0.27) # 縦軸を描く z <- floor(frq$counts) z2 <- 1:10*10^(log.min <- min(z)-2) if ((n <- max(z)-log.min) > 0) { for (i in 1:n) { z2 <- c(z2, z2*10^i) } } log.z2 <- log10(z2) axis(2, at=log.z2, labels=z2, pos=-4.4) par(old) 2つの軸を描く部分が冗長なので，一つの関数にまとめるとよい。このようなヒストグラムは見たことない。 -- &new{2005-11-25 (金) 00:24:01}; #ref(pdf2.png) #comment **補正AICについて >[[山形]] (2005-11-10 (木) 23:16:59)~ ~ モデル選択をする場合、AICを基準にモデル選択をするのは知ってますが、~ 扱っているサンプルが少ないため、AICではモデルを過大評価してしまうことが~ あると聞きました。~ その場合、補正したAIC(c-AIC)を使うといいということも聞きましたが、この~ 補正AICをRで行なうことは出来ますか？探してみたのですが、力不足で探せませんでした。よろしくお願いします。~ // -[[stepAIC():http://finzi.psych.upenn.edu/R/library/MASS/html/stepAIC.html]]の引数にkというのがあります。AICと補正AICの式を比べて適当に設定すればいいのでは。 -- [[aki]] &new{2005-11-17 (木) 01:54:44}; #comment **時系列データの生成について >[[yuta]] (2005-11-09 (水) 15:41:12)~ ~ 表題の件で質問です。~ ~ 以下に例示するような、サンプリングのインターバルが不定で欠損値を含む時系列データを処理する方法はないでしょうか？~ なお、行列testmatの１列目は日付のシリアル値（windowsで使用する、1900/1/1＝1とする整数値）です。単純にこの列を列ラベルにしたりしてts()で処理することはできるのですが、インターバルについて思うように処理できていないように思います。 test<-c({38482,NA,NA,NA,NA,NA, 38531,NA,0.992,NA,NA,NA, 38552,NA,0.988,0.994,0.990,NA, 38562,NA,0.993,0.995,0.992,NA, 38576,NA,0.989,0.990,0.994,NA, 38595,NA,0.994,0.996,0.990,NA, 38608,NA,0.931,0.992,0.990,NA, 38621,0.985,0.956,0.991,0.991,NA, 38637,0.984,0.670,0.980,0.976,NA, 38649,0.989,0.879,0.972,0.983,NA }) testmat <- matrix(test,c(10,6),byrow=T) ts(testmat) 使用バージョンはR2.2.0（闇国際化）、OSはWindows2000です。 また、時系列オブジェクトの生成には基本パッケージのts()を使用しています。 宜しくお願いします。~ // -何で，c({データ...}) なんだろう。私の環境ではエラーになるんだけど。 -- &new{2005-11-09 (水) 18:16:34}; -tseriesのirtsではダメですか？ -- &new{2005-11-09 (水) 19:16:33}; -ありがとうございます。当方の力不足で探し出せていませんでした。helpを読んでトライしてみます。 -- [[yuta]] &new{2005-11-09 (水) 19:52:57}; #comment
