初心者のための R および RjpWiki に関する質問コーナー

新規投稿はできません。

encodingをグローバルに設定したい†

akira (2007-12-31 (月) 22:26:19)

WinXPとubuntuを並用していています。ubuntuでSHIFT_JISを指定せずに読み込んで、grepするとエラーになってしまいます。
仮にhogeの中身を

test\n
テスト\n
test\n

とします。

事前にencodingがわかっていれば、

> con <- file("hoge.tsv", open="r", encoding="SHIF_JIS")
> x <- readLines(con)
> close(con)
> grep("hoge", x) <- SHIFT_JISを指定しないと In grep("hoge", x) :  入力文字列 2 はこのロケールでは不適切です となります
1 2

としています。

これを、UTF-8とSHIFT_JISを区別せずに使いたいのですが、グローバルに設定を変える方法、もしくは検索キーワードを教えてほしいです。

localeのLC_CTYPEを変更すれば良いと思うのですが...

> sessionInfo()
R version 2.6.1 (2007-11-26)
i486-pc-linux-gnu

locale:
LC_CTYPE=ja_JP.UTF-8;LC_NUMERIC=C;LC_TIME=ja_JP.UTF-8;LC_COLLATE=ja_JP.UTF-8;
LC_MONETARY=ja_JP.UTF-8;LC_MESSAGES=ja_JP.UTF-8;LC_PAPER=ja_JP.UTF-8;
LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=ja_JP.UTF-8;
LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

loaded via a namespace (and not attached):
[1] rcompgen_0.1-17
↑

計数データのPoissonGLMによるモデルからカウントされる0の度数を求めるコマンドについて†

idaten (2007-12-27 (木) 04:25:41)

0が多い計数データ、solder.balanceというデータをGLMでPoisson回帰しました。

> P=glm(formula = skips ~ ., family = poisson, data = solder.balance2)

そのモデルからカウントされる0の数を調べたいときに、色々調べていましたら、

"the obserbved zero counts are compared to expected number of zero counts for the likelihood model:"
> round( sum( dpois( 0,fitted(P) ) ) )

として得られると書いてありました。
fitted(P)は推定されたPoisson分布の平均パラメータλですから、このコマンドの意味は各λで0をとる確率の和ということになると思うのですが、これでなぜ0をカウントすることができるのでしょうか？対数尤度関数として考えるのでしょうか？
そもそも訳が間違っていて、このコマンドは0を得るのでしょうか？
もし情報不足、不適な記述がありましたら、申し訳ございません。指摘していただけると助かります。

↑

マッピングのためのgridの設定†

R 初心者 (2007-12-23 (日) 22:13:19)

こんにちは。
空間統計学のためにRを使っている者ですが、次の質問に対して皆さまのヘルプを頂けると嬉しいです。

アメリカ(127.5W-67.7W, 24N-52N)領域を4x5のgridに区切り、それぞれの値をgridごとに計算し、その値の空間時間分布を解析しようとしています。
すなわち、gridの中央点は070W42N,070W46N,...となります。
しかし、この地理情報をshape fileとして与えなければ絵を描くことができません。
そのようなgrid情報はArcIGSで書かなければならないのでしょうか。
もしそうであれば、ArcGISを一度も使ったことがないものに対して、簡単なアドヴァイスを頂けると嬉しいです。
どうかよろしくお願いします。

sp-grid-20080112.png
↑

正規分布とロジスティック分布の密度関数の重ね描きは？†

idaten (2007-12-21 (金) 04:26:30)

　正規分布N(0,1)とロジスティック分布の密度関数の違いを図で確認するために、この2つを一つの図で重ねがきしたいのですが、どのようにしたらいいのでしょうか？それぞれの密度関数はかけても、x軸の範囲が違ったりするので、一つの図できれいに表現することができません。どなたかコマンドを教えていただけると助かります。よろしくお願いします。

↑

あるベクトルから条件に合致した要素を取り出す際に、条件が複数の場合はどうすればいい？†

AobaAoba (2007-12-20 (木) 18:24:29)

どうぞよろしくお願いします(WinXPで、バージョン2.5.1使用です)。

下記のxのベクトルのうち、yのベクトルの要素に含まれるものだけを取り出して、zというベクトルを作りたいのです。

x <- c(2,2,3,6,2,0,6,7,9)
y <- c(2,7)

そこで、単純には以下のようにすれば良いことは初心者ながらにわかります。

z <- x[x == 2 | x == 7]

しかし、yがもっと長いベクトルの場合はどうすれば良いのでしょうか?
私が思いつけるのは「forループを使って、yの要素を1つずつ取り出してxにあるかどうかチェックし、ベクトルを結合していく」という以下のような方法です。

z <- c()
for(i in 1:length(y)){
z <- c(z,x[x == y[i]])
}

もっと簡単な方法があればよろしくお願いいたします。

↑

データフレーム名の自動生成or変換†

うち (2007-12-20 (木) 17:13:04)

irisのデータを例にして説明します。
種名ごとに抽出したデータフレームを自動で生成し，そのデータフレーム名を種名と対応させたいと考えております。
下記のような方法を試してみましたが，当然ながら全くダメでした。データフレーム名の自動生成方法，またはデータフレーム名変換方法を御教示いただければ幸いです。

種名 <- levels(iris$Species) #種名を抜き出す
種数 <- length(種名) #種数を取得
for(i in (1:種数) ) {
  paste(種名[i]) <- subset(iris, Species==種名[i]) 
}
#　「関数 "paste<-" を見つけることができませんでした」とエラーが出る。
for(i in (1:種数) ) {
  df <- subset(iris, Species==種名[i])
  names(df) <- paste(種名[i]) 
}
#　これだとデータフレーム名はそのままで，列名が変わってしまう。
↑

ARIMAモデルの次数選択†

kuri (2007-12-20 (木) 02:14:02)

ARIMAモデルの次数選択をAICを選択基準として求めるプログラム(すべての組合せを行い、最小のAICのときの次数を返す)を作ってみたのですが、途中でエラーが出て停止してしまいます。
エラーが出ても無視して、全通りを実行することはできますか？

↑

微分できない関数のoptim()での最適化†

Katoh (2007-12-18 (火) 12:37:52)

はじめまして。
optim()関数について質問です。

こちらの「汎用最適化関数 optim() の使用法」ページによると、
＞微分できないような関数では"Nelder-Mead" 法 や"SANN" 法 しか選択肢は無い。
とありますが、
実際に準ニュートン法や共役勾配法でも微分できない関数の最適（と思わしき）値を見つけることができます。
どのようにして微分値を求めているのでしょうか？それとも求めずに何か他の手を使って最適化を行っているのでしょうか？

よろしくお願いします。

↑

barplotで、特定のbarだけをかえたいのですが†

[[<ふ>]] (2007-12-15 (土) 22:10:24)

たとえば、
x <- matrix(seq(10,20,2),6,2)
のようなデータを用意して、

barplot(t(x))
とすると、各barは、同じ二色（defaultでは、glay scaleですが）で塗り分けられます。

これを、（例えば）右から二本目のbarの色を、"lightblue","green"に変更するには、どのようにしたらいいでしょうか。

（bar内が分割されてない場合は、そのbarだけ色を変えることはできています。）

よろしくお願いいたします。

↑

ftableで作成したクロス集計で合計がゼロの行を消したい†

Transalp (2007-12-14 (金) 18:32:49)

お世話になります。
要素が４つ以上あるデータフレームの値をtableで集計し、ftableしてクロス集計表を作りましたが、行合計がゼロのデータが集計表に表示され、困っております。

hoge.txtの中身
所属	名前	月度	出席
営業部	営業太郎	01月	はい
総務部	総務花子	02月	はい
営業部	営業次郎	03月	はい
購買部	購買好男	03月	はい
経理部	経理京子	04月	はい
技術部	技術冴子	04月	はい
営業部	営業太郎	04月	はい

Rコンソールから

df <- read.delim("hoge.txt", header=TRUE)
tbl <- table(df[,1], df[,2], df[,3])
ftbl <- ftable(tbl)

ftable実施後

                01月 02月 03月 04月
                                   
営業部 営業次郎     0    0    1    0
       営業太郎     1    0    0    1
       技術冴子     0    0    0    0
       経理京子     0    0    0    0
       購買好男     0    0    0    0
       総務花子     0    0    0    0
技術部 営業次郎     0    0    0    0
       営業太郎     0    0    0    0
       技術冴子     0    0    0    1
       経理京子     0    0    0    0
       購買好男     0    0    0    0
       総務花子     0    0    0    0
経理部 営業次郎     ・・・以下略

営業部には営業太郎と営業次郎しかいないはずなので、

tbl <- tbl[rowSums(tbl) > 0,, drop=FALSE]

としましたが「次元が違う」とエラーになってしまいます。そもそも、こういうやり方ができるのかどうかをここ２〜３日調査しているのですが、ftableオブジェクト(?)に関する情報があまり見つけられなくてはまっております。

この「行合計がゼロ」のftableを削除する方法をご教授下さい。
実行環境はWindowsXP SP2 R-2.6.1 です。

よろしくお願い申し上げます。

↑

ケースの選択方法†

mikasa (2007-12-13 (木) 20:10:39)

SPSSからの転向に向けて、Ｒコマンダーを使って勉強中です。
(1) データセットから特定のデータだけを分析対象にする方法（SPSSで言うところのケースの選択）ですが、「データセットの部分集合を抽出」で新たなデータセットを作る以外に方法はないでしょうか？
(2) データセットを分割して一回の指示で分割したデータセットそれぞれの分析をすることは可能でしょうか？（SPSSで言うところの、ケースの分割）
よろしくお願いします。

↑

値を文字列として加えるには？†

AobaAoba (2007-12-11 (火) 17:03:50)

エクセルを卒業してRの世界に足を踏み入れたばかりのものです。
さきほどからこのサイトやらなにやらを3時間ばかり彷徨っているのですが、それでもわからず、たまらず質問させて頂きます。

やりたいことはもう少し複雑なのですが、以下のことができれば後は自分でできるはずなのです。


以下のように変数に数字を入れて、それぞれの値を代入したいのです。
x_1 <- 1
x_2 <- 2
x_3 <- 3
. . . x_10 <- 10


そこで次のようなことを書けば良いと思ったんですが、うまくいきません。

for (i in 1:10){
x_i  <- i
}


x_iのところのiの前後に何か入れるんじゃないかと思うのですが（エクセルのVBAではそうだったので）、どうすれば良いのでしょうか？
お手数おかけしてしまって本当にすみませんが、どなたかどうぞお願いいたします。

↑

Bray-Curtis Distanceの関数について†

R.,M. (2007-11-22 (木) 10:56:38)

初めて質問させていただきます。
昨日、クラスター解析の準備としてBray-Curtis Distanceの計算を試みたのですが以下のエラーメッセージが出て計算できません。

