R および RjpWiki に関する質問コーナー

過去の記事のアーカイブ

Ｑ＆Ａ (初級者コース)/13 の目次

場合分けに応じた新しい列を追加†

asap (2012-01-12 (木) 18:36:36)

今ある行列irisに、n列の場合分けに応じた新しい列を追加したいのですが、どのようにすればよいでしょうか
新しい列の変数には次のような値を入れたいのですが

if (iris$Species[i] == "setosa") {
  iris$Sp[i] <- 1
} else {
  iris$Sp[i] <- 0
}
↑

geoRパッケージの標本バリオグラムについて†

まち (2012-01-11 (水) 21:49:52)

geoRパッケージのvariogで標本バリオグラムの計算を行うと最大距離を持つポイントが省かれてしまいます。
同じデータでgstatパッケージvariogramを使うと省かれません。
パラメータの指定が悪いのでしょうか?
自己相関が無くなる距離なので推定には問題ないと思いますが気になります。
よろしくお願いします。

library(geoR)
library(gstat)
# データの読み込み
data(ca20)
#   geoR
vcloud <- variog(ca20)
(vcloud)
plot(vcloud)
#   gstat
ca20gstat <- cbind(ca20$coords[, "east"], ca20$coords[, "north"], ca20$data)
colnames(ca20gstat) <- c("east", "north", "data")
ca20gstat <- as.data.frame(ca20gstat)
vfgstat  <- variogram(data ~ 1, loc = ~east+north, data = ca20gstat,
                      cutoff = vcloud$max.dist, width = vcloud$max.dist / 13)
(vfgstat)
plot(vfgstat)

下記の 13 行目で geoR の場合 $n が 12 点でなく 11 点になってしまいます。
geoRの場合

　　$n     $u       $v
1   543   43.77376  55.50552
2  1648   131.3213  78.09254
3  2364   218.8688  97.05838
4  2246   306.4163  117.1554
5  2531   393.9638  129.0646
6  2152   481.5114  146.4889
7  1708   569.0589  151.9713
8  1309   656.6064  153.6482
9   694   744.1539  144.7442
10  343   831.7014  146.8878
11  144    919.249  131.8958
12   59   1006.796  139.2288
13   11   1094.344  139.2727

gstatの場合

     np       dist     gamma
1   543   60.35719  55.50552
2  1648  132.58886  78.09254
3  2364  220.90992  97.05838
4  2246  305.65494 117.15539
5  2531  391.80321 129.06460
6  2152  482.02453 146.48885
7  1708  568.46999 151.97131
8  1309  652.62969 153.64820
9   694  739.59074 144.74424
10  343  826.14460 146.88776
11  144  913.61963 131.89583
12   59  995.10190 139.22881
13   12 1082.94416 165.16667

使用環境

R version 2.12.1 (2010-12-16)
Platform: i386-pc-mingw32/i386 (32-bit)
gstat_1.0-8
geoR_1.6-35
↑

R-2.14以降、Tcl/Tkのスライダーを使うとグラフが一部消える†

hiro (2012-01-10 (火) 16:39:38)

R-2.14.1 for Windows(32bit)を使っています．
R-2.13.2までは、以下のサンプルでスライダーを動かすと、グラフはヒストグラムとカーネル密度共に滑らかに表示されます．
しかし、2.14以降、スライダーを動かしている間は一部のグラフが表示されなくなりました．この例ではカーネル密度が時々消えます．
グラフが複雑化すると、目立って見苦しくなります．
複数のWindowsPC（XP/Win7）で試しましたが、100％再現します．
もし解決策を御存知の方がいらっしゃれば、ご教示頂けないでしょうか．
宜しくお願いします．

library(tcltk)
draw.graph <- function(...) {
  upper <- as.double(tclvalue(Slider))
  hist(x, freq = FALSE, ylim = c(0, upper),
       col = "cyan", border = "white")
  lines(den, lwd = 3)
}
x <- rnorm(200)
den <- density(x)
tt <- tktoplevel()
Slider <- tclVar("0.5")
slider <- tkscale(tt, from = 0.1, to = 1,
            resolution = 0.01, showvalue = TRUE,
            variable = Slider, orient = "horizontal",
            command = draw.graph)
tkgrid(tklabel(tt, text = "upper"), sticky = "w")
tkgrid(slider)
tkfocus(tt)
↑

メモリの割り当てについて†

kate (2012-01-06 (金) 17:47:23)

メモリーのエラーがでます。

エラー：  サイズ 176.4 Mb のベクトルを割り当てることができません 
> memory.limit()
[1] 2047
> memory.size()
[1] 817.03

上記のエラーを疑問に思うのですが、2047MBの割り当てに対して、現在817MB使用で、あと約1200MB使えるという認識は違うのでしょうか？解決策を教えていただけないでしょうか
Win vista 32bit ,R2.14.1です。

↑

gstatライブラリーのvariogramについて†

Yasu (2012-01-05 (木) 14:41:54)

データとして（東経，北緯，高さ，値）の３次元空間の座標値と属性値からなるものを解析対象としています。
このデータに対して，東経，北緯，高さから計算される２点間の距離を尺度としてバリオグラムを計算したいのですが，どのようにすれば良いのか分かりません。
東経，北緯から計算される距離を用いたバリオグラムの計算はできました。

↑

繰り返し処理について†

manabe (2012-01-05 (木) 14:11:49)

month <- c("1月", "2月", "3月", "4月", "5月", "6月",
           "7月", "8月", "9月", "10月", "11月", "12月")

全ての『月』の文字を取り除くにはどのようにすれば良いのでしょうか？
逆に、1:12を使って上記のベクトルを作るにはどのようにすれば良いのでしょうか？

↑

変数を使って，levels()などの関数を使う方法†

ぞうに (2012-01-05 (木) 01:54:36)

いくつかのデータセットがあり各変数のlevels()を求めたいのですが，for文を使ってはできないのでしょうか？よろしくお願いします．
例えば，以下のデータセットがあります

sex    <- c("F", "F", "M", "M", "M")
height <- c(158, 162, 177, 173, 166)
weight <- c(50, 50, 60, 60, 70)
( x    <- data.frame(SEX=sex, HEIGHT=height, WEIGHT=weight) )

求めたいことは以下の事なのですが，

levels(factor(x$SEX))
levels(factor(x$HEIGHT))
levels(factor(x$WEIGHT))

イメージとしては以下のような形で求めたいのですが，うまくいきません．for文を使ってできないでしょうか？

colname <- names(x)
for(i in 1:3) {
  levels(factor(x$colname[i]))
}

以下，参考にしましたが解決できませんでした
http://q.hatena.ne.jp/1316156944

  • x[, colname[i]]とすればいいのです。もっといえばx[, i]でOK。そもそもforを使わずにlapply(lapply(x,factor),levels)でもOKでは？ -- 2012-01-05 (木) 02:03:17
  • 解決しました．特に後者が参考になりました．どうもありがとうございました -- ぞうに 2012-01-05 (木) 03:10:50
  • lapply を二重に使うのは初めてみました。sapply(x, function(z) levels(factor(z))) の方がわかりやすいかな。 -- 2012-01-05 (木) 10:23:55
↑

function定義と代入の違い†

Aizawa (2012-01-03 (火) 00:30:20)

functionの使用方法が理解できず悩んでいます。
abc <- function()(print("ABC"))　と関数定義をすれば、> abcと入力した時点でprint("ABC")が実行されて"ABC"と返ってきそうな気がします。しかし、実際には、function()(print("ABC"))、と返ってきます。これだと、abcにfunction以下の文字式が代入されているだけに想えます。print文が関数として登録されているのではないようでしょうか。何か、関数定義に関して基本的なことを勘違いしているように想われます。全く初歩的な質問で恐縮ですが、どなたか疑問を解いて頂けると嬉しいです。

↑

なぜか発生した誤差で、四捨五入が妨げられる†

MasHARA (2011-12-30 (金) 20:55:44)

Rのroundで小数第2位を丸めるために、round(1.05,1) を実行すると、日本の四捨五入とは異なり、1.1ではなく、1.0になります。
これは、切り上げ後の値の末尾の数字が偶数になるようにする、「JIS Z 8401」の2.c)で決められている丸め方のようですが、Rのプログラミングの練習を兼ねて、自分で四捨五入のための関数を作ることにしました。 正の数のみを扱うものとして、以下のような関数を作りました。

> gonyu <- function(x, digits = 0) {
>	     ceiling(x * 10 ^ digits) * 10 ^ (-digits) - 
+           round(ceiling(x * 10 ^ digits) *
+           10 ^ (-digits) - x, digits)
>	}

作業環境は、Windows XP SP3で、version 2.13.0のRを使用しています。
この関数を試したところ、ほとんどの場合は問題なく四捨五入になりますが、時々、以下のように四捨五入にならない場合が発生しました。
例えば、関数の定義後、次の式を評価すると、

> gonyu(1.45, 1)

返ってきた値は、1.4 でした。

> gonyu(2.45,1)

であれば、ちゃんと 2.5 が帰ってきます。そこで、問題の発生した

> x <- 1.45
> digits <- 1

としてから、関数の中を一つ一つ見ていきました、すると、

> ceiling(x * 10 ^ digits) * 10 ^ (-digits) - x

が 0.05 を返すものの、

> round(ceiling(x * 10 ^ digits) * 10 ^ (-digits) - x, digits)

は、なんと 0.1 を返したのです。

> round(ceiling(1.45 * 10 ^ 1) * 10 ^ (-1) - 1.45, 1)

も、なぜか0.1を返してきます。

> round(0.05, 1)

を実行すれば、ちゃんと、0 (ゼロ)を返してきます。

> ceiling(1.45 * 10 ^ 1) * 10 ^ (-1) - 1.45

は 0.05 を返してくるが、表示されていない桁で誤差を含むのでは、と思い、

> ceiling(1.45 * 10 ^ 1) * 10 ^ (-1) - 1.45 - 0.05

を実行すると、返された値は、4.163336e-17 でした。
後ろにある引き算、2つをカッコでまとめて、

> ceiling(1.45 * 10 ^ 1) * 10 ^ (-1) - (1.45 + 0.05)

とすると、返される数字は、ちゃんと 0 です。
一体、この表示されない誤差は生じたり、生じなかったりするのでしょう。御存じの方が居られれば、御教授のほど、御願い申し上げます。

↑

scatterplot3dについて†

takayuki (2011-12-27 (火) 20:39:17)

３次元の散布図を書くのに，scatterplot3dを使用しています。Z軸についてグラフを回転させたいのですがそのようなことは可能でしょうか。もし可能であればご教授お願い致します。

↑

スクリプトからのコンソール画面消去は実行することが出来ますか？†

Montecarlo (2011-12-26 (月) 19:35:10)

現在では以下のマウス操作またはショートカットにてコンソール画面を消去することが出来る事は確認しています。

「編集」→「コンソール画面を消去」
または
「Ctrl」＋「L」

この作業をスクリプトの命令から再現することは出来ますでしょうか。
ご存じの方がいらっしゃいましたら、よろしくお願い致します。

↑

例えば dgamma の lse を求める†

(2011-12-22 (木) 17:40:53)

回帰解析をしたいのですが、plot(x,y) を見ていたら、ガンマ分布の pdf を当てはめてみたくなりました。２つパラメータがありますが、それぞれを最小二乗法で求めるにはどうしたらいいのですか？

さらに、例えば勝手に自分が作った関数のパラメータを使うことも可能ですか？でっちあげの例ですが、