> dist.BC(x)
 以下にエラー dimnames(x) <- dimnames : 
   'dimnames' はリストでなくてはなりません

そこでdimnames(x)を実行すると”行の番号”と”変数名”が正常に出てきます。ちなみにdist.BC関数を含むpackage"clusterSim"中にある他の関数（例えばdist.SM）では問題なく計算できています。また、packageのヘルプを実行したところ以下の警告が出ています。

> help(package="clusterSim")
Warning message:
In readLines(f) :
   'C:/PROGRA~1/R/R-26~1.1RC/library/clusterSim/INDEX' で不完全な最終行が見つかりました

データは観測点における各生物分類群の現存量（対数変換済み）でデータの規模は173測点、52分類群のデータとなっています。他の関数（ユークリッド距離など）では問題なく計算できていますので、取り込んだデータの形式に問題は無いと思うのですが・・・

原因と修正方法をお教えていただけませんか。

↑

T,FとTRUE,FALSE†

初心者 (2007-11-16 (金) 20:50:15)

TRUE, FALSEは上書きは不可能なんですが，それらのエイリアスであるT，Fは上書き可能で，いろいろな値が代入できてしまいます(こないだこれでドツボにはまりました)．

プログラミングの際にはT，Fは使わずにTRUE，FALSEを使うのが好ましいということなんでしょうか

↑

R 言語でベクトル解析†

Ｅｕｌｅｒ (2007-11-16 (金) 20:31:13)

　てな、Ｗｅｂぺーじないかな？

↑

boxとプロットの間の隙間†

初心者 (2007-11-16 (金) 16:16:59)

いつもお世話になっています．

>theta <- seq(0,2*pi,0.01)
>plot(theta,sin(theta),xlim=c(0,2*pi));grid()

このようにプロットした際，boxの大きさはxlimで指定した範囲+アルファに設定され，プロットしたデータとboxの間に隙間ができてしまいます(つまり，boxの両端がx=0,x=2*piになるようにしたい)．この隙間を調整するにはどのようにすれば良いのでしょうか．

↑

インストーラなしのWin版Rのインストール†

内部統制 (2007-11-14 (水) 10:30:43)

　Win版Rをインストーラを使わないでインストール方法はあるでしょうか？

↑

２次輸送問題†

OR mania (2007-11-13 (火) 11:14:42)

　Ｒには２次輸送問題(QTP:Quadratic Transportation Problem)を解く関数はあるのでしょうか？

　こちらの「数理計画」のページを見ましたが見当たりませんでした。

↑

barplotのlegendについて†

K (2007-11-12 (月) 14:05:32)

barplotのlegendについて教えて下さい。
グラフをプレゼン用に使おうと、par(ps=22)として文字サイズを大きくしました。

par(ps=22)
par(mar=c(5,5,3,3))
barplot(sum, names.arg=breaks2,
        col = c("darkred", "darkgreen", "lightblue", "pink"),
        xlab="胸高直径階(cm)",  ylab="立木本数（本/ha）", main="")
legend("topright", c("アカマツ", "ヒノキ", "スギ", "その他" ),
       fill=c("darkred", "darkgreen", "lightblue", "pink"))

そして、以下のようなヒストグラムを得ました。

histgram_k2.png

この凡例中の文字列の前の四角のサイズを大きくしたいと思い、legendのマニュアルを見たのですが方法が分かりませんでした。どなたか方法をご教授いただければ幸いです。よろしくお願いします。

↑

Mac OS X 10.5にRのインストールができません†

suzuk (2007-11-11 (日) 11:37:40)

（Rがインストールされていない）Mac OS X 10.5 (Leopard)にRのインストールを試みましたが、以下のメッセージが表示されてしまいます。 「このソフトウェアのより新しいバージョンがシステムにインストールされています。この古いバージョンをインストールする必要はありません。」 R-2.6.0、R-2.5.0、R-2.4.1の何れのバージョンにおいても同様の現象が認められます。アドバイスをよろしくお願いします。

↑

キュー、スタックの実現†

kbi (2007-11-11 (日) 09:38:12)

google等の検索エンジンとサイト内検索を試したのですが
有効な情報が見つからなかったので質問させてください。
Rにおいて、基本的なデータ型としてキューやスタックは
用意されていないのでしょうか。
待ち行列などのモデルをシミュレーションしてみたいと思っています。
調べるのに有効なキーワードだけでも教えていただけると幸いです。

↑

パッケージのインストールができない†

tomo (2007-11-09 (金) 20:27:05)

Rの初心者です。　R2.4.1を使用しています。
パッケージのインストールができません。ダウンロードはしていますが，解凍ができていないようです。

> utils:::menuInstallPkgs()
URL 'http://prs.ism.ac.jp/bioc/1.9/bioc/bin/windows/contrib/2.4/Heatplus_1.4.0.zip' を試しています
Content type 'application/zip' length 436881 bytes
開かれた URL
downloaded 426Kb
~
以下にエラーzip.unpack(pkg, tmpDir) : 
ファイル 'C:/Program  Files/R/R-2.4.1/library/file57d3669a/Heatplus/chtml/Heatplus.chm' を開くことができません

違うパソコンからだと，同じ操作でインストールできました。パソコンの設定が悪いのでしょうか。 よろしくお願いいたします。

↑

ウィンドウサイズの拡大†

クラスタ (2007-11-07 (水) 20:30:56)

hclust()を使い樹形図を作成しました。
しかし、データ数が600個あり文字が重なって読めません。
ウィンドウサイズを拡大しようにも画面サイズより広がりませんでした。
なんとか読めるように文字が重ならない大きさまで拡大したいです。
ご教示お願いします。

  • 初級Ｑ＆Ａ アーカイブ(6)の『クラスター分析におけるフォントサイズ変更方法について』を見た上で書いてる？ -- 2007-11-07 (水) 21:51:46
  • 読んでませんでした。同じ問題を抱えた人の質問でしたが、フォントサイズではなくウィンドウサイズの拡張限界の変更方法は、ないでしょうか？ -- クラスタ 2007-11-07 (水) 23:16:38
  • 600個のデータを一度に見よう（一部分だけ見てもそれはそれで意味不明になる可能性大）というのが無茶なのでは。example(dendrogram) で見られるような文字グラフィックスで眺めてみたらどうでしょう。 -- 2007-11-08 (木) 00:12:34
  • こんな表示方法もあるんですね、ウィンドウの拡張は諦めて、この方法で試行してみます。ありがとうございます。 -- クラスタ 2007-11-08 (木) 00:33:27
  • これと同じ質問にここ数週間内に答えたんだけど，元記事が消去されたのか，なくなっていますね。
    それはさておき，画像を画面一杯に限って重ならないようにフォントを小さくしても，今度は，フォントが読めなくなるだけ。アラース。
    画面に描かれたものを見るだけでよいという局面は少ないわけで，何らかの形でプリントするのが最終形?とすると，なにも，グラフのサイズはコンピュータの画面のドット数には限定しなくても良いということ。
    いったん画像ファイルに保存すれば，それをそのまま印刷するには巨大すぎても，そのファイルを画像処理プログラム（フォトショップとかその他なんでも）で読み込めば，部分を拡大してみることもできるわけ。普通の大きさのフォントを使っても，重ならないように大きな描画領域をとればよいだけ。
    つまり，横幅（縦幅も）を十分大きく取って描画し，それをファイルに落とすようにしておけばよいだけ。
    以下を実行して，できるファイルをご覧下さいませ。 -- 2007-11-08 (木) 00:54:04 
    x <- matrix(rnorm(1200),600)
pdf("cluster.pdf", height=10, width=150)
plot(hclust(dist(x)), hang=-1)
dev.off()
  • おお、できました。これを求めていたんです。お二方ありがとうございました。 -- クラスタ 2007-11-08 (木) 01:40:35
↑

関数 cut() の引数 right について†

実験万 (2007-11-01 (木) 14:40:22)

ある数値データベクトルの度数分布を求めるときに、階級の端点に一致する数値がある場合、その数値を、下側と上側、どちらの階級に含めるべきか、という問題を検討しています。そこで、次のような実験をしてみました。

まず、0.1 から 9.9 まで、0.1 刻み、99 個の数値からなるデータベクトルを用意します。このデータベクトルについて、最初の階級の下限 0、最後の階級の上限 10、階級幅 1 の条件で度数分布を求めます。

実験 (1)
階級の端点に一致する数値がある場合、その数値を「下側」の階級に含めることにすれば、次のような結果を得ます。

命令:
x <- seq( 0.1, 9.9, 0.1 )
table( cut( x, breaks=seq(0,10,1), right=TRUE ) )
結果:
(0,1]  (1,2]  (2,3]  (3,4]  (4,5]  (5,6]  (6,7]  (7,8]  (8,9] (9,10]
   10     10      9     11     10     10     10     10     10      9


実験 (2)
階級の端点に一致する数値がある場合、その数値を「上側」の階級に含めることにすれば、次のような結果を得ます。

命令:
x <- seq( 0.1, 9.9, 0.1 )
table( cut( x, breaks=seq(0,10,1), right=FALSE ) )
結果:
[0,1)  [1,2)  [2,3)  [3,4)  [4,5)  [5,6)  [6,7)  [7,8)  [8,9) [9,10)
    9     10     10     10     10     10     10     10     10     10


それぞれ、最初と最後の階級の度数をみると、期待したとおりの結果になっています。

しかし、ひとつ問題が見つかりました。それは、実験 (1) の結果で、(2,3] の度数が 9、(3,4] の度数が 11 であることです。本当は、それぞれ 10 であるはずだと思うのです。このような結果になる原因について、ご教示いただけませんでしょうか。

↑

0.999999999†

ビットマップ (2007-11-01 (木) 12:18:08)

R version 2.5.1です。

set.seed(100)
vo <- round( rnorm(40,100,5), 1)


table(diff(sort(vo)))

とすると

                0 0.0999999999999943  0.100000000000009  0.200000000000003 
                4                  8                  5                  1 
0.299999999999997  0.300000000000011  0.399999999999991  0.400000000000006 
                5                  1                  1                  1 
              0.5  0.599999999999994  0.700000000000003                  1 
                7                  1                  1                  1 
 1.20000000000000              2.800                3.5 
                1                  1                  1


という結果が得られました。
0.0999999999999943 と 0.100000000000009　はどちらも 0.1 と扱ってほしいのですが、そうするにはどうしたらよいですか？

↑

グラフの軸を太くしたい†

surg (2007-10-28 (日) 07:40:31)

Rで作成したグラフを、Microsoft Wordの書類に挿入しようと思っております．WindowsとMacintoshの双方で利用したいため，グラフを大きめのbitmapで作成し，Word上で縮小表示しようとしております．

png("test.png", width=324*4, height=324*4)
par(cex=4, lwd=4)
hist(rnorm(1000))
dev.off()

上記のようにしてWord上で15%程度に縮小して印刷すると，グラフの軸だけがかすれてしまいます．軸を太くしたいのですが，何か方法はありますでしょうか？ help(par) は読んだのですが，それらしい項目は見当たりません．

↑

PLSとCCAは同等か†

yishii (2007-10-24 (水) 23:04:38)

初心者なので教えていただきたいことがあり、投稿させていただきます。あるケモメトリックスの入門書に、PLS(Partial least squres)法と正準相関解析(Canonical correlation alalysis,CCA)が同等だと書いてあったのですが、それは正しいのでしょうか。もしそうだとすれば、この両者の関係を論じた文献があれば、ご教示いただければ幸いです。よろしくお願いします。

↑

R commanderでの生存時間解析(cmprsk)における多変量解析(crr)について教えてください†

takeiteasy (2007-10-24 (水) 00:22:52)

先日、Rcommanderでのcmprskを教えていただいたものです。おかげさまで、cumincに関しては、Rcommanderで解析可能になりました。ありがとうございます。たびたび申し訳ありませんが、cmprskにおいての多変量解析であるcrr(competing risks regression analysis)をRcommanderに機能実装しようと試みましたが、うまくいきません。
Rcommanderハンドブックの6章のCox回帰を参考にやってみましたが、根本的な理解が不足していて、できません。
申し訳ありません。

↑

e1071パッケージのsvmメソッドの損失関数について†

codedir (2007-10-19 (金) 12:51:21)

e1071パッケージのsvmメソッドの回帰における損失関数はε-インセンシティブ損失関数で線形だと思うのですが、この損失関数を線形だけでなく２次関数に変更したりできるように、プログラムを触りたいんです
損失関数を定義しているプログラムを探したんですが、見当がつきません。プログラム内のepsilonが出てくる部分を見てみても、epsilonを使って計算しているようには見えません。
よろしくお願いします

e1071:::svm.default
e1071:::svm.formula
↑

Rcommanderでcmprskはできますか？†

takeiteasy (2007-10-17 (水) 23:46:49)

医学データの解析に使用したいのですが、
いままで、JMPなどで統計解析していましたが、
cmprskはJMPなどでは解析できません。
コマンド入力などは素人でなかなか難しいです。
Rcommanderなどを用いてcmprskは解析できるでしょうか？

 ### MyProgram.R の中身
 Mycuminc <- function(){
  initializeDialog(title="cumulative incidence function")
  timeBox   <- variableListBox(top, Numeric(), title="時間変数")
  eventBox  <- variableListBox(top, Numeric(), title="イベント/打ち切り変数")
  groupBox  <- variableListBox(top, Factors(), title="群変数")
  Var1      <- tclVar("0")
  Var1Entry <- tkentry(top, width="6", textvariable=Var1)
  onOK <- function(){
    XXX <- getSelection(timeBox)
    YYY <- getSelection(eventBox)
    ZZZ <- as.factor(getSelection(groupBox))
    AAA <- as.numeric(tclvalue(Var1))
    if (length(XXX) != 1 || length(YYY) != 1 || length(ZZZ) != 1) {
      errorCondition(recall=Mycuminc, message="変数を1つずつ選択してください。")
      return()
    }
    closeDialog()
    library(cmprsk)
    if (is.na(AAA)) {
      command <- paste("cuminc(", 
        ActiveDataSet(), "$", XXX, ", ",
        ActiveDataSet(), "$", YYY, ", ",
        ActiveDataSet(), "$", ZZZ, ", ",
        "cencode=0)", sep="")
    }
    else {
      command <- paste("cuminc(", 
        ActiveDataSet(), "$", XXX, ", ",
        ActiveDataSet(), "$", YYY, ", ",
        ActiveDataSet(), "$", ZZZ, ", ",
        "cencode=",           AAA, ")", sep="")
    }
    doItAndPrint(command)
    tkfocus(CommanderWindow())
  }
  OKCancelHelp(helpSubject="cuminc")
  tkgrid(getFrame(timeBox), getFrame(eventBox), getFrame(groupBox), sticky="w")
  tkgrid(tklabel(top, text=gettextRcmdr("打ち切り")), Var1Entry, sticky="e")
  tkgrid(buttonsFrame, columnspan=2, sticky="w")
  tkgrid.configure(Var1Entry, sticky="w")
  dialogSuffix(rows=2, columns=1)
}

Mycifplot <- function(){
  initializeDialog(title="cumulative incidence function")
  timeBox   <- variableListBox(top, Numeric(), title="時間変数")
  eventBox  <- variableListBox(top, Numeric(), title="イベント/打ち切り変数")
  groupBox  <- variableListBox(top, Factors(), title="群変数")
  Var1      <- tclVar("0")
  Var1Entry <- tkentry(top, width="6", textvariable=Var1)
  onOK <- function(){
    XXX <- getSelection(timeBox)
    YYY <- getSelection(eventBox)
    ZZZ <- as.factor(getSelection(groupBox))
    AAA <- as.numeric(tclvalue(Var1))
    if (length(XXX) != 1 || length(YYY) != 1 || length(ZZZ) != 1) {
      errorCondition(recall=Mycuminc, message="変数を1つずつ選択してください。")
      return()
    }
    closeDialog()
    library(cmprsk)
    if (is.na(AAA)) {
      command <- paste("tmp <- cuminc(", 
        ActiveDataSet(), "$", XXX, ", ",
        ActiveDataSet(), "$", YYY, ", ",
        ActiveDataSet(), "$", ZZZ, ", ",
        "cencode=0); plot(tmp)", sep="")
    }
    else {
      command <- paste("tmp <- cuminc(", 
        ActiveDataSet(), "$", XXX, ", ",
        ActiveDataSet(), "$", YYY, ", ",
        ActiveDataSet(), "$", ZZZ, ", ",
        "cencode=",           AAA, "); plot(tmp)", sep="")
    }
    doItAndPrint(command)
    tkfocus(CommanderWindow())
  }
  OKCancelHelp(helpSubject="cuminc")
  tkgrid(getFrame(timeBox), getFrame(eventBox), getFrame(groupBox), sticky="w")
  tkgrid(tklabel(top, text=gettextRcmdr("打ち切り")), Var1Entry, sticky="e")
  tkgrid(buttonsFrame, columnspan=2, sticky="w")
  tkgrid.configure(Var1Entry, sticky="w")
  dialogSuffix(rows=2, columns=1)
}

MySetup <- function(){
  initializeDialog(title="cmprsk のセットアップ")
  onOK <- function(){
    closeDialog()
    install.packages('cmprsk')
  }
  OKCancelHelp(helpSubject="install.packages")
  tkgrid(buttonsFrame, columnspan=1, sticky="w")
  dialogSuffix(rows=1, columns=1)
}

次に、同フォルダにある Rcmdr-menus.txt の一番最後に以下を追記して下さい．

### Rcmdr-menus.txt の一番最後に追記
### メニュー追加・ここから
   menu  MyMenu   topMenu  ""                       ""         ""  ""
   item  topMenu  cascade  "メニュー"               MyMenu     ""  ""
   item  MyMenu   command  "CIF（計算）"            Mycuminc   ""  ""
   item  MyMenu   command  "CIF（プロット）"        Mycifplot  ""  ""
   item  MyMenu   command  "cmprsk のセットアップ"  MySetup    ""  ""
### メニュー追加・ここまで

まずは R commander を起動し，「メニュー」→「cmprsk のセットアップ」を行った後，
サンプルデータを R Commander でアクティブにしてから，「メニュー」→「CIF（計算）」などを実行してください．

### サンプルデータ
x <- data.frame(time=c(1:10,1:10), censor=floor(3*runif(20)), group=c(rep("A",10),rep("B",10)))

R Commander への機能追加については 説明スライドの最後の方 や R Commander本の6章を参照ください．ちなみに，コマンドでやられる場合は・・・

x <- data.frame(time=c(1:10,1:10), censor=floor(3*runif(20)), group=c(rep("A",10),rep("B",10)))
library(cmprsk)
( tmp <- cuminc(x$time, x$censor, x$group) )
plot(tmp)
cif.gif

です．-- 舟尾 2007-10-18 (木) 02:02:03

↑

音声データをスペクトル解析したい†

trantila (2007-10-16 (火) 18:28:30)

音声データ（WAV形式です。）をスペクトル解析したいのですが、どのようにしたら良いのでしょうか？
ここのサイトのスペクトル解析のところは見たのですがいまひとつよく分かりませんでした。
Rでは音声データを解析対象として扱えないのでしょうか？
もし、そうであれば音声データを時系列データに変換する方法が必要なのでしょうか？

よろしくお願いします。

↑

RESTFul な Web Service　で Rを利用†

tantan (2007-10-15 (月) 11:48:37)

RWeb なような対話型の web service ではなくて、CGI の Get や Put などで、入力パラメータを指定して、サーバサイトで R を実行し、戻り値を XML 形式にして出力するような Web アプリってないでしょうか？

↑

RからWSDLを読み込みたい†

akira (2007-10-14 (日) 17:17:22)

なかまさんのRからのBBDJアクセスメモを試したところ、エラーがでてしまいました。proxyは変更していません（<-これか？）。実のところ、SSOAPがよくわかりません。

> library("SSOAP")
 要求されたパッケージ XML をロード中です 
 要求されたパッケージ RCurl をロード中です 
> ddbj <- processWSDL ("http://xml.nig.ac.jp/wsdl/GetEntry.wsdl")
> iface <- genSOAPClientInterface(def=ddbj)
 以下にエラー parse(text = txt) :   予想外の 入力  です  以下中: 
 "       server = .defaultServer, .convert = .operation@returnValue, .opts = list(), ..., 
         nameSpaces =  &#65533;"

とエラーになってしまいました。
関数を読み込むと、