y <- function(x,a,b){ log(a*x + x^-(a/(2*b)) }

として、x と y が与えられた時に、a と b を推定したいということです。

よろしくお願いします。

↑

type="b"のplotと整合するlegend()の使い方†

T (2011-12-21 (水) 20:02:30)

超初心者質問で大変恐縮なのですが、ヘルプを読んでも解決せず、質問します。
例えば、

matplot(outer(1:5, 1:3), type="b", pch=1:3, lty=1:3)
legend("topleft", LETTERS[1:3], pch=1:3, lty=1:3, col=1:3, seg.len=4)

とすると、左上に凡例が表示されますが、線にシンボルが上書きされた type="o" のような表示になります。
グラフは type="b" で描いていますので、凡例も同様にしたいのですが、legend() のヘルプを読む限りは、type="b" と同じ表現になるオプションはないようです。
もちろん、legend(..., type="b") とすると、エラーになります。
上記 2 行の出力画像を貼った方がよければ貼ります。
ヘルプの見落としのような気もいたしますが、よろしくお願いします。

↑

cat関数の特殊文字について†

Montecarlo (2011-12-21 (水) 18:10:04)

文字列を表示する「cat」関数があると思いますが、大文字小文字のアルファベットを調べたところ、以下３個の文字が何を表しているかGoogle等で調べてみても分かりませんでした。
「\b」「\r」「\v」
実行コード

cat("\a") #音が鳴る
cat("\b") #不明
cat("\f") #フィード文字（半角スペース？）
cat("\n") #改行
cat("\r") #不明
cat("\t") #タブ
cat("\v") #不明（半角スペース？）

この３個特殊文字の意味をご回答頂ければ幸いです。よろしくお願い致します。

↑

結果の解釈†

いぬ (2011-12-21 (水) 10:18:21)

成績が勉強時間の2乗に比例しているのか３乗に比例するのか調べたいです。
式としては
ｙ＝a0+a1*x+a2*x^2
ｙ＝a0+a1*x+a2*x^2+a3*x^3
という式です。

手元の資料からは決定係数と調整済みの決定係数が１に近い方がいい、特に調整済みの決定係数が近いほうがいいと読み取れます。
また、ｔ検定、F検定、ＡＩＣも見たほうがいいと思うのです。

どれを１番優先してみればいいのでしょうか？
それとも違うものを見たほうがいいのでしょうか？
どなたか教えてください。

あとＲでanovaの結果の見方がわかりませんので教えていただきたいです。

↑

RDA解析の解釈†

ほし (2011-12-21 (水) 00:13:46)

初めて投稿させていただきます。
パッケージVeganのRDAを用いた解析・解釈に苦戦しております。
20以上のデータを用いて解析を行い、plotで表にしました。
しかし、表には4つ程のデータ名しか表示されません。
また、summaryで詳細な情報を見ても、表示されたデータ名についてしか書かれていません。
これは、表示されなかったデータは取るに足らないデータであったということなのでしょうか？
初歩的な質問で恐縮ですが、どなたかアドバイスの程よろしくお願いいたします。

↑

deprecatedとは？†

S (2011-12-20 (火) 18:14:37)

表を作って平均をとろうとすると、次のような警告がでます。

mean(<data.frame>) is deprecated

R2.13では出なかったのですが、R2.14で出るようになったと思います。これは何でしょうか？

↑

一般化線形モデルの推定結果の図示（offset項を入れた時）†

カス (2011-12-20 (火) 14:03:20)

offset項を入れた時の推定結果を図示するには、どのようにすればよいのでしょうか？
二項分布でリンク関数を"cloglog"にした場合と、ポアソン分布（リンク関数はデフォルトの"log"）の場合についてご教授いただけないでしょうか？
例えば、二項分布でoffsetなしの場合、

bi <- glm(y ~ x, family = binomial(link = "cloglog"))
lines(x, 1-exp(-exp(b + a*x)))

とすれば良いと思うのですが、二項分布でoffsetありの場合

bioff <- glm(x ~ y + offset(log(z)),
             family = binomial(link ="cloglog"))

この推定結果の図示はどのようにすればよいでしょうか？
また、同様にポアソン分布でオフセットなしの場合、

po <- glm(y ~ x, family = poisson)
lines(x, exp(b+ax))

とすればよいと思うのですが、

pooff <- glm(y ~ x + offset(log(z)), family = poisson)

の場合はどうでしょうか？ よろしくお願い致します。

↑

read.csvで、header=TRUE, row.names=1として読み込んだとき†

matak (2011-12-15 (木) 12:15:41)

CSVファイルを例えば、

> data <- read.csv("read_sample.csv", header=TRUE, row.names=1)

という風に読み込んだ時、CSVファイルの一番左上のセルだけ読み込まれないですよねぇ。
結果、読み込んだ行列で何らかの演算を行って、書き込んだCSVファイルも一番左上のセルが空白になってしまいます。
どうすれば一番左上のセルも読み込み&書き込みできるでしょうか。
読み込むCSVファイルの1行目および1列目は文字列で、第2行第2列目以降から数値データなので、header=TRUE, row.names=1で読み込んだ方が都合がいいのですが。

↑

行列の行に対してcumsumを適用することは出来るでしょうか？†

Montecarlo (2011-12-12 (月) 21:25:13)

初の投稿で、なにぶん至らない点があるかと思いますがよろしくお願いします。
使用環境は以下です。

R version 2.14.0 (2011-10-31)
Platform: i386-pc-mingw32/i386 (32-bit)
 
locale:
[1] LC_COLLATE=Japanese_Japan.932  LC_CTYPE=Japanese_Japan.932
[3] LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C
[5] LC_TIME=Japanese_Japan.932

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

質問になりますが、行列の行に対してcumsumを適用することは出来るでしょうか？

> x <- c(1, 1, 1, 0, 0, 0, 0, 1, 1, 1,
         1, 1, 0, 1, 1, 1, 1, 0, 1, 0,
         1, 0, 0, 1, 0, 1, 1, 0, 1, 0,
         0, 0, 0, 1, 1, 1, 0, 1, 1, 0,
         0, 1, 1, 0, 0, 1, 1, 0, 0, 1,
         0, 1, 0, 1, 1, 0, 1, 0, 0, 1,
         1, 0, 0, 0, 0, 1, 0, 0, 1, 0,
         1, 1, 0, 0, 1, 0, 0, 0, 1, 1,
         1, 0, 1, 0, 0, 1, 1, 0, 0, 1,
         1, 0, 0, 0, 0, 1, 1, 1, 0, 0)
> y <- matrix(x, nrow=10, ncol=10, byrow=TRUE)

を実行すると、私の意図通りに10*10の行列ができあがります。
ここからが詰まっているところです。
この「y」の行に対するcumsumを適用したいと思っています。
列に対するcumsumは

> apply(y, 2, cumsum)

によって計算出来ていることは確認しています。
ですが

> apply(y, 1, cumsum)

ではうまくいかないのです。（結果はapply(y, 2, cumsum)に似ていますが違った答えが返ってきます。）
質問事項としてまとめると以下の点となります。
１．行に対してcumsumを適用することは可能か（cumファミリー全般）
２．列はうまくいくのに行でうまくいかない事からすると、cumsumの挙動はどのようになっているのか
以上２点について分かる方がいらっしゃいましたらお願い致します。

↑

readBin関数のサイズ指定†

Alin (2011-12-12 (月) 17:05:20)

C言語のint型、4バイト長のバイナリデータ(sample.bin)をreadBin関数を使って読み込もうとしているのですが、

> conn <- file("sample.bin", "rb")
> readBin(conn, integer(), n=100*100, size=4)

のようにして読み込むと正常な値が読み込まれず(100*100はデータの要素数です）、

> readBin(conn, integer(), n=100*100, size=2)

のようにして、sizeを２に指定すると、正常な値が読み込まれます。
これはRで扱うデータ型（？）の基本的な概念が、C言語等と異なっているからなのでしょうか？

↑

ccfでp-value†

vivaTS (2011-12-12 (月) 13:32:48)

R64 2.14.0 (Lion)を使用しています。
ccfで計算される２変数の時系列相関のp値, 95%信頼区間を求めるには，どのような方法があるのでしょうか。

↑

プロキシサーバーの設定について†

波 (2011-12-06 (火) 17:13:05)

お世話になります。
CentOSにR 2.13.1をインストールし、プロキシサーバー経由でパッケージのアップデートをしようとしています。
Rのインストールディレクトリに以下のコマンドを記述した.Rprofileを保存しているのですが、update.packeges()で「HTTPステータスは407 Proxy Authentication Requiredです」といわれて接続できません。

options(CRAN="http://cran.md.tsukuba.ac.jp/bin/windows/contrib/2.13")
Sys.setenv("http_proxy"="http://(アドレス):(ポート番号)")

また、.Rprofileを削除してRの起動直後の画面に上記のコマンドを入力しても接続できませんでした。
Windows XPのR 2.13.1ではマイドキュメントに上記の.Rprofileを保存することでアップデートが可能なことを確認しています。
Linuxの場合はWindowsと異なる設定が必要なのでしょうか？対応策をご存知であれば教えてください。
宜しくお願いします。

↑

文字列ベクトルからの共通要素の抽出†

K (2011-12-06 (火) 10:29:04)

こんにちは。WinXP R 2.13.2を使用しています。
文字列ベクトルA,B,C・・・があり、これらに共通する要素を抽出したいと考えています。
例えば以下のような場合は"aaa"が共通しているので、"aaa"を抽出したいということです。

A <- c("aaa","aaaa","aaaaa")
B <- c("bbb","aaa","bbbbb")
C <- c("aaa","ccc","cccc")

実際には文字列ベクトルとその要素は数百個あります。apply()を使えばスマートに記述できそうな気はするのですが、どうにもできませんでした。
どなたか御助言を宜しくお願いします。

↑

ヘッダに括弧が含まれるCSVから正しくヘッダを取得したい†

ひろ (2011-12-05 (月) 12:55:27)

ヘッダに括弧が含まれるCSVから正しくヘッダを取得したいのですが，方法がわからずに困っております．
下記例ではファイル中のヘッダは「ほげ（ほげ）」なのですが，read.csvで読み込むと「ほげ．ほげ．」と括弧→．となり困っております．
もし改善方法をご存知でしたらご教示いただければ幸いです．
環境：win7 64bit + R 2.14.0 64bit
調べたこと：?read.csvや?read.tableの結果からparenthes,bracket(丸かっこ)の単語がないか．googleで本サイト(site:)指定/指定なしに対して，「括弧 R言語 ヘッダ名or列名」で調べてみました．
testCSV.csvの中身は下記通りです．

ほげ（ほげ）
1
2
3


> read.csv("input/testCSV.csv", header=TRUE)
  ほげ.ほげ.
1          1
2          2
3          3
↑

column の名前を入力せずに lm()†

笹川よつ斗 (2011-12-05 (月) 12:42:13)

例えば、

x <- data.frame( matrix(rnorm(500), ncol=10) )
names(x) <- c( sample(LETTERS[1:26], 9, FALSE), 'y' )

というデータがあります。
column の名前の 'y' 以外を入力せずに lm() を使うことは可能でしょうか？
大きなデータフレームから column をサンプルして繰り返し解析する、ということを試しています。 for loop の中で応用することになります。　今のところ以下のようにやっていますが、もう少し簡潔な方法があるものかと思い投稿しました。

f <- paste( 'y ~', paste(names(x)[1:9], collapse='+') )
lm( as.formula(f), data=x )
> x <- data.frame(matrix(rnorm(10000),ncol=100))
> Formula <- unique(sapply(seq(10^4),FUN=function(i) paste("X100 ~ ",
             paste(sort(sample(names(x)[-100],9)),collapse="+"))))
> res <- lapply(seq(10^4), FUN=function(i) lm(Formula[i],data=x))
> Formula[9999]; res[[9999]]
[1] "X100 ~  X20+X25+X30+X41+X55+X6+X68+X80+X99"

Call:
lm(formula = Formula[i], data = x) 

Coefficients:
(Intercept)          X20          X25          X30          X41          X55  
   0.054830     0.020898     0.036853     0.025012    -0.007347    -0.138181  
         X6          X68          X80          X99  
   0.013766     0.071437    -0.010897     0.083487
↑

ライブラリーscatterplot3について†

斉藤正 (2011-12-02 (金) 12:26:18)

これまで scatterplot3 を使って、計算結果を3次元表示していました。Rのバージョンを2.14.0に上げたら、

> library(scatterplot3d)
以下にエラー library(scatterplot3d) : 
  '‘scatterplot3d’' という名前のパッケージはありません
　　　：  関数 "scatterplot3d" を見つけることができませんでした

というエラーメッセージがでて、図が描けなくなりました。対処法をお教えください。

↑

Rserveを使ってload()をする回数を抑える方法について†

@garuby (2011-11-30 (水) 18:58:49)


すでに重回帰モデルやSVMモデルを作成しており、このモデルを使って新しいサンプルの目的変数を予測することを現在Rscriptで行なっています。この方法の場合、毎回重回帰モデルなどを保存したファイルをloadして、目的変数を予測させているので、いつも計算するのに時間がかかります。

計算する量は、サンプル数が１つしかない場合もありますし、10万を超える場合もあります。複数の人が利用するため頻度も多いです。

そこで、Rserveを使って最初にモデルをloadして予測ジョブごとでわざわざloadすることを避けられないか検討しています。ヘルプを見ますと、セッション番号を共有すれば何とかなりそうではないかと思っています。が、まったく分からず立ち往生しています。

このようなことはRserveで実現できるのでしょうか？あるいは他の方法で解決するのでしょうか？

↑

Rで出した結果の見方について†

さ (2011-11-29 (火) 19:09:19)

・Rでlm関数を使って分析をしました。summaryで結果を出したのですがEstimateの所の数字をどのようにとらえていいのかがわかりません。＊の数も同じです。この数字や＊は大きかったり、多いほうが相関関係があるのでしょうか？
・もうひとつstep関数も用いました。これはlm関数よりも詳しい値が出ると考えてよいのでしょうか？
・最後に題とは少し変わるのですが、目的変数が見付からないときにはなしでそれぞれの相関をもとめる（lm関数とstep関数を用いて）ことは可能でしょうか？例えば県ごとの脳梗塞の死亡率と癌の死亡率と糖尿病の死亡率しか出ていない場合です。県を目的変数にしようとしましたが、数字でないせいかエラーになりました。

どなたか教えてください。よろしくお願いします。

↑

非線形最適化について†

ｙｙｙ (2011-11-29 (火) 15:30:47)

お世話になります．

例えば「3*x^2-4*x*y+3*x+6*y^2+8*y」のような2変数の最適化を行うにはどうすればよいのでしょうか．
uniroot()は1変数でしか使えませんよね？
どなたかご教示下さい．宜しくお願いします．

↑

集合演算について†

zoo (2011-11-29 (火) 11:55:26)

お世話になります。

・集合xと集合yがある
・サイズは x > y で、yはxに全て含まれる
・xに含まれるyをyの順番を維持したままxから抜き出したい
という処理をしようとしています。
たとえば

x <- c("A","B","C","D","E")
y <- c("D","E","B")

x[x %in% y]

このように入力すると「"B" "D" "E"」が得られますが、yの順番である「"D" "E""B"」として出力するにはどうすればよいでしょうか？
どなたかご教示ください。宜しくお願いします。

↑

Mac OS X上のpng()で文字化け†

m (2011-11-29 (火) 02:23:48)

Mac OS X（10.6.8）でR2.14.0を使用しています。plotがデフォルトの画面上（Quartz）では問題無いのですが、png()では文字化けしてしまいます。

png("hoge.png", width=640, height=640, unit="px")
plot(1,1)
text(1,1,"ほげ")
dev.off()

こうすると「ほげ」部分が□□になってしまうのですが、何か見落としているためでしょうか？
ターミナルから「R」と入力して起動した場合も、Rのアプリケーション上の「Rコンソール」でも同じ結果になります。
よろしくお願いいたします。

↑

データの型を指定して出力†

shannon (2011-11-28 (月) 00:13:55)

既出でしたらすいません。
Rではデータ型がinteger(整数),numerical(実数),complex(複素数),character(文字)に分類されるとのことですが、
C言語などで使われるfloat(32bit),double(64bit)でビット幅を固定して、データを書き込むことはできないのでしょうか？
writeBin関数ではそのようなオプションがないように思うのですが、
書きこむ前のデータのビット幅を何らかの方法で固定して、その後writeBin関数で書き込む、という流れになるのでしょうか？
分かりにくい質問かもしれませんが、よろしくお願い致します。

↑

並列処理したい関数における外部変数の扱い†

QDY (2011-11-24 (木) 22:47:47)

最新バージョンのRには並列処理用のパッケージ parallel が同梱されています。特にマルチコアのパソコンでは便利だと思いますが、並列処理したい関数が(その関数にとっての)外部変数を含む場合はエラーになるようです。おそらく子プロセスが親プロセスの実行環境を知らないからだと思いますが(?). 関数の実行環境(外部変数)をすべての子プロセスで共有する簡単なおまじないがあるのでしょうか。（Linux, R2.14.0）

> library(parallel)
> cl <- makeCluster(rep("localhost", 4), type="SOCK")
> y <- 3
> tmp <- function(x) x+y
> clusterApply(cl, 1:4, tmp)
 以下にエラー checkForRemoteErrors(val) : 
  4 nodes produced errors; first error:  オブジェクト 'y' がありません
> stopCluster(cl)
↑

コマンドのコンソール出力を抑えるには†

バイソン (2011-11-22 (火) 13:20:05)

Linux(CentOS5.5)上で、R version 2.11.1を使用しています。

a = 1
cat(a,"\n")

上のような内容のRファイル(test.r)を作成して、

> R --vanilla --quiet < test.r

のようにして実行すると、コンソール上には実行コマンドとcat関数による出力の両方が表示されてしまいます。

=====コンソール表示======
> a = 1
> cat(a)
1
=========================

cat関数で出力する値だけをコンソール画面に表示させたいのですが、R --helpを見てもそのようなオプションが無いように思います。
sink関数を使って、出力先を別ファイルに変更させて確認するという方法もあるとは思うのですが、 出力させたい部分で一回一回sink関数で囲うのが面倒なので、
できればRコマンドのコンソール出力を抑えて、cat関数等で出力させる値(上の例では1)のみをコンソールに表示できないものかと考えているのですが、良い方法はありますでしょうか？

↑

persp()で軸のタイトルの文字サイズの変更†

まあくん (2011-11-22 (火) 09:05:29)

タイトルの通りなんですが、関数persp()で軸のタイトル（xlab,ylab,zlab）が、角度によってお互いが重なってしまい、見にくくなってしまうので、文字サイズの変更もしくは、表示位置の微調整ができないでしょうか？

↑

plyrのロードができない†

小田ちん (2011-11-19 (土) 00:51:34)


ggplot2を使いたいのですが、plyrのロードで失敗します。
MacのR64を使っています。どうすれば使えるようになるのでしょうか？
よろしくお願いいたします。

> library(ggplot2)
要求されたパッケージ reshape をロード中です 
要求されたパッケージ plyr をロード中です 
Error in dyn.load(file, DLLpath = DLLpath, ...) : 
  共有ライブラリ '/Library/Frameworks/R.framework/Versions/2.14/
    Resources/library/plyr/>libs/x86_64/plyr.so' を読み込めません： 
 dlopen(/Library/Frameworks/R.framework/Versions/2.14/Resources/
    library/plyr/libs/x86_64/plyr.so, 6): Library not loaded:
    /Library/Frameworks/R.framework/Versions/2.13/Resources/lib/libR.dylib
 Referenced from: /Library/Frameworks/R.framework/Versions/
    2.14/Resources/library/plyr/libs/x86_64/plyr.so
 Reason: image not found 
エラー：  パッケージ '‘plyr’' をロードできませんでした
↑

Boxプロットで、項目軸を数値にする方法†

T.I. (2011-11-16 (水) 07:11:07)

Rの初心者です。グラフの描き方に苦戦しております。
Boxプロットで、項目軸（X軸）を数値軸にしてグラフを描きたいのですが、方法をご存知の方がおられましたら、教えていただけないでしょうか？

↑

行列の割合表示†

yasushi (2011-11-15 (火) 18:56:52)

Ｒ初心者です。
次の２つの実行結果は同じになると思ったのですが，なりません（２番目の数値がおかしい）。

> x <- matrix(c(1:4), nrow=2, ncol=2, byrow=T)
> (x/apply(x, 1, sum))

y <- matrix(c(1:4), nrow=2, ncol=2, byrow=F)
t(y/apply(y, 2, sum))

２番目のものを，colSumsにしてみても同様の結果でした。R 2.14.0でもR 2.2.0でも同様の結果でした（Window XP上で）。

どういうことなのでしょうか？

↑

ARIMAモデルでの解析結果を詳細に表示したい場合…†

前田 (2011-11-14 (月) 22:45:36)

R初心者のものです。
大学の卒業研究でポンドの解析を行っています。
ARIMAモデルでの解析結果は出せるのですが
解析結果の詳細を表示させたいのですが
どなたかお知恵を貸していただけませんでしょうか？
パワースペクトルなどの結果も見れたらと思います。
宜しくお願い致します。

↑

2階の導関数のもとめかた†

みのむし (2011-10-27 (木) 13:26:24)

Rの初心者です。簡単な質問で大変申し訳ないのです。
ある関数の、2階の導関数をもとめ、それを他の計算に使いたいと思っています。

f <- deriv(~x^2*y, c("x", "y"), function(x, y) { }, hessian=TRUE)

と書くと、4種類の2階の導関数を計算してくれるのはわかるのですが、個別の2階の導関数をどう選択して良いのかわかりません。たとえば、この場合、xとyで1回ずつ微分したものはどう選択すればよいのでしょう？

attr[f, "hessian"][x, y]

こんな風にかいてみたのですが、だめなようです。ご存知の方教えていただけませんか？

↑

連続の検出†

nkoji (2011-10-08 (土) 14:45:51)

ベクトルのある値の連続と各連続における順番を取得したいです
例えば，c("a","a","a","b","c","a","b","a","a")に対してc(1,2,3,0,0,1,0,1,2)を得たいのです。今はこんなコードでやっています。

A <- c("a","a","a","b","c","a","b","a","a")
B <- rep(0,length(A))
A.str <- paste(A, collapse="")
#開始位置と長さの取得
C <- gregexpr("a+", A.str)
#C.mat[1,] 開始位置 C.mat[2, ] 長さ
C.mat <- rbind(C[[1]], attr(C[[1]], "match.length"))
B[A=="a"] <- unlist(apply(C.mat, 2, function(x){c(1:x[2])}))

ベクトルの要素を一つずつ比較するよりは速いと思うのですが，上よりもっと速くなるコードがあれば教えてほしいです。

↑

read.tableで読み込んだ行列がis.finiteでFALSE??†

もうり (2011-10-06 (木) 08:38:06)

> A <- read.table("matrix.dat")
> B <- matrix(A[1,], 3, 3)
> kappa(B, exact=TRUE)

としたところ，

Error in svd(z, nu = 0, nv = 0) : infinite or missing values in 'x'

というエラーが出力されました．
matrix.dat は１０００行9列スペース区切りで実数を並べたテキストファイルです．
そこで，svd 関数の定義を見たところ，

if (any(!is.finite(x))) 
   stop("infinite or missing values in 'x'")

との記述がありました．そこで確認のため，

> any(!is.finite(B))

としたところ結果はTRUE．
ここまでは納得できなくもないのですが，

> kappa(A, exact=TRUE)

はエラーが出ず，しかも，

> any(!is.finite(A))

の結果はTRUEです．

> A <- matrix(scan("matrix.dat"), ncol=9)

と読み込むことによって，

> any(!is.finite(A))
> any(!is.finite(B))

が両方FALSEとなったため一応解決はしましたが，すっきりしません．
原因に心当たりがある方がいらっしゃったらお教えいただければと思います．
最後に実行環境について記述しておきます．

$ uname -sr
Linux 2.6.18-92.el5
$ R --version
R version 2.10.0 (2009-10-26)
↑

LSMEANを求めながら，そのデータのsdを求める†

怪獣 (2011-09-01 (木) 12:21:54)

よろしくお願いします。初級Q&A（１０）のところで，LSMEANについてのやりとりがされていました。私のデータでLSMEANを得ることはできたのですが，グループの変数が３つあり，多重比較をしなければいけません。教えて頂きたいことは，LSMEANをベースにした3つの群のそれぞれのsdの求め方と二次データを使った多重比較の方法です。被験者数はすべて同じ数です。どうぞよろしくお願いいたします。

↑

csvからxmlに変換したときの日本語文字化けについて†

さとう (2011-08-29 (月) 18:16:59)

Rでcsv形式のデータからXMLへ変換しています。日本語のタグやデータも扱いたいのですが，以下のコードを実行すればわかるように，日本語が文字化けします。
(タグは文字化けしませんが，データの値は文字化けします。)
これを防ぐ方法はあるでしょうか？

library(XML)
data <- c("変数1の値です。","value of variable 2.")
xml <- xmlTree()
xml$addTag('データセット', close=FALSE)
xml$addTag("データ", close=FALSE)
xml$addTag("変数1", data[1])
xml$addTag("変数2", data[2])
xml$closeTag();xml$closeTag();xml$closeTag()
cat(saveXML(xml))
↑

座標値の等間隔プロット†

R初心者 (2011-07-25 (月) 17:59:00)

今大学でRを使って研究をしているものですが初心者なので苦戦しています

今、等間隔でない座標値がいくつかあって線で結ばれたデータがあります
このデータを等間隔に（例えば２ｍ間隔とか５ｍ間隔とか）プロットさせるようなプログラミングを作成したいのですがどのようにやったらいいのでしょうか
簡単な質問でしたらすみません。何分初心者なもので・・・
ご指導よろしくお願いします。

↑

Excelから１セルの読込†

kero_10625 (2011-06-22 (水) 12:26:37)

節電対策などで、建物のエネルギー解析が必要になっています。ビル管理システムにExcelで記録された数百から一万程度の多数のファイルの処理を行っております。

Excelの単一のセル、１列あるいは１行だけを読み込む必要が出てきました。

RODBCを用いて、sqlQuery を使うとデータフレームで値を返すためか、空欄のセルがあるとエラーが発生します。

xlsReadWriteパッケージのread.xls()関数を使うと可能なのですが、xlsx形式に未対応であり基本的にはシェアウエアであり将来の対応に不安が残ります。

建築技術者にRを使ってもらおうと考えていますので、Pearlなどあまり追加のソフトのインストールは避けたいです。

何か良い方法をご存じでないでしょうか。

↑

Rで他言語の展開†

ktr (2011-06-21 (火) 21:30:21)

Rで，ある関数を展開し，中身を見ると.CとしてC言語で書かれた関数が使用されていました．Rの関数は関数名をタイプすれば中身を見ることが出来ますが，このように他言語で書かれた関数をRで展開することは可能なのでしょうか?

↑

t.test()の機能†

竹澤 (2011-05-31 (火) 16:27:12)

　t.test()の機能を確認するために以下のようなRプログラムを試しました。

function() {
  nt <- 10000
  pt <- 0
  for(kk in 1:nt){
    set.seed(12+kk*144)
    xx1 <- abs(rnorm(50, mean = 0, sd = 1))
    xx2 <- abs(rnorm(50, mean = 0, sd = 1))
    xx <- c(xx1, - xx2)
    ttest1 <- t.test(xx, mu = 0)
    pp <- ttest1$p.value
    if (pp < 0.05){
      pt <- pt +1
    }
  }
  pt <- pt/nt
  print(pt)
}

正規分布に従うデータを作成する方法がややひねくれています。結果は、0.05に近くなるはずなのに0.0018になりました。正規分布に従うデータを作成するときに普通の方法を使うと0.0517になりました。このような違いが生じるのはなぜでしょうか。

↑

順序ロジットモデルの適合度について†

miya (2011-05-02 (月) 10:11:04)

polr関数を用いて，順序ロジットモデルを推計しています．
SUMMARY()を利用すると，推計結果が表示されます．
その際，AIC指標が表示されますが，尤度比は表示されません．
初期尤度，最終尤度，ならびに尤度比を求めたいのですが，polr関数にどのようなオプションを付ければよいのか，教えていただけると助かります．

もしかすると，定義式にしたがって，自分で計算式を作成して計算するしか方法がないのでしょうか．

↑

散布図のプロットの中で条件によって色を変える方法†

ウルトラ初心者 (2011-04-18 (月) 15:33:39)

BioconductorのDESeqというパッッケージを使ってます。ほとんどコピペ程度でしか使えてません。

plot(res$baseMeanA, res$baseMeanB, pch=20, cex=.1, col = ifelse( res$padj < .1, "red", "black" ))

というコマンドで散布図がマニュアルに載っていて、そのままコピペしてかけていました。
resという表のpadjカラムの値がが0.1より小さければ赤、それ以外を黒という風にプロットされています。

さて、別のパッケージBaySeqで統計解析しか結果もこの散布図に反映させたいと考えました。
そこで、一つカラムを増やしてanalysisとなづけ、DESeqで有意なデータを1、BaySeqで有意なデータを2、DESeqでもBaySeqでも有意なデータを3、それ以外を0というように入力しました。

そこで、analysisカラムをみて1ならば赤、2ならば青、3ならば緑、0ならば黒というように上のplotコマンドを改良したいのです。
たぶん、ifを重ねて使っていけばいけそうなのですが、どのようにググってよい物かわかりません。

全くの丸投げ質問で恐縮ですが、よろしくお願いします。

↑

マルチコアCPUでの、Rの効率的な利用について†

さと (2011-04-17 (日) 05:21:20)

現在、クアッドコアのCPUを搭載したWINDOWSマシンにRをインストールして利用しています。
そのうえで、for文を使って繰り返し計算を行っていますが、非常に時間がかかっております。

Rを立ち上げて計算させても一つのコアしか使わないので、マルチコアの資源を有効に使えておりません。
仕方がないので、Rを4つ一度に立ち上げて、for文の繰り返しを, 1:25, 26:50, 51:75, 76:100というように分けて実行しています。
これで計算時間は短くなるのですが、もっとスマートにやる方法をご教授いただければ幸いです。

UNIX環境で使える multicore というパッケージは見つけたのですが、WINDOWSで使える方法を見つけられませんでした。

こちらの環境ですが、
OS: Windows 7 Professional, Service Pack 1
プロセッサ: Intel(R) Core(TM) i7 CPU M620 @ 2.67GHz
実装メモリ(RAM):　8.00GB
システムの種類：64ビット　オペレーティングシステム

なにとぞよろしくお願いいたします。

↑

Rで作図をする際に、レイヤー表示のように手前に持ってくるキャラクターの指定はできないのでしょうか？†

Tetsu (2011-01-27 (木) 22:59:49)

現在、ts.plotを使用して、時系列データの作図をしております。
その際に、１つの点を黒丸、もう一方の点を白抜き三角でプロットを行っております。
三角が手前に来るようにデータ系列の指定の順番を変えたり、par(plot=new)を使用して重ね書きをしてみたのですがうまくいきません。
結局、三角の中身が白で塗りつぶされていないのが問題だと思うのですが、pchで指定可能な白三角と黒丸を一通り試したところ、全て同様の図が表示されてしまいました。
このような問題を解決する方法は無いのでしょうか?

あと、一点質問があります。
作図をする際に見易くするために、軸に主目盛りと副目盛りをつけようとしております。
このような場合はaxis()を2回使用して、目盛りの長さが異なる軸を重ねて描く以外に方法はないのでしょうか？
教えていただけないでしょうか?
よろしくお願いいたします。

↑

Rコマンダーがインストールできない†

Natsu (2011-01-26 (水) 00:02:56)

Rをインストールして、Rコマンダーを使いたいのですが、パッケージからRコマンダーを選ぶと、以下のように出てしまいます。

install.packages(NULL, .libPaths()[1L], dependencies = NA, type = type) 中で警告がありました： 
 'lib = "C:/PROGRA~1/R/R-2.12.1/library"' は書き込み可能ではありません  
以下にエラー sprintf(msg, userdir) : 
  <83>ｼ<e3>Ν縺吶ｋ縺溘ａ縺ｫ蛟倶ｺｺ逧<84>縺ｪ繝ｩ繧､繝悶Λ繝ｪ 
'%s' 
繧剃ｽ懊ｊ縺溘＞縺ｮ縺ｧ縺吶°<ef>ｼ<9f> に不正なマルチバイト文字があります

どのようにすればインストールできるのでしょうか？
初歩的な質問で申し訳ないのですが、教えてください。

↑

GAM(一般化加法モデル)における変数選択について†

S.K (2011-01-21 (金) 23:19:26)

GAMを用いてデータ解析をしているのですが、変数選択の方法が分からなくて
困っています。web上にはstep.gamというものがあるそうなのですが、どのパッケージに入っているかもわからない状態です。何かご意見を頂けないでしょうか。よろしくお願い致します。

↑

コネクションエラー†

R初心者 (2011-01-10 (月) 09:58:15)

R初心者です。
JavaからRserveを利用して、Rを起動しようとしてます。
そこである程度同時実行をさせて、最終的に200件ほどの要求をRserveになげるようにしています。そこで、問題なのですが、Rserveからコネクションエラーが数十件返ってきます。そのコネクションエラー数は毎回毎回まちまちで数も特定できません。私が調べた限りでは、Rserveでコネクションの制限を設定している部分がないことから現状なぜこのようになるか不明な状態です。
そこでもしこのような現状に対する情報をお持ちでしたら、教えていただけないでしょうか。よろしくお願いいたします。

↑

パネルデータ分析　固定効果†

H.A (2010-12-13 (月) 21:07:38)

はじめまして
大学で企業財務を勉強しているものです。

さて、初歩的な質問で恐縮なのですが、Rでパネル分析を行った際に、企業固定効果の係数を出力することはできないのでしょうか？
各企業ごとの企業固定効果を分析上使用する必要があるため、それらをぜひとも抽出したいのです。

もしコードなどをご存じのかたがいましたらご教授お願いします。

↑

interp1?†

XC (2010-12-11 (土) 12:04:28)

はじめまして
Rを調べ出しの初心者です

MATLAB系のスクリプトをRで動くようにする変換プログラムを作っているのですが、元のスクリプト群で多用されている、スカラを返す補間関数がマニュアルを探しても見付かりませんでした。
例えば、エンジン回転数ー発生トルクの2次元のマップ（回転数1000毎とか）、rpm_tqがあったとして、任意の、例えば2850回転の時の発生トルクであれば、
interp1()を使用して、
tq = interp1(rpm_tq(:,1), rpm_tq(:,2), 2850)
（rpm_tq(:,1)、 rpm_tq(:,2)は、回転数とトルクのベクタになります）
の様になっているので、Rで
tq <- interp1(rpm_tq[1], rpm_tq[2], 2850)の様に書けると嬉しいのですが
皆様のご教示をいただけましたら幸甚です。

↑

ARIMAモデルの定数項について†

arima (2010-10-31 (日) 15:09:12)

直線的なトレンドのある時系列に対して

arima(時系列,order=(1,1,0))

とすると、ＡＲ部分の係数が出力されるのみで定数項が出力されません。

これにpredictを施しても、定数項がないため、将来の値が一定となってしまいます。

定数項を出力させるためにはどうすればよいのでしょうか。

宜しくお願い致します。

↑

splancsのkhatについて†

風の歌を聴け (2010-10-30 (土) 10:23:53)

いつも参考にさせていただいております。
splancsのK関数を求めるkhat関数について質問させていただきます。
ランダム分布のsimulation envelopeの上下限を求めるシミュレーション回数(nsim）はどのように決めるのでしょうか？19回や99回などが多く使われているようですが、回数はどのように決めるのでしょうか？また、回数が多いほうが上下限の幅が狭くなるのでしょうか？谷村晋さんの新刊「空間地理データ分析」（共立出版）を拝見しましたがnsim=99の例がのっていました。
よろしくお願いいたします。

↑

クラスター分析での項目名†

ST (2010-10-29 (金) 23:46:16)

クラスター分析での項目名についてです。

Rコマンダーを使用してデータを読み込んでいます。
クラスター分析を行い、樹形図を表示させた場合、各項目が項目名にならず、1からの連番になってしまいます。（項目名がAA,BB,CCでも1,2,3になるということです。）

どのようにすれば項目名が表示されるようになるのでしょうか？
よろしくお願いします。

↑

リッジ回帰について†

TKです (2010-09-16 (木) 17:54:37)

リッジトレースをRで実行したいのですが、どのようにすればいいですか？

プログラムを作るしかないですか？

リッジ回帰のコマンドは一応あるのですが、実行してはみたものの、あまり使い道がないと思うのですが。

↑

パッケージFEARのdeaでNDRSを実行したいのですが。†

atsuo (2010-07-26 (月) 18:01:11)

�FEARのdeaでは、CCRとBCCとNIRSの結果は出ます。さらに、NDRSの結果を出すことは可能でしょうか？
�deaの分析を行うためにFEARのdeaを改良したいと考えた場合、元のソースプログラムを手に入れる方法はありますか。（もしかしたら、基本的な質問かもしれません。申し訳ありません。Rを最近始めたばかりです。なんとかヒントをいただけませんでしょうか。よろしくお願いします。）

↑

線型等式制約Ax=bにおける最適化問題を解く方法†

kazuaki (2010-06-03 (木) 01:11:15)

目的関数 : x[1]^2+x[2]^2+...+x[n]^2
制約条件 : Ax=b , x>=c
xとcはn次元ベクトル,bはm次元ベクトル,行列Aはn*m行列。
xが決定変数。

上記の最適化問題（最小化問題）を解きたいのですが、適切な手法やパッケージが特定できずにおります。

↑

２元配置分散分析の場合の等分散性の検定†

shigeru (2010-04-09 (金) 06:35:07)

２元配置分散分析を行う前の等分散性の検定についてお伺いします。たとえば、要因A(2水準)と要因B(3水準)で繰り返しのある２元配置実験データが
y A B
y1 A1 B1
y2 A1 B2
y3 A1 B3
・・・
のような形で得られているときに、
bartlett.test(y ~ A * B)
とすると、自由度はdf=1となり、また
bartlett.test(y ~ B * A)
とすると、df=2となり、最初に出てくる要因しか評価されていないようです。AとBの組み合わせ要因（この場合はdf=2×3-1=5）に対して等分散性の検定をするにはどうすればよいのでしょうか。AとBの組み合わせからCという要因を新たにつくって対処することもできるでしょうが、その場合でも少し面倒になるような気がします。簡単に対処する方法はないのでしょうか。ご教示頂ければ大変助かります。

↑

plot()出力の部分拡大†

ohgi (2010-03-20 (土) 10:46:20)

仰木と申します．
Rで時系列グラフをplot()を使って描画していますが，実験最中に取れたデータ（約数十万行分）を描画して，必要そうな部分，注目したい部分を拡大して描画したいとき，Scilabではマウス右クリックでエリア指定の拡大ができるのですが，Scilabではあいにく今回の実験データが最大メモリ許容量を超えてしまうため，Rで同じことができないかと思っています．Rでデータの読み込みと全体描画はできるのですが，グラフ画面上で拡大したい部分をマウス操作で拡大するようなことはできないでしょうか？

もちろん，描画する範囲をスクリプトで描けば再描画することはわかっていますが，次々にその時刻を計算して拡大し，閲覧する作業が必要なので出来ればグラフ画面上で出来ないか，と思っています．

「できない」という答えでもよいので教えていただけませんでしょうか．

↑

Help†

あんどう (2010-03-14 (日) 20:57:00)

Help にある .rds, .rdb, .rdx はどのように開くのですか。

↑

R CMD BATCH < test.R > test.Rout†

boston (2010-03-08 (月) 05:47:02)

test.R は

cat("Hello World")

だが、test.Routに出力されない

R --vanilla < test.R > test.Rout

は出力される BATCHの書き方に問題がありますか

↑

パッケージSOMの仕様†

初心者 (2010-02-15 (月) 02:52:47)

パッケージSOMの、関数somの学習回数は何回なんでしょうか？

↑

lm について†

あんどう (2010-02-07 (日) 12:28:45)

lmの処理を調べようとtest_lma.Rを作成し実行すと、match.call()で

以下にエラー match.call() : 
  'match.call' がその中から呼び出されたクロージャを見付けることができません。

となります。 プログラムの冒頭は以下のとおりで

formula <- y ~ 1 + x
xx <- c(1.386294, 2.079442 ,2.772589, 3.465736)
yy <- c(0.06254298 ,0.71444544, 1.27103337, 1.28339677)
ww <- c( 1.0000000 ,0.5000000, 0.3333333, 0.2500000)

data <- data.frame(x = xx, y = yy, w = ww)
subset <- 0
weights <- ww
na.action <- 0 
method <- "qr"
model <- TRUE 
x <- FALSE
y <- FALSE
qr <- TRUE
singular.ok <- TRUE 
contrasts <- NULL
offset <- 0  

ret.x <- x
ret.y <- y
cl <- match.call()
mf <- match.call(expand.dots = FALSE)

これ以降はlmのままです。

目的は余分な処理を削除したいのです。

↑

関数の追加†

あんどう (2010-02-06 (土) 08:50:29)

libraryに関数を追加する方法を教えてください。
新たな関数を作り既存の関数をリンクさせたいのですが。
.Rをつくり cl <- match.call()を組み込んだのですがエラーになります。新たな関数をlibraryに追加できればエラーがなくなると考えています。

↑

多項ロジット†

こうや (2010-01-18 (月) 11:13:26)

それぞれの説明変数に対する回帰係数が応答変数のカテゴリーごとに出ると思うんですが、それらの係数がどのカテゴリーと対応しているかがよく分かりません。
教えていただけないでしょうか。

↑

グンベル分布推定†

吉川 (2010-01-14 (木) 17:37:24)

Rで極値統計をやっております。
extRemes, ismev を使ってGEVの推定を行っています。
gev.fit(x)を使って推定を行うと、
グンベル分布・ワイブル分布・フレシェ分布の中から推定されます。
これを、グンベル分布のみで推定したい場合にはどのようにしたらよいでしょうか。
ご存知の方がいらっしゃいましたら、よろしくお願いいたします！

↑

hclustでのward計算の原理†

genri (2010-01-03 (日) 13:51:31)

どう計算しているのか原理をお教えください。

この関数はdist構造のデータを使って計算します、たとえば距離マトリクスhogeにたいして
hclust( as.dist(hoge), method="ward")
で結果を返します。

Ward'sの方法では、クラスター間の距離はクラスターの統合によって増加した分散によって定義されています。そこで、距離マトリクスを算出した最初のデータがないと距離が算出できないのではないかと思うのです。

どうやって距離マトリクスからwardの距離を算出しているのでしょう？

↑

2次元正規分布関数のパラメータ推定について†

zhang (2009-12-07 (月) 10:25:06)

おはようございます。張＠弘前大学です。R初心者で何か間違ったらご指摘ください。

今研究で，リンゴの形状推定を行っています。リンゴのお尻のくぼみを三次元計測して，2次元正規分布関数として推定しようと考えています。

R初心者で，色々調べたら，kde2d()やbked2d()でできるらしいが，しかし，最終的に出力してほしいのは中心値（μ1，μ2）と標準偏差(σ1,σ2）です。これでできるでしょうか？

ほかに何か方法がありませんか？宜しくお願いします。

↑

コクランのC検定†

Shigeru (2009-11-21 (土) 17:17:18)

等分散性の検定にコクランのC検定を使った文献をよくみかけますが、どのような方法なのでしょうか。またRでは実行可能でしょうか。お教えいただけるとありがたいです。

↑

平滑化された値の出力について†

b.b. (2009-11-04 (水) 10:50:45)

質問です。
スプライン平滑化を行ってグラフを得たのですが、
xがある値の時の平滑化されたyの値を出力するには
どうしたらいいのでしょうか？

↑

constrOptimのhessianは表示可能でしょうか？†

koba (2009-09-03 (木) 12:51:05)

線形不等式制約付きの最適化関数constrOptimにてhessianでヘッセ行列を表示しようとしてもエラーが出てしまいます。
そもそも存在しないのでしょうか？

↑

Map (mapply)から関数のリストを返すにあたって†

dtak (2009-08-13 (木) 07:06:37)

はじめまして、dtak と申します。
R (version 2.9.1 i386-apple-darwin8.11.1) でのMap関数(mapplyのラッパー関数)の振る舞いについて少し気になることがありましたので、投稿いたします。

Mapを用いて関数のリストを作ろうとしたのですが、Mapの引数fの内部で作られるクロージャでのみ評価されるfの引数が、最後の要素の値で上書きされてしまいました(なんだかややこしいですね、以下の単純化したコードをみてください)。

> Map (function (x) x(), Map (function (y) {function () y}, c(1, 2, 3))) ## list(1, 2, 3) が返ることを期待
[[1]]
[1] 3

[[2]]
[1] 3

[[3]]
[1] 3
>

あらかじめ評価してしまうことで回避することはできるのですが、

> Map (function (x) x(), Map (function (y) {y;function () y}, c(1, 2, 3)))
[[1]]
[1] 1

[[2]]
[1] 2

[[3]]
[1] 3

>

(遅延評価が行われるのであっても)少し奇妙な印象を受けました。これはこのような挙動をするものと考えるべきなのでしょうか、それとも不具合と考えるべきなのでしょうか？
ご存知の方がおられましたら、よろしくお願いいたします。

↑

順序ロジスティック回帰について†

fall (2009-07-07 (火) 13:28:37)

　パッケージMASSのplorを用いて順序ロジスティックを行おうと考えております。
　ここで私が困っていることが以下のような内容であります。
　
（目的変数にhelpでは”factor”を指定しなければならないと書かれており、数値を目的変数に入れるとエラーになってしまいます。目的変数を１，２，３など順序のある数値として解析をおこなうためのパッケージやplorの設定方法があるのであれば、教えて頂きたいと思います。）
　
　どうぞよろしくお願いいたします。

↑

非線形連立方程式の解†

かわお (2009-06-23 (火) 23:14:25)

非線形の連立方程式を解くプログラミングを教えて下さい。
例えば、
x^2+y^2=5, 2×x^2+y=4の解x,yを解く方法が知りたいです。
よろしくお願いします。

↑

Rでエラーログ解析†

たけぼう (2009-06-14 (日) 01:53:04)

R超初心者です。
ある半導体関連装置のサポートエンジニアをやっております。担当している装置はPC上で実行されるプログラムに従って動作するようなシステムです。

担当している装置が出力する日付、時間を含んだエラーログを見る機会が多く、その解析にRを利用できないかと考え入門しました。　担当している装置はログファイルの出力が「コマンド実行履歴」と「エラーログ」に分かれて記録されているので、「エラーが発生したとき、その装置は何をしていたか？」を調べるのに毎回エクセルで悪戦苦闘していましたが、この部分をRで一発で視覚的に表示できないかな？、と考えています。ログのフォーマットは次の通りです(csv)。

[エラーログ]
YYYY/MM/DD, HH:MM:SS, ALARM1
YYYY/MM/DD, HH:MM:SS, ALARM2

          :

[コマンド実行履歴]
YYYY/MM/DD, HH:MM:SS, COMMAND1
YYYY/MM/DD, HH:MM:SS, COMMAND2

          :

「上記二つの時系列ファイルから、両者の因果関係を視覚的に表示したい」
というのがやりたいことなのですが、Rでできそうでしょうか？
ヒント・アイディアだけでも良いので、何か閃いた方はご教示下さいませ。

↑

作図の繰り返しのコマンド†

路地裏の少年 (2009-05-25 (月) 22:57:16)

低水準作図の繰り返しコマンドについて教えてください。

データは以下のようになっています。
緯度　経度　年
35　　156　2000
38　　165　2000
37　　150　2001
38　　148　2001・・・・

年毎の緯度経度のプロット図（20年分、各図の題名にはその年を充てたいと思っています）を書かせたいと思っています。
年ごとに抽出すれば、書けるのですが、20回繰り返せねばなりません。
初歩的な質問で恐縮ですが、よろしくお願いします。

↑

ポアソン分布のロジスティック回帰†

ま (2009-04-16 (木) 19:13:38)

非線形最小二乗法を実施する関数として，nlsがありますが，このポアソン分布版はありませんでしょうか？教えていただければ幸いです。

↑

splancsのzoomについて†

K.K (2009-01-27 (火) 16:15:38)

クラスター分析で、plotされた樹形図の細かい統合部分をsplancsのzoomで見ようとしましたが、樹形図が表示されないまま拡大されてしまいます。
plotされた樹形図をzoomで細かく表示させる方法はありませんか？

↑

順序ロジット　関数 polr() について†

y.t (2009-01-15 (木) 18:10:28)

MASSパッケージにある順序ロジットの関数 polr() に 定数項を入れたいのですがどうすればよいでしょうか？

↑

BRugsでのAIC，DICの計算について†

yuki (2008-12-16 (火) 01:49:28)

WinBUGSではDICは簡単に計算できるコマンドがあるようですが，BRugsでもそのようなコマンドはあるのでしょうか？

↑

R-2.8.0 x64 Linux上でのビルドについて†

大津起夫 (2008-12-05 (金) 17:33:15)

R-2.8.0 を CentOS5.2 x64上でコンパイルし、make check したところ, tests/Examples/methods-Ex.Rの実行でエラーが生じます。 173行目あたりの setMehods("f","B0", function(x) c(x@b0^2, callNextMethod()))) 実行時にsprintf実行例外が発生してしまいます。これは他の版などで生じていないのでしょうか？　また、すでに知られたバグでしょうか？

エラーログ

> f <- function(x) class(x)
> 
> setMethod("f", "B0", function(x) c(x@b0^2, callNextMethod()))
 以下にエラー sprintf(gettext(fmt, domain = domain), ...) : 
   引数が少なすぎます 
Calls: setMethod -> message -> gettextf -> sprintf
 実行が停止されました

OSは gcc 4.1.2環境ですが、コンパイルはIntel C++ 11.0 +Intel ifort 11.0+MKL 10.0でおこなっています。
下記が、configure.siteの指定、および ./configure の引数です(インストールマニュアルにある　B.Ripley先生となかまさんの指定にしたがっていますが、config.siteのパラメータは -mpを-ieee-fpにしました）。

上記の命令は２回目以降はエラーを起こしません。上記の箇所の直前に try(setMethods(....))　）のように同様の命令を挿入すると、テストをパスします。
また,CentOS5.2のパッケージ compat-libstdc++33-3.2.3-61 から/usr/lib64/libstdc++.so.5をインストールしています。

--- config.site ----
CFLAFGS="-g -O2 -wd188 -ip -std=c99"
F77=ifort
FFLAGS="-ieee-fp -g -O3"
CXX=icpc
CXXFLAGS="-g -O3 -ieee-fp"
FC=ifort
FCFLAGS="-g -O3 -ieee-fp"
ICC_LIBS=/opt/intel/Compiler/11.0/074/lib/intel64
IFC_LIBS=/opt/intel/Compiler/11.0/074/lib/intel64
LDFLAGS="-L$ICC_LIBS -L$IFC_LIBS -L/usr/local/lib64"
SHLIB_CXXLD=icpc

----- configure 実行引数 -----
OMP_NUM_THREADS=8

MKL_LIB_PATH=/opt/intel/Compiler/11.0/074/mkl/lib/em64t

MKL="   -L${MKL_LIB_PATH}                               \
            -Wl,--start-group                               \
                    ${MKL_LIB_PATH}/libmkl_intel_lp64.a        \
                    ${MKL_LIB_PATH}/libmkl_intel_thread.a     \
                    ${MKL_LIB_PATH}/libmkl_core.a           \
            -Wl,--end-group                                 \
            -liomp5 -lguide -lpthread -lgomp"

./configure --with-lapack="$MKL" --with-blas="$MKL"

リストの終わり

↑

線形モデルでの変量効果（ランダム効果）・混合効果の指定†

Rcommander利用者 (2008-09-02 (火) 10:45:00)

線形モデルでの変量効果（ランダム効果）・混合効果の解析方法を教えてください。変数の指定方法などわかりません。よろしくおねがいします。

↑

JACKNIFE　という統計用語はどんなことですか？†

むかしちゃん (2008-06-25 (水) 22:20:15)

統計用語にJACKNIFE　という物騒な用語をインターネットで見つけましたが
どんなところでどんな風に使うのですか？

↑

組み込みデータセットを参照している関数の検索法？†

気になった人 (2008-05-20 (火) 09:32:44)

Rの魅力の一つは、各関数が比較的詳しい参考実行コードを持っており、関数の使い方がわかることです。もう一つ豊富な組み込みデータセットがあり、関数参考実行コードはしばしばそれを用いることで、単に参考のための参考に止まらない迫真性を与えることができます。

ところで逆ができるでしょうか。つまり、ある組み込みデータセット（例えば co2）を例示用コードで参照している関数を一覧することです。

例えば Linux なら端末から次のような命令を実行すればそれらしいことができますが、Rには元々そうした機能は無い？

$ grep -r "co2" /usr/lib/R/* | grep "/help"
↑

CIでのtick marks†

NP (2008-04-22 (火) 14:27:46)

Rを用いてCumulative incidence curveを書くところまではできたのですが、tick marksを打ち切りとなったものの観察期間を示すためにつけろといわれました。どのようにすればできるでしょうか。教えてください。よろしくお願いいたします。

↑

ACDモデル†

ＵＯ (2008-04-09 (水) 11:48:43)

皆さんに相談があります。
RでAutoregressive Conditional Durationモデルを使いたいのですが。
可能でしょうか？
ご存知の方がいらっしゃっいましたら、教えてください。

↑

rownameの付け方†

rowhelp (2008-04-07 (月) 10:57:23)

はじめまして、現在、Rを使用し遺伝子解析を始めたばかりです。
データフレーム内のある特定の列に含まれる文字列をrownameに変換したいと思います。

どのような作業を行うのが一番よろしいでしょうか。
ご教授をお願いいたします。

↑

凡例中のラインの角度の変更方法†

yuta (2008-04-03 (木) 16:43:46)

こんにちは。いつも参考にさせていただいています。
以下に示すような折れ線グラフと縦線プロット(plot()のtype="h")の混在したグラフがあります。このグラフに凡例を付けたいのですが、図にもあるように、legend()を普通に使うだけではラインはすべて横線になってしまいます。

(plotのスクリプト省略)
par(xpd=NA)
legend(locator(1),c("a-pro","a-height","b-pro",
  "b-height"),col=c("black","black","Royalblue","Royalblue"),
  lwd=c(2,3,1,1.5),pch=c(1,NA,1,NA))

凡例中で縦線を使う方法はあるのでしょうか。お尋ねいたします。

#ref(): File not found: "QA.png" at page "Ｑ＆Ａ (初級者コース)/13"

↑

特定のクラスターに属するデータだけをランダムに抽出するには？†

たいち (2008-02-12 (火) 18:20:13)

こんにちは。こんなデータがあるとします（実際はもっと長く，同じschoolに属するデータの個数は一定ではありません）。

   school y x
 1      1 4 1
 2      1 3 2
 3      1 2 3
 4      5 4 2
 5      5 5 3
 6      5 3 4
 7      8 4 3
 8      8 4 3
 9      8 1 2

ここから，同じschoolに属する人のxとyだけをランダムに抽出したいのですが，どのようにすればいいのでしょうか？
よろしくお願いします。

↑

ルート権限のないシステムにパッケージを個人的にインストールするには？†

使用歴だけは長い人 (2008-02-06 (水) 14:10:12)

ルート権限の無いシステムに個人的にパッケージをインストールしようとしています。ほとんどのパッケージは問題なく個人ディレクトリにインストールできるのですが、(これまた個人的にインストールした)外部ライブラリへのラッパーであるパッケージをインストールしようとすると、(当然でしょうが)　/usr/local/lib を参照しようとしてエラーになります。こうしたとき必要なライブラリ位置をRに教えてあげるにはどうすればよいのでしょうか。ちなみに問題のパッケージは gsl でこれは GNU GSL ライブラリーへのラッパーです。OS は　Linux, R 2.6.1 を使用しています。GSL ライブラリは個人ディレクトリにインストール済みです。管理者に GSL をインストールしてもらえば済むことでしょうが、同様の問題が今後も起きそうで、何とか個人的にインストール出きると便利なのですが。うまい方法があれば教えてください。パッケージのインストールにはいつも R CMD INSTALL xxx.tar.gz としています。

↑

パワーアナリシス†

大学生K (2008-01-24 (木) 04:45:22)

Rでパワーアナリシスをできる検定方法は、
t-test,anova,propのほかにはないのでしょうか？
個人的にノンパラメトリックを含め、いろいろな検定手法に対する、パワーアナリシスを行う必要がありまして・・・ご教授ください。

↑

Cコードで三角関数を使用すると０が返る件†

田島 (2008-01-18 (金) 22:06:14)

はじめまして。田島と申します。突然のお尋ねで恐縮ですがよろしくお願いいたします。
WindowsXpにてR-2.6.1＋MinGWにて下記のCコードをコンパイルしました。コンパイルに関するエラーは出ませんので、includeの問題はなさそうですが、seq[i]という入力ベクトルにどんな値が入っていても返り値であるsfというベクトルのすべての要素に０が返ります。ためしに、すべての返り値にsin(1.0)という固定値を代入するようにしても０が返るはずはないのですが、すべて０と返ります。私だけの環境で生じている可能性もあろうかと思いますので、初心者Q&Aがふさわしいかもしれませんが、なぜか書き込みができませんでしたので、こちらに投稿させていただいております。

その他の四則演算をするコードでは、期待した結果が返ります。

あまりコンパイルするような作業の経験はないのですが、なんかのパスが悪いのかと思いますが、お恥ずかしいことにかなり検索などもしましたがどこに問題があるか見当がつかない状況です。三角関数を使いたいだけということですので、四則演算で計算するアルゴリズムをつかって計算すればいいのですが、本来できることができていないのであれば、非常に気持ち悪いことですし、今後のために理解を深めたいということもありまして、お尋ねしているしだいです。
もし、ほかにお尋ねすべき適切なセッションがあれば、ご指導いただければ幸いです。どうぞよろしくお願いいたします。
test1.c

#include <math.h> 
void test1(double *seq, int *seqn, double *sf)
{
	int i;
	for(i=0;i<=*seqn-1;i++){
	sf[i] = sin(seq[i]);
	}
}

やったこと

Rcmd SHLIB test1.c

R上でやったこと

dyn.load("test1.dll")
a<-seq(0,pi,length=10)
.C("test1",as.double(a),length(a),sf=double(length(a)))$sf
↑

対応のある多重比較†

こはく (2008-01-12 (土) 22:14:03)

繰り返しのない2元配置（乱塊法または反復測定一元配置（対応のある一元配置？））で分析できるデータがあります。

これをボンフェローニの補正をおこなって多重比較するために、たとえば

a <- rep(1:3, each=4)
b <- rep(1:4, 3)
x <- c(7, 8, 5, 8, 1, 14, 16, 11, 7, 7, 6, 8)
pairwise.t.test(x,a,p.adjust.method="bonferroni")
pairwise.t.test(x,b,p.adjust.method="bonferroni")

ってかんがえたのですが、Ｒをはじめたばかりでよく分かりません。
ご意見いただけるとうれしいです。

↑

3次元グラフでのz軸方向の回転†

atte (2007-12-01 (土) 11:05:31)

scatterplot3dのグラフを回転させたいのですが、angleではxy方向しか回転できません。xz、yz方向に回転する方法はないのでしょうか。

↑

display関数†

似非R使い (2007-11-30 (金) 16:50:43)

初めて投稿させていただきます。
Rを用いた統計の教科書の中で、display()とあるのですが、その通り書き込んでもエラーが出ます。display()は削除されたか、別の関数に置き換わったのでしょうか？

↑

barplot(x,names=label,horiz=T)の長いラベルを横書きに†

SAStoR (2007-11-19 (月) 16:23:14)

初めて投稿させていただきます．

barplot()を使って棒グラフを作成したいのですが，長いラベル（全角１６文字程度）がついています．dotchartのようにラベルを横に表示させるにはどうすればよいでしょうか．棒グラフの場合，horiz=Tを指定してもラベルは垂直方向になってしまうため，全て表示がされず困っています．
”barplot(),text(),names,長いラベル”の組み合わせで検索してみましたが，うまく情報を見つけられませんでした．
どうぞよろしくお願いいたします．

A1 <- c(9,8,19,8,2,23,2,10)
A1label<-c("AAAAAAA","BBBBBBBBBBBBBBBBBBBBBBB","CCCCCCCCCCCCCC",
           "DDDDDDDDDDDD","EEEEEEEEEEEEEEEEEE","FFFFFFFFFFFFFFFFFF",
           "GGGGGGGGG","HHHH")
barplot(A1,names=A1label,horiz=T)
↑

Rで重回帰　再投稿†

Saito (2007-11-13 (火) 15:55:26)

皆様、的確なご指摘ありがとうございました。
それぞれのコメントにお答えする前に、まずマルチポストに関してお詫び申し上げなければなりません。その件に関しては非常に軽率であり、マナーに欠けた行為であったと反省しております。統計学関連なんでも掲示板に投稿したあとに自分の質問がその掲示板に不適切だと判断してこちらに投稿したのですが、もう少しインターネット上のルールを勉強してから投稿するべきでした。さらに言葉を重ねますと、統計学関連なんでもありの掲示板に書き込みをしたあと確認を怠り、こちらの掲示板で皆様の非難の声を聞いてから統計学なんでもありの掲示板にコメントがなされているのを知った次第です。皆様に不快な思いをさせてしまい申し訳ありませんでした。
　また、中級者用掲示板にこのような初心者の質問をしてしまったこともお詫びいたします。私の全くの勘違いなのですが、昨日の時点では初級者用掲示板がスパム防止のため封鎖した、というようなことが書いてあったので。しかし今はそんな文はどこにも記載されておらず、完全に私の勘違いでした。皆様にご迷惑をかけてしまい、申し訳ありませんでした。お詫びの言葉にもなっていないかもしれませんが、なにとぞご容赦ください。

＞簡単でいいですから疑問点が再現できるような例（もしくは組込みデータを使った例）をつけてください。
質問する側として当然の配慮が欠けていました。申し訳ありません。私が扱っているデータは以下のような感じです。すみません何度も試したのですが、うまく表示できませんでした。

> x
   StandardLength Feed Area
1             100   10    N
2             100   11    N
3             100   12    N
4             100   13    N
5             100   14    N
6             100   15    N
7             100   16    N
8             100   17    N
9             100   18    N
10            100   19    N
11            100   20    N
12            100   10    N
13            150   21    W
14            150   22    W
15            150   23    W
16            150   24    W
17            150   25    W
18            150   26    W
19            150   27    W
20            150   28    W
21            150   29    W
22            200   31    S
23            200   32    S
24            200   33    S
25            200   34    S
26            200   35    S
27            200   36    S
28            200   37    S
29            200   38    S
30            200   39    S
> result=lm(StandardLength~Feed+Area)　 #StandardLengthを独立変数、
                                         FeedとAreaを従属変数として重回帰
> summary(result)
Call:
lm(formula = StandardLength ~ Feed + Area)

Residuals:
       Min         1Q     Median         3Q        Max 
-4.422e-14 -1.798e-15  1.179e-15  4.160e-15  9.431e-15 

Coefficients:
             Estimate Std. Error   t value Pr(>|t|)    
(Intercept) 1.000e+02  9.435e-15 1.060e+16  < 2e-16 ***
Feed        2.209e-15  6.171e-16 3.579e+00  0.00139 ** 
AreaS       1.000e+02  1.332e-14 7.506e+15  < 2e-16 *** #問題2の該当箇所。
AreaW       5.000e+01  7.750e-15 6.452e+15  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 9.815e-15 on 26 degrees of freedom
Multiple R-Squared:     1,      Adjusted R-squared:     1 
F-statistic: 1.791e+32 on 3 and 26 DF,  p-value: < 2.2e-16 

> sd(Feed)　　#問題3の該当箇所。
               今となっては恥ずかしいかぎりです。
[1] 9.123986
> sd(AreaS)
 以下にエラーsd(AreaS) :  オブジェクト "AreaS" は存在しません
             #問題1の該当箇所。
              sd()が計算されないため、
              前に投稿したやり方では標準化できませんでした。

＞result() とは？ ＞summary()のことだと，思いますよ（^_^）
ご指摘のとおりです。Summary(result)の間違いでした。申し訳ありません。

＞問題３はあなたの理解通りですが、なぜ両者が同じだと思ったのですか(そもそも sd() で何を計算したのかそこがわからぬ)
　単純に標準誤差という言葉で表されるので同じ意味なのかな、と。
＞稲垣宣生著「数理統計学改訂版」定理12.5（正規誤差仮定の下での回帰係数の分布）を見てください。
　ご紹介ありがとうございます。また私の書いた問題3の中で、「その違いを述べている参考書は見当たりませんでした」とありますが、これは統計学的にではなく、Rの参考書の話です。誤解があったようで失礼致しました。

＞sd() で，従属変数の標準偏差を調べたのでは?推測に過ぎませんが。
推測のとおりです。説明不足、配慮不足でした。申し訳ありません。

＞問題1は，そもそも，変数の単位に明確な意味があるケースなのですから，標準偏回帰係数を算出する必要性が薄い気がします。
ご指摘のとおりです。カテゴリカルデータの傾きを知ったところで意味がないのは承知しているのですが、以前読んだ論文にカテゴリカルデータにも関わらず、GAM Coefficientというものを算出しているものがあったので…。現在本当にその意味（Coefficientを回帰係数と訳すこと）であっているのか読み直しているところです。
＞他のソフトウェアでは，おそらく，カテゴリカルデータを0, 1で表現して，連続変量と同様の算出法をしていると思いますので，ご確認して，ご報告頂けると幸いです。
＞問題2は，? contrasts を読むとわかるかと思います。
＞問題3は，すでにご紹介された教科書を参照なさると解決するかと思います。
丁寧なご回答ありがとうございました。自分の勉強不足を実感しました。少し勉強してみます。稚拙なコメントしか返せず申し訳ありません。

＞推測の必要はないので以下のように答えておきます。…
わかりやすい説明ありがとうございました。返す言葉もありません。もう一度勉強させてください。

この掲示板をマナーを守って利用している方々にご迷惑をおかけしたことを重ねてお詫び申し上げます。

↑

Ｒで重回帰†

Saito (2007-11-12 (月) 15:39:20)

初めて投稿させていただきます。

現在、統計ソフトRを使って卒業論文用のデータを解析しているのですが、重回帰をしているときに3つほど問題が出てきてしまいました。もしどなたか解決法をご存知でしたら、ご教授のほどよろしくお願いします。

問題1.カテゴリカル型データの標準化偏回帰係数が計算されない。

私は連続型とカテゴリカル型を同時に重回帰にかけているのですが、連続型のほうは
res <- lm(y~x1+x2+x3)
sdd <- c(0,sd(x1),sd(x2),sd(x3))
stb <- coef(res)*sdd/sd(y)
のような形でできたのですが、カテゴリカル型データはできません。しかし種々の参考書を見ると、カテゴリカル型データでも標準化偏回帰係数が出力されており、途方に暮れております。

問題2.カテゴリカル型データで一つの項目に複数の場合があるとき、うち一つが表示されない。

正確には表示されないことが問題なのではなくて、計算できないことが問題です。
例えば、
Yield=fertileのような関係を調べたいときに、カテゴリカル型データであるfertileの中にfertile1,fertile2,fertile3があったとします。これをそのまま重回帰にかけると、fertile2,fertile3の偏回帰係数やp値は表示されるのですが、fertile1は表示されません。Minitabの場合でしたら普通、fertile1の偏回帰係数はfertile1,fertile2,fertile3を足して、0になるように計算できる（らしい）のですが、Ｒはどうもそうではないような気がします。おそらく基準点のとり方が違うと思うのですが、どなたかＲでの基準点をご存知でしたら教えていただけないでしょうか。

問題3.重回帰を行った際に、result()で出力される結果に含まれるStd.Errorと、sd()で出力される結果が違う。

おそらく前者が傾きに対する標準誤差の計算結果で、後者がデータそのものの標準誤差ではないか、という予測はしているのですが私の調べた限り、その違いを述べている参考書は見当たりませんでした。データがカテゴリカル型であった場合には、sd()では出力されないのに、Std.Errorでは出力されているので、傾きに対する標準誤差なのかな、と。せめて傾きに対する標準誤差の計算方法がわかれば検算できるのですが、それもわからなかったため質問させていただきました。

以上の3つです。もしどなたか解決法をご存知でしたら、ご教授のほどよろしくお願いします。

↑

ヒストグラム†

星野 (2007-11-01 (木) 14:42:36)

ヒストグラムのソースなのですが、

# 2 群のヒストグラム
hist2 <- function(     x1,        # 第一群のデータ
                       x2,        # 第二群のデータ
                       brks=NULL, # 階級分割点
                       ...)       # barplot に引き渡す任意の引数
{
       if (is.null(brks)) {       # 階級分割点が与えられないときには，適切に設定
               brks <- hist(c(x1, x2), right=FALSE, plot=FALSE)$breaks
       }
       c1 <- hist(x1, breaks=brks, right=FALSE, plot=FALSE)$counts # 度数1
       c2 <- hist(x2, breaks=brks, right=FALSE, plot=FALSE)$counts # 度数2
       barplot(rbind(c1, c2), beside=TRUE, space=c(0, 0.2),        # 棒を並べて描く
               names.arg=brks[-length(c1)],                        # 横軸の目盛りラベル等
               axisnames=TRUE, axis.lty=1, ...)
}

というソースを用いているのですが、このソースだとx1、x2が数値でないといけないというエラーがでてきてしまいます。また、Ｒに関係する本でヒストグラムについてみても同じように数値の例ばかりです...。これだと、１０００個のデータをヒストグラムにする場合に入力が無理ではないでしょうか？データフレームからベクトルにして、ヒストグラムを描くような方法はないのでしょうか？よろしくお願いします。

↑

gstatパッケージでのクリギングについて†

青木.学生 (2007-10-31 (水) 14:17:21)

gstatでクリギングを行いたいのですが、予測グリッドデータの作成方法がわからなくて難航しています。
gstatマニュアルを読んで「RのGRASSパッケージにおいて国際横メルカトール・グリッドに変換した」ということはわかったのですが、具体的な方法がわからずにいます。
どなたかご存知の方がいらっしゃればご教示願えないでしょうか？
ちなみに現在The R BookについているR Ver1.9.0を使っています。

↑

ARFIMAXについて†

緑茶 (2007-10-16 (火) 15:41:17)

RでARFIMAXも推定することはできるのでしょうか？
出来るのであれば、よかったらやり方を教えてください。

↑

Rでの GPU 利用†

等々力渓谷 (2007-09-24 (月) 04:43:47)

　ＣＰＵよりも浮動小数点計算が高速で並列処理に向いていると言われる、最近のGPU（Graphics Processing Unit）ですが*1、Ｒでもこれを非グラフィックの数値計算に利用するプロジェクトはあるのでしょうか？

↑

power.prop.testの使い方（比較する2群のサンプル数が異なる場合）†

堀内 (2007-09-02 (日) 11:46:25)

実験計画を作成するにあたり、検出力80％以上となるようにサンプル数を設定したいと考えておりますが、その際比較する2群のサンプル数を同数ではなく、1：9もしくは2：8程度としたいと考えております。
他のソフト（SAS、Power and precisionなど）ではこのような設定での算定ができるようですが、power.prop.testに関するマニュアル等を見る限りではやり方が見当たりません。
対応方法をご存知の方がいらっしゃればご教示願えないでしょうか？
（特別な引数設定の仕方がある！他の関数でできる！アドインを入れてこのようにコーディングすればできる！など）

↑

mgcvライブラリのgam関数で平滑化スプラインが扱えるか？†

Blueblink (2007-08-22 (水) 01:18:17)

スプラインを用いたセミパラメトリック回帰をRで行いたいと考えております。まず，mgcvライブラリのgam関数を用いることを試みました。

こちらやこちらの文献では，gam関数のsオブジェクトを，平滑化スプラインを表すために用いています。例えば次のスクリプトで，変量Ozoneを3つの平滑化スプライン関数の和で説明する，加法モデルを求めることができるものとしています。

>library(mgcv); data(airquality)
>airq.gam<-gam(Ozone~s(Solar.R)+s(Wind)+s(Temp),data=airquality)
>summary(airq.gam)

ところが，help(gam)で見られるヘルプによると，mgcvのgamにおけるノンパラメトリック項（smooth term）には，罰則付きの回帰スプライン（Penalized Regression Spline）もしくは自由度を固定した回帰スプラインを用いると書いてあります。私の理解では，平滑化スプラインは共変量の数と同数の節点を用いますが，罰則付き回帰スプラインは計算負荷低減のため節点の数を抑えるもので，両者は異なります。

私は罰則付きの回帰スプラインの理論が良くわからないので，できればスプラインとして平滑化スプラインを用いたいのですが，本当にmgcvライブラリのgam関数で平滑化スプラインを用いることができるのでしょうか？

ご存知の方がいらっしゃいましたら，ご教授頂ければ幸いです。
どうぞよろしくお願いいたします。

↑

大きな行列の特異値分解・固有値分解†

有吉 (2007-07-11 (水) 17:15:20)

１万×３千〜４千ぐらいのサイズの行列（各要素は非負の実数）の特異値分解をしたいと思っています。
WindowsXPでR2.4を使って、
options("object.size"=160e+006)
memory.limit( size = 3072)
と設定してsvd()したのですが、メモリが足りなくなってできませんでした。
それで、転置行列とかけ合わせて３千〜４千×３千〜４千の対称行列にして固有値分解することにしたら、eigen()はできました。でも、そこでメモリ（ヒープ？）が足りなくなって、その後は行列演算ができなくなりました。
そこでいったんquit()してRを起動しなおして、行列演算を１つして、quit()して、起動しなおして、...の繰り返しになって作業が進みません。
上記のようなサイズの行列の特異値分解や固有値分解をして、さらに行列演算を続けていくにはどうしたらいいか教えていただけないでしょうか？
＃ 64bitOSで、64bitモードでmakeしたRを使えば、それで解決？

↑

qcc で X-Bar 管理図を作成するためのサブグルーピング†

岩田 (2007-07-08 (日) 18:33:35)

はじめまして。最近 R を使い始めた、岩田と申します。RjpWiki の皆様の努力のおかげで、R に関する情報が手軽に手にはいることを、感謝しています。

さて、qcc を使って X-Bar 管理図を作成したいと思っています。元となるデータは Excel で一つの列に納めて作成し、csv として書き出したものを、read.csv で R に読み込ませました。

X-Bar 管理図を作成するには、このデータをサブグループに分けなくてはなりませんが、qcc にそのような機能はありますか？　それとも、事前に複数列にまたがるテーブルに整形しなくてはならなかったのですか？　たとえばサブグループサイズが５の場合、R または Excel 上で簡単に複数列にまたがるテーブルに整形しなおす方法はありますか？　皆様のお知恵を拝借できればと思います。

↑

コンター図作成時のエラーについて†

Ovation (2007-06-12 (火) 11:11:51)

Rでwavelet変換をやっていますが、コンター図を描こうとしたら、

> filled.contour(Time,Frequency,t(z),nlevel=50,color.palette=topo.colors)
以下にエラーdiff(x) : 既定引数の再帰的な参照です

というエラーが出ました。
原因が分かる方、ご教授いただけると幸いです。

↑

NAN値を含むデータをcプログラムに渡すには†

kd (2007-06-12 (火) 10:17:56)

NAN値もcプログラムに渡す方法はあるでしょうか？
試しに下記のcプログラムを書いて実験を行いますと，.C() のところでエラーを出すようで，
NAN値を含むデータはcプログラムに渡らないようです．

#include <R.h>
void countnan(double *x, int *n, int *count)
 {
   int i;
   for (i=0,*count=0;i<*n;i++) 
     if (ISNAN(x[i])) (*count)++;
 }


> countnan <- function(x){ 
 n=length(x);
 .C("countnan",arg1=as.double(x),
  arg2=numeric(n),arg3=numeric(1))$arg3
}
> x<-c(1, NA , 3 , 4)
> countnan(x)
 以下にエラーcountnan(x) :  外部関数の呼び出し(引数 1) 中に NA/NaN/Inf があります
↑

R と Windows CSS†

パックマン (2007-06-07 (木) 16:28:32)

R で Windows CSSを利用している例はあるのでしょうか？

↑

dendrogramの解析†

とも (2007-06-06 (水) 16:49:44)

クラスタリング結果のデンドログラムを解析して、ノードごとの高さやリーフの情報を得る関数を書いています。
ノードやリーフの情報をベクトルとして取得する方法を教えてください。例えば、

("leaf","leaf","node","node")

のような形で返してほしいです。
dendrapplyという関数を参考にして書いたコードを示しました。
ノードやリーフの情報を再帰的に取得して表示することはできます。
dendrogramでRjpWikiやGoogleを検索してみましたが分かりませんでした。
引数の参照渡しの方法が分かるとできると思います。
dendrapplyはもとのデンドログラムのグラフ構造と同じ構造を返してくるので使えませんでした。

DNapplyTest <- function(d) {
	if (!is.leaf(d)) {
		for (j in seq_along(d)) Recall(d[[j]])
		cat("node:",attr(d,"height"),"?n")
	}else{
		cat("leaf:",labels(d),"?n")
	}
}

> DNapplyTest(as.dendrogram(hclust(dist(USArrests))))
leaf:Florida
leaf:North Carolina
node:38.52791
...
node:87.32634
node:168.6114
node:293.6228
↑

wavelets パッケージ（不等間隔 x_i vs. y_i のデータ系列）†

kd (2007-06-02 (土) 21:07:02)

wavelets パッケージを調べてみたのですが，「不等間隔な x_i とそれに対する y_i 」というデータ系列を与えて処理できるものが見当たらないようです．　このようなパッケージはどこかにあるでしょうか？

↑

glmmML 実行時の警告 (info = ) について†

たく (2007-04-25 (水) 11:33:54)

glmmML でモデルを作成した際、info = といった警告が出ることがあります。
モデルは作成されている用なのですが、このモデルは採用してもいいのでしょうか？
また、この際の AIC なども信頼できるのでしょうか?
どなたかご教授いただけますでしょうか?

↑

重回帰分析について†

たつや (2007-01-25 (木) 16:34:54)

重回帰分析して出た結果が英語で意味が読み取れませんでした。簡単な例を挙げて重回帰分析を行いましたので解説をお願いします。
質問はまず、この２つの違いについて。

F=lm(Z~X)
G=lm(Z~X-1)

つぎにsummaryで表示された内容がわかりません。できれば詳しい説明をお願いします。 自分で行う回帰分析、もしくはExcelの回帰分析は多少理解しているつもりなので、そっちと関連させて解説してもらえれば幸いです。

Z=c(1,2,3,4,5,6,7,8,9,10)
X=c(1,3,4,2,5,6,7,8,10,10)
F=lm(Z~X)
G=lm(Z~X-1)
summary(F)
Call:
lm(formula = Z ~ X)
Residuals:
Min     1Q       Median  3Q      Max
1.11283 -0.47400 0.09181 0.27600 1.80531
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.35841  0.59297    0.604   0.562
          X 0.91814  0.09329    9.842   9.56e-06 ***
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.887 on 8 degrees of freedom
Multiple R-Squared: 0.9237, Adjusted R-squared: 0.9142
F-statistic: 96.86 on 1 and 8 DF, p-value: 9.561e-06
summary(G)
Call:
lm(formula = Z ~ X - 1)
Residuals:
Min    1Q      Median 3Q     Max
0.9035 -0.5006 0.1770 0.2494 2.0644
Coefficients:
  Estimate Std. Error t value Pr(>|t|)
X 0.96782  0.04255    22.75   2.91e-09 ***
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8552 on 9 degrees of freedom
Multiple R-Squared: 0.9829, Adjusted R-squared: 0.981
F-statistic: 517.5 on 1 and 9 DF, p-value: 2.909e-09
  • Excelで実行して，内容を比較してみたら？ -- 2007-01-25 (木) 16:48:51
  • 何がわからないのか？というか何が知りたいのか？ -- 2007-01-25 (木) 17:34:36
  • 何でも聞けばよいというものではない。 -- 2007-01-25 (木) 17:35:28
  • ＞ 結果が英語で意味が読み取れません --- なのでしょうね。 -- 2007-01-25 (木) 17:38:19
  • 「重回帰分析」でキーワード検索、既に幾つかあるR本を見る、それでもわからなければ、改めて聞く。 -- 2007-01-25 (木) 17:42:33
  • Estimate Std. Error t value Pr(>|t|) の意味を教えてほしいです。 -- 2007-01-25 (木) 18:10:30
  • あと、F=lm(Z~X) とG=lm(Z~X-1)の違いが意味がわからないからできない…。 -- 2007-01-25 (木) 18:12:25
  • すみません。お願いします。 -- 2007-01-25 (木) 18:13:38
  • いっこめは係数の推定値、標準誤差、t統計量、p値。二個目は切片項を含むか含まないか。excelの結果と比べれば判ると思うんだけど。 -- 2007-01-25 (木) 18:19:16
  • １つ目の方は比較すればわかりました。すみません。 -- たつや 2007-01-25 (木) 18:29:11
  • 返答ありがとうございました。だいたい理解できました。あとは自分で考えます。 -- たつや 2007-01-25 (木) 18:32:45
↑

回帰係数の当てはめ†

佐藤 (2007-01-15 (月) 16:57:21)

> attach(learning)
> learning.glm<-glm(def~lastprice+volume+size+rtn+next1rtn+lsize+epr1
                    +greps1+dy+bpr+wcr+rrl+rrs+cfpr+roe+spr+ebitda+igr
                    +sgrowth+esr1+lev+cgrowth+vol+roa1,family=binomial)
> summary(learning.glm)
> detach(learning)
> test<-read.csv("property200006.csv",header=TRUE)
> attach(test)

ある（learning）データから回帰分析して得られた回帰係数を用いて、違う（test）データを回帰分析しようと考えています。
現状は回帰係数を自分で入力しているのですが、自分で入力することなく回帰分析することはできないでしょうか。
よろしくお願いします。

↑

重回帰分析のステップワイズ変数選択でのエラー†

GEO (2007-01-06 (土) 18:34:20)

はじめまして。
結果率（果実/花）を従属変数、各花形質を独立変数としてロジスティック重回帰を行おうとしております。この際に、関数step()を用いて、変数選択しようと試みたのですが、以下のようなエラーが出てうまくいきません。どのような原因でこのようなエラーが出るのか、また、このまま解析を進めても良いのか、教えていただけないでしょうか？

> h0 <- read.delim("c:/2000R2.txt")
> fs <- h0$fru/h0$flo
> result <- glm(fs~h0$ftl+h0$ptl+h0$ptw+h0$sta+h0$pis+h0$blt+h0$wit, weights=h0$flo)
> result2 <- step(result)
Start:  AIC= -2819.3 
 fs ~ h0$ftl + h0$ptl + h0$ptw + h0$sta + h0$pis + h0$blt + h0$wit 

         Df Deviance      AIC
- h0$wit  1    22.05 -2821.19
- h0$ptw  1    22.05 -2821.16
- h0$sta  1    22.06 -2821.13
- h0$blt  1    22.12 -2820.58
- h0$ptl  1    22.14 -2820.44
<none>         22.04 -2819.30
- h0$pis  1    22.33 -2818.70
- h0$ftl  1    23.20 -2811.30
Error in step(result) : number of rows in use has changed: remove missing values?

また、remove missing values? という問いかけが出てきますが、これに対しての答え方が良くわかりません。。どのようにすれば、欠損値を除いて、解析を先に進めることができるのでしょうか？
よろしくお願いいたします。

↑

Rでの単位根検定†

秋田県民Ａ (2006-12-17 (日) 17:43:33)

最近時系列分析を行うためにRを使い始めた初心者です。
まず定常性を確認するためにadf.test()関数を用いているのですが、
?定数項あり、トレンドあり　?定数項あり、トレンドなし
?定数項なし、トレンドあり　?定数項なし、トレンドなし
などの各状況に応じた使い方が分かりません。そのままadf.test()を使った場合にはどの場合の検定になるのでしょうか？
どなたかお分かりになる方がいらっしゃいましたら、教えてください。

↑

intel mac でのRのインストール、起動†

やまだ (2006-11-21 (火) 11:14:18)

初めて書き込ませていただきます。どうぞよろしくお願いします。
今までPowerBookG4でMac版を動かして使っていたのですが、今回先月発売になったCore2DuoのMacBook Proを新しく購入したのでインストールしてみたところ、起動時に以下のコメントが赤字で表示されるようになりました。

「2006-11-21 10:48:31.215 R[261] CFLog (21): dyld returns 2 when trying to load /Users/（ユーザ名）
/Library/ScriptingAdditions/YouHelper.osax/Contents/MacOS/YouHelper」

このときインストールしていたのがバージョン2.3.1でして、CRANを調べたらユニバーサル版の最新版として2.4がありましたので、それを再度インストールしてみましたが状況は同様で毎回起動時に表示されます。一応試しに使ってみましたが、問題なく動いてはいるようなんですが。。。ちょっと気持ち悪いので。
もしこの状況に見覚えのある方等おられましたらコメントいただけるとありがたいです。 どうぞよろしくお願いします。

↑

wavelet†

lostway (2006-11-12 (日) 01:58:50)

Rでwavelet解析やっている方いますでしょうか？

↑

cgiからRの起動†

イシ (2006-11-10 (金) 15:33:30)

cgiスクリプトからRを起動し、グラフを作成したいのですが、Rのjpegデバイスを開く段階でエラーが出てしまいます。対処法ありますか？
R-Version 2.3.1

#### エラーメッセージ #####

Error in X11(paste("jpeg::", quality, ":", filename, sep = ""), width,  : 
        unable to start device JPEG
In addition: Warning message:
unable to open connection to X11 display ''
Execution halted
↑

glmの説明変数に割合は使えるのか？†

服部 (2006-11-02 (木) 14:06:43)

ある目的変数(比率）を説明する変数を特定するときに、glmを用いて解析をしたいと思っています。
その際説明変数に割合の変数含めたいのですが、どのようにしたらよいのでしょう？
現在以下のようにglm関数を使っています、訂正すべき場所があったらそこも指摘していただけると幸いです。

> a<-adult/colonysize
> yy<-glm(cbind(x,y)~a+b+c+d,binomial)
↑

R2.4.0からの，無意味な正規表現におけるWarning†

アール (2006-10-27 (金) 16:26:24)

R 2.4.0 からは "?.foo" のような無意味なエスケープシークエンスに warning を出すようになりました。それはそれでいいのですが，これを考えていてちょっとした疑問が出てきました。
"foo.bar.baz" のような正規表現で，最初のピリオドは文字通りピリオド（他の言語の正規表現では"?."）二番目のピリオドは何でもよい一文字を表す正規表現"."なんかは，書きようがないんじゃないかと。fixed 引数も，全てに適用されるものなので無力だしね。

> sub("foo?.bar.baz", "match", c("foo.bar.baz", "foo.bar,baz", "foo,bar.baz", "foo,bar,baz"))
[1] "match" "match" "match" "match"
Warning messages:
1: '?.' is an unrecognized escape in a character string 
2: '?.' is an unrecognized escape in a character string 
> sub("foo.bar.baz", "match", c("foo.bar.baz", "foo.bar,baz", "foo,bar.baz", "foo,bar,baz"))
[1] "match" "match" "match" "match"
↑

閉曲線で囲まれた領域を塗りつぶす関数†

青木繁伸 (2006-10-19 (木) 15:41:02)

お絵かきソフトでよくある「ペイント缶」の機能を果たす関数がないのかなあと。
曲線にこだわるわけではない（曲線だって，グラフに描くときには折れ線）ので，polygon関数でも良いわけだが，点を与えてポリゴンを描いて内側を塗りつぶすのではなくて，経過は問わず既にできている閉局面で囲まれている領域を塗りつぶす関数が欲しいということ。

#ref(): File not found: "ex.png" at page "Ｑ＆Ａ (初級者コース)/13"

↑

メモリについて†

しも (2006-10-02 (月) 21:40:16)

マニュアル等で調べても、解決法が見つかりませんでしたので質問します。
関数を作成し、計算をさせています。計算途中で、以下の警告が出て計算がストップします。

以下にエラーodbcQuery(channel, query, rows_at_time) : 
Calloc がメモリー (263168 of 1) を割りあてられませんでした
追加情報: Warning message:
Reached total allocation of 1015Mb: see help(memory.size)

使用するメモリーサイズを"--max-mem-size=2G"を使用して拡大しても、計算が止まる状況です。
関数内では「rm()」を使用し、こまめにオブジェクトは削除しています。
Windows タスク マネージャのPF使用量を見るとrm()を使用していても使用量がどんどん積み重なり、ある一定値で計算が止まっり、計算が止まっても使用量は減りません。使用量を減らす唯一の方法はq()でRを閉じる方法しかありませんでした。

メモリーを上手に割り当てるようにするにはどうしたらよろしいでしょうか？

↑

optim()でエラーが出ます†

SS (2006-09-28 (木) 14:53:42)

2パラメータの最尤推定を行う際、optim()を使って尤度関数を最大化する値を探そうとしたら、「以下にエラーchol(M) : 次数 1 の主対角行列が正定値ではありません」というエラーメッセージを得ました。nlm()でも同様のエラーメッセージを得ました。

optimの場合、optionとして"Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN"の全てを試しましたが、同様のエラーメッセージが出ました。

2パラメータなので、尤度関数をプロットしてみた感じだと、山頂はある様子でした。初期値も、プロットした尤度関数の山頂付近に設定しました。

このエラーメッセージは、どういう状況のときに出力されるのか、ご存知の方がいらっしゃったら、教えていただけるとありがたいです。

↑

legendにグラデーション†

nonami (2006-09-14 (木) 19:45:14)

散布図を作っています。
x, yのデータ以外に第3のデータがありまして、それが0〜100まであります。
rainbowをの一部を使って、0なら青、100なら赤その間は緑から黄色、オレンジというようになんとなくサーモグラフィっぽい感じで色をつけています。
plot自体はうまく行ったと思うのですが、legendを付けようとしたときに
はたと困ってしまいました。青→緑→黄色→オレンジ→赤の順に数値が高くなるよというlegendは作成できるのでしょうか?

↑

安定分布の擬似乱数†

TH (2006-08-13 (日) 18:39:13)

Rで安定分布に従う擬似乱数を発生するパッケージ等をご存知の方はご教示ください。

↑

tune.svm を実行すると .Random.seed のエラーがでて実行できない†

はるか (2011-11-11 (金) 11:56:20)

R初心者です。先日は疑問に対して丁寧にお答えいただいてありがとうございました。

その後、tune.svmを使おうと先日からお世話になっているサイトで丁寧に解説があったので、下記のように同様に実行してみたところ

gammaRange <- 10^(-5:5)
costRange <- 10^(-2:2)
t <- tune.svm(Species ~ ., data = iris, gamma=gammaRange,
cost=costRange,tunecontrol = tune.control(sampling="cross", cross=8))
以下にエラー sample(n) : 
  .Random.seed は整数ベクトルではなくタイプ 'list' で

と出てしまい実行することができませんでした。.Random.seedの設定？など特にしていないのですが、何か特別な設定等必要なのでしょか？webを検索したのですが、同じような点で悩んでいる例がなかったので、質問させていただきました。
環境は、OS WindowsXp, R version 2.12.1です。お手数ですが、アドバイスいただけると大変助かります。よろしくお願い致します。

↑

作業ディレクトリを、スクリプトが置かれてるディレクトリに変更†

mat (2011-11-10 (木) 13:03:16)

あるスクリプトを開いてそれを実行する際に、作業ディレクトリを(スクリプト内で動的に)そのスクリプトファイルが置かれている場所に変更することは出来ないのでしょうか？

↑

図の上での計測について†

KT9 (2011-11-09 (水) 09:45:58)

医学系の研究者です。
JPEG形式の超音波画像があります。これをRで取り込んで、ある一点の(x,y)座標、2点間の距離、面積などを計測する方法はあるのでしょうか。
　 (x,y)さえ打ち出せれば、後は自分で別個に組んだプログラムで、距離、面積、角度などは簡単に測れると考えています。
　以前他の研究者がMathLabでは、JPEGで取り込んだ画像上で、様々な処理を行っているのを見ました。
　宜しくお願い致します。

↑

CSVファイルからのデータのインポートができない†

diereinevernunft (2011-11-03 (木) 00:17:54)

R初心者です。
Rを2.13.0から最新バージョン(2.14.0)に更新し、Rコマンダーも新たに1.7-0をインストールし直したのですが、データのインポートでCSVファイルからデータを読み込もうとすると、

Dataset <- 
read.table("C:/Users/Documents/ワークショップ/Rデータ処理用/（ここで強制改行）
修論素材/加工後/R用/111024/tsushin_growth.csv",
header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE)

が実行されたあと、メッセージ欄に

[2] エラー: NA

と表示され、データがインポートできません。
ちなみに、CSVファイルをtxtファイルに変えるとインポートできるようなので、現在はそれで対応しています。
googleで検索してみたのですが、現状では有効な対策を見つけることができませんでした。

OS Windows7(64bit)
R version 2.14.0
Rcmdr version 1.7-0
↑

svmの結果, Total Accurancyがでない†

はるか (2011-11-02 (水) 15:38:48)

R初心者です。svmに興味があって、googleを検索したところ、丁寧な解説があったのでこのWebに沿って、同じようにやると

Total Accuracy:94
Single Accuracies:
90 96 96

と同じように結果が再現され正答率がでたところまではよかったのですが、いざ自分の持っているデータでやって結果を見ようとすると...

model <- svm(judge~., data=dataset, corss=3)
summary(model)

Parameters:
   SVM-Type:  eps-regression 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.01190476 
    epsilon:  0.1 

Number of Support Vectors:  51

3-fold cross-validation on training data:

Total Mean Squared Error: 0.6049732 
Squared Correlation Coefficient: 0.06871624 
Mean Squared Errors:
 0.4767778 0.5592772 0.7788645

と出てしまい、Accurancyが出てきません。
何がいけないのか分からず、困っています。
OS WindowsXp, R version 2.12.1でやっています。

扱ったデータから解を求められないということなのでしょうか？
それとも何か別に問題点があるのでしょうか。
周りにアドバイスを求める人がおらず困っています。
お手数ですが、アドバイスいただけると大変助かります。よろしくお願い致します。

↑

導関数の求め方について†

みのむし (2011-10-30 (日) 16:58:40)

R初心者です。導関数を求めるのに、deriv()を使っています。直接、カッコ内に関数を打ち込まないといけないようなのですが、そのほかの方法はないのでしょうか？

たとえば、

> Df <- deriv(~x^2, "x", func=TRUE)
> Df(2)
[1] 4
attr(,"gradient")
     x
[1,] 4

と答えがえられるのですが、カッコ内に打ち込む関数x^2を、最初にfと定義し、それを代入するとうまくいきません。

> f <- expression(x^2)
> Df <- deriv(f, "x", func=TRUE)
> Df(2)
expression(x^2)
attr(,"gradient")
     x
[1,] 0

ご存じの方がいらしゃれば、教えていただけないでしょうか？お忙しいところ申し訳ありません。よろしくお願いいたします。

↑

plotの軸目盛りの指定について†

ringori (2011-10-27 (木) 01:48:51)

お忙しいところ失礼します。plotを使って折れ線グラフを作成しているのですが、x軸の目盛りの内容の一部を文字列で指定したく思っています。少し調べた所

test <- c(0.8699, 0.9029, 0.9133, 0.9014, 0.9160)
plot(test, type="l", lwd=2, xlab="test", ylab="AUC", xaxt="n")
mtext(1:4, 1, 1, at=1:4)
mtext("Q", 1, 1, at=5)

として近いものは出来たのですが、目盛りの線が無くなってしまうのと、もっとスマートな方法があればご教授いただければと思います。 もし、過去の投稿等見逃していましたら申し訳ありません。環境は

sessionInfo()
R version 2.10.1 (2009-12-14) 
i386-pc-mingw32

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

↑

grepでの検索について†

Toy (2011-10-26 (水) 14:12:37)

お世話になります。
Windows 7 で R version 2.13.0 を使用しています。
様々な記入形態がある文字列について抽出処理を行おうとしています。
下記のような文字ベクトルについて

x <- c("マツホ−ム", "タケハウス", "ウメコウムテン",
		"カブ）マツホ−ム", "タケハウスカブシキガイシャ", "ウメコウムテンカブ")

次のような例では、正しく結果が表示されます。

> grep("_*ウメコウムテン_*", x, value = TRUE)
[1] "ウメコウムテン"     "ウメコウムテンカブ"

ところが、次の例では

> grep("_*タケハウス_*", x, value = TRUE)
character(0)

と全く結果が得られず、次の例では、

> grep("_*マツホ−ム_*", x, value = TRUE)
[1] "マツホ−ム"         "ウメコウムテン"     "カブ）マツホ−ム"   "ウメコウムテンカブ"

と余計なものまで出力されてしまいます。
どうしてこのような結果になるのか、お教えいただければ幸いです。

↑

CSVを読み込んだあとで、ある特定の文字列の行数まで読み飛ばす方法†

nya (2011-10-23 (日) 01:21:25)

read.csv でカンマ区切りのCSVファイルを読み込んだ後、ある特定の文字列から始まる（たとえば"hoge test"）行数まで読み飛ばす方法はありませんか？
この行数はファイルによって異なるため、read.csvのskip=xxxでは対応できず、困っております。
ご教授いただければ幸いです。

filetest.csv
↑

複数のブロック行列†

複数のブロック行列 (2011-10-18 (火) 22:35:10)

A = (3,7,4,6,...,)のような列ベクトルがあり、I[N*N]を単位行列としたとき、B = (I[3*3], I[7*7], I[4*4], I[6*6], ...)のような、Aの要素を次元とするブロック対角行列Bを作成したいのですが、方法が分からず苦慮しております。
ご教授願えませんでしょうか。

↑

ｔ検定に関して†

tau (2011-10-17 (月) 22:09:13)

R 2.13.1を使用しております。

下記のようにt.testを行いました。

> group1
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
[1,]   15   12   12   14   14   12   16   11   19    18    15    11    10     9     7    15
[2,]   14   16   12    9   14   16   14   11   11     6     7    11    16    12     9     9
[3,]   16   10   11   16   20   20   16   18   17    16    12     8    14    12    17    10
> group2
    [,1] [,2] [,3] [,4]
[1,]   18   19   12   16
> group3
    [,1] [,2] [,3] [,4]
[1,]   14   17   17   15
> group4
    [,1] [,2] [,3] [,4]
[1,]   13   19   11   17
> t.test(group1,group2,var.equal=F)
       Welch Two Sample t-test
data:  group1 and group2 
t = -1.9197, df = 3.667, p-value = 0.1337
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
-7.810966  1.560966 
sample estimates:
mean of x mean of y 
  13.125    16.250 
> t.test(group1,group3,var.equal=F)
       Welch Two Sample t-test
data:  group1 and group3 
t = -2.9049, df = 6.241, p-value = 0.02597
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:~
-4.8156168 -0.4343832 
sample estimates:
mean of x mean of y 
  13.125    15.750 
> t.test(group1,group4,var.equal=F)
       Welch Two Sample t-test
data:  group1 and group4 
t = -0.9899, df = 3.474, p-value = 0.3861
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:~
-7.463495  3.713495 
sample estimates:
mean of x mean of y 
  13.125    15.000


group1とgroup2,3,4の平均値を比較していますが、group2とgroup4の平均値の中間の値を取るgroup3でのみP値が小さくなるのはなぜでしょうか？
ご教示頂ければ幸いです。

↑

複数回答の処理†

pp (2011-10-13 (木) 15:34:41)

こんにちは。
複数回答のアンケートを作って読み込んだところ、下記のように表示されます。（データフレーム）

> data

1     にんじん,たまねぎ
2     にんじん
3     にんじん,たまねぎ,じゃがいも

これをtable(data)とすると上記がそのまま出てしまいますが、

にんじん　3
たまねぎ   2
じゃがいも 1

とするような方法はないでしょうか。（文字をカンマで分割したい）

初心者で恐縮ですが、ご回答いただけますと幸いです。

↑

入力された変数名を文字列に変換†

Minako (2011-10-12 (水) 08:31:41)

はじめまして。R初心者（31歳前田敦子似）です。
データフレームのhogeという変数名を入力すると、その変数のサマリーや検定結果などをresult.txtという名前のファイルに保存する関数を作ろうと思ってます。計算結果の前に、どの変数について調べたか記載したいので、
あるデータフレームに年齢が入力されたageという変数名があるときに、ageage(x)という関数

ageage(x)<-function(x) {
m<-sprintf("title\t%s\n",x)
cat(m,file="result.txt",append=T)
}

を作って、result.txtに、title age　と表示させたいです。
ところがこれでageage(age)を実行すると

title  50
title  54
title  76

などとageが年齢を含んだベクトルとして理解されて、全部記載されてしまいます。関数のxに入力された変数名を文字列として取得するためにはどのようにすればよいでしょうか？？

↑

１秒以下の扱いについて†

tx (2011-10-11 (火) 11:00:45)

プレートリーダーが出力する以下のようなフォーマットのエクセルファイルから、Rへのデータの読み込みを試みておりますが、

Plate	Repeat	Well	Type	Time	Luciferase (CPS)	Time	 Luciferase (CPS)
1	1	A01	M	00:00:06.80	53530	00:02:38.05	47900
1	1	A02	M	00:00:07.07	107230	00:02:38.31	94150
1	1	A03	M	00:00:07.34	103510	00:02:38.58	95070


取り込むと、一秒以下の部分が削られてしまいました。
読み込みにはこちらのサイトに書かれていた方法を用いました。

> library(RODBC)
> sheet <- odbcConnect("Excel Files")
> tab <- sqlQuery(sheet, "select * from [List ; Plates 1 - 1$]")
> t <- as.numeric(difftime(tab[,7], tab[,5], units="secs"))
> t[1]
[1] 152


コンマ秒の部分を扱う方法を教えていただけると嬉しいのですが、
どうぞよろしくお願い致します。

↑

時間差の計算について†

tx (2011-10-11 (火) 07:18:37)

POSIXct形式の時刻の差を使って乗除計算をしたいのですが、defftimeの値を使っての割り算をRは認めてくれません。"difftime"オブジェクトを数値の配列に変換する方法、または他の方法での時間差の計算法がありましたら御教授いただけませんでしょうか。どうぞよろしくお願い致します。

↑

GIF形式での保存†

tetora (2011-10-10 (月) 00:22:25)

既出でしたらすみません。Rの出力をGIF形式で保存したいのですが、デフォルトでGIF出力する関数が見当たりませんでした。こちらのページを参照すると、パッケージでもGIF形式での保存に対応したものはないようなのですが、どうしてもGIF形式で保存したい場合は、一度別の形式(jpeg,pngなど)で保存した後、外部のプログラムで形式変換するしか方法はないのでしょうか？

↑

image関数の分解能†

まつだ (2011-10-07 (金) 15:32:59)

image関数を使って、1または0のバイナリデータの表示を試みているのですが、データ量が多い時にうまく表示できずに困っています。

例えば、以下のようなデータを作ります。

x <- c(rep(0, 30000), 1, rep(0, 19999))

これをimage関数を使って表記するために以下を実行します。

m <- matrix(x, nrow=50000, ncol=2)
par(mfrow=c(2,1)) # plot関数との表示比較のため
image(m, col=c("white", "black"))
plot(1:nrow(m), m[,1], type="l")

plot関数では、50000データの中に1つだけある1のデータを表記できているようですが、image関数ではそれが消えてしまいます。

このようなときに、皆さんはどのように対処なされていますでしょうか。
ご教授いただければ幸いです。

ちなみに、以下くらいのデータであれば問題はないようです。

x <- c(rep(0, 300), 1, rep(0, 199))
m <- matrix(x, nrow=500, ncol=2)
par(mfrow=c(2, 1))
image(m, col=c("white", "black"))
plot(1:nrow(m), m[,1], type="l")

どうぞよろしくお願いします

↑

2要素ずつのマッチング (forループを使わないで）†

shannon (2011-10-05 (水) 18:38:21)

質問のタイトルが不適当かもしれませんが、よろしくお願いします。

a <- c(0, 0, 1, 1, 0, 0, 0, 1, 0, 0)

上のような"0"と"1"が並んでいるベクトル(a)に対して、

a2 <- c(0, 1, 1, 0, 1, 0)

"0"が2つ以上ならんでいる部分を削除して、"1"の前後にそれぞれ"0"が１つだけ入るように（a2）したいのですが、forループを使って

a2 <- 0
for (i in 1:9) {
   if (a[i]==0 && a[i+1]==0) next
   a2 <- c(a2, a[i+1])
}

のようにすると、ベクトルの要素数が膨大になった場合にかなり時間がかかってしまいます。
そこで0が２つ並んでいる要素番号（上記のaでいえば、要素番号1,5,6,9)に、FALSEを当てはめたような、論理ベクトル(a3)を作成して、下のように計算したいのですが、方法が分からず悩んでいます。forループを使わずに、a3のようなベクトルを作るにはどうしたらよいのでしょうか。もしくは別の方法でa2を作りだす良い方法はないでしょうか。分かりにくい質問で申し訳ないですが、ご教授お願いします。

> a3 <- c(FALSE,TRUE,TRUE,TRUE,FALSE,FALSE,TRUE,TRUE,FALSE,TRUE)
> a2 <- a[a3]
> a2
[1] 0 1 1 0 1 0
↑

軸の入れ替えについて†

afromonkey (2011-10-05 (水) 16:26:51)

統計ソフトRで、横軸をy軸、縦軸をx軸にして線グラフやヒストグラムを描くには、どうすればよいのでしょうか。
barplot関数の場合は引数に horiz = T を与えればよいということが分かったのですが、 plot や hist などの他の関数についてはどのように設定すればよいでしょうか。
よろしくお願いします。

↑

R ver 2.13.1以降　Rコマンダーで不正なマルチバイト文字エラー†

３６５ (2011-10-04 (火) 22:28:06)

Rの2.13.1以降でRコマンダーで２バイト文字のファイル名などを使用すると不正なマルチバイト文字があるというエラーが出ます。
Rconsoleでは当該ファイル名は認識され、データもロードできますので、Rコマンダーの方の問題かと思うのですが原因がよくわかりません。過去に同様な投稿も見つけられませんでした。
また、R2.13.0以前では同様なエラーは出ませんでした。
OSはwindowsXP sp3、RcmdrはそれぞれのRバージョンで更新しています。
バージョンアップが出来ずに困っています。どなたかお心当たりのある方、ご教授下さい。

↑

パッケージの中身を書き換える†

西田 (2011-10-01 (土) 15:48:10)

あるパッケージに含まれる動作を拡張して使いたいと思っています。
今まではlibrary()で読み込んだ後、関数等を書き換えていたのですが、今回拡張したい部分がどうやらC言語で書かれているようで（.Call）書き換え方がわかりません。
パッケージの中身を直接いじることができれば一度で済みますし、余計な不具合の原因にもならないと思い質問させて頂きました。どうかご教授ください。

↑

類似度によって欠損値補完†

tak (2011-09-30 (金) 16:41:40)

下のような、縦方向にケース番号、横方向にデータ番号を配置した欠損値ありの行列(1)があるとします。(実際はケース数、データ数の種類がもっと多いものを使います。)

> data
       data1 data2 data3 data4  data5
case1     NA    10   100  1000  10000
case2      2    NA   200    NA  20000
case3      3    30   300  3000  30000
case4     NA    40    NA  4000     NA
case5      5    NA   500  5000  50000
case6      6    60   600    NA     NA
case7      7    70   700  7000  70000
case8     NA    80   800  8000     NA
case9      9    90    NA  9000  90000
case10    10    NA  1000 10000 100000

上の行列に対して、ケース毎の類似度を計算した行列(2)が下のようになります。(例えば1行2列は、case1のcase2との類似度です。『case1の、case2との類似度』と『case2の、case1との類似度』は同様に計算しているので対称行列になっています。)

> similarity
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,] 0.00 0.16 0.46 0.50 0.77 0.84 1.38 1.41 1.61  1.73
 [2,] 0.16 0.00 0.20 0.00 0.60 0.67 1.00 0.67 1.17  1.61
 [3,] 0.46 0.20 0.00 0.17 0.46 0.63 1.05 1.00 1.42  1.61
 [4,] 0.50 0.00 0.17 0.00 0.11 0.25 0.50 0.67 0.84  0.67
 [5,] 0.77 0.60 0.46 0.11 0.00 0.17 0.46 0.47 0.80  1.15
 [6,] 0.84 0.67 0.63 0.25 0.17 0.00 0.21 0.33 0.53  0.67
 [7,] 1.38 1.00 1.05 0.50 0.46 0.21 0.00 0.20 0.47  0.69
 [8,] 1.41 0.67 1.00 0.67 0.47 0.33 0.20 0.00 0.17  0.31
 [9,] 1.61 1.17 1.42 0.84 0.80 0.53 0.47 0.17 0.00  0.20
[10,] 1.73 1.61 1.61 0.67 1.15 0.67 0.69 0.31 0.20  0.00

そこで、行列(1)に存在する欠損値を類似度によって補完する処理を作成したいのですが、どのようにすればよいのか思いつかず困っています。

処理の内容は、欠損値補完対象のケース番号と類似度の高い上位k個のケース(但し、kは任意で、補完対象のデータと同じ番号のデータが欠損値であるケースは除く)の、同じ番号のデータの平均値を入れるようにします。

例えば(k=3の場合)、case1のdata1が欠損値になっていますが、case1と類似度の高い上位3ケースは、case10, case9, case8です。しかし、case8のdata1は欠損値なので、次に類似度の高いcase7を使って、case10, case9, case7のそれぞれのdata1である10,9,7の平均をとって、8.66を補完値とします。この処理を全ての欠損値に対して行います。

説明がわかりづらくて申し訳ありませんが、宜しくお願いします。

↑

指定行へのジャンプ†

ohno (2011-09-30 (金) 15:37:38)

C言語などで使うgoto文のように、指定行へジャンプしたい時に使う関数はあるのでしょうか？
下の行へジャンプする場合は、以下のようにジャンプ先の命令文の前までをwhileなどの条件式で囲って途中でbreakすれば、一応はジャンプできますが、
この方法では上の行へジャンプする時が複雑になるので、なにか良い方法がないかと考えています。

while (1) {
    ...
    if (1) break
    ...
}
...　#ここへジャンプ

良い方法がありましたらご教授お願いいたします。

↑

文字列の何番目がマッチするかを調べる†

do-san (2011-09-29 (木) 15:30:24)

例えば、文字列"abcde"があったとして、文字パターン"d"にマッチするのが何番目かを調べるにはどうすればよいのでしょうか？（答えが4になる。）
文字列"abcde"を一文字ずつ分解・ベクトル化して、charmatchで調べる、というようなことをすればいいのでしょうか？

↑

名前のマッチの不具合†

BCD (2011-09-27 (火) 18:38:03)

時折、名前のマッチを行うと不具合が生じます。これは何故でしょうか? 今回起こった不具合までの経過は以下のようになります。どなたか解決法をご存知の方がいらっしゃれば教えていただけないでしょうか？

数値と値に対する名前が入っているベクトルA、ベクトルBがあります。

> AB <- c(names(A), names(B))   # 名前を取り出して一つのベクトルにします。
> CD <- AB[duplicated(AB)]   # 重複している名前を一つにする。
ここで以下のような行列xがあります。
V1には名前が入っていて、ベクトルAとベクトルBの名前は
この中に全て含まれています。
> x
      V1   V2
1  名前1   値
2  名前2   値
3  名前3   値
:
> x[x$V1 == CD, ] # CDの要素と行列xの名前とが
                    一致する行を取り出します。

これらを実行しても最後に名前が一致する行が取り出されないのです。この時エラーが出るのですが「長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません 」と表示されます。

↑

ベクトルからのマッチした要素の除去†

BCD (2011-09-27 (火) 17:58:16)

c("水", "金", "地", "火", "木", ......)　のようなベクトルがあります。
ここから　c("水", "地", ....)　を除いたベクトルを作りたいのですが、ifやforを使わずに簡単に行う方法はありますか？ご存知の方がいらっしゃいましたらどうか教えてください。

↑

任意の形状でプロット†

宗二 (2011-09-26 (月) 11:25:26)

plot関数を使用する際に、グラフィックスパラメータ"pch"で点の形を選択できますが、用意されている形状以外で、任意に点の形を変えることはできるのでしょうか。
具体的には、一点ごとに角度と長さの異なる矢印をプロットしたいのですが、良い方法がありますでしょうか。

↑

マルコフ連鎖†

square (2011-09-25 (日) 18:59:28)

マルコフ連鎖とは、乱数発生の方法みたいですが、sample関数による乱数発生とは何が違うのでしょうか？
どんな分布もsampleで発生できると思うのですが。

↑

閾値を求める方法について†

ランゲル・ハンス (2011-09-21 (水) 08:49:14)

いつも掲示板を参考にしております。
これまでExcelで計算していたのですが、Ｒでもぜひ計算できるようにしたいと思います。
ベースラインy=cと回帰直線y=ax+bの交点から閾値を求めたいと思います。
次のようなデータ（dataXY)があります。

X	Y
4	1.1
4.1	0.6
4.2	0
4.3	0.9
4.4	0.9
4.5	0
4.6	0
4.7	0
4.8	0.7
4.9	0.9
5	0
5.1	1
5.2	1.1
5.3	0
5.4	0.4
5.5	1.5
5.6	0.9
5.7	1.5
5.8	2.1
5.9	2.6
6	3.2
6.1	4.3
6.2	6.1
6.3	8.3
6.4	10
6.5	11.9
6.6	15.3

data <- read.table("dataXY.txt", header=T)

(1)データの1行目から16行目までのYの値の平均cを求める。

data1 <- data[1:16,]
c <- mean(data1$Y)
c

(2)データの20行目から27行目までのX,Yの値を一次回帰してy=ax+bのaとbを求める。

data2 <- data[20:27,]
L <- lm(Y~X, data= data2)
L

(3)y=cとy=ax+bの交点から閾値（X座標）を求める。
(4)最後に作図する。

plot(data)
abline(h=c, col=2)
abline(L, col=4)

(3)をご教示いただけないでしょうか？
また、グラフ中に交点の座標を表示したいと思います。
どうぞよろしくお願いいたします。

↑

readShapeLines関数とプロット範囲†

宗二 (2011-09-20 (火) 15:57:37)

maptoolsパッケージの"readShapeLines"関数を使って、GISデータをプロットしようとしているのですが、plot()関数内の引数xlim,ylimを使って範囲指定しても、実際のプロット範囲（緯度経度）がずれてしまいます。

練習用にshapefileライブラリにあるtokyok.zipを使って、以下のプログラムでプロットしたのですが、

library(maptools)
map <- readShapeLines("tokyok.shp")

lon1 <- 138   ### 経度138E-141E
lon2 <- 141
lat1 <- 35    ### 緯度35N-35.5N   
lat2 <- 35.5

par(usr=c(lon1, lon2, lat1, lat2))
plot(map, xlim=c(lon1, lon2), ylim=c(lat1, lat2))

作図領域を確認すると

> par("usr")
[1] 138.00000 141.00000  33.90625   36.59375

というように、（緯度の）指定範囲がずれてしまいます。
これは何が原因なのでしょうか？また、修正方法はあるのでしょうか？

↑

複数ページのpdfの作成†

まあくん (2011-09-20 (火) 11:35:53)

pdfを作成するときに、複数ページのpdfの作成はできないでしょうか？。複数のグラフを１度に作成していますが、現在は各ページごとにpdfファイルを作成し、手動で１つのpdfファイルにしています。

↑

発症率の検定？†

tau (2011-09-19 (月) 22:56:26)

いつも、色々と教えていただきましてありがとうございます。

暴露Yが疾患Xの発症に及ぼす影響を調べたいと考えています。
調査の対象はA県、B県の県民です。両県の県民共に暴露を受けており、暴露前後での疾患の発症数が判明しています。
A県では暴露前1ヶ月で100例、暴露後1ヶ月で110例の疾患Xの発症があり、B県ではそれぞれ、80例、120例の発症がありました。A県よりB県で暴露による疾患の発症率が増加したことを証明したいと考えています。
１．A県とB県の人口が分かれば単位人口当たり(/1000人など)の発症率を計算してそれを比較と言うことになるかと思われますが、この場合はどのような検定が用いられるのでしょうか？
２．また、A県とB県の人口が十分に大きいと言うことは分かるものの、詳細な人口が分からない場合、上記の仮説を証明する方法はありますでしょうか？(暴露前後での発症数の比を取ってその比を検定する方法はありますでしょうか)
ご教示いただければ幸いです。

↑

多目的最適化で特殊制限†

西田 (2011-09-15 (木) 00:08:11)

複数の入力値と制約条件から一つの出力が最適になるように分散共分散行列利用の進化的戦略cmaesというパッケージを使ったプログラムを作ろうと頑張っています。

制約条件である説明変数の下限上限がlower(0,0,0),upper(10,10,10)のように定数であればなんとか動いてくれるようになったのですが、下限上限が変動する、具体的にはx1+x2+x3=10のような条件が追加されたときにきちんと動かず困っています。

現在はcmaes内部で、新しく子集合を作る際に

x <- xmean + σ # 親集合の重心周辺に点を取っているのだと思う(σは変数)
x <- 10 * ( x / sum(x) )      #xの合計が5なら5で割ってその後10倍
...##色々計算　よくなっていたら更新
...##くり返し

として無理やりx1+x2+x3=10を実現しています。が、これできちんと最適化してくれているのか不安です。できればcmaes関数を変更することなく、つまりその外部でx1+x2+x3=10という条件を満たしたいと思っています。 何かいい方法はないでしょうか。(cmaes関数自体にコダワリはないです)

  • 最適化する関数内のx3を10-(x1+x2)に置き換えるとかではダメなのですか？ -- Iona 2011-10-01 (土) 03:30:25
  • わかりづらい質問にお答えいただきありがとうございます。どうやら関数内ではxをリストで持ち一括実行しているようで、動きとしてｘ１が変わったとしてもｘ３はそれを受けつけず一つ前の状態のｘ１で実行されてしまうようです。渡されたものをそのまま使っていましたがもっと一般的なRパッケージを用いてやってみようと思います。ありがとうございました -- 西田 2011-10-02 (日) 22:58:22
↑

一般化加法モデル(gam)の回帰式の出力の仕方†

Saito (2011-09-14 (水) 16:58:49)

いつもお世話になっております。ネットや過去ログで調べましたが、見つからなかったので質問させてください。mgcvパッケージの中に、一般化加法モデルが扱えるgamという関数があります。これで回帰式の推定を行うとそれらしい曲線が引けるのですが、その曲線の関数形が知りたいのです。以下にサンプルを示します。

> n <- 30
> x <- seq(-n, n)
> y <- 0.04*x^3+0.03*x^2+0.02*x+rnorm(2*n+1, 0, 1) # 真の関係
> plot(x, y)
> 
> res <- lm(y ~ x) # 直線回帰
> func <- function(x){
+ coef(res)[1] + coef(res)[2]*x # これが直線回帰の回帰式
+ }
> points(x, func(x=x))
> 
> res2 <- gam(y~s(x, k=3)) # 加法モデル（
                           # デフォルトは薄板スプライン）による回帰
> points(x, predict(res2), type="l", col=2, lwd=3)
                           # predictによる予測で何故か直線回帰と一致
> coef(res2) #回帰係数が何か出力されてはいるがヘルプを見ても分からない
  (Intercept)        s(x).1        s(x).2 
 9.309025e+00 -1.199402e-09  3.932092e+02 
> 
> res2 <- gam(y~s(x, k=4)) #節点の数を増やすと曲線となる
> points(x, predict(res2), type="l", col=2, lwd=3)
> coef(res2) 
(Intercept)      s(x).1      s(x).2      s(x).3 
   9.309025 -843.683785 -129.386077 1162.199292

パッケージを扱うのに問題はありませんが、結局、どのような回帰式になったのかを知りたいのです。ヘルプを見れば書いてあるはずなのですが、私には見つけられませんでした。どこかの本に、回帰式の求め方が書いてある、という情報だけでも構いません。
どなたか、ご教授頂けると幸いです。どうぞよろしくお願い致します。

↑

ablineで横軸の指定した範囲に直線を描く†

くま (2011-09-14 (水) 15:50:47)

lmで得られた回帰直線式をグラフに描くのに、ablineを使っていますが、これを使うと、直線は横軸の全領域にわたって引れてしまいます。指定した横軸の範囲にのみその直線を描くにはどうすればよいでしょうか。たとえば、グラフの横軸の範囲のxlim=c(0,100)のときに、回帰直線式を描く範囲をc(20,80)にすることは可能でしょうか。

↑

行列の任意の列に含まれる欠損値NAの割合算出†

(2011-09-08 (木) 13:46:05)

行列の任意の(全ての)列に、欠損値NAがどれだけあるか、割合を算出したいのですが、どうすればよいでしょうか。全ての列に対して、それぞれNAの割合を求めたいです。

例えば、100行10列の行列において、3列目にNAが40個あれば、NAの割合は40% ← この40(または0.4)を求めたい

宜しくお願いします。

↑

理論的な標本分布について†

manabu (2011-09-06 (火) 12:03:05)

本で勉強していて、以下のような問題がありました。

サンプルサイズをn=1、3、5、10と変化させた時に、標本分布の形状がどのように変わるか調べてみましょう。
※母集団分布は標準正規分布N(0,1^2)

答えは以下のように書いてありました。

curve(dnorm(x, sd=sqrt(1/10)), -3, 3)
curve(dnorm(x, sd=sqrt(1/5)), -3, 3, add=TRUE)
curve(dnorm(x, sd=sqrt(1/3)), -3, 3, add=TRUE)
curve(dnorm(x, sd=sqrt(1/1)), -3, 3, add=TRUE)

母集団の標準偏差や分散が1である時に標本の偏差の二乗の合計が1と言えるのは何故でしょうか？

↑

SPSSのWEIGHTと同等の機能†

じゃんぽけ (2011-09-05 (月) 10:41:06)

Rを利用してアンケート分析を実施したいのですが、
SPSSのWEIGHT(重み付け）の様な機能はどうすれば再現できるでしょうか？

↑

image における時系列軸の反転†

てるよ (2011-09-01 (木) 11:59:52)

投稿の前にテキスト整形のルールを読んでも、ほとんど意味がわかりませんでした。お見苦しい投稿になっていましたら、どうぞお許しください。

image() における時系列軸の反転について教えていただけたら幸いです。

"2008/12/19"から"2009/3/24"の96日間に、等間隔でデータを記録したデータがあります。等間隔とは、1時間に1点で1日に24点を記録したデータで、要素数が94日×24点の2304です。このデータは1または0のバイナリデータです。

このデータで1がどのような時間帯に現れるのかを調べるために、縦軸（y軸）に日付、横軸（x軸）に時刻（0時から24時）をとってプロットした図を書くために以下のコマンドを実行しました。

st <- c("2008/12/19", "2009/3/24")
st2 <- strptime(st, "%Y/%m/%d", tz="")
y.days <- seq(st2[1],st2[2], by="1 day")
x.hours <- seq(0, 24, by=1)

m <- matrix(sample(c(0, 1), 24*length(y.days), rep=T),
            nrow=24, ncol=length(y.days))

image(x.hours, y.days, m, col=c("white", "black"), yaxt="n",
      xaxp=c(0,24,4), xlab="hours", ylab=NA)

r <- as.POSIXct(round(range(y.days), "days"))
axis.POSIXct(2, at=seq(r[1],r[2], by="20 day"), 
             format="%Y/%m/%d", las=2)

これでとりあえずグラフは書けたのですが、y軸の時間軸を反転させて、日付の進行を上から下に向かって進むように変更したいと思ったのですが、これがなかなかうまくいきません。

四苦八苦しているうちに、

image(...,  ylim=c(as.numeric(y.days[96]), as.numeric(y.days[1])))

を書くことでプロット自体の反転はできたのですが、今度は

r <- as.POSIXct(round(range(y.days), "days"))
axis.POSIXct(2, at=seq(r[1],r[2], by="20 day"), 
             format="%Y/%m/%d", las=2)

を実行してもy軸に日付が表示されなくなってしまいました。

y軸の時系列軸を反転させた後も、日付を表示するうまい方法がありましたら教えていただけると幸いです。

↑

Linuxでのメモリ不足解消†

BCD (2011-08-29 (月) 15:17:00)

初歩的なことかもしれないので恐縮なのですが、教えてください。
Linuxのスパコンを使ってRのheatmap.2関数で作図をしているのですが 'サイズ 1.4 Gb のベクトルを割り当てることができません'というエラーがでました。調べてみたところ、メモリが足りないということでそれまでに使ってきた変数を削除したりgc()をしてメモリの節約もしたりしました。windowsではメモリの制限を解除するようなコマンドがあるのは分かりましたが、Linuxではどのようにすれば良いのか分かりませんでした。この問題を解決する方法を教えてください

  • それって32bitのRではないでしょうか？ 本当に64bitのRですか？ -- 2011-09-08 (木) 18:23:10
  • 確かに64bitのRだと思います -- BCD 2011-09-17 (土) 20:47:22
↑

データフレームの抽出と並び替え†

BCD (2011-08-28 (日) 12:16:08)

>x
	1	2	3	4	.	.
A	20	50	60	29    .       .
B	49	74	90	38    .       .
C	57	49	20	39    .       .
>Y
[1] "C"	"A"	.....

のようなデータフレーム'x'(行、列ともに20000行ある)とxの全行名の一部がある順番に従って格納されているベクトル'Y'(C, A,......7000個。要素の順序は重要)が入っているベクトルがあります。このベクトル'Y'を参照して行名から'x’を抽出&並べ替える方法を探しています。

現在、私は以下のスクリプトでこの作業を実行しています。しかし、データが膨大で時間がかかってしまいます。この作業が早くなる方法をご存知の方がいらっしゃればどうか御教示願います。

result <- NULL
for (i in 1:length(Y)) {
    result <- rbind(result, x[rownames(x) == Y[i], ])
}
↑

一様性の検定について†

tau (2011-08-27 (土) 10:46:47)

一様性の検定に関して質問致します。
発症が希な疾患を有する患者が一年間に何人入院したかを2008-2011年の4年間で比較したいと考えています。2008年から2011年までそれぞれ1,2,1,7例ずつ入院があり、2011年のみ入院が多い事を証明しようと思います。
カイ二乗検定を用いて下記のように検定を行いましたが、度数が少ないためか「不正確かもしれません」のコメントがでます。

> chisq.test(c(1,2,1,7))
Chi-squared test for given probabilities
data:  c(1, 2, 1, 7)
X-squared = 9, df = 3, p-value = 0.02929
警告メッセージ： 
In chisq.test(c(1, 2, 1, 7)) :  カイ自乗近似は不正確かもしれません

分割表でカイ二乗検定を行うとき同様、度数が少ない場合は、この検定法を用いない方が良いのでしょうか？
不適当な場合、Fisherの正確検定を用いるのかと考え、調べてみましたが、同方法による一様性の検定の方法がわかりませんでした。
可能であればFisherの正確検定による一様性の検定の方法を、不適当であれば他の検定方法をご教示頂ければ幸いです。

↑

質問紙尺度データのクロス集計の値をパーセントになおす方法†

Wombat (2011-08-25 (木) 09:16:29)

私の書き方が不明瞭だったようなのでもう一度投稿させて下さい。以下のデータは５が各質問に対して「強くそう思う」で１が「まったくそう思わない」になる質問紙尺度と呼ばれるものになります。

     ID class Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q11 Q15
1     1     a  5  4  4  4  4  2  4   4   5
2     2     a  4  4  3  3  4  4  3   4   2
3     3     a  4  4  3  5  4  4  3   4   5
4     4     a  4  4  1  3  3  4  5   3   5
5     5     a  5  4  1  5  5  3  2   4   4
6     6     a  4  4  4  4  4  2  1   2   1
         　　　-- 省略--
94   94     d  4  4  3  4  4  2  3   4   4
95   95     d  4  3  4  3  3  4  4   4   4
96   96     d  4  4  1  4  3  3  2   4   3
97   97     d  4  5  2  5  5  5  4   4   4
98   98     d  5  4  4  4  3  4  5   4   4
99   99     d  4  4  2  4  4  3  3   4   4
100 100     d  4  4  2  3  4  4  5   4   4
101 101     d  5  4  2  4  4  5  4   4   1
102 102     d  5  4  4  4  3  5  2   4   2

これを各クラスごとに項目別にsapply(d01[,3:11],function(x)　table(d01$class,x))でクロス集計を行なう（河童さんのおかげです。）と以下のようになります。

$Q1
   x
     3  4  5
  a  2 15  8
  b  0 21  6
  c  2 19  4
  d  0 14 11

$Q2
   x
     2  3  4  5
  a  1  5 19  0
  b  1  6 20  0
  c  0  8 14  3
  d  0  4 17  4
 -- 以下省略--

これにより、Ｑ１はclass aは強くそう思うと答えた回答者が８人いることがわかりますが、クラスの数が同じではないので、８人ではなくて、32%というふうにパーセントで報告するのが普通になると思うのですが、いちいち手で計算するのは効率が悪いので、この８という頻度を表す数字をパーセントになおす方法をご存知でしたら教えてください。よろしくお願いします。

↑

クラスと項目別のクロス集計の値をパーセントになおす方法†

Wombat (2011-08-24 (水) 19:26:14)

河童さんのおかげでクラスと項目別のクロス集計ができるようになりましたが、論文に書くときにはクラス間で人数の差があるので、パーセントに変換して書きます。１つ１つ計算してもいいのですが、もし、その値を簡単にパーセント値に変換する方法がありましたら教えてください。よろしくお願いします。OSはWindows XP, RのVersionは2.13.1です。

$Q1

  x
    3  4  5
 a  2 15  8
 b  0 19  3

$Q2

  x
    2  3  4
 a  1  5 19
 b  1  5 16
↑

一度にクラスごとに各項目を集計（和を求める）する方法†

Wombat (2011-08-24 (水) 01:52:47)

次のようなデータがあります。これはある英語教授法を４つのクラスでおこなって、英語力のどの部分が伸びたかを学生に尋ねたアンケートの結果を表したものです。複数回答項目ですので、各項目ごとに変数を準備して、選択したら１、しなければ０で処理しています。このデータの集計でクラスごとに各項目の和を一度に求めるにはどのようにしたら良いでしょうか。このクラスごとにどの項目が多いかを調べらばいいので、クロス集計でなく、単に和を求めるだけで良いと思います。sapply(dat01[1:25,3:9],sum)でclass aにおける各項目の集計は求めることができますが、何度も[ ]の中の行の値を変えるのはあまり効率の良い方法とは言えないと思います（４回くらいやれよと言われそうですが）ので、よろしくお願いします。ＯＳはWindo coef(res2)

(Intercept)      s(x).1      s(x).2      s(x).3 
   9.309025 -843.683785 -129.386077 1162.199292

パッケージを扱うのに問題はありませんが、結局、どのような回帰式になったのかを知りたいのです。ヘルプを見れば書いてあるはずなのですが、私には見つけられませんでした。どこかの本に、回帰式の求め方が書いてある、という情報だけでも構いません。
どなたか、ご教授頂けると幸いです。どうぞよろしくお願い致します。

↑

グラフの細かいx軸範囲制御†

BCD (2011-08-22 (月) 19:09:35)

plot関数で折れ線グラフを作る時にxlimでグラフの中のx軸の範囲を決めることができますが、その中のある領域だけそのグラフ内で広げるようにすることはできますか？
例えば xlim = C(0, 10000)　のグラフの場合、グラフの左端が０で右端が10000になり中心が5000になります。これを操作してグラフのx軸の左端から1/4までが0~4500の範囲、1/4から3/4までが4500~5500の範囲、3/4から右端までが5500~10000の範囲と部分的にx軸を広げたいのです。

↑

read.tableでファイルを読み込む時に、先頭が数字である列ラベルに付加される'X'を取り除きたい†

BCD (2011-08-22 (月) 13:13:44)

read.tableでファイルを読みこむと列ラベルの先頭が数字の場合、その前にXが付加されてしまいます。Xが付加されないようにするにはどうしたらよいですか？(以前、Rjpwikiに投稿があったと思いますが、見つけ出すことができませんでした)

ex)
     0002jh 0043sss
  1         2           4

のようなtxtデータをread.tableで読み込むと

    X0002jh X0043sss
   1       2        4

となります。これを解消したいです。

↑

複素数を含む積分(変数は実数)†

square (2011-08-21 (日) 17:47:20)

下のような複素数の積分はできませんか？
f=function(x) exp(i*x)
integrate(f,0,1)
integrateでは実数限定のようですが。

↑

data.frameをtransactionに変換する際の問題†

アレックス (2011-08-17 (水) 15:58:26)

R 2.13.1を利用しています。

要素1,要素2,要素3
りんご,バナナ,みかん
バナナ,なし,
パイン,もも,みかん

という中身のファイル(data.txt)を以下のコマンドで読み込み

data1 <- read.table("C:/lab/kankou.txt",sep=",",
                    header=TRUE)
　要素1  要素2  要素3
1 りんご バナナ みかん
2 バナナ なし
3 パイン もも   みかん

というdata.frame形式のデータを読み込むことはできたのですが、これをtransactionに変換する際に

data.tran <- as(data1, "transactions")
data.frame <- as(data.tran, "data.frame")

としてdata.frameを開いてみると

1{要素1=りんご,要素2=バナナ,要素3=みかん}
2{要素1=バナナ,要素2=なし,要素3=}
3{要素1=パイン,要素2=もも,要素3=みかん}

となってしまいます。これを

1{りんご,バナナ,みかん}
2{バナナ,なし,}
3{パイン,もも,みかん}

となるようなtransactionを生成する方法はないでしょうか。

↑

Plmをインストールする方法†

z (2011-08-12 (金) 18:52:21)

Plm Pacakgeを用い、下記の分析を行いたいのですが、うまくインストールできません。Rは2.13.1を利用しており、スペックは満たしているようです。
このぱっページをインストールする方法か、類似の分析方法をご教示頂けましたら、幸いです。自身でも、plmを使わずに分析を行うコードを色々と試行錯誤しているのですが、なかなか解決できていません。宜しくお願い致します。
the Fixed effects model (within),
• the pooling model (pooling),
• the first-difference model (fd),
• the between model (between),
• the error components model (random).

↑

混合モデルの重回帰†

酔鯨 (2011-08-12 (金) 17:05:34)

使用しているOS=Windows　XP　Rのバージョン=2.10.0
独立変数が、定量データと定性データが存在する混合モデルを重回帰分析をしました。定性データは（０，１）の2値データです。この場合は、重回帰データとして同じ扱いが出来ると市販本に書かれていたので、AIC最小モデルを変数減増法と群馬大青木先生の総当り法で求めました。結果は同じになりました。定性データは、最小モデルには残りませんでした。
単回帰の結果では、定性データでモデルを別々にしなければならないという結果が出たのですが、重回帰の方法に誤りは無いでしょうか?

↑

RMA回帰の予測区間について†

ryuji (2011-08-11 (木) 23:37:05)

現在、ある昆虫の体の一部から体長を推定するために、これまで採集した標本をもとにRMA回帰を求めました。そこまではRを使ってできたのですが、実はこの昆虫の死骸が多く見つかるため、一体これらの死骸の体長がどのくらい大きさであったかを推定したいと思っております。そこでRMA回帰に代入して体長を求める前に、予め予測区間内に入るかどうか、グラフで表示しておきたいのですが、その方法が分かりません。超初心者で申し訳ないのですが、ご教授願えないでしょうか？よろしくお願いいたします。

↑

cannot coerce type 'closure' to vector of type 'any'　というメッセージ†

tavia (2011-08-11 (木) 20:46:50)

お世話になります。

Rを使って、クラスター分析をしようとしています。
データの読み込みまではできるのですが、その後
d <- as.dist(data)
を実行すると、

以下にエラー as.vector(x, mode) : 
 cannot coerce type 'closure' to vector of type 'any'

というメッセージが出ます。
このようなメッセージが出たら、どういうところをチェックすればよいでしょうか。
教えていただければ幸いです。

Rのバージョンは2.10.1
OSはWindows Vista　Business　S.P.2
です。
よろしくお願いいたします。

↑

configure で失敗してしまいます。†

nob (2011-08-11 (木) 11:31:15)

Rそのものは超初心者ですが、結構前から使っています。
Plamo Linux 4.7.3でのつまずきです。
いつもパッケージをCRANからダウンロードしてconfigure、make、make install でシステムをインストールしています。が、plamo Linuxのバージョンを4.7.3にしたら、
checking for xmkmf... /usr/X11R7/bin/xmkmf
configure: WARNING: I could not determine FPICFLAGS.
configure: error: See the file doc/html/R-admin.html for more information.
で止まってしまいます。2.13ばかりでなく古いバージョンも同様です。こちらのHPを見たらconfigureで失敗した場合不足しているものをconfig.logで確認して再度configureせよとのことだったのですが見ても何が不足なのかよくわかりません。
どなたか同じような件で解決された方がおりましたら、どのように解決されたか教えていただけないでしょうか？
よろしくお願い致します。

↑

グループごとに各質問項目の答えの頻度を計算する方法†

Wombat (2011-08-10 (水) 12:18:35)

Ｒ超初心者です。御回答よろしくお願いします。以下のようなデータがあります。

     ID class Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q11 Q15
1     1     a  5  4  4  4  4  2  4   4   5
2     2     a  4  4  3  3  4  4  3   4   2
3     3     a  4  4  3  5  4  4  3   4   5
4     4     a  4  4  1  3  3  4  5   3   5
5     5     a  5  4  1  5  5  3  2   4   4
6     6     a  4  4  4  4  4  2  1   2   1
7     7     a  5  4  2  2  4  3  4   4   2
8     8     a  4  4  2  5  4  3  3   4   3
9     9     a  4  4  2  3  3  4  3   5   5
10   10     a  5  4  4  4  5  3  3   4   4
11   11     a  4  3  2  3  2  3  1   3   3
12   12     a  3  2  3  4  3  5  4   3   2
13   13     a  5  4  2  4  4  4  5   4   4
14   14     a  3  3  4  3  4  3  1   3   2
15   15     a  4  4  3  4  4  4  4   4   3
16   16     a  5  3  2  4  4  3  3   5   5
17   17     a  4  4  2  4  3  3  4   3   3
18   18     a  4  4  2  2  4  3  2   3   4
19   19     a  5  4  4  4  4  2  5   4   5
20   20     a  4  4  3  4  4  5  4   4   5
21   21     a  4  4  1  4  3  4  4   4   2
22   22     a  4  3  2  3  3  5  3   3   3
23   23     a  4  3  3  3  4  3  3   4   3
24   24     a  5  4  1  4  4  4  4   4   2
25   25     a  4  4  3  4  4  5  4   3   3
26   26     b  4  4  1  4  4  5  3   5   5
27   27     b  5  4  3  4  3  3  3   4   3
28   28     b  4  3  2  4  3  5  4   3   2
29   29     b  4  4  2  3  3  4  4   4   2
30   30     b  4  4  2  3  4  4  3   4   3
31   31     b  4  4  4  3  4  4  2   4   4
32   32     b  4  4  4  4  3  4  3   4   3
33   33     b  4  4  4  3  3  4  4   3   3
34   34     b  4  4  5  3  3  4  3   3   3
35   35     b  5  4  3  4  4  3  3   4   4
36   36     b  4  3  1  3  3  4  5   3   1
37   37     b  4  4  5  4  3  4  4   4   2
38   38     b  4  4  5  2  2  4  2   3   3
39   39     b  4  4  4  4  3  4  2   4   3
40   40     b  4  4  3  4  4  3  3   4   4
41   41     b  4  3  2  3  3  4  4   4   2
42   42     b  5  3  3  4  3  4  5   3   2
43   43     b  4  3  2  3  4  5  5   4   5
44   44     b  4  2  3  3  3  4  3   3   3
45   45     b  4  4  3  5  3  3  5   4   1
46   46     b  4  4  3  3  5  4  5   4   4
47   47     b  4  4  3  4  3  4  4   4   2

Ｑ1からＱ１１の数値はリッカートスケール（５が強く思う〜１全くそう思わない）です。この表ではグループ（Class)はaとｂしかありませんが、実際のデータではdまであります。例えば、aクラスにおけるQ1における数値の頻度は> table(dat01[1:25,3])で求めることができことは「Ｒによるやさしい統計書」インターネットでわかりました。しかし、この方法では時間がかかってしまうので、これを一度にグループことの各質問項目におけるリッカートスケールの数値の頻度を求めるにはどのような関数を使えば良いのか教えてください。

↑

grepで検索元の単語に「ー」が含まれていても正しく検索できる方法について†

ひろ (2011-08-09 (火) 18:14:41)

grepコマンドの実行で，検索元(パラメータ1)の単語の中に「ー」が含まれていると下記通りエラーメッセージが表示されて失敗だと失敗します．他の「日」「月」「年」といった言葉であれば問題なく検索できました．

検索元の単語に「ー」が含まれていても正しく検索できる方法はございますでしょうか？

grep("リー","あ","い","う","え")

以下にエラー grep("リー", "あ", "い", "う", "え") : 
  "リー" は不正な正則表現です, 理由は 'Missing ']''
↑

Rコマンダーで日本語が使えません†

nappa (2011-08-04 (木) 12:50:45)

下記のように日本語があると「不正なマルチバイト文字があります」と表示され、うまく動作しません。コメントなのですがダメです。以前は大丈夫だったのですが、エンコードの仕様が変わったためでしょうか？
対策がありましたら、教えてください。

#サンプル
x <- 2
x

Ｒ　Version:  2.13.1
ＲコマンダーVersion:  1.6-4
↑

パッケージ"adapt"がインストールできません†

しげ (2011-08-03 (水) 16:04:22)

パッケージ"adapt"をインストースしようとすると

In getDependencies(pkgs, dependencies, available, lib) :
  package ‘adapt’ is not available (for R version 2.13.1

と警告メッセージが出てインストールできないのですが、どのようにしたらインストールできるようになりますか？使用環境は

> sessionInfo()
R version 2.13.1 (2011-07-08)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Japanese_Japan.932
    LC_CTYPE=Japanese_Japan.932   
[3] LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C                  
[5] LC_TIME=Japanese_Japan.932    

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

other attached packages:
[1] rgl_0.92.798

loaded via a namespace (and not attached):
[1] tools_2.13.1

です。かなり初歩的なことだと思いますが、教えてもらえると嬉しいです。よろしくお願いします！

↑

小数点のデータをrankでソートできますか？†

ひろ (2011-07-29 (金) 16:01:26)

小数点(ex 4.564e-05, 4.234e-02 etc.)のデータセットはrankでソートできないのでしょうか？これをソートすると正しく順序づけられなく困っております．並び替える方法はございますでしょうか？お知恵をお貸いただけますでしょうか？

■ コマンド
res.sorted <- res[rank(res$D),]
データセットresの中のD(eviance)の列は小数点データで，これをrankコマンドで昇順に並べた添え字リストを得て，resを並び替えた結果をres.sortedに入れようとしています．

↑

ピーク面積の求め方†

どぜう (2011-07-26 (火) 10:41:16)

　お世話になっております。過去ログにも無いようですので、質問させてください。
　とあるスペクトルデータがあり、ピーク面積を求めたいのです。前処理でバックグランドデータを引いたデータを用意するまではよかったのですが、ピークのあるチャンネル番号や、そのピークが囲む面積（ピクセル数？）の求め方が分かりません。
　よろしくお願い申し上げます。

↑

差の検定について†

Tau (2011-07-21 (木) 22:41:07)

統計手法の選択についてご教示ください。
ある疾患が、1週目に4例、2週目に5例、3週目に6例、4週目に17例生じたとします。1週目から3週目までの疾患発症数の平均と4週目の疾患発症数の差を検定する場合は、t.testあるいはノンパラメトリックな手法(U.testなど)を用いて良いものでしょうか？(つまり、比較したい2群があって一方は複数のcaseがあり平均値が算出されるのに、他方は一つのcaseしかない場合です)
ご教示いただければ幸いです。

↑

propensity score matchingについて†

sui? (2011-07-14 (木) 15:45:22)

propensity scoreを用いたマッチングの手法について質問させて頂きます。
Rはver.2.13.1、Matching 4.7-14を用いています。
datasetはMaching package内の「lalonde」を用いています。
以下の通り、propensity scoreを算出します。

> library(Matching)
> data(lalonde)
> Y <- lalonde$re78
> Tr <- lalonde$treat
> glm1 <- glm(Tr~age+educ+black+hisp+married+nodegr+
              re74+re75, family=binomial, data=lalonde)

glm1$fittedがpropensity scoreです。 これを用いてMatchingを行います。

> rr1 <- Match(Y = Y, Tr = Tr, X = glm1$fitted)
> summary(rr1)

Estimate...  2624.3 
AI SE......  802.19 
T-stat.....  3.2714 
p.val......  0.0010702 

Original number of observations..............  445
Original number of treated obs...............  185 
Matched number of observations...............  185 
Matched number of observations  (unweighted).  344

サマリーを見ると185例でMatchingが行われたようですが、Matchig後のデータセットを表示したいと考えております。お分かりの方がいらっしゃいましたら教示頂けませんでしょうか？
よろしくお願い申し上げます。

↑

行列データの修正†

内藤 (2011-07-14 (木) 10:28:10)

以下の行列データがあります（本物はもっとデータが大きい）

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    0    0    0    1    0    0    0    0    0     0
 [4,]    0    0    0  999    0    0    0    0    0     0
 [5,]    0    0    0  999    0    0    0    0    0     0
 [6,]    0    1    0  999    0    0    0    0    0     0
 [7,]    0    1    0  999    0    0    0    0    0     0
 [8,]    0    1    0  999    0    0    0    0    0     0
 [9,]    0    1    0  999    0    0    0    0    0     0
[10,]    0  999    0  999    0    0    0    0    0     0
[11,]    0  999    0  999    0    0    0    0    0     0
[12,]    1  999    0  999    0    0    0    0    0     0
[13,]    1  999    0  999    0    0    0    0    0     0

999以降の列データを2を8個,3を10個というようにプログラムで変更したいのですが，上手くいきません。
初歩的な質問で申し訳ありませんが，上手いやり方を教えてください。

↑

latticeパッケージにおけるxyplotの重ね描きの方法†

Ryohei (2011-07-12 (火) 22:12:18)

以下のxxというデータセットをID別にTIME2対VAR1とTIME2対VAR2の折れ線グラフを重ね描きしたいと考えております。

#1に示した図をアウトプットのイメージとして、#2の2種類のグラフを1枚で重ね描きをして変数VAR1とVAR2に対する時間(TIME2)の推移を同時に表現したいと考えております。xyplotではどのようにプログラムを書けばよいのでしょうか？
今回は便宜上2行2列のグラフでしたが、実際には4行4列、5行5列といったグラフを描こうと考えております。
xyplotの特長である、軸ラベルが両端・最上下段のみに表示されすっきりする、IDごとのグラフエリアが大きくなって見やすくなるといった点を活かしたいと思っております。
どうぞよろしくお願いいたします。

#1
N <- 4 ; TIME <- seq(6)
TIME2 <- rep(TIME, N)
ID <- rep(seq(N), each=length(TIME))
VAR1 <- rep(seq(3), length(TIME)/3*N)
VAR2 <- round(runif(N*length(TIME), min=1, max=9))
(xx <-  data.frame(ID, TIME2, VAR1, VAR2))
layout(matrix(1:4, byrow=TRUE, ncol=2))
old <- par(mar=c(2, 2.7, 1.3, 1), mgp=c(1.6, 0.7, 0))
for (i in 1:N) {
	plot(xx[xx[,"ID"]==i, "TIME2"], xx[xx[, "ID"]==i,
            "VAR1"], xlim=c(0, 7), ylim=c(0, 10),
	     xlab="", ylab="", main=sprintf("ID%02d", i),
            col="red", type="b")
	points(xx[xx[, "ID"]==i, "TIME2"],
              xx[xx[, "ID"]==i, "VAR2"],
	xlim=c(0, 7), ylim=c(0, 10), col="blue", type="b")
	}
par(old)
layout(1)
#2
#library(lattice)
#xyplot(VAR1 ~ TIME2 | ID, data=xx, ylim=c(0, 10),
        type="b", col="red")
#xyplot(VAR2 ~ TIME2 | ID, data=xx, ylim=c(0, 10),
        type="b", col="blue")
↑

package:rattleがインストールできない†

sakura (2011-07-09 (土) 07:57:43)

PCのシステムは、
Microsoft Windows XP
Professional
Version 2002
Service Pack 3
です。
R version 2.130(2011-04-013)で、package:rattleをロードしてdeta mining を学ぼうとしているのですが、パスが通らず？次のエラーが出てしまいます。
どう対処したら良いのか、分かりません。ご教示頂ければ、幸いです。お願いします。

Learn more about GTK+ at http://www.gtk.org
If the package still does not load, please ensure
that GTK+ is installed
and that it is on your PATH environment variable
IN ANY CASE, RESTART R BEFORE TRYING TO LOAD THE
PACKAGE AGAIN
 以下にエラー as.GType(type) : Cannot convert
              RGtkBuilder to GType
 追加情報:   警告メッセージ： 
In unzip(path, exdir = .windows_gtk_path) :
   zip ファイルから抽出中に書き込みエラーが生じました
↑

Rexcel:RthroughExcelWorksheetsが無い†

Ackin (2011-07-01 (金) 14:02:41)

Office2007、WindowsXPの職場PCにRAndFriendsSetup2130V3.1-15-1.exeをインストールしてRexcelを学び始めました。
インストール時は管理者権限で行いましたが、使用時に権限はありません。C:\Program Files\RExcelにインストール出来たようです。
職場PCはマイドキュメントをサーバにあるP:としてあり、getwd()すると"P:/"となります。

【症状】
ExcelメニューからRexcel>RthroughExcelWorksheetsとしたいのですが見当たりません。
書籍「ExcelでR自由自在」P.34図3.1によるとDemo WorksheetsとAbout Rexcelの間にあるようです。
C: P: を検索しましたがBookFilesTOC.xlsmは見当たりません。

どこからかBookFilesTOC.xlsmを探してきて、パスを通さなければ使えないでしょうか？

↑

scatterplot3dの軸調整について†

ランゲルハンス (2011-06-22 (水) 14:29:19)

いつも掲示板を参考にしております。
三次元散布図を描きたいと思います。
下記のプログラムでは、ｘ、ｙ、ｚ軸の目盛は違っても長さが同じで立方体の中に各点が散布されます。
それぞれの軸を実際の長さの比（10：20：5 =2:4:1）にした直方体の中に点を散布する方法をご教示いただけないでしょうか？
見た目の錯覚を防ぎたいと思います。
よろしくお願いいたします。

library(scatterplot3d)
x <- runif(100, 0, 10)
y <- runif(100, 0, 20)
z <- runif(100, 0, 5)
scatterplot3d(x, y, z, pch=20)
↑

ヒストグラムの横軸に度数、縦軸に階級を取る方法†

buta (2011-06-21 (火) 13:07:20)

ヒストグラムの横軸に度数、縦軸に階級を取るにはどのようにすればよいでしょうか。よろしくお願いします。

↑

pasteでベクトル内の文字列を短いコマンドで連結する方法†

poro (2011-06-18 (土) 04:00:33)

str <- ""
for (i in c("a", "b", "c")) {
    str <- paste(str, i, sep="")
}

としていますが、ベクトルをpasteコマンドにそのまま渡して文字を連結してくれないかと期待するのですができないものでしょうか？

↑

breaksが小数となるヒストグラムの密度表示が正しくない†

石ころ太郎 (2011-06-15 (水) 19:19:09)

お世話になります。

hist(c(0.1*1:10), probability=TRUE)

を実行したときにdensityが0.1となるべきところが、1.0となってしまいます。いろいろ確認したところ、breaksの閾値が小数となるときにあらわれる現象と見受けられます。
解決法はございますでしょうか？

↑

四分位数にバグ？†

三十路 (2011-06-14 (火) 21:23:44)

c(1,2,3,4,5,100,1001,1002,1003,1004,1005)
をsummaryで見ると、

  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.0     3.5   100.0   466.4  1002.0  1005.0

という結果が表示されるのですが、この3rd Qu.は間違ってますよね。
c(1,2,3,4,5,10,101,102,103,104,105)
ならば正しく、

  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.00    3.50   10.00   49.09  102.50  105.00

が返されるのですが。

↑

RColorBrewerが使えません†

HS (2011-06-14 (火) 18:15:53)

R (2.13.0, x86_64-pc-mingw32/x64) 上で RColorBrewer (1.0-2) を使いたいのですが，

> library(RColorBrewer)
エラー：  パッケージ 'RColorBrewer' は R 2.10.0 以前に造らました。新しくインストールして下さい

となってしまいます．64bit環境のためなのでしょうか？
解決方法等ご存知でしたらご教示ください．

↑

lme()で平方和を出力する方法†

TK (2011-06-14 (火) 12:45:30)

線形混合モデルを使った時に、複数の固定効果の相対的な重要性を検討するために、分散分析表を出力しようとしています。分散分析表は、自由度、平方和(Sum of Square)、平均平方(Mean Square)、F値、P値からなっていることが多いので、このような標準的な分散分析表を執筆中の論文に含めたいと考えています。標準的な表を載せる方が読者も理解しやすく、査読者からも受け入れられやすいと考えるからです。
しかし、lme()を使ってもlmer()を使っても、このような出力ができません。Rで呼び出せるOrthodontのデータを例にすると、以下のような作業です。
使用環境は、Windows XP、R2.13.0、nlme_3.1-100、lme4_0.999375-39 です。

library(nlme)
res.lme <- lme(distance ~ age + Sex,
               random= ~ 1|Subject, data=Orthodont)
anova(res.lme)

library(lme4)
res.lmer <- lmer(distance ~ age + Sex + (1|Subject),
                 data=Orthodont)
anova(res.lmer)

lme()の場合、分子自由度、分母自由度、F値、P値が出力されますが、平方和と平均平方が出力されません。
lmer()の場合、分子自由度、平方和、平方平均、F値が出力されますが、P値が出力されません。lmer()がなぜP値を出力しないかは、https://stat.ethz.ch/pipermail/r-help/2006-May/094765.htmlにDouglas Batesさんの解説があるのを見つけました。混合モデルにおける分母自由度とP値の計算は複雑な問題なので、lmer()では出力しないということのようです。ということはlme()が出力するP値も何らかの問題があるものなのでしょうが、lme()のモデルの平方和が分かれば、分散分析表を完成させられるので、今回の論文ではlme()に基づこうと考えました。
lmer()が返す平方和を用いようかと考えたのですが、lme()とlmer()とでは変数によってF値が微妙に異なっていますので、平方和も異なっていそうです。どうしたらいいのでしょうか。

↑

一般プロクラステス分析†

タカ (2011-06-13 (月) 17:27:24)

いつもお世話になっております。
GPA（一般プロクラステス分析）について質問です。

   Acid Strange Hard Acid1 Strange2 Hard2
C1  2.0       2  2.0     2        1   2.0
C2  4.0       2  1.0     4        3   4.0
C3  7.0       1  2.5    11        1   7.0
C4  4.5       2  1.5     1        2   2.8
C5  5.0       3  4.0     5        6   4.0

以上のデータ（test）に対して、

library(FactoMineR)
res.GPA <- GPA(test, group=c(3, 3))

としたのですが、「以下にエラー if (a >= 0) prob <- pgamma(rvstd - (-2/a), shape = (4/a^2), scale = (a/2), : TRUE/FALSE が必要なところが欠損値です」とエラーが表示されてしまいます。 御存じの方いらっしゃいましたら、ご教授の程よろしくお願いいたします。
使用環境は、WindowsXP、R 2.13.0です。

↑

パッケージneuralのインストールについて†

M_Saito (2011-06-11 (土) 16:12:29)

初めて投稿します。この投稿欄にあるhiro?氏と同じような状況ですが、libraryからneuralのパッケージが削除されているようです。そこでhiro?氏の記事を参考にインストールを試みましたが駄目です。当方はOS WindowsXp, R2.12.2です。ダウンロードしたneuralアーカイブは.nueral1.4.2.1.tar.gzですので解凍後ファイルの中身を確かめ再度そのファイルをzipに変換してインストールを試みました。しかし次のようなメッセージが出てしまい、完了しません。

> utils:::menuInstallLocal()
 以下にエラー gzfile(file, "r") :  コネクションを開くことができません 
 追加情報:   警告メッセージ： 
In gzfile(file, "r") :
  圧縮されたファイル 'neural/DESCRIPTION'
  を開くことができません,
  理由は 'No such file or directory' です 
> install.packages("C/neural.zip",repose=NULL)
ファイル名から 'repos = NULL' を推測
 以下にエラー zip.unpack(pkg, tmpDir) : 
   zip ファイル 'C/neural.zip' が見付かりません 
> install.packages("C/neural.zip",contriburl=NULL)
 以下にエラー zip.unpack(pkg, tmpDir) : 
   zip ファイル 'C/neural.zip' が見付かりません

ファイルは作業デレクトリに置いてあります。どこが間違っているのか、又その解決策をご教授願いたいのです。下らない質問で済みません、宜しく御願い致します。

↑

ロジスティック回帰の曲線のグラフを描きたい†

うー (2011-06-02 (木) 17:03:52)

ロジスティック回帰をしています。
量xと反応の割合pについて、glm関数でロジスティック回帰式を出しました。
xとpについてプロットを行い、この図に回帰曲線を加えたいと思います。

単回帰であればabline関数とlm関数を用いてコマンド一つで回帰直線を描くことができますが、ロジスティック回帰についてはそれが見当たりません。

いくつか本やウェブページを調べたのですが、独立変数を非常に細かくとってpredictとpointsで作図するとか、glmで求めた予測値についてlines関数で補間するなど、ちょっと回りくどい感じがします。

curve関数で書いてしまおうかとも思ったのですが、高水準作図関数なので実データとの作図の順番が逆になってしまいます。

何かいい方法はないでしょうか？　自分で関数を書く以外ないでしょうか？

↑

テーブル出力時の列名のズレ†

ぺーぺー (2011-06-01 (水) 21:35:13)

相関行列をテキストファイルに出力したいのですが、どうしても列名（1行目のデータ）が1列左へズレてしまいます。

例として、以下のコードで再現可能です。

a<-iris[,1:4]
b<-cor(a)
write.table(b,"C:/b.txt",append=F,
            quote=F,col.names=T,row.names=T)


コンソール上でbを出力すると、1行1列目は空白としてきちんと表現されますが、テキストファイルに出力すると、1行1列目に、本来1行2列目の"Sepal.Length"がズレこんできます。

現状、テキストファイルの頭にスペースを追加して読み込んで対応しているのですが、気持ちが悪いのでコード上の問題があるならば解決したいです。

宜しくお願い致します。

↑

pairs()で対数軸にする方法†

はと (2011-05-29 (日) 16:21:03)

R2.13.0を利用しております。
pairs()で対数軸を利用する方法について質問があります。
plot()では、plot(data-frame,log="xy")で対数軸に出来るので、pairs(data-frame,log="xy")としたところ、一応対数軸になったものの、「"log"はグラフィックスパラメータではありません」とメッセージが出るとともに、図内のラベル（データフレームのheader）が表示されなくなってしまいます。
dataのlogをとることも考えましたが、相関係数は元のdataで計算したいのです。
よろしくお願い致します。

↑

subset関数でワイルドカードは使えますか†

たろ (2011-05-23 (月) 17:46:02)

初心者です。データフレームの中の変数 variable (数字と文字列の混ざった値です)が 20, 29, 2A, 2B の値のデータを選び出したいのですが、すべてのケースを列挙するのではなく、

subset(data, data["variable"] == "2*")

のような書き方で一遍に取り出すことは可能でしょうか。

↑

対応のあるデータのグラフ作成†

ようやくRが楽しくなってきたとこ (2011-05-19 (木) 18:49:52)

http://software.ssri.co.jp/statweb2/column/column0708.html に示されている、トレーニング前と後の違いを見るような、対応のあるデータを一覧できるグラフを作成したいのですが、Rではどうやればよいか教えていただければ幸いです。

↑

nnetに教師用データを追加する方法†

R始めて1ヶ月 (2011-05-07 (土) 12:28:53)

ニューラルネットワークのパッケージnnetについて質問です。
すでにデータAを学習したネットワークに対して、さらに別のデータBを学習させるにはどうすれば良いのでしょうか？先にAとBを結合して読ませるには時間がかかるために困っております。
別のパッケージならできるという情報もお待ちしております。

↑

eval(parse(text ="文字列"))のエラー†

質問者 (2011-05-06 (金) 11:07:52)

お世話になります。
fv1からfv25という変数にそれぞれread.csv()で数値データ（30KB~90KBなどまちまち）を代入しました。
これらのデータを順番に同じ処理をしたく以下のようなコードを書きました。

for (i in 1:25) {
	eval(parse(text = paste("fv <- fv", i)))
                      # 変数fvにi番目のfvを代入
	#処理
}

コードを実行すると

以下にエラー parse(text = paste("fv <- fv", i)) : 
   <text>:1:10:  予想外の 数値定数 です  
 1: fv <- fv 1 
             ^

というエラーが出てきてしまいます。
この原因が分かる方がいましたら、どうかご教授願います。

環境：
R version 2.12.1 (2010-12-16)
Platform: x86_64-pc-mingw32/x64 (64-bit)
↑

3次元空間の図形の色付け†

初心者KT (2011-05-04 (水) 09:19:54)

三次元の座標を持つ点AからEが存在します
下記の様な三角形が三つあります。
各三角形は連続して隣接しています。
下記の三角形の表と裏を別の色で塗ることは可能でしょうか。
またZの値によって、gradationがかかる様な塗り方も可能でしょうか。
実際には、100個以上の三角形の組み合わせによる構造物を扱う予定です。
三角形１
Point X (mm) Y (mm) Z (mm)
A -14.54 8.39 60.82
B -12.36 10.37 61.72
C -9.4 7.88 61.2

三角形２
A -14.54 8.39 60.82
D -15.77 5.74 60.3
E -10.51 6.07 61.2

三角形３
A -14.54 8.39 60.82
E -10.51 6.07 61.2
C -9.4 7.88 61.2

↑

パッケージが更新できない†

質問者 (2011-05-03 (火) 21:29:33)

パッケージを更新しようとすると以下のような警告出ます。

警告：  パッケージ 'urca' の前のインストールを取り除くことが出来ませんでした  
パッケージ 'zoo' は無事に開封され、MD5 サムもチェックされました 
ダウンロードされたパッケージは、以下にあります 
C:\Users\***\AppData\Local\Temp\RtmpVFYyw4\downloaded_packages 
install.packages(update[instlib == l, "Package"], l,
contriburl = contriburl,  中で警告がありました： 
 'lib = "C:/PROGRA~1/R/R-212~1.1/library"' は
 書き込み可能ではありません  
以下にエラー install.packages(update[instlib == l, "Package"], 
l, contriburl = contriburl,  : 
パッケージをインストール出来ませんでした


解決策はありますでしょうか。ちなみに作業状態ですが、Rを終了する時にワークスペースは保存しないを選んで終了し、Rを再起動後すぐにパッケージの更新を実行しました。

↑

パッケージのバージョン†

質問者 (2011-05-03 (火) 13:43:25)

あるパッケージをインストールした後、そのバージョンの確認が必要になったのですが、確認方法はありますでしょうか。

↑

ESSで日本語を表示する方法？†

Toy (2011-05-02 (月) 15:14:07)

お世話になります。
初歩的な質問で恐縮です。
ESSで日本語を表示する方法について教えてください。
まず.emacsに、次のように記述しました。

;;essのロード
(load "c:/Meadow/site-lisp/ess/Lisp/ess-site.el")

そしてess-site.elの299行目を、次のように書き換えました。

(setq-default inferior-R-program-name "c:/r/bin/i386/Rterm")

その上でMeadowを起動し、Alt+x rと入力するとRは起動します。
しかし、起動した後の初期画面には日本語は表示されず、作成されたプロットでも日本語は文字化けしてしまいます。
どのように設定すれば、日本語が正しく表示されるのか、御教示いただければ幸いです。
当方の環境は、以下のとおりです。

Editor:Meadow 3.01-dev
OS:Windows XP
R:Version 2-13-0
↑

Rコマンダーでクラスター分析以降の選択ができない†

R初心者 (2011-05-02 (月) 11:14:32)

R、統計ともに超のつく初心者です。
RおよびRコマンダーを用いて、クラスター分析を行い、樹形図を描きたいと思っております。
Rコマンダーのインストールはできたのですが、クラスター分析を選択した後、次の選択欄が全てグレイ表示になってしまい選択できません。
Rコマンダー入門によりますと、「現在の状況において利用できないものはグレイで表示され，選択できないようになっている．」との記載があり、現在、それらの項目が使用できない状況にあることはわかるのですが、他に何が足りないのかがわかりません。
クラスター解析をする際にインストールが必要なパッケージ等があるのでしょうか。
使用OSはWindowsXP, Rのバージョンは2.13.0,データはテキスト形式で読み込みは可能でした。
どうかご教授お願いいたします。

↑

数量化?�爐里笋衒�擇咼┘蕁爾砲弔い†

cyako (2011-05-01 (日) 18:44:20)

R初心者ですが、心理学論文を書くのにコレスポンデンス分析、もしくは数量化?�爐鰺僂い燭い塙佑┐討い泙后? そこで、解説サイトや本を参考にしつつ、下記�のスクリプトで実行しました。
(windows7、Rのバージョンは2.13.0、データはexcelで作成したカテゴリーデータです。)

&#65533;1つ目に使用したスクリプト
> df <- read.table("clipboard", header=TRUE)
> df
> library(MASS)
> corresp(df, nf=2)

このスクリプトで1つ目のデータはなんとか結果？を出すことができました。
しかし、clipboardの内容(別のカテゴリーデータ)を変えて、同じスクリプトで実行してみても下記�のようなエラーが出て結果が出ません。

&#65533;2つ目のエラー
> df <- read.table("clipboard", header=TRUE)
> df
> library(MASS)
> corresp(df, nf=2)
以下にエラー corresp.matrix(as.matrix(x), ...) :
empty row or column in table←これが表示がされます

このようなエラーに対してどのように対処すればいいのでしょうか？
もしくは、そもそもスクリプトが間違っているのでしょうか？
初歩的な質問で大変申し訳ないのですが、ご指導いただきたく存じます。

↑

ショートカットキーの変更†

hiro (2011-05-01 (日) 11:46:01)

Mac OS X 上で R 2.13.0 を使用しています。以前は標準エディタ上で ctrl+H で直前の1文字を消去できていたように思うのですが、このバージョンではヘルプが表示されてしまいます。ctrl+H の機能をもとへ戻す方法があれば教えてください。

↑

数量化I類における再計算について†

R初心者 (2011-04-30 (土) 18:30:19)

仮にRでDBからデータフレームに対して数量化I類を行って、stepを噛ましてAIC最小のモデルを選ぶとします。また、summaryで各カテゴリ変数に対しての有意な結果が分かっているとします。
この時点で、「有意じゃない」カテゴリ内の変数に関して除去しつつ、lmで再計算を行うにはどのようにすれば良いのでしょうか。オリジナルのデータ、あるいはデータフレームの内容を変更するしか手立てが無いのでしょうか。

↑

3次元配列を１次元の添え字順に二次元配列として表示する†

質問者 (2011-04-30 (土) 09:54:22)

Aployに次のような三次元配列を定義します。

	Apoly <- array(c(1,-0.5,0.3,  0,0.2,0.1,
                      0,-0.2,0.7,  1,0.5,-0.3),
                      c(3,2,2))

Aployを１次元の添え字ごとに2次元配列として表示するために次のようなコマンドを打ち込みました。

	Apoly[1,,]
	Apoly[2,,]
	Apoly[3,,]

結果は次の通りです。

> Apoly[1,,]
     [,1] [,2]
[1,]    1    0
[2,]    0    1
> Apoly[2,,]
     [,1] [,2]
[1,] -0.5 -0.2
[2,]  0.2  0.5
> Apoly[3,,]
     [,1] [,2]
[1,]  0.3  0.7
[2,]  0.1 -0.3

これらを一つのコマンドで実行する方法はないでしょうか。例えばApoly[1:3,,]を試しましたが、Apolyの表示方法と同じの3次元の添え字順に2次元配列として表示され、うまくいきません。

> Apoly[1:3,,]
, , 1

     [,1] [,2]
[1,]  1.0  0.0
[2,] -0.5  0.2
[3,]  0.3  0.1

, , 2

     [,1] [,2]
[1,]  0.0  1.0
[2,] -0.2  0.5
[3,]  0.7 -0.3
↑

数量化I類での基準変数の指定の仕方†

R初心者 (2011-04-30 (土) 05:53:05)

DBにRに接続しながらlm関数で数量化I類を試しています。
Rで数量化I類を行うと、キチンと回帰結果を返してくれるのですが、一方、基準となるカテゴリカルデータをこちらから指定したいのですが、その方法が分からないのです。
Interceptや、基準となるカテゴリー内のデータ名はlm関数の引数に与える事が可能なのでしょうか。

↑

Rcmdrのエラー†

miyucka (2011-04-30 (土) 00:25:29)

こんばんわ．いつもお世話になっています
環境はMac OSX　10.4.11です．
Rcmdrを使いたくて
library(Rcmdr)
と入力すると

Error in structure(.External("dotTclObjv",
objv, PACKAGE = "tcltk"), class = "tclObj") : 
[tcl] invalid command name "font".
Error :  .onAttach は 'attachNamespace' で
失敗しました 
エラー：  'Rcmdr' に対するパッケージ
もしくは名前空間のロードが失敗しました

となってしまいます．

パッケージのアップデートをCRANのバイナリもソースも両方やってみました
（違いは理解してませんがとりあえず．．．）
BioConductorもやってみましたが変化無しです．
おかしなことをしたのかと思ってもとからもっていたRを消して
もう一度インストールして
library(Rcmdr)と入力しても

要求されたパッケージ tcltk をロード中です 
Tcl/Tkインターフェースのロード中  終了済 
要求されたパッケージ car をロード中です 
要求されたパッケージ MASS をロード中です 
要求されたパッケージ nnet をロード中です 
要求されたパッケージ survival をロード中です 
要求されたパッケージ splines をロード中です

という記述はでるものの，結果は変わりません

他のQ&Aのページでは勝手に治ったという結果でバグだったのか？
という記述で終ってました．
これはどこをどうしたらいいのでしょうか？？？
ご指導頂ければ嬉しいです．よろしくお願い致します．

↑

効果的なグラフ作成について†

初心者 (2011-04-27 (水) 13:21:40)

立体的なグラフの表現方法について
初心者KT? (2011-01-31 (月) 22:24:19)

R初心者です。

自分は数学系、物理系、工学系の専門教育は受けたことがなく、医学系の仕事をしています。しかし今回幾つかの解析や、グラフによる表記を行う必要性に迫られています。

３０個の点が、xyzを与えられています。
１． ３０個の点が、“最も平均的に通る平面”（平面Aとする）を設定しました
（以前のこのコーナーで質問に答えていただき、平面Aは、ax+by+cz+d=0 と表現した時、optim関数でa = 0.09792275, b = 0.37648968 c = -0.32677544, d = 1となりました）。

２． 平面Aの上に８１個の点の様な構造物がドーム状にかぶさっています（曲面Bとします）。
� 30個の点と平面A
� 30個の点と平面A,および平面B
� �平面Aと平面B
のそれぞれの3種類の組み合わせを、3次元的にグラフで綺麗に表すには、どの様な方法があるでしょうか。
LatticeからCloundやwireframeなどを使おうとしましたが、なかなかうまくいきません。よろしくお願いします。
理想的には、平面Bは高さによって（Zの値によって）、色が変化するように表せたら最高です。平面Bに関する理想的な色の使い方参考例として、PDFで添付します。
３０個の点
X Y Z
56.8458 56.3972 95.2184
58.1262 55.2657 96.9028
61.0080 53.8846 97.6925
63.7117 52.5191 98.1147
65.9630 51.1603 98.0884
68.0598 49.8240 97.4626
70.3595 48.6390 96.2383
73.0480 47.8221 94.6007
76.0860 47.6060 92.8562
79.2539 48.1658 91.3380
82.2386 49.5663 90.3180
84.7249 51.7420 89.9494
86.4678 54.5090 90.2453
87.3370 57.6021 91.0912
87.3282 60.7262 92.2840
86.5478 63.6101 93.5838
85.1780 66.0516 94.7681
83.4311 67.9458 95.6757
81.5042 69.2908 96.2331
79.5429 70.1705 96.4590
77.6191 70.7182 96.4472
75.7289 71.0677 96.3317
73.8088 71.3037 96.2427
71.7678 71.4224 96.2643
69.5302 71.3165 96.4051
67.0781 70.7936 96.5927
64.4841 69.6327 96.6992
61.9234 67.6774 96.5962
59.6531 64.9476 96.2303
57.9525 61.7413 95.6913
曲面Bの構成要素
X Y Z
66.5214 52.3693 100.1890
69.0172 55.2805 98.6580
71.1465 57.6134 98.0792
73.4542 60.2416 96.9707
74.2816 61.4516 95.2830
75.8534 63.4528 93.5104
78.1356 65.5883 94.6485
79.3441 66.4813 96.3968
81.8057 68.7355 97.8613
69.4359 50.9051 99.5420
71.1025 54.5844 99.2688
71.9337 56.6633 98.3382
72.8589 58.7782 97.9511
73.5783 60.6862 96.7916
74.3058 62.7413 95.2100
74.8671 64.3807 93.8145
76.4086 67.2760 95.2155
78.2037 70.3165 97.9255
72.8598 48.6686 97.5556
73.0553 53.8495 97.0382
73.1433 56.6261 96.5826
73.3167 59.6391 96.9136
73.3231 61.7117 95.9461
73.3418 63.9996 95.0022
73.3420 65.9646 94.0229
73.5266 67.9228 94.9986
73.9282 71.3346 97.5466
76.6417 48.3653 95.4308
75.5219 52.4555 96.4391
74.9038 54.7740 97.3086
74.0257 57.8910 97.6356
73.1674 60.8503 97.5045
72.5119 63.0188 96.9356
71.6926 65.7293 96.2245
71.1617 67.6709 96.7142
70.0804 71.6452 97.8153
81.0857 49.4645 93.7371
77.7238 54.0901 94.3492
75.9176 56.5834 95.7828
74.4364 58.6303 97.2458
73.4941 59.9329 98.2465
72.1696 61.7559 98.5740
71.0152 63.3431 98.6402
69.0276 66.0708 98.0351
66.2486 69.8901 97.9546
85.2473 52.2150 93.0255
81.2928 54.9562 92.2683
79.3713 56.2273 93.0238
76.9566 57.8602 93.3170
75.3852 58.7892 95.9743
72.8658 60.3924 98.1361
70.5020 62.0146 97.9852
67.0597 64.3694 97.9051
62.4689 67.4921 98.1250
88.0380 56.8771 94.4088
84.4938 57.8743 92.7572
81.0642 58.7234 92.3931
78.9147 59.1500 93.2910
75.6953 59.7277 95.2876
72.8867 60.1105 98.3210
68.8648 61.0533 98.4587
65.0531 61.9484 98.5723
59.3881 63.4470 96.9480
87.7385 61.3387 95.9346
84.2589 61.3849 93.4009
80.7455 61.1783 92.9229
78.1388 60.8384 94.1021
74.5410 60.2960 96.3323
70.9528 59.7412 98.6698
67.4871 59.4730 98.7284
63.3083 59.2362 98.0871
58.8063 58.9720 97.4714
85.0653 65.4926 97.8556
83.9668 65.5278 94.3416
80.8889 64.4879 91.7567
77.4482 62.5331 93.9201
75.0270 61.1638 95.4027
73.0548 59.9064 97.5156
70.2444 58.5210 97.9352
66.1804 56.5065 98.6132
62.4675 54.5575 99.9249

↑

R2.12.2が64-bit版Windowsでないとアンインストールできないとは？†

NT (2011-04-23 (土) 16:50:39)

R2.12.2をアンインストールしようとしたら「このプログラムは64-bit版Windows上のみでアンインストール可能です。」という表示がでてきます。
使用OSはWindowsXP Proffesional Ver.2002 Servic Pack3です。
私としてはこれまでRがバージョンアップされる度に、古いバージョンと入れ替えてきており、問題は生じておりませんでした。

↑

PHP の decbin に相当する関数あるでしょうか？†

tadashi (2011-04-21 (木) 13:24:10)

PHP には、10進法を2進法に変換するdecbin とい関数があります。Rでは同様の関数はあるでしょうか？　つくるのは簡単で、方法もわかり、つくってもいますので、作り方を知りたいわけではないです。このライブララリを読むと、こんな関数があるよというのがごく当たり前にあるのでしたら、ぜひ教えてください。

↑

glmのpredic(…, se.fit)で出力される$se.fitについて†

まつよ (2011-04-20 (水) 12:36:01)

下のようなデータフレーム"data"がを用いて、glm関数による回帰分析を行った後に、predict関数のse.fit=Tを指定したときに出力される$se.fitについて教えてください。

> data
   x  y
1   1  2
2   2  3
3   3  5
4   4  9
5   5  4
6   6  4
7   7 10
8   8 12
9   9 15
10 10 18
11 11 19
12 12 19
13 13 20
14 14 22
15 15 25
16 16 29
17 17 30
18 18 35
19 19 36
20 20 35

> m=glm(y~x, data,
  family=gaussian(link=identity))
> predict(m, newdata=data.frame(x=x),
  se.fit=T)

$fit
        1          2          3          4          5          6 
-0.3428571  1.5458647  3.4345865  5.3233083  7.2120301  9.1007519 
        7          8          9         10         11         12 
10.9894737 12.8781955 14.7669173 16.6556391 18.5443609 20.4330827
       13         14         15         16         17         18 
22.3218045 24.2105263 26.0992481 27.9879699 29.8766917 31.7654135 
       19         20 
33.6541353 35.5428571 
$se.fit
       1         2         3         4         5         6 
0.9877626 0.9129465 0.8408734 0.7723116 0.7082816 0.6501237 
       7         8         9        10        11        12 
0.5995491 0.5586215 0.5295821 0.5144480 0.5144480 0.5295821 
      13        14        15        16        17        18 
0.5586215 0.5995491 0.6501237 0.7082816 0.7723116 0.8408734 
      19        20 
0.9129465 0.9877626 
$residual.scale
[1] 2.292081

おそら$se.fitの値は標準誤差（SE）なのだと思いますが、この値から標準誤差（SD）の値を求めることは可能なのでしょうか？例えば、rnorm()などの関数を使ってデータを再現し、ブートストラップ法でデータの信頼区間を求めるにはSEよりもSDがわかった方が便利な気がするのですが。

Googleで検索して少し調べてみるとSE=SD/sqrt(N)の関係があるそうですが、ここのときのNの値には何を使えばいいのでしょうか？N=nrow(data)でいいのでしょうか？

ちなみにsessionInfo()で出力される使用環境は以下の通りです。

> sessionInfo() 
R version 2.12.1 (2010-12-16)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Japanese_Japan.932 
    LC_CTYPE=Japanese_Japan.932   
[3] LC_MONETARY=Japanese_Japan.932
    LC_NUMERIC=C                  
[5] LC_TIME=Japanese_Japan.932    

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

loaded via a namespace (and not attached):
[1] MASS_7.3-9

se-fit.png

↑

Coxとforest plot†

michi (2011-04-18 (月) 19:03:14)

R 2.12.2を使用しています（mac）．
&ref(sorafenib SHARP trial NEJM.pdf)のfigure 3のような，Cox proportional hazard modelの結果をforest plotにする方法（パッケージ）をご存知でしょうか？

↑

お礼†

kazusata (2011-04-18 (月) 18:39:14)

ご回答ありがとうございます。
close.screen( ) をこまめに入れてみます。

↑

screen(n, new=FALSE)が効かない。†

kazusata (2011-04-17 (日) 21:52:27)

windows7で、R ver.2.12.0を便利に使わせて頂いています。
以下の使用で、思いどおりに行かず、ヘルプをお願い致します。
jpegデバイスに複数のグラフを書く中で、screen(n, new=FALSE)で、自由にscreenを移動したいのですが、以下の手順で、前に描いたグラフパーツが消えてしまい、new=FALSEが効かないのですが、何か他に設定が必要なのでしょうか？

jpeg(file.name, quality = 50)
split.screen(c(8,4) )
screen(1)
plot(x,y,type="n")
screen(2)
......
screen(1, new = FALSE)
liens(x,y) 
........
dev.off()

これで、jpegを見ると、座標軸が消えている。

↑

Kaplan-Meierのグラフについて†

hana (2011-04-13 (水) 17:21:10)

Mac OS X 10.6.7でR ver. 2.12.2を使用しています初心者です．
Kaplan-Meierの生存曲線はこちらの既出の方法で描けたのですが，たとえば再発率の解析などでy軸の0から曲線を描くにはどうすれば良いでしょうか（生存曲線の上下反転のような）？
またnumber at riskをグラフに盛り込むにはどうすればよいでしょうか？
よろしくご教示ください．

↑

R.12.2.2のインストール後、install.packages(　)ができなくなった†

sizu (2011-04-08 (金) 14:16:51)

よろしくお願いします。
R.12.2.2のインストール後、install.packages(　)で、下記に貼り付けたエラーメッセージが出て失敗するようになり困っています。

utils:::menuInstallPkgs()でも同じで、「以下にエラーメッセージを･･･」の後と同じメッセージがでます。

OSはWindowsXP (SP3) 、インストール方法は、full installを選択しました。
ミラーサイトは、Tsukuba, hyogo どちらも試しました。
また、私のPCでは、下記メッセージの「圧縮されたファイル･･･」以下に記載されている、「Rimp70CFbVb]の下には「downloaded_packages」というフォルダが入っており中身は空で、メッセージに書かれてあるよう名前のフォルダは見あたりません。

R2.12.1までは、普通に動いていました。バージョンアップしてからも、前に作成したワークスペース（.rdata）から起動した場合に限って、普通に動いていたこともありました。
現在では、R2.12.1、R2.12.2 いずれをインストールしなおしても、できませんし、以前保存したワークスペースから起動しても同じです。
インストール方法等については、Google などでいろいろ検索し、試行錯誤をしましたが、解決できず困っております。よろしくお願いします。

以下にエラーメッセージを貼り付けます。

> install.packages("Rcmdr")
 パッケージを ‘C:\Documents and Settings\
                Administrator\Local
                Settings\Application Data\
                R-core/R/win-library/2.12
   中にインストールします 
 (‘lib’ が指定されていないので) 
 以下にエラー gzfile(file, mode) :
 コネクションを開くことができません 
 追加情報:   警告メッセージ： 
In gzfile(file, mode) :
  圧縮されたファイル 'C:\DOCUME~1\ADMINI~1\LOCALS~1\
                       Temp\RtmpNubfVC/libloc_C
%3a%5cDocuments%20and%20Settings%5cAdministrator
%5cLocal%20Settings%5cApplication%20Data%5cR-core%2fR%2fwin-
library%2f2.12Version,Priority,Depends,Imports,
LinkingTo,Suggests,Enhances,OS_type,License,Archs,Built.rds'
を開くことができません, 
理由は 'No such file or directory' です
↑

クリップボードの読み込みについて†

miyuka (2011-04-07 (木) 22:14:08)

R初心者です．おねがいします．
環境はMac OS X 10.4.11．R version 2.10.1 です．
クリップボードからデータを読み込みたいのですが

x <- read.delim("clipboard")
以下にエラー file(file, "rt") :  コネクションを開くことができません :
追加情報:   警告メッセージ： 
In file(file, "rt") :
クリップボードを開くことができないか，中身がありません

となります．

ファイルのアクセス権の問題かと思いアクセス権をユーティリティから修復したのですがうまくいきません．
ディレクトリはファイルのあるところ，デスクトップ，ホームで試しましたが全て同じ結果です．
何がおかしいか教えて頂ければ幸いです．

↑

このグラフの呼び名†

K (2011-04-01 (金) 03:58:55)

こちらの論文（PDF）の図4のようなグラフはなんと呼ぶのでしょうか？　また、このグラフをRで作成するパッケージはありますでしょうか？
名前だけでも分かれば検索できるようになると思うので、どうぞよろしくお願いします。

↑

CSV fileのimport について†

初心者KT (2011-03-31 (木) 12:02:24)

初心者的質問です。
　CSV fileからデータをimportすると、時々データ行列の一番下に XXX levels: X, X, X, X, と表示が出ます（Xは数字）。このXXX levels・・・の表示を伴ったデータ行列になってしまった場合は、足し算や引き算などの行列の単純な計算を行う時に”因子に対しては無意味です“との表示が出て、計算が出来なくなってしまいます。
　またこの様な時は通常のplotも不可能となります。
　また行列の値を[ ] を用いて取り出しても、やはり値の下に、XXXlevels・・・の表示があります。
　CSV fileから、一部分のみコピーして、他のExcel fileにペーストしてからCSV形式に保存し、それをRにimportすると、XXX levelsの表示が消えて、行列計算が可能となる場合もあれば、そのままLevels・・・・表示が出て、単純な行列計算が不可能な場合もあります。
このXXX levelsの表示は何でしょうか？
また消す方法はあるのでしょうか？
宜しくお願いします。

下記にその表示をコピーします。

> AllL[2,3]
[1] 87.5242
47 Levels:  84.7952 86.1777 86.4596 86.4788
  86.5069 86.6655 ... Z
> AllL[3,3]
[1] 88.8609
47 Levels:  84.7952 86.1777 86.4596 86.4788
  86.5069 86.6655 ... Z
> AllL[2,3]- AllL[3,3]
[1] NA
 警告メッセージ： 
In Ops.factor(AllL[2, 3], AllL[3, 3]) :
 -  因子に対しては無意味です
↑

latticeのマルチプロットでパネルごとに注釈を付ける方法†

Toy (2011-03-29 (火) 16:08:34)

お世話になります。
latticeのマルチプロットで、gridにより各パネルごとに注釈を付けたいと考えております。
しかし、注釈をベクトルで値を変えながら付けていきたいのですが、上手くいきません。
どなたか、解決策を御存知の方、御教授いただければ幸いです。

library(lattice)
library(grid)
levels <- levels(as.factor(ethanol$C))
xyplot(NOx ~ E|C, data=ethanol,
	panel = function(x, y) {
	panel.xyplot(x, y)
	panel.loess(x, y)
       # 各パネルごとにCの水準を示したい。
	grid.text(label=paste("C=", levels, sep=""), 
	   x=unit(0.05, "npc"), y=unit(0.95, "npc"),
          just="left")})
↑

一般化加法モデル関数gamでの平滑化について†

hevel (2011-03-28 (月) 16:35:42)

gamを使って一般化加法モデル解析をしています。
gamの中で平滑化スプライン関数sを使って、“s(x, df=4, spar=0)”のように等価自由度を変更設定できると言うことなのですが、設定できません。
つまり、“df=某”という式を入れると、「ExtractVars 中のモデル式が不正です 」と言うメッセージが出ます。
バグでしょうか？それとも何か別のパッケージを入れておかなければならないのでしょうか？

↑

Silhouette plotの個別データラベル表示†

ken (2011-03-18 (金) 21:41:05)

はじめまして。
お世話になります。
私はpackage(cluster)内の関数pamを用い、クラスター分析を行っています。
たとえば、

> plot(pam(data, 3), ask=TRUE)

で、Silhouette plotやClusplotを作成した場合、グラフ上に個々のデータラベル（個体番号）を表示させるにはどうしたらよいのでしょうか。
上の例で、Silhouette plotやClusplot上で3クラスターに分けることができても、どの個体がどのクラスターに属しているのかが分からず困っています。
大変初歩的な質問で恐縮ですが，御教示いただければ幸いです．

↑

GLMのモデル選択について†

taipapa (2011-03-14 (月) 17:40:04)

お世話になります．
長文になりますが，お許しください．
ある疾患104人の治療効果の判定をロジスティック回帰モデルで行っています．
応答変数はOutcome は良と悪の2値で，説明変数は，9つのカテゴリー変数（a to i），6つの連続変数（J to O）からなり，eventであるOutcomeが悪は25例でした．
まず，下記の用にフルモデルを作成しました．

> FullModel <- glm(Outcome ~ a+b+c+d+e+f+g+h+i+J+K+L+M+N+O, family=binomial, data=mydata)

モデル選択はstepで行ったのですが，得られたモデルは，以下の通りです．

summary(stepModel)
Call:
glm(formula = Outcome ~ a+b+c+d+ J + K,
  family = binomial, data = mydata)
Deviance Residuals: 
    Min        1Q    Median        3Q       Max  
-1.53608  -0.65721  -0.36439  -0.09402   2.83204  
Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -6.7502     1.8614  -3.626 0.000287 ***
a                      0.8752     0.6204   1.411 0.158350    
b                     -1.2636     0.5959  -2.120 0.033974 *  
c                      1.4079     0.8331   1.690 0.091028 .  
d                      0.8834     0.6045   1.461 0.143914    
T1                    3.6828     0.9318   3.952 7.74e-05 ***
T2                    -0.8963     0.4643  -1.930 0.053555 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’
  0.05 ‘.’ 0.1 ‘ ’ 1 
(Dispersion parameter for binomial family taken to be 1)
   Null deviance: 114.717  on 103  degrees of freedom
Residual deviance:  81.622  on  97  degrees of freedom
AIC: 95.622

しかし，イベントが25しかないのに，説明変数が6つもあり，良いモデルとは言えず，含まれる変数も臨床的に意味のなさそうなものが入っており，困りました．
そこで，遺伝的アルゴリズムを用いてAIC/BICがもっとも小さいモデルを探してくれるglmulti (http://cran.r-project.org/web/packages/glmulti/index.html)と言うパッケージを用いてみました．
やり方やオプションは，付属のマニュアルにしたがいました．基準はデフォルトのAICcです．

GeneticModel <- glmulti(Outcome ~ a+b+c+d+e+f+g+h+i+J+K+L+M+N+O, family=binomial, data=mydata,
method="g", report = FALSE, marginality = TRUE, deltaB = 0, deltaM = 0.01, conseq = 6,
sexrate = 0.15, imm = 0.2, level = 1)

これで得られたモデルは4つの説明変数を有し，かつ，臨床的にも納得のいくものでした．
さらに，このモデルは，anovaを用いた尤度比検定で上記のstepで得られたモデルより有意に優れていました．
glmのsummaryは以下のとおりです．

> summary(GeneticModel.glm)
Call:
glm(formula = Outcome ~ a+b+J+K,
  family = binomial, data = mydata)
Deviance Residuals: 
   Min       1Q   Median       3Q      Max  
-1.6986  -0.5952  -0.4402  -0.1794   2.3542  
Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -5.3365     1.5652  -3.410 0.000651 ***
a                   -1.1344     0.5692  -1.993 0.046268 *  
b                    1.1870     0.7798   1.522 0.127945    
J                    3.4808     0.8931   3.897 9.72e-05 ***
K                   -0.8159     0.4454  -1.832 0.066993 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’
  0.05 ‘.’ 0.1 ‘ ’ 1 
(Dispersion parameter for binomial family taken to be 1)
   Null deviance: 114.717  on 103  degrees of freedom
Residual deviance:  85.879  on  99  degrees of freedom
AIC: 95.879
Number of Fisher Scoring iterations: 5

さらに，このモデルは，anovaを用いた尤度比検定で上記のstepで得られたモデルより有意に優れていました．
これで決まりかと思ったのですが，説明- pam(iris[,1:4],3)

変数bのp値が0.05をかなり超えてます（0.128）．~

またオッズ比を見ると，

exp(cbind(OR=coef(GeneticModel.glm),
   confint(GeneticModel.glm)))
Waiting for profiling to be done...
                    OR                2.5 %         97.5 %
(Intercept)    0.00481246   0.0001591826    0.08237885
a                   0.32160075   0.1006644837    0.96164214
b                   3.27732742   0.8119375358   18.53407469
J                  32.48671705   6.5308248197   224.16425104
K                   0.44222572   0.1700117573   0.98476474

やはり，説明変数bのオッズ比の95%信頼区間は0.81 to 18.5で1を跨いでいます．大部分は1より大ですが．

お聞きしたいのは，このようにinformation criteriaが低いという基準でモデルを選択したところ，一部の変数が有意でない場合どうしたら良いかということです．
基準がAICcなので，有意差検定は気にしないで，このモデルで解析を進めてよいものでしょうか?
primitiveな質問で恐縮ですが，ご教示いただければ幸いです．