> processWSDL
function (fileName = "KEGG.wsdl", handlers = WSDLParseHandlers(fileName), nameSpaces = character()) {
    library(XML)
（中略）
    types = processWSDLTypes(doc[["types"]], doc)
（以下、省略）

と、"processWSDLTypes"を見つけられず、トレースできませんでした。 で、exampleを確認すると、

>  tmp = processWSDL(system.file("examples", "KEGG.wsdl", package = "SSOAP"))
Warning messages:
1: In function (node)  :
  skipping import node with no schemaLocation attribute
2: In function (node)  :
  skipping import node with no schemaLocation attribute

となり、KEGG.wsdlの読み込みでエラーになったようで、

> tmp = processWSDL("http://soap.genome.jp/KEGG.wsdl")
Warning messages:
1: In function (node)  :
  skipping import node with no schemaLocation attribute
2: In function (node)  :
  skipping import node with no schemaLocation attribute

と、大元にアクセスしてもだめでした。
一つのスレで2つの疑問が生じていますが、

ddbj <- processWSDL ("http://xml.nig.ac.jp/wsdl/GetEntry.wsdl")
tmp = processWSDL("http://soap.genome.jp/KEGG.wsdl")

のエラーを回避したいです。

> sessionInfo()
R version 2.6.0 (2007-10-03) 
i486-pc-linux-gnu 
locale: LC_CTYPE=ja_JP.UTF-8;LC_NUMERIC=C;LC_TIME=ja_JP.UTF-8;
LC_COLLATE=ja_JP.UTF-8;LC_MONETARY=ja_JP.UTF-8;LC_MESSAGES=ja_JP.UTF-8;
LC_PAPER=ja_JP.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;
LC_MEASUREMENT=ja_JP.UTF-8;LC_IDENTIFICATION=C
attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     
other attached packages:
[1] SSOAP_0.4-6 RCurl_0.8-3 XML_1.93-2 
loaded via a namespace (and not attached):
[1] rcompgen_0.1-15

です。

  • SSOAPの0.4-3をお使いください. 直しかけましたが, 挫折しました... -- なかま 2007-10-14 (日) 20:52:56
  • ありがとうございます。DDBJにアクセスできました！
    ただ、SSOAPの0.4-3でも、example(processWSDL)はこけますね。これはexampleにあるkegg.wsdlに"schemaLocation attribute"がないっていうエラーなんでしょうか？ SAXとか、DOMとか、よくわからないんです...スミマセン
    RSiteSearch("processWSDL")すると、"http://www.ncbi.nlm.nih.gov/entrez/eutils/soap/eutils.wsdl"でもこけてるみたいですね。 
    > ncbi <- processWSDL("http://www.ncbi.nlm.nih.gov/entrez/eutils/soap/eutils.wsdl") 
 以下にエラー ans[[ns]] <<- o :  一つ未満の要素を選択しようとしました 
 追加情報:   50 件以上の警告がありました (警告を見るには warnings() を使って下さい) 
> warnings()
 警告メッセージ： 
1: In FUN(X[[2L]], ...) :
  Failed to handle node HasLinkOut of type attribute&simpleType in processWSDLType
2: In FUN(X[[2L]], ...) :
  Failed to handle node HasNeighbor of type attribute&simpleType in processWSDLType
3: In FUN(X[[1L]], ...) :
(以下、略)
    うーん、DDBJだけじゃなく、KEGG APIもeUtilsも使いたかったのですが... -- akira 2007-10-14 (日) 21:18:15
  • 最新のSSOAP_0.4-6でも,options(useFancyQuotes=FALSE)で`ddbj'なら動きます(dQuoteがエンコードにより異なる仕様になったのが原因でした). ちょっと, これに取り組む余力はありませんです. -- なかま 2007-10-15 (月) 11:18:18
  • ありがとうございます。今はWinXPで使っています。SSOAP_0.4-4ならoption指定無しでddbjは動くようです。 -- akira 2007-10-15 (月) 12:30:17
↑

グラフのフォントサイズを設定したい†

田中 (2007-10-14 (日) 15:41:02)

R2.6.0，Windows XPを使っております。
barplot()を使って棒グラフを描き，PowerPointのスライドに貼り付けてプレゼンをしようと考えています。
PowerPointに貼り付けてみると，文字が少し小さく，フロアからは見えにくいようです。
x軸，y軸のラベルや図のタイトルのフォントサイズを変更することはできるでしょうか。
どなたか対処法をご存知の方がいらっしゃいましたら，お返事をどうぞよろしくお願い致します。

↑

線分の描き方について†

真矢 (2007-10-11 (木) 18:27:55)

座標（X,Y）を始点とした傾きｍ、長さLの線分の書き方を教えてください。
R2.4.1を使っています。

↑

ESSからRcmdrを日本語で呼び出したい†

＜ふ＞ (2007-10-02 (火) 21:57:19)

（連続で失礼します。）

essをMeadowで使ってます。
（WindowsXP、R2.5.1、Rcmdr1.3-5、ess-5.3.4、Meadow3）
このessのR processからlibrary(Rcmdr)を実行すると、langが渡らないのか、起動したRcmdrのメニューが日本語になりません。

そこで、Rのインストールディレクトリの /etc/Rconsole を開いて、

## Language for messages
language = ja
のように、強制的にjaにしてみました。

これで、起動するRcmdrのメニューなどは、めでたく日本語になったのですが、「データセットの表示」「データセットの編集」では、日本語が表示されません。
やはり、Tcl/Tkにちゃんとlangを伝えないとだめなようです。

Meadowの中のなにかを設定すればいいのではないかと思おうのですが、よくわかりません。

どなたかアドバイスをお願いいたします。

↑

datasetにhelpをつける方法を知りたいのですが†

＜ふ＞ (2007-10-02 (火) 14:58:10)

RcmdrでR2.5.1を使っています。

データの集計は、Excelでおこなって、それをRのDatasetとして読み込ませたあと、数値を因子に変換したり、あれこれ修正をして、.Rdaで保存し、それをもとに分析を行っています。

この.Rdaでsaveされたデータセットに、helpをつける方法がわかりません。

Packageで提供されているDatasetには、helpがついていて、便利なものですあから（というよりも不可欠..）、一度整形したデータを使いまわすためにデータの定義など、きちんとしたドキュメントをデータにつけておきたいと思っております。

ご教示ください。

↑

折れ線プロットは？(plot)†

kd (2007-09-25 (火) 03:46:53)

折れ線プロット（すなわちポイント点は無印）を描こうと思っております．
（Mac版 R-2.5.1.dmg を使っております．）
しかし下記を試しても折れ線プロットを描いてくれませんでした．何か方法はあるでしょうか？

plot(x$V2~x$V1,pch=-1,lty=1) ---> 丸印でポイント点をプロットするだけ

plot(x$V2~x$V1,pch="",lty=1) ---> 無印，無線

↑

igraphについて†

あひる (2007-09-20 (木) 23:42:45)

現在、大学でigraphを用いた授業を受けております。家でもじっくり課題ができるようにと早速Ｒをダウンロードしたのですが、どうしたらigraphが使えるようになるのかがわかりません。Googleで検索したページでigraphのダウンロードはできたようなのですが、Ｒでlibrary(igraph)と打ってもエラーになってしまいます。どうしたらigraphが使えるようになるのでしょうか？どなたかご教授頂けたら幸いです。

↑

R Consoleの表示容量について†

カイロ (2007-09-20 (木) 22:25:29)

Rを利用して3日目です．医療データを用いたROC曲線を用いて予後推定のcut-off値を求める作業をしていますが，データの取得で躓いています．RjpWiki内，Google，等を検索しても分かりませんでしたので，こちらでお伺い申し上げます．大変初歩的な質問と思いますが，ご教示いただけましたら幸いです．ROC曲線作成には青木先生のライブラリからsourceを使用させて頂いております"http://aoki2.si.gunma-u.ac.jp/R/ROC.html"．

図のように，R Console上で数千行ある計算結果が得られるのですが，スクロールバーを一番上まで持っていっても全ての結果を表示できません．GUIプリファレンスでbuffer bytes等を増やしてもうまくゆきません．R Console上で表示，あるいは別ファイルに書き出し，などして全ての計算結果を取得する方法はどうすれば宜しいでしょうか．

宜しくお願いいたします．

当方の環境：
R version 2.5.1 (2007-06-27)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;
LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
attached base packages:
[1] "stats" "graphics" "grDevices" "utils" "datasets"
[6] "methods" "base"

consolegamen.png


↑

nlme ライブラリ，lme のデフォルト†

キュロシ (2007-09-20 (木) 08:32:11)

Mixed Effects Models in S and S-PLUS を勉強中のキュロシです．
作者（Pinheiro & Bates）に訊きなさい，という種の質問かもしれません．

> library(nlme)

として nlme パッケージを使います．いくつかのデータセット (groupedData) がありますが，ここでは Orange を使って説明します．

> f1 <- lme(Orange)

とすると，以下の結果が得られます．

> f1
Linear mixed-effects model fit by REML
Data: Orange 
Log-restricted-likelihood: -139.9061
Fixed: circumference ~ age 
(Intercept)         age 
17.3996502   0.1067703 
Random effects:
Formula: ~age | Tree
Structure: General positive-definite
           StdDev      Corr  
(Intercept)  1.89732541 (Intr)
age          0.02402639 -1    
Residual    10.00622121       
Number of Observations: 35
Number of Groups: 5

つまり，lme の argument である，fixed, random　を data を指定することによって，はしょっています． とすると，

> f2<-lme(fixed=circumference ~ age, data=Orange, random=~age|Tree)

とすることで同じ結果が得られそうなものですが，得られません． f1 のモデルが 収束しない f2 と何が違うのかがわかりません．何かヒントになりそうなことがありましたら，是非とも教えてください．よろしくお願いいたします．

↑

回帰分析についての質問†

netbean (2007-09-18 (火) 16:46:17)

Rの初心者です。他ツールで解析経験がありますが、Ｒはまだ数時間の実践しかやっていません。
step関数を用いてステップワイズ変数選択を実施する際、下記エラーが出ました、

> step(RegModel.2)
Start:  AIC= 189.44
。。。
object$residuals + predict(object) 中で警告がありました：
         長いオブジェクトの長さが
         短いオブジェクトの長さの倍数になっていません
以下にエラーlm.fit(x[, jj, drop = FALSE], y, offset = offset) : 
        矛盾した次元です

回帰のモデル精度が非常に低く、stepwiseを期待したいところが、ここで固まって進めない状況です。どなたか対処法をご存知でしたら、教えてください。ありがとうございます。

↑

判別分析の変数選択†

mayu (2007-09-15 (土) 22:49:29)

discの線形判別関数のソースを使い、判別分析をしています。
自分が選んだ変数だけで解析したい場合どうすればいいのでしょうか？
変数選択が出来るソースを探していたのですが見当たらず質問させて頂きました。
自分でプログラムを組めばいいのかもしれませんが、プログラムが分からず、
勉強不足ですみませんが、よろしければ回答お願いします。

↑

数値データ、非数値データが混在するデータのカウント方法について†

terada (2007-09-14 (金) 21:20:20)

こんにちは。数値データ、非数値データが混在するデータのカウント方法について質問させてください。

> dat
    　　記号1 　 　記号2 　　  記号3 　     記号4
90         B1         A1          A2           A2
91         B1         A1          A1           A1
92         B1         A1          A1             
93         B1         A1          A1           A1
94         B1         B1          B3            K
95                                              K                               
96         B0         B0          B1           A1
97         B0                              
98         B1         A1                    
99          3                              
100                    3                              
101         2                              
102         2                              
103         1

以上のようなデータ(14行4列のデータで、空いているところはnullです。)がある場合に、このデータフレームのすべての要素、B0、B1、B3、A1、.........　、1、2、3がそれぞれ何個ずつあるかカウントし、かつこのデータフレームにはない要素、B2、A3は0であるとカウントさせ、その数値データをベクトルの形で得たいのですが、何か方法があれば教えてください。table関数を使って色々やってみたのですがデータフレーム単位の結果がどうしても出せません。
よろしくお願いします。

すみません。質問の仕方が不適切でした。訂正します。 以下のような結果が得たいという意図で質問しました。

A1  A2  A3  B0  B1  B2  B3   K  1  2  3
11   2   0   3   7   0   1   2  1  2  2

再度アドバイスをいただければ幸いです。よろしくお願いします。

↑

ベクトルから特定の要素を削除するには？†

青葉 (2007-09-13 (木) 13:40:42)

昨日から始めた初心者です。壁にごつんごつん当たりながら進んでます。

xから、yの中にある要素を数も考慮して取り除きたいのですが、どうすればよいでしょうか？
x <- c(1,1,1,2,2,2,3,3,3)
y <- c(1,2,2)
この2つから、(1,1,2,3,3,3)というベクトルを作りたいのです。
また、xから要素の種類だけ特定して取り除くことが可能でしょうか？
例えばx <- c(1,1,1,2,2,2,3,3,3)から2という数字だけを指定して、すべての2を取り除く方法です。
検索の仕方が悪いのかもしれませんが、苦戦しております。どなたかどうぞよろしくお願いします。

↑

ベクトルの要素の種類数を知るには？†

青葉 (2007-09-13 (木) 00:13:47)

今日使い始めたばかりの初心者です。
ベクトルの要素の数を知るにはlengthで良いことは分かったのですが、種類数を知る方法が見つからず1時間ほど苦戦しています。
例えば、x <- c(1, 1, 2, 2, 2, 3)だったら、要素の数は6ですが、要素の種類数は、1と2と3の3つです。この3という値を知りたいのです。
どなたかどうぞよろしくお願いします。

↑

tcltk のロードができません†

yosito (2007-09-11 (火) 18:21:39)

R 初心者です。環境は MacOS X (10.4) にて、R.app を使っています。
tcltk は標準で package に含まれているのですが、いざロードしようとしますと、
以下のようなメッセージが出てロードできませんでした。

要求されたパッケージ tcltk をロード中です
Tcl/Tk インターフェース ロード中... Error in fun(...) : couldn't connect to display ":0"
Error :  .onLoad は 'tcltk' のための 'loadNamespace' で失敗しました
エラー： パッケージ 'tcltk' をロードできませんでした

どなたか対処法をご存知でしたら、教えてください

↑

nlme ライブラリ内の plot のマージンを変更する方法†

キュロシ (2007-09-11 (火) 04:42:46)

Pinheiro & Bates の "Mixed-Effect Models in S and S-PLUS" という本で勉強しています．

library(nlme)

とした後

plot(Dialyzer,outer=~QB)

を実行すると，下図のようなきれいな図が描けます．

groupedData1.jpg


ひとつだけ大きな外れ値があったとします．すなわち

Dialyzer$rate[1] <- 500

です．

ここで，先ほどと同じ図を描こうとすると

plot(Dialyzer,outer=~QB)

マージンが妙に大きい奇妙なプロットになってしまいます．

groupedData2.jpg


このマージンを自由に変更したいのですが，いまひとつやり方がわかりませんでした．
ヘルプファイルは

?plot.nffGroupedData(nlme)

だと思うのですが，該当するヘルプがわかりませんでした．
このマージンはどのように変更できるのでしょうか？

↑

特定の行と列を削除した行列の作り方†

ＫＫ (2007-09-10 (月) 11:57:13)

ある行列Ａの特定の行および列を削除した行列Ｂを作りたい場合、どのようにすればよいでしょうか。
よろしくお願いいたします。

↑

日付＋時刻の時系列のプロットでの扱い†

たけ (2007-09-01 (土) 15:54:42)

たとえば2006/01/01 00:00:00から2007/09/01 15:48:32までの時系列で欠損値-999、行の欠落はないデータを扱う際、解析自体は適当に行番号でできるのですが、plotする際に、x軸に日付をどころどころ表示したプロットを出力したいです。ts()でも年+月までであればそれっぽいものができますが、この日付+時刻の状態を扱ってプロットにもっていけるパッケージやtsオブジェクトの作成方法をご存知の方がいましたらお教いただけますでしょうか。環境は2.51 w32です。
どうぞよろしくお願いします。

↑

2つの因子の重複†

KEN (2007-08-30 (木) 17:53:43)

2つのデータ（data1とdata2）の各因子（1〜6）の一致をみたいと思っています。一致したデータ数（今回の場合は３行目の１つ）とその割合（今回の場合は1/5）を計算するにはどうすればよろしいでしょうか。Excelに限界を感じ先週よりRの勉強を始めました。超初心者ですが、よろしくお願いいたします。

	data1	data2
1	A	B
2	B	A
3	B	B
4	B	A
5	A	B

使用環境：R version 2.4.1 (2006-12-18) i386-pc-mingw32

↑

heatmapのlegend†

みゅー (2007-08-29 (水) 11:02:21)

かなり調べたつもりなのですが、heatmapのlegendのつけ方がわかりません。heatmapのカラー表示がどのレベルに対応するのかを明示したいので、方法が書いてあるページなどをお教えいただけないでしょうか。

↑

Rのコンパイラー†

あつ (2007-08-28 (火) 18:31:23)

http://www.cs.uiowa.edu/~luke/R/compiler/
RのコンパイラーをWindowsのマシンに入れようとするのですがうまく行きません。
tar.gz形式なので、一度展開してzip形式にしてGUI画面から読ませると、昔のRしか動かないよ。。と言った感じのメッセージが出ます。
仕方なく

R CMD INSTALL  compiler_0.1.-3.tar.gz

としても、テンポラリーの位置が合わないのか展開されせん。tar,gz,perlは、入っています。
どうやったら、これがR-2.5.1で動くか、どなたかお教え願えませんでしょうか？

↑

変数の交換†

とも (2007-08-27 (月) 17:29:39)

中間変数を使わずに内容を交換することは可能ですか。例としてa1とa2の内容を交換する方法をお示しします。

a1 <- c(1,2,3)
a2 <- c(2,3,4)
z <- a1
a1 <- a2
a2 <- z
↑

関数内の変数が展開されない†

kutakuta (2007-08-24 (金) 20:46:58)

よろしくお願いします。 R version 2.5.1 (2007-06-27) powerpc-apple-darwin8.9.1 を使用しています。 ggplotというライブラリをつかって箱ヒゲ図+jitterplotを作ろうとしています。 行列"data"の中に"class"というカラムと"colx"(x=1-N)というカラムがあって、classごとの箱ヒゲ図を作るのが最終目的です。
例： http://bg9.imslab.co.jp/Rhelp/R-2.4.0/src/library/ggplot/man/images/ggjitter_004.png

p <- ggplot(data,aesthetic=list(y = class, x = col1))
ggjitter(ggplot(p))

とするとうまく目的の図が得られるのですが、

i <- 1
p <- ggplot(data,aesthetic=list(x = class, y = paste("col",i,sep=""))
ggjitter(ggplot(p))

とするとうまくいきません。 summary(p)で変数の中身をみてみると、xの値が望みの"col1"ではなく"paste("col",i,sep="")"というように変数展開されていない事が分かりました。 いろいろ試してみたのですが、どうしてもこの問題が解決できません。 どなたかヒントをいただけないでしょうか。

↑

Rコマンダーの日本語表示について†

J (2007-08-23 (木) 20:55:38)

Mac OS 10.4.10 上で R ver.2.5.1 とRcmdr 1.3-5 を使っています。
Rコマンダーの日本語メニュー表示がやや読みにくいのですが，これはfontの設定等で改善できるものでしょうか？

↑

ファイル出力加工†

Bun (2007-08-23 (木) 20:21:21)

Ｒ初心者です。分かりづらい標記でしたら、ご指摘の程よろしくお願いいたします。

以下のプログラムで、クロス集計の結果を外部出力しました。

x<-c(1,2,1,2,1,2)
y<-c(1,1,1,2,2,2)
z<-data.frame(x,y)
z[,1] <- factor(z[,1], levels=1:2, labels=c("男", "女"))
XTBL<-table(z[,1],z[,2])
write.table(XTBL, file="C:/XTBL.txt", append=F, quote=F, col.names=T)

下記は、出力結果

   1 2
男 2 1
女 1 2

上記の結果ファイルをエクセルで読み込み、作業したいのですが、Ｒからの出力のデータファイル時に、男女の上に変数名を付けて、出力することは可能でしょうか？

↑

Rコマンドファイル内からディレクトリ作成は出来ますか？†

みゅー (2007-08-23 (木) 16:21:27)

R初心者につき言葉足らずな質問になっていたらごめんなさい。ご指摘ください。

大量のファイルが生成されるため種類ごとに出力を分けたいと思いました。
シェルコマンドのmkdirに対応するようなRコマンドはありますか？
みつけられませんでした。
(ls()やrmは使えるようですね。)

単に、./以下にhogeディレクトリを作成していない状態のまま"./hoge/filename.pdf"とやってももちろん./hogeは自動生成しませんでした。
先に./hogeを用意しておかないとだめ、
ディレクトリ生成はできないということでしょうか。

↑

変数名の修正†

Bun (2007-08-23 (木) 15:22:23)

既に変数名を含んだデータセットについて、一部の変数名を新しい変数名にしたい場合、どのような関数、プログラムを利用すればよいのでしょうか。

↑

関数実行時の入力文字列の取得方法†

r_user (2007-08-18 (土) 10:14:36)

たとえば、f(x)と入力して関数fを実行すると、関数にはオブジェクトxの中身が渡されますが、
中身だけでなく、文字列xも取得して利用できるような関数を定義したいと考えています
引用符を付けてf("x")のようにするのではなく、普通にf(x)と入力して、
xの中身に加えて、"x"という文字列も取得して利用できないものでしょうか。

↑

条件によるIDの振り分けについて†

suzuki (2007-08-17 (金) 17:17:46)

はじめまして。こんにちは。
データによるIDの振り分けについて質問させてください。

DATA1	 DATA2	 DATA3
1	 30	 1
1	 30	 1
1	 30	 1
1	 30	 2
1	 30	 2
1	 31	 1
1	 31	 1
1	 31	 1
1	 31	 1
2	 30	 1
2	 30	 2
2	 30	 2
2	 30	 2
2	 30	 3
2	 30	 3


以上のようなデータがあって、以下のようにDATA1、DATA2、DATA3をキーとして全ての値が同じであれば、同じIDとするように
IDを振り分けていきたいのですが、各IDの行数がランダムであるためうまく処理できる方法が見つかりません。
ある程度、手作業でできなくもないのですが、20万行ほどのデータであるため、何とか自動化できないものかと考えています。
アドバイスよろしくお願いします。

DATA1	 DATA2	 DATA3	 ID
1	 30	 1	 1
1	 30	 1	 1
1	 30	 1	 1
1	 30	 2	 2
1	 30	 2	 2
1	 31	 1	 3
1	 31	 1	 3
1	 31	 1	 3
1	 31	 1	 3
2	 30	 1	 4
2	 30	 2	 5
2	 30	 2	 5
2	 30	 2	 5
2	 30	 3	 6
2	 30	 3	 6
↑

Rcmdr QQプロット†

KT (2007-08-15 (水) 18:25:40)

Rcmdrで2つのデータをQQプロットに同時にのせたいと思います。可能でしょうか？

↑

RのグラフをPower point上で変更†

sh (2007-08-08 (水) 01:40:14)

Rのグラフ（軸ラベルのフォント、プロット記号）をPower point上で変更することは可能でしょうか？

↑

同じ要素を持つ行を選択するには†

ななし (2007-08-08 (水) 01:14:02)

次のようなデータがあったとき

> x <- matrix(c("a","b","c","c","e","d","a","c","c","c","c","c"), c(4,3), byrow=T)
> x <- data.frame(x)
> x
  X1 X2 X3
1  a  b  c
2  c  e  d
3  a  c  c
4  c  c  c

X1が等しくかつX3が等しい行を選択したい（この場合は1行と3行）のですがどのようにすればよいのでしょうか？

>  which( x[,1]==x[,1] & x[,3]==x[,3])
[1] 1 2 3 4

や

>  x[,4] <- 1:nrow(x) #自分と同じ行を除けばいけるかも…
>  which( x[,4]!=x[,4] & x[,1]==x[,1] & x[,3]==x[,3])
integer(0)

や

> which(duplicated(x[,1]) & duplicated(x[,3]))
[1] 3 4

では（当然ながら）意図する結果を得られませんでした。

単純にforで１行ずつ回すしか無いのでしょうか？

↑

データハンドリング　縦横変換†

蔵正 (2007-08-07 (火) 09:56:23)

Aのデータを読み込み、Bのようにデータ変換したいのですが、どうすれば良いでしょうか、ご教示いただければ幸いです。

A----------------------------------------------------------------
PATNO VISIT SUBSCALE QUESTION RESP
1103 Baseline Total Score 99 0.571428571
1103 Month 1 Total Score 99 1.214285714
1103 Month 6 Total Score 99 1.4
1103 Month 9 Total Score 99 0.547619048
1104 Baseline Total Score 99 2.048780488
1104 Month 1 Total Score 99 1.095238095
1104 Month 3 Total Score 99 1
1104 Month 6 Total Score 99 0.904761905
1104 Month 9 Total Score 99 0.880952381
1105 Baseline Total Score 99 1.571428571
1105 Month 1 Total Score 99 1.666666667
1105 Month 3 Total Score 99 0.976190476
1105 Month 6 Total Score 99 0.928571429
1105 Month 9 Total Score 99 1.023809524
　　:

B----------------------------------------------------------------
PATNO Baseline Month 1 Month 3 Month 6 Month 9
1103 0.571428571 1.214285714 1.142857143 1.4 0.547619048
1104 2.048780488 1.095238095 1 0.904761905 0.880952381
1105 1.571428571 1.666666667 0.976190476 0.928571429 1.023809524
　　:

↑

バックスラッシュを置換できますか？†

ぼう (2007-08-06 (月) 16:51:00)

winのファイルパス名が入ったテキストをいじる際、バックスラッシュを置き換えたいと思います。以下のようにすると「??」を置換することはできますが、エスケープシーケンスを減らしても、「?」1文字だけを置換することができません。私の正規表現の理解が不十分なのかもしれませんが、どなたかご教示いただければ幸いです。環境はWinXp R2.51です。

> x <- "aa?bb??cc"
> gsub("????", "/", x)
[1] "aa?bb/cc"    #本当は"aa/bb??cc"としたい
↑

中央値でソートして箱ヒゲ図を表示するには†

sh (2007-08-04 (土) 13:43:20)

グループの中央値でソートしての箱ヒゲ図を表示するにはどうすればよろしいでしょうか？

boxplot(count ~ spray, data = InsectSprays, col = "lightgray") # 1
boxplot(decrease ~ treatment, data = OrchardSprays, log = "y", col = "bisque") # 2
# 1を2のように表示したい
↑

降順のソート関数†

* (2007-07-31 (火) 09:10:01)

降順のソートの関数はあるのでしょうか？
探しているんですが昇順しか見当たりません。
探し不足かもしれませんが、よければ回答お願いします。

↑

関数の出力を表示させない†

RYUI (2007-07-30 (月) 16:37:13)

Rのバージョンは2.5.1です．
非計量多次元尺度構成法を行っているのですが，そこで次元数PとストレスSの折れ線グラフを描く関数を作ってみました．

library(MASS) # isoMDS関数を使用するため
#類似度データ
ruijiData <- data.frame(
  c(0,1,5,3),
  c(1,0,2,4),
  c(5,2,0,5),
  c(3,4,5,0)
)
#次元数とストレスのプロット関数
plot.mds.stress <- function(data) {
  stress <- NULL
  for(i in 1:c(nrow(data)-1)) {
    mds <- invisible(isoMDS(dist(data),k=i))
    ps <- c(i,mds$stress)
    stress <- rbind(stress, ps)
  }
  plot(stress,main="ストレスのプロット",xlab="次元数 P",ylab="ストレス S",type="o")
}

plot.mds.stress(ruijiData) を実行すると，図と以下の（isoMDS関数による）出力が得られます．

initial  value 13.343941 
iter   5 value 0.313496
final  value 0.006694 
converged
initial  value 0.000000 
final  value 0.000000 
converged
initial  value 0.000000 
final  value 0.000000 
converged

ここで，上記の出力が表示されないようにしたいのですが，invisible関数ではできませんでした．
どのようにしたら，出力を消すことができるのでしょうか？ よろしくお願いいたします．

↑

3次元グラフで座標値に文字列ラベルをつける方法†

niko (2007-07-27 (金) 15:41:34)

scatterplot3dを使って3次元グラフを作成しています．
グラフ中に表示される座標値に文字列ラベルをつける方法を探しています．
2次元グラフの場合は

plot(ai.i$points, type="n")
text(ai.i$points, labels=rownames(ai.i$points))

としていたのですが，3次元グラフ（scatterplot3d)での方法を見つけられずにいます．
どなたかご教示いただけないでしょうか．

↑

barplotの凡例の順序について†

K (2007-07-24 (火) 16:15:10)

以下のようなスクリプトで、下図のヒストグラムを作りました。

  x <-    c( 4, 6, 8,10,12,14,16,18,20,22,24,26,28,30,32)    #直径階
remain<-  c( 8, 3, 7, 5, 4, 3, 5, 3, 4, 0, 2, 2, 3, 1, 1)      #スギ間伐前本数
fell  <-  c( 5, 2, 5, 3, 3, 2, 4, 2, 3, 0, 1, 1, 2, 1, 1)      #スギ間伐後本数
dead  <-  c( 3, 1, 1, 1, 2, 3, 1, 1, 1, 0, 1, 1, 2, 1, 0)      #スギ(Cryptomeria japonica D.Don)枯損木
DBH_dist <- function(x, remain, fell, dead, title, xlabel, ylabel){
  par(ps=20)
  par(mar=c(5,5,3,3))
  dummy <- x * 0
  sum   <-  matrix(rbind(remain, fell, dummy),nrow=3, length(x))
  yrange <- c(-max(dead)-10, max(remain + fell) + 10)
  par(new=F)
  barplot(-dead, names.arg=x,  ylim=yrange, col = "gray", axis.lty = 1, xlab=xlabel, ylab=ylabel, main=title)
  par(new=T)
  barplot(sum, ylim = yrange,  col = c("black","white","gray"),
           axes=F, ann=F, legend=c( "残存木","伐採木", "枯損木" ))
}
title = "テスト用のダミーデータ"
xlabel="DBH"; ylabel="直径階別本数分布(本/plot)"
DBH_dist(x, remain, fell, dead, title, xlabel, ylabel)


Kuni_070724.png

できれば凡例を上から順に「伐採木」「残存木」「枯損木」のように配置したいのですが、
legendのヘルプを見てもよく分かりません。どのようにすればよいのでしょうか。よろしくお願いします。

↑

モード/最頻値†

とくめい (2007-07-23 (月) 17:43:33)

Rはsummary関数でmedianとmeanを出してくれますが、modeは出してくれません。
modeを計算する関数はあるのでしょうか？

↑

Rmapのインストール†

内田　直樹 (2007-07-20 (金) 17:25:44)

Rmap_1.1.0zipをインストールしましたがRのversion2.4.1では”有効なパッケージではありません”となってしまいます。ためしにversion1.9にインストールするとlibraryコマンドでロードできます。version2.4.1が日本語表示なのでこれを使いたいのですが、現在のversionにインストールする方法はあるでしょうか？

↑

barplotによるヒストグラムの作成方法について†

K (2007-07-19 (木) 17:07:04)

Rで以下のようなヒストグラムを作りたいのですが、思うようにいかずに困っています。

Kuni_070719.png

barplotを使ってやろうとしているのですが、helpをみても見当が付きません。
barの幅を0.4にして、間伐前と後のbarの位置をそれぞれ元の位置から
±0.25にしてやればできそうな気がするのですが、barplotのパラメータに
width = rep(0.4, length(x)) と書き加えても何も変わらないし、offsetは縦方向へのシフトだし...。
こんなイレギュラーなグラフは無理かなぁ、と思う一方、Rならきっとできるに違いないとは思うのですが。
どなたかお知恵を貸して下さい。よろしくお願いします。
以下が現状のスクリプトです。

# barplot03.R
#直径階別本数分布を、間伐前と間伐後を上向きに並べて、枯損木を下向きにして表示したいのだが...
x <-  c(4,6,8,10,12,14,16,18,20,22,24,26,28,30,32)         #直径階
N0 <- c(3,0,1, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0)    #スギ枯損木
L0 <- c(2,0,3, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1)   #広葉樹枯損木
N1 <- c(8,0,7, 5, 4, 0, 0, 3, 0, 1, 0, 2, 0, 0, 1)         #スギ間伐前本数
L1 <- c(2,0,4, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1)         #広葉樹間伐前本数
N2 <- c(6,0,5, 3, 3, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0)         #スギ間伐後本数
L2 <- c(1,0,3, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1)         #広葉樹間伐後本数
sum0 <- N0+L0; sum1 <- N1+L1; sum2 <- N2+L2
before <- matrix(rbind(N1, L1),2, length(x))
after  <- matrix(rbind(N2, L2),2, length(x))
koson <- matrix(rbind(N0, L0),2, length(x))
par(ps=20)
par(mar=c(5,5,3,3))
xlabel="DBH"; ylabel="直径階別本数分布(本/plot)"
title = "テスト"
yrange = c(-max(sum0)-1, max(sum1)+1)
par(new=F)
barplot(before, width = rep(0.5, length(x)), space = NULL, names.arg=x, ylim = yrange,
       axis.lty=1, xlab="", ylab="", main="")
par(new=T)
barplot(-koson, names.arg=x,  ylim=yrange,
      axes = F,  xlab=xlabel, ylab=ylabel, main=title)
↑

データフレームの列の（範囲）指定†

PC初心者 (2007-07-19 (木) 04:00:45)

初心者です．（WindowsXP，R-2.5.1）
以下のサンプルデータ("sitsumon")があるとします．

q1 <- c("F","M","M")
q2_1 <- c(60,70,75)
q2_2 <- c(70,75,70)
q2_3 <- c(90,95,80)
q3 <- c(1,0,0)
q4 <- c(3,4,5)
sitsumon <- data.frame(q1=q1, q2_1=q2_1, q2_2=q2_2, q2_3=q2_3, q3=q3, q4=q4)

（＊ここで，実際には，より多くのデータ数（行数）と質問の項目数（列数）があります）

実際に分析する際は，"q2_1"と"q2_2"，"q2_3"，"q3"のみを使用するので，扱いやすいように

UFO <- sitsumon[,2:5]

とすればよいのですが，実際には，列数が多い(100くらい)ので，"q2_1"，"q2_2"，"q2_3"，"q3"が何列目か非常にわかりにくいです．
また，抜き出す列数も実際には，もっと多いですので，

UFO <- sitsumon[ ,c("q2_1","q2_2","q2_3","q3")]

というのも避けたいです．

したがって，

UFO <- sitsumon[,"q2_1":"q3"]

のような感じで，2列目("q2_1")〜5列目("q3")を抜き出したいのですが，できませんでした．

そこで，皆様にご教授お願いしたいのが，「列名」からそれが何列目になるのかを数値で返す関数や，もしくは他のよい解決方法はございますでしょうか？
（できればOS(プラットフォーム)などに依存しないようなやり方がよいです）
愚直な質問ではございますが，ご教授くだされば幸いです．

↑

ベイズ−ＢＵＧＳでのＭ−Ｈアルゴリズムの実行†

azu (2007-07-18 (水) 23:31:51)

ベイズでM-Hアルゴリズムを実行したい場合、
ＢＵＧＳで動かすことは可能なのでしょうか？
ＢＵＧＳの名前から考えるとできなさそうなのですが・・・・
M-Hアルゴリズムを動かすソフトウェアがあれば、
教えてください。
お願いします。

↑

spgrass6のインストールが出来ません_ubuntu7.04†

しまだ (2007-07-18 (水) 20:40:33)

どなたか教えてください。
linux超初心者です。grassとRを使いたくて、ubuntu-7.04をセットアップしました。
ubuntu:7.04
grassのバージョン：6.2.1
Rのバージョン：2.4.1
R version 2.4.1 (2006-12-18)~ i486-pc-linux-gnu
ところが、install.packages("spgrass6")を実行すると、
下記のようなメッセージがでてspgrass6とrgdalをインストールすることができません。
spやmaptoolsは、r-base-devを入れたら、無事インストールできました。
どなたかアドバイスをいただけないでしょうか。

install.packages("spgrass6") 中で警告がありました：
   argument 'lib' is missing: using /usr/local/lib/R/site-library
--- Please select a CRAN mirror for use in this session ---
Loading Tcl/Tk interface ... done
* Installing *source* package 'spgrass6' ...
** R
** inst
** preparing package for lazy loading
Loading required package: sp
Loading required package: maptools
Loading required package: foreign
Error: package 'rgdal' required by 'spgrass6' could not be found
Execution halted
ERROR: lazy loading failed for package 'spgrass6'
** Removing '/usr/local/lib/R/site-library/spgrass6'
The downloaded packages are in
/tmp/RtmpQLX6fo/downloaded_packages
Warning message:
installation of package 'spgrass6' had non-zero exit status 
   in: install.packages ("spgrass6")
↑

boxplotの目盛りと外れ値†

sh (2007-07-15 (日) 11:10:40)

boxplot(c(0.1, 1:100, 1000), log="y")　で、boxplotのy軸目盛りを「1e-01,1e+00,1e+01,1e+02,1e+03」ではなく「0.1,1,10,100,1000」と表示する方法を教えてください。また、○印（外れ値または異常値。ここでは1000）を決定するアルゴリズムを教えてください（標準偏差の何倍以上？）

↑

バイナリデータの扱い方について†

YURI (2007-07-14 (土) 14:43:14)

初心者です。バイナリデータの扱い方が良く分りません。
バイナリファイルから読み込んだデータを4バイト単位にまとめるためにmatrixでncol=4を指定してみました。

> f <- file('/tmp/a.out', open='rb')
> a <- readBin(f, what='raw', n=40)
> a
[1] 7f 45 4c 46 01 01 01 00 00 00 00 00 00 00 00 00 02 00 03 00 01 00 00 00 90
[26] 82 04 08 34 00 00 00 50 1d 00 00 00 00 00 00
> matrix(a, ncol=4)
> dump("a", file="/tmp/a.R")
> system('cat /tmp/a.R')
"a" <-
c(7f, 45, 4c, 46, 01, 01, 01, 00, 00, 00, 00, 00, 00, 00, 00, 
00, 02, 00, 03, 00, 01, 00, 00, 00, 90, 82, 04, 08, 34, 00, 00, 
00, 50, 1d, 00, 00, 00, 00, 00, 00)
> source('/tmp/a.R')
Error in parse(file, n = -1, NULL, "?") : syntax error on line 2
> dput(a)
c(7f, 45, 4c, 46, 01, 01, 01, 00, 00, 00, 00, 00, 00, 00, 00, 
00, 02, 00, 03, 00, 01, 00, 00, 00, 90, 82, 04, 08, 34, 00, 00, 
00, 50, 1d, 00, 00, 00, 00, 00, 00)

matrixの結果は何も表示されません。what="raw"で読み込んだデータには matrixは使えない(1要素のサイズが分らないので?)ということでしょうか。
また、dumpしたものがsourceで読み込めないのも、こういうものなのでしょうか?
結局、"rawデータ"というものは、どのように使えるデータなのかを知りたいです。

↑

散布図を右90度回転†

ssk (2007-07-13 (金) 07:02:13)

plotで右90度回転させた散布図を書くためにはどうすればよろしいでしょうか？
boxplotやbarplotのようにオプションを指定して簡単にできないものでしょうか？

boxplot(1:10, horizontal=T)
barplot(1:10, horiz=T)
plot(1:10, horiz=T) # Warning messages:

環境：R version 2.5.0 (2007-04-23) i386-apple-darwin8.9.1
参考：「2変量散布図とヒストグラム図」「クラスター図を横向きで描画したい」

↑

メタファイルからビットマップファイルへの変更って？†

チョメスケ (2007-07-12 (木) 15:47:06)

最近，Rの勉強を始めたばかりの素人です。質問内容はごく簡単な内容だと思いますが，お分かりになる方がいらっしゃれば，よろしくご指導のほどお願いいたします。

質問内容：MCMCpackのMCMCregressを実行しました。独立変数が4つあるのですが，最初にグラフ化される（全2ページ中の1ページ目）のが切片を含め，x1,x2の3つの時系列プロットとposteriorが表示されます。ここで，次のページへ移行する前に，1ページ目の画像を保存で，どうしてもメタファイル形式しか選ぶことができません（2ページ目は保存形式が多数あるので，問題なありません）。WSwordに張り付けることを考えると，ビットマップ形式の方がストレスがなくて良いと考えているのですが，どのようにすれば1ページ目の画像をビットマップ形式での保存ができるのでしょうか？
いろいろと試行錯誤してみたのですが，できませんでした。もし，できないとしてもメタファイルで保存した画像をビットマップファイルで保存する方法を知っていらっしゃいましたら，ぜひご教唆いただけないでしょうか？
よろしくお願いいたします。

↑

barplotの横軸目盛り†

K (2007-07-12 (木) 06:37:14)

　barplotについての質問です。

x <-    c(4,6,8,10,12,14,16,18,20,22,24,26,28,30,32)
Grp1 <- c(8,0,7, 5, 4, 0, 0, 3, 0, 1, 0, 2, 0, 0, 1)
Grp2 <- c(2,0,4, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1)
Grp3 <- c(3,0,2, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0)
y <- matrix(rbind(Grp1, Grp2, Grp3),3,length(x))
barplot(y, xlab="DBH")

として、下図を作図しました。このbarplotの横軸にベクトルxの値を目盛りとして付けたいのですが、どのようにしたらよいのでしょうか。よろしくお願いします。

Kuni_barplot.png
↑

時系列データ処理の高速化の仕方について†

蓮 (2007-07-11 (水) 22:33:30)

アーカイブの（６）にも似たような内容がありまして恐縮ですが、どうしても自分の力で解決できないため、お力添えいただければと思います。

データ：２分間隔で測定した１年分のデータ、50箇所分。
目的　：週毎に平均した１日の変動の値を算出する。
（0:00の週の平均値、0:02の週の平均値…23:58の週の平均値の720のデータからなるセットを52週分算出する。）

#結果収納用データフレームを作成
y <- x[1,] #xは測定値が入っている262800行50列の行列

#52週間分平均データ計算
for (i in 1:52){
 for (j in 1:720){
  x <- rbind(y,
       apply(x[seq((i-1)*720*7+j, 720*7*i, by=720),], 2, mean))
 }
}

これだと、１回計算するのにかなりの時間（６〜７時間）がかかってしまうなので、なんとか高速化できないかと思っています。
行列を作ってみようともしたのですが、結局、平均するデータの行を指定する行列を繰り返し発生させるプログラムを書くことしか思いつきませんでした。

作業を効率的に行う方法をご教示いただけませんでしょうか。
よろしくお願い申し上げます。

↑

データフレームからルックアップ†

でみあん (2007-07-11 (水) 10:01:43)

いつまで経っても，編集されないので，代理で編集しておく
以下のような二つのデータフレームがある場合、

df1<- data.frame(cbind(c("a","a","a","b","b"),c("abc","def","ghi","jkl","mno")))
df2<- data.frame(cbind(c("a","b"),c(3,1)))

df2の各レコードに対して、それぞれ一列目の値をキーに一定の基準（対応する最後のの行、など）に基づき、df1のレコードを選択してMergeする良い方法はありませんでしょうか？現在は、for　ループを用いていますが実行にかかる時間を短縮する良い方法があればと悩んでいます。

最終的には、

  X1 X2  X3
1  a  3 ghi
2  b  1 mno

のようなデータフレームを作成したいのです。

↑

文字列を分解してベクトルに変換する方法†

学生 (2007-07-10 (火) 18:01:03)

はじめまして。文字列を分割してベクトルに変換したいのですが、strsplit関数で

str <- "1.1 2.3 3.4"
str.v <- strsplit(  "1.1 2.3 3.4", " " )

などを自分では試してみたいのですが、うまくいきません。どのようにすればよいでしょうか？
よろしくお願いします。

↑

平方ユークリッド距離†

KO (2007-07-10 (火) 16:27:11)

非類似度として平方ユークリッド距離を使用したいのですが、helpで調べましたが、ユークリッド距離はありますが、平方ユークリッド距離については記載がありません。
計算可能でしょうか？

↑

アポストロフィーを含む文字列の読み込み†

K (2007-07-10 (火) 09:42:04)

read.table("ImportFile", header=F, sep=";")のようにファイルを読み込む際、以下の例のようにデータの一部にアポストロフィーが入っていると、次にアポストロフィーが出てくる場所までがひとかたまりに読み込まれてしまいます。

2094;0;MACY'S;3;5;7
2095;0;MACY'S;3;4;8

read.tableのオプションなどでアポストロフィーを無視するような設定はできますか？

↑

Rcmd SHLIBでDLLファイルを作成する際にR2.5だとうまくいかない？†

たけ (2007-07-05 (木) 17:45:38)

仮に解決済みなので恐縮ですが、Rの他言語利用の、windowsの場合 にしたがって、dllファイルを作成してみたところ、Rの2.5.1もしくは2.5.0では以下のエラーになり・・・

D:?test?r>Rcmd SHLIB test.c
windres -I C:/PROGRA~1/R/R-25~1.1/include  -i test_res.rc -o test_res.o
c:?MinGW?bin?windres.exe: unknown format type `C:/PROGRA~1/R/R-25~1.1/include'
c:?MinGW?bin?windres.exe: supported formats: rc res coff
make: *** [test_res.o] Error 1

Rの2.3.1では以下のようにうまく行きました。

D:?test?r>Rcmd SHLIB test.c
latex: not found
gcc  -shared -s  -o test.dll test.def test.o -LC:/PROGRA~1/R/R-23~1.1/bin   -lR

これはバグでしょうか？私の環境はWinXPです。ちなみにtest.cとは上記サイトの例のmysum.c と同じです。ここで適切か分かりませんが、結構悩んだのでご報告まで。

↑

任意の相関係数を持つポアソン乱数列を作成したいのですが†

taro (2007-07-02 (月) 00:51:12)

与えられたlambdaを持つポアソン乱数列に対して、任意の相関係数を持つポアソン乱数列を生成したいと思っております。
こちらのサイトの「初級Ｑ＆Ａ アーカイブ(1) - 層別の相関値」や、下記のサイトを拝見して、任意の相関係数を持つ乱数列を生成する関数 gendat2() と、ポアソン分布への適合度検定関数 poissondist() を知りました。

そこで、gendat2() 関数の最初のデータ作成部分を変更して、任意の相関係数を持つポアソン乱数列を生成するようにしてみて、その結果が果たしてポアソン分布に従っているのかを poissondist() 関数で確認しようとしましたが、うまくいきませんでした。

質問させていただきたいのは、
1. あるlambdaを持つポアソン乱数列に対して、任意の相関係数を持つポアソン乱数列の作成方法と、
2. poissondist() 関数の出力結果で、lambda が設定した値とかなりずれてしまっている原因について
の2点です。

下記に実行結果を記載させていただきましたが、gendat2() 関数を真似して実行した変換では、元の2つのポアソン乱数列を、異なるlambdaを持つ2つのポアソン乱数列へとまとめて変換してしまうこととなってしまいました。(確かに、変換後の2つのポアソン乱数列の相関係数は指定した通りなのですが、あるlambdaを持つポアソン乱数列に対して、指定した相関係数を持つポアソン乱数列を作成する、という本来の目的とは違ってしまいました)

また、下記以外にも、いろいろなlambdaを指定してrpois()関数で生成したポアソン乱数列を、poissondist()関数で検定してみたのですが、出力結果のlambdaが、設定したlambdaとかなり異なってしまっていて、原因が分からず困っております。

どうか、上記2点に関しまして、原因と対策をご教授いただけませんでしょうか。よろしくお願いします。
長文失礼しました。

nc <- 1000 # 乱数列の個数
r  <- 0.8  # 相関係数
z <- matrix(rpois(2*nc, lambda=1/10), ncol=2)
# 最初に生成したポアソン乱数列の確認
poissondist(hist(z[,1])$counts)
   $result
          n   lambda
   1000.000    0.419
   $table
       x   d            p            e
   c-0 0 899 6.577042e-01 6.577042e+02
   c-1 1   0 2.755781e-01 2.755781e+02
   c-2 2   0 5.773360e-02 5.773360e+01
   c-3 3   0 8.063460e-03 8.063460e+00
   c-4 4  98 8.446474e-04 8.446474e-01
   c-5 5   0 7.078145e-05 7.078145e-02
   c-6 6   0 4.942905e-06 4.942905e-03
   c-7 7   0 2.958682e-07 2.958682e-04
   c-8 8   0 1.549610e-08 1.549610e-05
   c-9 9   3 7.528501e-10 7.528501e-07
   $result2
    chi sq.     d.f.  P value
   2993.188    2.000    0.000
   Warning messages:
   1:  長いオブジェクトの長さが
            短いオブジェクトの長さの倍数になっていません in: x - e
   2:  長いオブジェクトの長さが
            短いオブジェクトの長さの倍数になっていません in: (x - e)^2/e
res <- eigen(r2 <- cor(z))
coeff <-  solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=T))*res$vectors)
z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff
x <- z %*% chol(matrix(c(1, r, r, 1), ncol=2))
# 変換して生成したポアソン乱数列の確認
poissondist(hist(x[,1])$counts)
   $result
          n   lambda
   1000.000    0.817
   $table
         x   d            p            e
   c-0   0 820 4.417549e-01 4.417549e+02
   c-1   1   0 3.609138e-01 3.609138e+02
   c-2   2   0 1.474333e-01 1.474333e+02
   c-3   3   0 4.015100e-02 4.015100e+01
   c-4   4 161 8.200841e-03 8.200841e+00
   c-5   5   0 1.340017e-03 1.340017e+00
   c-6   6   0 1.824657e-04 1.824657e-01
   c-7   7   0 2.129635e-05 2.129635e-02
   c-8   8   2 2.174890e-06 2.174890e-03
   c-9   9  16 1.974317e-07 1.974317e-04
   c-10 10   0 1.613017e-08 1.613017e-05
   c-11 11   0 1.198032e-09 1.198032e-06
   c-12 12   0 8.156599e-11 8.156599e-08
   c-13 13   1 5.442424e-12 5.442424e-09
   $result2
    chi sq.     d.f.  P value
   5244.693    4.000    0.000
   Warning messages:
   1:  長いオブジェクトの長さが
            短いオブジェクトの長さの倍数になっていません in: x - e
   2:  長いオブジェクトの長さが
            短いオブジェクトの長さの倍数になっていません in: (x - e)^2/e

下記の関数を参考にさせていただきました。

## http://cat.zero.ad.jp/~zak52549/R/haebara/chap04.txt
## 任意の相関係数を作成する関数を定義 gendat2(標本サイズ,相関係数)
gendat2 <- function(nc, r)
{
    # 仮のデータ行列を作る。この時点では変数間の相関は近似的に0。
    z <- matrix(rnorm(2*nc), ncol=2)
    # 主成分分析を行い，主成分得点を求める。この時点で変数間の相関は完全に0。
    res <- eigen(r2 <- cor(z))
    coeff <-  solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=T))*res$vectors)
    z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff
    # コレスキー分解の結果をもとに，指定された相関係数行列 を持つように主成分得点を変換。
    z %*% chol(matrix(c(1, r, r, 1), ncol=2))
}

## 青木先生のサイト http://aoki2.si.gunma-u.ac.jp/R/poissondist.html
## ポアソン分布への適合度の検定
poissondist <- function(d)                              # 度数分布ベクトル
{
        k <- length(d)                                  # 階級数
        n <- sum(d)                                     # データ数
        k1 <- k-1
        x <- 0:k1
        lambda <- sum(d*x)/sum(d)                       # 平均値（=λ）
        result <- c("n"=n, "lambda"=lambda)
        p <- exp(-lambda)*lambda^x/factorial(x)         # 確率
        p[k] <- 1-sum(p[-k])                            # 最後の確率は合計が 1 になるように
        e <- n*p                                        # 期待値
        table <- data.frame(x, d, p, e)                 # 結果をデータフレームにまとめる
        rownames(table) <- paste("c-", x, sep="")
        while (e[1] < 1) {                              # 1 未満のカテゴリーの併合
                d[2] <- d[2]+d[1]
                e[2] <- e[2]+e[1]
                x <- d[-1]
                e <- e[-1]
                k <- k-1
        }
        while (e[k] < 1) {                              # 1 未満のカテゴリーの併合
                d[k-1] <- d[k-1]+d[k]
                e[k-1] <- e[k-1]+e[k]
                x <- d[-k]
                e <- e[-k]
                k <- k-1
        }
        chisq <- sum((x-e)^2/e)                         # カイ二乗統計量
        k <- k-2                                        # 自由度
        p <- pchisq(chisq, k, lower.tail=FALSE)         # P 値
        result2 <- c("chi sq."=chisq, "d.f."=k, "P value"=p)
        return(list(result=result, table=table, result2=result2))
}
添付ファイル: fileiris78975839789.png 2814件 [詳細] fileggplot98.png 2357件 [詳細] filesp-grid-20080112.png 3090件 [詳細] filebarplot2.png 2259件 [詳細] filescatterplot3d-example.png 2349件 [詳細] fileKuni_070724.png 2227件 [詳細] filelegend_K1.png 1349件 [詳細] filecif.gif 3203件 [詳細] fileKuni_barplot.png 2277件 [詳細] fileconsolegamen.png 1584件 [詳細] filehistgram_k2.png 2993件 [詳細] fileggplot99.png 2505件 [詳細] filegroupedData1.jpg 2319件 [詳細] filenlme-plot.png 2408件 [詳細] fileplot999.png 2421件 [詳細] fileKuni_070719.png 2523件 [詳細] filegroupedData2.jpg 2792件 [詳細]
トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-03-25 (土) 11:19:17