新規投稿はできません。

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

過去の記事のアーカイブ



緯度経度から距離の換算

きむ (2011-01-07 (金) 18:12:19)

いつも参考にさせていただいています。
現在、緯度経度で座標が指定されている地域標準メッシュ(3次メッシュ)の座標をある地点を基準としたx,y座標に換算しようとしています。
そこで、RjpWikiの、「RでGPS」を参考にlibrary(sp)のspDistsN1という関数を用いて距離を出そうとしているですがうまくいきません。

library(sp)
ll <- matrix(c(34, 34, 138, 139), ncol=2)
km <- spDistsN1(ll, ll[1,], longlat=T)
sum(zapsmall(km))

とすると110.96kmとなりますが正確には92.38kmが正?また、

library(sp)
ll <- matrix(c(34, 35, 139, 139), ncol=2)
km <- spDistsN1(ll, ll[1,], longlat=T)
sum(zapsmall(km))

とすると結果が0になってしまいます。使用方法に誤りがあるのだと思いますがhelp等でも良くわかりません。すみませんが使用方法をご指導ください。
また、今のところ上記関数を用いて基準点からの距離で各メッシュの座標を設定しようとしていますが、直接緯度経度をx,y座標(平面直角座標系)に変換できる関数等がありましたら教えていただけると助かります。基本的な質問で恐縮ですがよろしくお願いいたします。

出力の一部をクリップボードにコピーして利用したい。

yamada (2011-01-07 (金) 13:46:05)

表のような結果をsink(),print()でテキストファイルに出力しています。実際にクリップボードに格納したいのは、一列目の行ラベルを除いた、jp,hk,usの列だけですが、テキストファイルで範囲選択をすると、行ラベルを含んでしまいます。現在は、エクセルでスペース区切りを指定して開いていますが、スペース区切り等開くまでの操作が手間です。そこで、行ラベルを出力しないか、あるいは、行ラベル付きのデータで行ラベル除いて出力する方法はないでしょうか。

               jp            uk            us
1    -3.030863e-02 -2.523841e-02 -4.387615e-03
2    -1.594305e-02 -1.853893e-02 -1.499784e-02
3     9.720749e-03  8.167201e-03  1.648989e-02
4    -1.138466e-02  5.133008e-03  1.188763e-03

read.csv で読み込んだデータの積集合をとりたい

tadashi (2011-01-05 (水) 15:08:35)

やりたいことは、1列しかないCSVファイルを2つ読み込んで、共通部分を出力したいのです。できるだけ簡単に行う方法はあるでしょうか?
read.cvs で読み込むと下記のようにうまくいきません。
CSVから、データの作り方がうまくないのだと思いますが、どのコマンドを使えばいいのかがわかりません。

> a <- read.csv('temp.1', header=FALSE)
> b <- read.csv('temp.2', header=FALSE)
> a
  V1
1  1
2  2
3  3
> b
  V1
1  3
2  4
3  5
> intersect(a, b)
data frame with 0 columns and 0 rows

パッケージVARの関数restrictの計算精度

松代 (2010-12-31 (金) 14:08:35)

グレンジャー因果性の検定量をrestrictで制約した残差と制約無し残差を使って計算すると、教科書の計算結果(Eview)やエクセルの計算結果と大幅な違いが出ました。原因は制約付残差の違いでした。restrictの計算精度は悪いのでしょうか。

初級Q&A アーカイブ(6)クラスター図を横向きの応用方法

hashi (2010-12-30 (木) 17:22:31)

いつもお世話になっております。
デンドログラムの回転方法ですが、青木先生著「Rによる統計解析」P230、図6.26を「初級Q&A アーカイブ(6)クラスター図を横向きで描画したい」の方法、horiz = T を前項のスクリプトのどの場所に挿入すれば回転したデンドログラムが表現できるのでしょうか。以下前項のスクリプトです。

set.seed(123)
x <- round(matrix(rnorm(100), ncol=5), 3)
d <- dist(x)
ans <- hclust(d^2, method="ward")
pdf("cluster.pdf", height=375/72, width=500/72)
plot(ans, hang=-1)
dev.off()

hang=-1の後ろに書き加え計算すると、

"horiz"はグラフィックスパラメータではありません

とメッセージが出ます。

mvpartパッケージで情報量エントロピーを使いたい

akira (2010-12-30 (木) 14:20:04)

いつもありがとうございます.
mvpartパッケージのrpart関数について質問です.
ヘルプを見ると、分岐指標に"gini"と"information"を選べるような記載がありますが、引数"parms"を変更しても結果が変わらないように思います.
同志社大の金先生のHPには、引数"split"とありますが、rpart関数、rpart.control関数は引数"split"を持たないようです.
一方、rpartのコード(40〜63行目ぐらいと思っていますが…)ではparmsが規定しているように見えます.
ご存知の方、いらっしゃいませんか?

$ operator

松代 (2010-12-28 (火) 15:54:34)

パッケージVARSを使うために解説書どおりに以下のコマンドを打ち込みましたが、最後コマンドの後に“以下にエラー Canada$e : $ operator is invalid for atomic vectors”という警告が出ました。対応策をご存じの方教えてください。

> library(vars)
> data(Canada)
> layout(matrix(1:4, nrow = 2, ncol = 2))
> plot.ts(Canada$e, main = "Employment", ylab = "", xlab = "")

Windows7+R64 2.12.1でパッケージadaptをインストールできません。

松代 (2010-12-25 (土) 00:07:03)

Windows7+R64 2.12.1でパッケージadaptをインストールできません。解決方法をご教示ください。

> install.packages("adapt")~
パッケージを ‘C:\Users\fujimot\Documents/R/win-library/2.12’ 中にインストールします 
(‘lib’ が指定されていないので) 
--- このセッションで使うために、CRANのミラーサイトを選んでください --- 
警告メッセージ: 
In getDependencies(pkgs, dependencies, available, lib) :
 package ‘adapt’ is not available

tapply関数について

ランゲル・ハンス (2010-12-21 (火) 09:00:38)

いつも掲示板を参考にさせていただいております。
以前も同じような質問をさせていただいたので、大変恐縮です。
NとMを10区間に区切って、その区間内に入るdの合計と平均を求めたいと思います。
下記の例(results1,results2)ではいくつかNAが出ます。NAを0として計算する方法を教えていただけないでしょうか?
最終的には区間にすべて色の入るlevelplotを作図したいと考えています。

N <- c(0, 1, 3, 7, 12, 20, 30, 45, 50, 81) 
M <- c(10, 20, 40, 60, 65, 75, 90, 95, 98, 100)
d <- runif(10, 0, 10)
data <- data.frame(cbind(N, M, d))

by1 <- cut(data$N, seq(0, 100, 10), right=TRUE, include.lowest=TRUE)
by2 <- cut(data$M, seq(0, 100, 10), right=TRUE, include.lowest=TRUE)

results1 <- cbind(tapply(data$d, list(by1, by2), sum))
results2 <- cbind(tapply(data$d, list(by1, by2), mean))
results1
results2

library(lattice)
levelplot(results1)
levelplot(results2)

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

特定の値が連続して出現する回数の最大値を計算

なつ (2010-12-20 (月) 22:13:28)

Rを使って,特定の値が連続して出現する回数の最大値を計算する方法を探しています。
例えば,特定の値=1としまして,

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

という2つのオブジェクトがあるとします。
この場合,aからは3,bからは5という値を計算したいということです。

何か良い方法がありましたら,教えて頂ければ幸いです。
初歩的な質問で申し訳ありませんが,よろしくお願いいたします。

少し高度なヒストグラム(か散布図)

yos (2010-12-19 (日) 20:39:01)

よろしくおねがいします.

1分おきで1日(0~1440分)に計測された,ある値があります.
1列目が,時間(0,1,2,...,1440)
2列目が,その時間における値(10, 40, ..., 50)
3列目が,グループ(A, B, C)
これを,
X軸を,例えば1時間おきの24分割
Y軸を,A, B, Cの3分割
で,値の平均値をプロットしたいのですが,色々調べてもどうしても方法が分かりません.
ご教授いただけないでしょうか

自由度調整済み決定相関係数の求め方

Kai (2010-12-12 (日) 15:31:30)

こんにちは。
lm()でデータAの一次回帰式を求め、predict()でデータBを回帰します。

A_res <- lm(formula=Volume ~ Girth, data=trees) # データAの回帰結果
B <- trees+rnorm(dim(trees)[1]*dim(trees)[2]) # データBの作成
B_res <- predict(A_res, B$Girth) # データBの回帰

B_resの回帰結果から自由度調整済み決定相関係数を求めたいのですが、どのようにすればよいのでしょうか?
ヘルプファイルを読む限りではse.fitという引数が関係ありそうな気がしたのですが、そこから先がよくわかりませんでした。

Rは2.11.1, OSはMac 10.5.8です。宜しくお願いします。

igraphのインストールエラー

kouka (2010-12-11 (土) 22:07:24)

R全くの初心者です。igraphのパッケージは下記のようにインストールされるのですが、

> utils:::menuInstallLocal()
パッケージ 'igraph' は無事に開封され、MD5 サムもチェックされました 

パッケージを読み込む際に、

> library(igraph)
 エラー:  パッケージ 'igraph' は 'arch=i386' に対してインストールされていません。

と表示されます。原因を教えていただけませんか?

rep.aovとSPSS

YK (2010-12-11 (土) 19:53:30)

はじめまして。
被験者間因子と被験者内因子を含む分散分析に関する質問です。
先日、中級Q&Aに投稿した内容ですが、初級Q&Aが適切ではないかというご指摘をいただき、こちらに再投稿します。
Rのaov、nlmeパッケージのlme、carパッケージのAnovaと岡田先生のrep.aovで同じ結果が得られるのですが、SPSSや桐木先生のWEBの結果と一致しません(SPSSと桐木先生の結果は一致します)。
下記にサンプルを示します。
P値はfactor1, factor2, timeの順に0.044, 0.441, 0.017となりますが、SPSSや桐木先生のWEBでは順に0.042, 0.397, 0.575となります。
ダミーの結果を加えて被験者内因子の水準を3にしてみても結果は一致しません。
他のサンプルデータでも結果が一致するときとしないときがあるようです。
SPSSはType IIIで計算していますが、R側でtype IIIを指定しても一致しません。
原因についてどなたかご教示いただけましたら幸いです。

ID <- factor(rep(c(1:32),2))
factor1 <- factor(c(0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 
1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 
0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 
0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0))
factor2 <- factor(c(0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 
0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 
0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 
0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0))
time <- factor(c(rep(1,32), rep(2,32)))
result <- c(102.56, 101.80, 123.72, 74.16, 34.12, 
61.52, 54.08, 115.08, 58.96, 60.24, 111.68, 122.20,
45.04, 104.84, 121.48, 125.20, 45.48, 217.00, 39.40,
157.48, 119.96, 65.80, 4.44, 39.40, 65.80, 90.20,
45.48, 117.72, 83.04, 156.40, 42.72, 110.92, 135.52,
22.60, 129.08, 127.76, 149.12, 69.92, 73.24, 97.12, 
114.92, 54.56, 170.08, 71.60, 179.96, 57.08, 199.88, 
126.48, 59.56, 120.40, 53.56, 127.12, 66.48, 75.24, 
268.64, 90.88, 102.48, 147.84, 116.84, 84.40, 147.84, 
138.76, 122.96, 120.04)
df <- data.frame(ID, factor1, factor2, time, result)

rep.aov(result ~ I(time) + B(factor1, factor2) + S(ID), df)

ネットワーク(sna)のグラフ描画(辺ごとの太さの調整)について

net (2010-12-10 (金) 19:24:00)

R初心者です.
ネットワーク分析をして,ネットワーク図を描くときに,辺の太さを辺ごとに調整できるとされています.太さを与えるデータはベクトルでも隣接行列でもよいという解説がありました.
たとえば,重み付きグラフの隣接行列(たとえばwg)でグラフを描き,その重みに応じて各辺の太さを変えたい場合,辺の太さを指定するedge.lwdをつかって具体的にどのように表現したらいいでしょうか.
ちなみに,wgは作業ディレクトリにあるcsvファイルで,
wg<-as.matrix(read.csv("XXXXX(ファイル名)".csv))
で既に定義してあるとします.
私は,
gplot(wg,
edge.lwd = wg)
として,グラフを描こうとしたら,エラーで「使われていない引数があります」として,edge.lwd = wgがでてしまいました.

具体的なスクリプトを教えていただけたら幸いです.

列数が不揃いなtxtファイルを読み込みたい

しょう (2010-12-10 (金) 02:48:02)

スペース区切りのテキストファイル(数百万行*9列)をRに読み込もうとしています。

ただ、元データの8列目が「備考」欄のため、scan(file="xxx.txt", sep="")を使うと備考のない部分は区切りのスペースと同じ扱いを受けてしまい、matrixを使って行列化するときにずれてしまいます。

read.tableを使うと「'289' 行目には,8 個の要素がありません」とエラーがでてしまい、読み込んでくれません。(288行目までは備考がなく、289行目で初登場のため)

read.delim("xxx.txt")を使うと、空白の備考欄に「NA」を入れてくれ、一番理想に近いのですが、1行目に備考がないため8列のフレームが作られてしまいます。

read.delimで列数を設定できればよいのですが、可能でしょうか?googleやRの書籍などで探したのですが見当たりませんでした。解決方法があればお教え願えると幸いです。よろしくお願いいたします。

R version=2.9.2, OS=Vista

R commanderのセットアップ

まさ (2010-12-08 (水) 17:45:28)

R commanderを動かそうとしている初心者です。MacOSX10.6、Rのバージョンは2.12.0です。

install.packages("Rcmdr", dependencies=TRUE)
の後に、

library(Rcmdr)
と入力すると以下のメッセージが来ました。

要求されたパッケージ tcltk をロード中です
Tcl/Tkインターフェースのロード中


となり、それから特に何も起こりません。 この後少し打ち込もうとすると、それ以上入力できず、終了せざるをえない状況です。良い解決策を教えていただけないでしょうか?またX11は立ち上がった状況でやってはいますが、解決できていません。よろしくご教授下さい。

plotの使い方

質問君 (2010-11-27 (土) 18:29:37)

plotの使いかたで質問です。

x <- c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
y <- c(5, 6, 7, 1, 3, 5, 6, 7, 1, 3)

これらのデータで、x[1]の時、y[1]をプロットし、x軸にそのままxの値を10〜100まで表示したいのですが、どうやればよいのでしょうか?
よろしくお願いします。

aggregate関数について

ランゲル・ハンス (2010-11-26 (金) 15:14:51)

aggregate関数について質問させていただきます。

以下のようなデータフレーム(data)があったとき、cut 関数で N を 10 区間に区切って、その区間における d の合計を求めたいと思います。ある区間において NA がある場合、NA を 0 として合計する方法を教えていただけないでしょうか?下記の例では (50, 60], (60, 70], (70, 80] の区間の合計を 0 としたいと思います。
どうぞよろしくお願いいたします。

N <- c(0, 1, 3, 7, 12, 20, 30, 35, 50, 91)
d <- runif(10, 0, 10)
data <- data.frame(cbind(N, d))
by1 <- cut(data$N, seq(0, 100, 10), right=TRUE, include.lowest=TRUE)
agg <- aggregate(data$d, by=list(by1), FUN="sum")
agg

改行、空行

超初心者 (2010-11-20 (土) 15:53:12)

空行の実行方法とプロンプト上で改行を行なう方法を教えてください!!
マニュアル等いろいろ探したんですが、全くわかりません。
初歩的すぎて申し訳ございませんが、よろしくお願いします。

日本語を含む固定長ファイルの読み込み

困り果ててます (2010-11-19 (金) 18:48:40)

日本語を含む固定長ファイルを読み込みたいのですが,うまくいきません。

123あ  456
321いう789

のような内容の test.txt というファイルがあり,1件目のデータをV1="123", V2="あ ", V3="456" と読みたいのですが,

read.fwf("test.txt", width=c(3, 4, 3))

だと意図したように読めません。

readChar("test.txt", c(3, 4, 3), useByte=TRUE)

としてみたのですが,1行目しか読めません。
解決法をご存知の方がいらっしゃいましたら是非 ご教示ください。
よろしくお願いします。

エクセルとの乱数の比較について

質問君 (2010-11-18 (木) 22:04:57)

エクセルでMM法を乱数で実行してるのですが、Rの乱数の方が良いと聞いたので、Rでも試したいのですが、エクセルとの擬似乱数比較について、Rが優れているなど分かりやすく、紹介されてるサイトや記事はないのでしょうか?たしかにエクセルよりも、乱数関数は多いです。性能的にはどうなのでしょうか?
探しても見当たりませんでした。もし、ご存知の方がいらっしゃれば、お願いします。

パッケージ ltm の grm 関数の出力結果の読み方

NT (2010-11-17 (水) 20:14:57)

3 段階評価 (0, 1, 2) の尺度の項目母数を grm 関数によって求めてみました。
出力結果の二行のみ示すと以下の通りとなります。

      Extrmt1	Extrmt2	Dscrmn
Item1  -1.719	 -0.126	 1.553

ここで Extrmt1 とは評価点 0 と 1 の項目曲線が交わる地点のθの値、Extrmt2 とは評価点 1 と 2 の項目曲線の交わる地点のθの値と解釈してよいのでしょうか。
ご教示願います。

ある座標と最も近い座標を特定する際の高速化の方法

Saito (2010-11-17 (水) 01:37:43)

いつもお世話になっております。
この手の話は既出かと思い検索しましたが、見つけることが出来なかったので質問させてください。

ある座標セットaを持っています。この座標セットには、座標とそれに付随する価が入っています。もう一つ別の座標セットbを持っています。こっちには、座標しか入っていません。今、座標bに最も近い座標aの持つ付随値を補間したいと考えています。下記にサンプルプログラムを示します。

> ###座標とそれに付属する価###
> a <- data.frame(
+         expand.grid(Long=seq(11, 13), 
+                      Lat=seq(20, 23)),
+         value=seq(1, 12)
+                 )
> ###マッチさせる数###
> num <- 1000
> 
> ###座標しか分かっていない###
> b <- data.frame(long=runif(num, 11, 13),
+                  lat=runif(num, 20, 23)
+                 )
> head(b)
      long      lat
1 12.00220 21.48183
2 11.14583 21.13535
3 12.31510 20.28334
4 11.54138 21.18872
5 11.66167 20.76173
6 12.52289 21.27551
> ###全ての組み合わせから最近接距離を持つ座標を特定###
> dis <- NULL
> for(i in 1 : nrow(b)) {
+ dis <- sqrt((b[i, 1] - a[, 1])^2 + (b[i, 2] - a[, 2])^2)
+ b[i, 3] <- a[which(dis==min(dis))[1], 3] #i番目のbがどのaに最も近いか探し出 してマッチング
+ }
> head(b)
      long      lat V3
1 12.00220 21.48183  5
2 11.14583 21.13535  4
3 12.31510 20.28334  2
4 11.54138 21.18872  5
5 11.66167 20.76173  5
6 12.52289 21.27551  6


しかし、実際にはaもbも超巨大で、この方法だとものすごく時間がかかることがあります。上記プログラムをもっと高速化するにはどうすればよいでしょうか。どなたかわかる方がいらっしゃいましたら、ご助言頂けると幸いです。

大量の従属変数について共分散分析を行う方法について

カワウソ (2010-11-15 (月) 19:11:55)

はじめまして。
n個の従属変数(y1,y2,…,yn)と一つの独立変数xに関して,if構文を使って以下のような作業を一度に行うにはどのようにプログラムを書いたらよいのでしょうか?

lm(y1 ~ x + g)
lm(y2 ~ x + g)
・
・
・
lm(yn ~ x + g)

データセットの削除方法

ダルビ (2010-11-10 (水) 21:24:20)

初めまして。
Rコマンダーで特定のデータセットを削除するにはどうすればよいのでしょうか?
いろいろなデータセットができてしまって困っています。
お手数ですが、よろしくお願います。

kmeans と pamのクラスタリング結果について

Jam (2010-11-08 (月) 17:07:02)

初めまして。
clusterライブラリで利用できるk-meansとpamの挙動について質問です。

2次元のxy座標で示せるデータをk-meansとpamでクラスタリングをしていて気になることがありました。それはpamの実行結果が毎回同じであるということです。
私の解釈として、どちらの手法も初期のランダムサンプルに結果が依存するため、クラスタリング結果は同じデータ・設定で実行したとしても毎回変わると思っています。確かにk-meansでは実行の度に結果は変わることは確かに確認できるのですが、pamではクラスタリング結果が毎回同じです。
クラスタリング結果が同じというのは、plotした図、result$centers (k-meansの場合)、result$medoids (pamの場合)を見て判断しています。データ数は1000弱です。

これはどういうことなのでしょうか?pamのランダムサンプルはクラスタリング結果に依存する、つまり結果は毎回多少でも異なってくると思うのです。k-menansの$cnetersは実行のたびに異なるが、pamの$medoidsは毎回同じ値になることは以下のプログラムで確認しました。

検証プログラム。

library(cluster)
# (x, y)形式のデータを読み込む
data <- read.table(var.in_file, header=F, sep="\t") 

################
# pam clustering          #
################
data.clust <- pam(data, 10)
data.clust$medoids

################
# k-means clustering     #
################
data.clust <- kmeans(data, 10)
data.clust$centers

LOESS平滑化曲線と信頼区間

MKI (2010-11-06 (土) 05:26:32)

散布図にLOESS()で平滑化曲線を加えたあと、その信頼区間を合わせてグラフに描画したいのですがそれらしい関数が見あたらず、もしどなたかご存じでしたら教えて下さい。できれば信頼区間を95%以外にも自分で変更できれば助かります。

リストの成分

nan (2010-10-31 (日) 18:35:42)

初めまして。

リストの成分にアクセスしたいのですが、例えば以下のようなデータについて

> test <- function(x) {y=log(x); z=sin(x); return(list(value=x, log=y, sin=z))}
> test(1:3)
$value
[1] 1 2 3
$log
[1] 0.0000000 0.6931472 1.0986123
$sin
[1] 0.8414710 0.9092974 0.1411200

$logの[0.0000000]や$sinの[0.9092974]にアクセスするにはどのようにすれば良いのでしょうか?

宜しくお願いします。

空の値を欠損値と指定する

のの (2010-10-31 (日) 04:07:34)

データフレームxの中の列名(変数名)aについて

> attach(x)
> x <- ifelse(a=="", NA, a)

とすると、空の値を欠損値NAと指定できます。
しかし、aがfactorだった場合には、factorの性質が失われてしまうようです。
何かもっとエレガントな方法はあるのでしょうか?

> x <- subset(x, a != "")
> table(a, useNA="always")

これも違うようです。

ブートストラップ法によるreliability

MKI (2010-10-29 (金) 22:29:49)

多変量解析を行い、意味のありそうな説明変数を取り出す作業をした後、それをブートストラップ法によってreliabilityを確認するよう求められました。Rでやってみたいと思いますが、「entry criteriaをp<=0.10にしてretention criteriaをp<0.05で確認する」ように指示されましたがstepAICでできるのでしょうか。方法としては
1.スプレッドシートからランダムにサンプリング(重複あり)
2.一般化線形モデルに当てはめ
3.stepAIC???で残った説明変数を抽出
4.1〜3を1000回繰り返す
といった感じになるかと思うのですが3.の部分が方法として正しいのかどうか分かりません。あるいは同じようなこと(多変量解析+ブートストラップ法)ができるパッケージはありますか。もしご存じでしたら教えて下さい。

abline(h = 5) との交点の X の値

hashi (2010-10-28 (木) 05:34:07)

Q&A (初級者コース)/11であった「任意の X 軸と曲線との交点(Y 値)の値の算出」の逆のことをしたいのです。
以下がデータです。

x <- c(2:10)
y <- c(0.084, 3.642, 12.472, 27.262, 47.035, 67.439, 83.974, 92.767, 100)
plot(x, y, xlim = c(2, 10), ylim = c(0, 100))
lines(spline(x, y, n = 20), col = 2)

このとき、abline(h = 5) を作図したときスプライン曲線との交点を読みたいのですが、どのようにしたらよいのでしょうか。
どうかよろしくご教授お願いいたします。

各列毎にある条件にあう要素の数をカウント

MKI (2010-10-27 (水) 13:17:49)

データフレームの各列毎にある条件(たとえば0.1以上など)を満たす要素をカウントするにはどうしたらよろしいでしょうか。apply関数で各列のmeanなどは簡単に得られますが、似たような方法でカウント数を得ることは可能でしょうか。

4つの集合のベン図

sh (2010-10-27 (水) 09:17:04)

4つの集合のベン図をRで描くにはどうすればよろしいでしょうか? [参考文献] http://ja.wikipedia.org/wiki/ベン図 多数の集合のベン図 http://www.ats.ucla.edu/stat/r/faq/venn.htm R FAQ: How can I generate a Venn diagram in R? vennDiagram Can't plot Venn diagram for more than 3 sets

行列操作について

ランゲル・ハンス (2010-10-26 (火) 10:29:21)

いつも掲示板を参考にさせていただいています。
さて、下記の行列mの各行の値を逆順にして、行列m2を作りたいと思います。
行列操作でmからm2を作る方法、あるいはdからm2を作る方法をご教示いただけないでしょうか?

d <- c(1, 2, 3, 4, 5, 6, 7, 8, 9)
m <- matrix(d, 3, 3, byrow=TRUE)
d2 <- c(3, 2, 1, 6, 5, 4, 9, 8, 7)
m2 <- matrix(d2, 3, 3, byrow=TRUE)

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

ハフ変換

yorudan (2010-10-25 (月) 17:22:10)

ハフ変換を行うためのパッケージはありますか!?

library(gregmisc)のバグ!?

のの (2010-10-24 (日) 03:42:03)

library(gregmisc)を導入して、excelシートからデータを直接取り込む方法を試しています。
いくつかおかしな挙動をみつけたので書き込みます。
こちらの環境は、MacOS10.6.4, GUI版R version 2.11.1 (2010-05-31) x86_64-apple-darwin9.8.0 です。
他の環境で、同じ状況が再現するのかも教えて欲しいです。

1. エクセルシートの中で、"'’"または"…"を使っていると次のエラーが出る。

Wide character in print at /Library/Frameworks/R.framework/Resources/library/gdata/perl/xls2csv.pl line 262.

調べた範囲では、たぶん内部処理でperlに文字列を渡すときに文字コードエラーを起こしているのだと思います。ただデータは正しく取り込まれるので実質的な問題はないように思います。
これがfileサンプルファイルです。
2. サンプル数が増えてしまう
こちらの方が問題が大きいです。Rはエラー表示しませんが、変なレコードが追加されて、データフレームに影響を与えます。
試行錯誤して、""',"この文字列を含むとおかしくなることは再現できました。
これがfileサンプルファイルです。
まだ、library(gregmisc)開発サイトの方は良く調べていないのでこれから調べて何か分かればこの記事を更新する予定です。

図の中に+/-

MKI (2010-10-22 (金) 13:29:37)

図(グラフ)の中に平均(変数)+/-標準偏差(変数)を書き込みたいのですがexpressionを使用してもうまくいきません。ひとつひとつ数字を書き込めば
可能ですが、変数などと組み合わせて表示する方法はないでしょうか。

scatterplot3dでの三次元プロットのサイズ設定について

しげゆき (2010-10-18 (月) 10:39:18)

XPでR2.11.1の初心者です。よろしくお願いします。

scatterplot3dでプロットの色、サイズを変えた図を作成していますが、色はデータどおりいきましたが、サイズがデータと食い違ってうまくいきませんでした。

以下データで左3列がXYZ軸、colorが色4種類、sizeが5段階です。

  sisu dens length color size
1 -115  1.4    9.8     A    5
2  -62  1.6    7.0     A    3
3  108  1.7  452.0     B    2
4  174  1.8  202.0     B    1
5  184  1.7  141.0     B    4
6   11  1.7   10.2     C    4
7   23  1.7   40.0     C    5
8  -14  1.8  232.0     D    5
9  127  1.7  132.0     D    3
scatterplot3d(x=d$sisu, y=d$dens, z=d$length,
color=c(2,3,4,5)[unclass(d$color)],
cex.symbols=c(1,2,3,4,5)[unclass(d$size)])

cexを書き出しましたが、データどおりでおかしくありませんでした。

> c(1,2,3,4,5)[unclass(d$size)]
[1] 5 3 2 1 4 4 5 5 3

どこか間違っているのでしょうか。

ベクトルの個数を返す関数の自作

のの (2010-10-14 (木) 11:05:44)

いつもお世話になります。 例えば、次のようなベクトルが二つあるとき、

x <- c(1:10)
y <- c(20:40)
test(x)とtest(y)の結果が1でtest(c(x,y))の結果が2となるような関数testを作りたいのですが、何か良い方法があれば是非教えて下さい。
test <- function(i) {
  text <- as.character(match.call()[2])
  if (grep("^c", text) != 1) {1}
  else {grep(",", text)}
}

のようになるのかと考えましたが上手くいきません。

周囲の値を用いた中心座標の平均値の算出

Saito (2010-10-13 (水) 13:00:02)

いつもお世話になっています。
過去ログ等検索しましたが、見つからなかったので質問させてください。
座標と座標に付与されたデータからなる3列のデータセットがあります。今、座標が細かすぎるので、もう少し粗い解像度に変換したいのですが、その変換プログラムが上手くいきません。以下に例を示します。

> ###仮想データセット###
> ###x, yが座標で、zが値###
> mat <- data.frame(expand.grid(x=1:6, y=1:4), z=1:24)
> mat
   x y  z
1  1 1  1
2  2 1  2
3  3 1  3
4  4 1  4
5  5 1  5
6  6 1  6
7  1 2  7
8  2 2  8
9  3 2  9
10 4 2 10
11 5 2 11
12 6 2 12
13 1 3 13
14 2 3 14
15 3 3 15
16 4 3 16
17 5 3 17
18 6 3 18
19 1 4 19
20 2 4 20
21 3 4 21
22 4 4 22
23 5 4 23
24 6 4 24
> ###y軸方向に足し合わせたときの行列を用意###
> mat2 <- data.frame(matrix(0, ncol=ncol(mat), nrow=nrow(mat)/2))
> 
> ###x軸方向にも足し合わせたときの行列を用意###
> mat3 <- data.frame(matrix(0, ncol=ncol(mat), nrow=nrow(mat2)/2))
> 
> ###y軸の値がx軸方向にいくつ詰まってるのか確認###
> a <- as.numeric(summary(as.factor(mat[, 2]))[1])
> 
> ###y軸方向に何度足せばよいのか確認###
> b <- length(levels(as.factor(mat[, 2])))
> 
> ###等比数列でy軸の値をx軸方向に詰まっている分だけ足し合わせ###
> for (i in 1 : (b/2)) {
+ mat2[((i-1)*a+1):(i*a), ] <- 
+     mat[((2*a)*(i-1)+1):((2*a)*(i-1)+a), ] + 
+     mat[((2*a)*(i-1)+(1+a)):((2*a)*(i-1)+(a*2)), ]
+ }
> mat2
   X1 X2 X3 
1   2  3  8
2   4  3 10
3   6  3 12
4   8  3 14
5  10  3 16
6  12  3 18
7   2  7 32
8   4  7 34
9   6  7 36
10  8  7 38
11 10  7 40
12 12  7 42
> ###同様にx軸についても###
> for(i in 1 : (nrow(mat2)/2)){
+ mat3[i, ] <- mat2[(2*(i-1)+1), ] + mat2[(2*i), ]
+ }
> mat3
  X1 X2 X3
1  6  6 18
2 14  6 26
3 22  6 34
4  6 14 66
5 14 14 74
6 22 14 82
> ###最終的に得たい座標セット###
> ###周囲4つの値の平均となってほしい###
> mat3/4
   X1  X2   X3
1 1.5 1.5  4.5
2 3.5 1.5  6.5
3 5.5 1.5  8.5
4 1.5 3.5 16.5
5 3.5 3.5 18.5
6 5.5 3.5 20.5
> 

一応、上記のへたくそなプログラムでもおおよそのやりたいことは出来ているのですが、実は、まだ上記に含めていない条件が二つあります。
一つは、周囲4つの平均でなく周囲9つの平均や16の平均(つまり正方形で扱いたい)と任意に変えたいときに、上記のプログラムでは一から書き直しになります。
そこを明示的に加えたいのですが、上手く書けませんでした(例ではそもそも24行しかないので、9や16は難しいですが…)。

もう一つの条件は、z軸にNAが含まれている座標があるのですが、そのときはNAを抜いた数で平均を返したいのです。
例えば、

mat[1, 3] <- NA

等とした時には、そのまま上記プログラムを走らせると最後までNA表記となり、やりたいことができません。
4で割るのではなく、そのときだけ上記例ですと3で割ってほしいのですが、複雑で出来ませんでした。
さらに、このデータ、実際は結構大きな行列(6万×3くらいです)ですので出来るだけベクトルで処理したいと思っています。
高速化した例ですと、大変助かります。

どなたか上記の条件で変換が出来る方がいらっしゃいましたら、ご教授頂けると幸いです。

ポリコリック相関係数行列の算出について

JJ (2010-10-12 (火) 20:17:51)

まったくの初級者です。
ポリコリック相関係数を用いた因子分析(最尤法、プロマックス)を3日間ほど悩んでいます。最終的に因子得点も算出したいと思っています。
現在までのところ、
factanal(covmat=cor(データ), factors=3,promax,regression)

という感じまでは来たのですが。。。数字は出ますが、これは一般の因子分析なのかと思っています。

ちなみにpolychor(x,y,)関数の場合、相関係数にはなので、相関係数行列にはなりませんよね。

まことに見当はずれかもしれませんが、誰か力を貸していただければありがたいです。
上記の corのままでは

library(psych)
library(polycor)
data(bfi)
dat <- bfi[1:17] # 17項目のデータ
pcr <- polychoric(dat)
pcr # polychoric関数はrhoやtauといった複数の値を返す。polychoric関数のヘルプ参照
ans <- factanal(covmat=pcr$rho, factors=2, rotation="promax")
ans

続いている数字の数を出力する方法

sayaka (2010-10-12 (火) 14:05:00)

以前、「T,Fを文字列として出力するには」でお世話になりました。
うまく説明するのが難しいのですが、下記の data のようなベクトルがあるとき、399,400,401・・のように数値が続いている場合に"何個続いているか"を下記の a のように出力したいと考えているのですが、どのようにすれば求めることが出来るでしょうか。
いろいろと試行錯誤したのですが、どうしても思いつきません。
data はあるデータ中のエラーデータの行番号を抽出したもので、エラーデータの開始位置と長さより削除するプログラムがあるために長さが必要です。
開始位置は下記のposiで求めました。
どうぞよろしくお願いします。

> data
 [1]  14  41  73 152 296 297 399 400 401 402 418 419 420 421 422 423 424 425 426 427 428 429
[23] 430 431 451 452 453 454 460 461 466 500 501 502 503 504 505 506 507 508
> posi <- data[diff(c(0,data)) != 1]
> posi
 [1]  14  41  73 152 296 399 418 451 460 466 500
> a
 [1]  1  1  1  1  2  4 14  4  2  1  9

RとWindows 7の相性

菊亭 (2010-10-06 (水) 22:12:34)

現在XPでRを使っております。そのうち7にアップグレードしょうかと画策中ですが、VistaとRの相性がかなり悪かった(というかVistaが酷かった)記憶があります。7とRには既知の問題はありますでしょうか?

plotirxのcolor2D.matplotの色表示について

ランゲル・ハンス (2010-10-02 (土) 09:51:51)

plotrixにあるcolor2D.matplot関数について質問させていただきます。

library(plotrix)
x <- c(1, 0, 10, 5, 3, 6, 4, 5, 4, 3, 8, 5, 1, 0, 7, 5)
data <- matrix(x, ncol=4)
color2D.matplot(data, c(1, 0), c(0, 1), c(0, 1), show.values=TRUE)

このmatirixの数値をカラー表示のマトリックスに変換して数値も表示したいと思います。
上記の例では変換できるのですが、マイナスがある場合のスケール変換の方法と色表示の方法をご教示いただけないでしょうか?
例えばyの場合について

y <- c(1, 0, 10, 5, -3, 6, -4, 5, 4, 3, 8, 5, 1, 0, 7, 5)

よろしくお願いします。

psych パッケージの factor.pa の結果表示

隣は何をする人ぞ (2010-10-01 (金) 21:34:02)

psych パッケージを使って因子分析を行ったとき、因子負荷量の大きい順にソートしたときとしなかったときで、因子負荷量の値が異なったものが表示されます。単に因子単位に因子負荷量の絶対値の大きい順に並べ替えて表示するだけだと思うのですが、なぜ因子負荷量の値が異なるのでしょうか。訳がわかりません。

library(psych)
data(bfi)
fa.parallel(bfi)
fpa.out <- factor.pa(bfi, nfactors=3, rotate="promax")
print(fpa.out) # 分析に使用した変数順そのまま
print(fpa.out, sort=TRUE) # 因子負荷量の大きい順に並べ替える

で分析を行いました。並び順は異なっても、各変数の因子負荷量は同じはずです。
まず、print(fpa.out) の結果です。いくつかの変数だけを選択して掲載します。

########## デフォルト(sort=FALSE)のとき
> print(fpa.out)
Factor Analysis using method =  pa
Call: factor.pa(r = bfi, nfactors = 3, rotate = "promax")
Unstandardized loadings based upon covariance matrix
            PA1   PA2   PA3     h2   u2     H2   U2
A1        -0.22  0.08  0.01 0.0600 0.94 0.0600 0.94
  :
E4         0.74 -0.06 -0.15 0.4823 0.52 0.4812 0.52
  :
gender     0.20  0.15  0.00 0.0490 0.95 0.0490 0.95
education -0.03 -0.04  0.10 0.0099 0.99 0.0099 0.99
age        0.04 -0.09  0.14 0.0398 0.96 0.0398 0.96

                PA1  PA2  PA3
SS loadings    3.30 2.66 2.15
Proportion Var 0.12 0.09 0.08
Cumulative Var 0.12 0.21 0.29

 Standardized loadings
          item   PA1   PA2   PA3     h2   u2
A1           1 -0.22  0.08  0.01 0.0600 0.94
  :
E4          14  0.73 -0.06 -0.15 0.4812 0.52
  :
gender      26  0.20  0.15  0.00 0.0490 0.95
education   27 -0.03 -0.04  0.10 0.0099 0.99
age         28  0.04 -0.09  0.14 0.0398 0.96

次に、print(fpa.out, sort=TRUE) の結果です。

########## sort=TRUE を指定したとき
> print(fpa.out, sort=TRUE)
Factor Analysis using method =  pa
Call: factor.pa(r = bfi, nfactors = 3, rotate = "promax")
Unstandardized loadings based upon covariance matrix
            PA1   PA2   PA3     h2   u2    H2   U2
E4         0.74 -0.06 -0.15 0.4823 0.94 0.339 0.66
  :
A1        -0.22  0.08  0.01 0.0600 0.73 0.076 0.92
gender     0.20  0.15  0.00 0.0490 0.58 0.078 0.92
  :
age        0.04 -0.09  0.14 0.0398 0.99 0.039 0.96
education -0.03 -0.04  0.10 0.0099 0.96 0.010 0.99

                PA1  PA2  PA3
SS loadings    3.30 2.66 2.15
Proportion Var 0.12 0.09 0.08
Cumulative Var 0.12 0.21 0.29

 Standardized loadings
          item   PA1   PA2   PA3    h2   u2
E4          14  0.62 -0.05 -0.12 0.339 0.66
  :
A1           1 -0.25  0.09  0.01 0.076 0.92
gender      26  0.25  0.19  0.00 0.078 0.92
  :
age         28  0.04 -0.09  0.13 0.039 0.96
education   27 -0.03 -0.04  0.10 0.010 0.99

Unstandardized loadings based upon covariance matrix については、並べ替えしてもしなくても、各変数の因子負荷量は同じです(これが当たり前だと思います)。
Standardized loadings については、並べ替えしたのとしないとでまるっきり違うものが表示されています。
なぜでしょう。
一応、psych クラスの print メソッド(psych:::print.psych) のソースもたどっては見たのですけど、明らかなバグというのではなく、書かれているプログラムがなぜそのようなパスをたどらなければならないのかがよくわかりませんでした(それ自身がバグと言うことなのかも知れませんが)。
多くの人が使っているパッケージなので、いまだにバグがあるとも思えませんが、不思議に思いましたので質問させて頂きます。

Rの実行を中断する方法

のの (2010-09-29 (水) 02:16:00)

MacでGUI版Rを使っています。version.string R version 2.11.1 (2010-05-31)
エディタにコマンドを連ねて、メニュー>編集>実行をする場合、途中でエラーがあっても、止まらずに最後まで流れてしまいます。
コマンドリストの任意の場所でRの実行を止める方法を探しています。
readline(),stopifnot()などを試しましたがだめでした。
今のところquit()を入れると、そこで、保存するかどうか聞いてくるので目的は達成できているのですが、何かもっと良い方法はないでしょうか?

パッケージのロードは起動後毎回必要ですか

shumei (2010-09-25 (土) 14:33:13)

Mac版R 2.11.1 をインストールした後、パッケージurcaをインストールしました。
urcaをロードし、パッケージの使用はできましたが、Rを起動するたびに未ロード状態にもどってしまいます。
何冊か書籍に当たってみましたが、該当する記述が見つからず、こちらに投稿させて頂きます。
パッケージのロード状態を維持するにはどうすればいいのでしょうか。それとも、これはRの仕様なのですか?
ご教授下さい。どうぞよろしくお願い致します。

ヒストグラムの階級分割

moyu (2010-09-21 (火) 14:54:22)

ヒストグラムで、ある一部の区間において細かく階級を分割したいのですが、どうしたらいいのか分かりません。ちなみに以下の操作をしました。

hangseng = read.csv("ASIA/hangseng/hangsengdaily-1986dec31~2010jun11.csv")
hangseng.ts = ts(rev(hangseng$Close),start=c(1986,12,31),frequency=248)
ts.plot(hangseng.ts)
ts.plot(diff(hangseng.ts))
hist(diff(hangseng.ts),breaks=30,col="magenta")

histの中のbreaksをどういじったらいいのでしょうか??

embedFontsで「・」が表示される

初級者です。 (2010-09-20 (月) 18:11:38)

CentOS 5.5 に Rをインストールし、以下の操作をしました。

pdf("test.pdf",family="Japan1")
plot(1:10,ylab="test test")
dev.off()
embedFonts("test.pdf")

すると、ylabに指定した"test test"文字列の間にある空白の位置に、「・」が表示されてしまいます。
正確には、"test ・test"という感じで、普通に半角スペースが表示されている上に、若干右側によって「・」が重なっているという感じです。
embedFonts()を呼ぶ前の状態では、「・」は表示されていません。
また、その状態で文章のプロパティから使用されているフォントを見ると、

「KozMinPro-Regular-Acro」
「KozMinPro-Regular-Acro.Bold」
「KozMinPro-Regular-Acro.BoldItalic」
「KozMinPro-Regular-Acro.Italic」
「Symbol」
「ZapfDingbats」

というフォントが使用されていると出ますが、embedFonts()を呼んだ後で同様にフォントを確認すると、「Sazanami-Gothic(埋め込みサブセット)」というフォントが使用されていると出ます。
何らかフォントの設定が足りないとは思い、いろいろ調べてみてはいるのですが、壁に当たった状態です。
同様の現象を解決された方がいらっしゃれば方法をご教示いただけないかと思い、こちらに投稿しました。
よろしくお願いいたします。

ヘルプファイルを表示エラー

てるてるぼうず (2010-09-19 (日) 23:22:48)

関数のヘルプファイルを表示したいと思い、以下を実行するのですが、
?…またはhelp(…)毎回以下のようなエラーが出ます

> ?glm
 警告メッセージ: 
In file.show(temp, title = gettextf("R Help on '%s'", topic), delete.file = TRUE) :
   file.show():ファイル 'C:\DOCUME~1\蜿、蟾晏忽蠢予LOCALS~1\Temp\RtmpIDu7yb\Rtxt678418be' は存在しません 
> help("glm")
 警告メッセージ: 
In file.show(temp, title = gettextf("R Help on '%s'", topic), delete.file = TRUE) :
   file.show():ファイル 'C:\DOCUME~1\蜿、蟾晏忽蠢予LOCALS~1\Temp\RtmpIDu7yb\Rtxt3d6c4ae1' は存在しません 

検索エンジン等で調べてはみましたが、未だ解決していません。
環境は以下の通りです。

Microsoft Windows XP Professional
Version 2002
Service Pack 3
Intel(R) Core(TM)2 Duo CPU
E8400 @ 3.00GHz
2.99 GHz、976 MB RAM

Rのバージョン
R version 2.11.1 (2010-05-31)

regexprのpatternの文字数によるエラー

sh (2010-09-16 (木) 11:50:21)

regexpr(pattern, text, perl = T)を実行したところ、以下のエラーが出力されました。patternの文字数(59291文字)を減らせば、正常に動作します。patternの文字数を減らすことなく、正常に動作させることは可能でしょうか?環境は、R version 2.11.1 (2010-05-31); x86_64-apple-darwin9.8.0です。

以下にエラー regexpr(pattern, text, perl = T) : 
追加情報:   警告メッセージ: 
In regexpr(pattern, text, perl = T) :
 PCREパターンのコンパイルエラー
	'regular expression is too large' 
	at ''

read.csv で大きな値を読み込み表示する

g (2010-09-14 (火) 11:54:44)

例えば,

A,B,C
3121234567890,10,20
3122345678901,20,30
3123456789012,30,40

という中身のcsvファイルa.csvをread.csvで読み込ませると

> read.csv("a.csv")
             A  B  C
1 3.121235e+12 10 20
2 3.122346e+12 20 30
3 3.123457e+12 30 40

となってしまいます.浮動小数点ではなくそのままの形で読み込ませるにはどうしたらいいでしょうか.

plot エラー xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' lengths differ

id (2010-09-12 (日) 21:34:03)

密度関数を直接積分することで、分布関数を求めてプロットしようとしているのですが、うまくいきません。
簡単に書きなおすと以下のようなコードで、実行するとエラーを吐きます。

> pdf <- function(x) exp(-x)
> cdf <- function(x) integrate(f, 0, x)
> plot(cdf)
 以下にエラー xy.coords(x, y, xlabel, ylabel, log) : 
  'x' and 'y' lengths differ

cdf(Inf) などとすると正常な答えを返してくれるので、関数として機能はしているようですがプロットしてくれません。
プロットするにはどうすればよいでしょうか?

ちなみに実際に考えている分布は beta prime distribution です。
この分布の cdf に出てくる 2F1 関数の扱いがわからず、上のような苦肉の策をしようと思っています。

pooled adjacent violator algorithmについて

sakura (2010-09-12 (日) 17:26:35)

Rで、pooled adjacent violator algorithm を使ったライブラリーには何があるのでしょうか?

for文内での連番オブジェクトの記述

mtanaka (2010-09-07 (火) 20:12:53)

for文で任意の関数に連番で命名したオブジェクトを入力する場合どのように記述すればよいのでしょうか?
for文内でのファイル名の記述様式がわかりません。

関数 func()
オブジェクト(入力データ): hoge01, hoge02, hoge03
オブジェクト(出力データ): piyo01, piyo02, piyo03
実行例
piyo01 <- func(hoge01)
piyo02 <- func(hoge02)
piyo03 <- func(hoge03)

この上記の実行例をfor文で記述したいです。
どなたか、ご助言いただける方がいらっしゃいましたら、宜しくお願いいたします。

パッケージに含まれているはずの関数が呼び出せない

rcddnsj (2010-09-06 (月) 23:54:37)

Windows7でR version 2.11.1を使用しております。

サポートベクターマシンを試してみたくて、パッケージe1071をインストールし、関数svmを実行するところまでは、たどり着きました。

しかし、e1071に含まれているはずの関数predict.svmで予測しようとすると、「関数 "predict.svm" を見つけることができませんでした」とのエラーでpredict.svmを呼び出せませんでした。

操作はRcmdrのスクリプトウィンドウから実行しています。

このように、パッケージに含まれているはずの関数が呼び出せない場合、どのような原因とあるいは、対処法が考えられるのでしょうか?

Rの再インストール、パッケージe1071の再インストールなどはもちろん試しています。

sessionInfo() の実行結果は下記の通りです。

R version 2.11.1 (2010-05-31) 
i386-pc-mingw32 

locale:
[1] LC_COLLATE=Japanese_Japan.932  LC_CTYPE=Japanese_Japan.932    LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C
[5] LC_TIME=Japanese_Japan.932

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

other attached packages:
[1] SparseM_0.86    e1071_1.5-24    class_7.3-2     Rcmdr_1.6-0     car_2.0-2       survival_2.35-8 nnet_7.3-1
[8] MASS_7.3-7

loaded via a namespace (and not attached):
[1] tools_2.11.1

どなたか、ご助言いただける方がいらっしゃいましたら、宜しくお願いいたします。

ある条件に対応する要素だけを置換したい

ito (2010-08-31 (火) 20:12:14)

あるベクトルについて、他のベクトルの中のNAに対応する要素だけをNAに変換したいのです。よい方法はないでしょうか。

x <- c(1,3,4,6,8,9)
y <- c(1,NA,2,NA,3,4)

変換後

> x
[1]  1 NA  4 NA  8  9

となるようにしたいのです。

ifelse(x[is.na(y)],NA)

とか考えたのですが、うまくいきません。
よろしくおねがいします。

pls回帰分析の有意検定について

初心者 (2010-08-30 (月) 16:56:25)

いつもお世話になっております。
 
pls回帰分析をしたいのですが,Rですと,ベータは求められても,そのベータが有意かの検定はできないのでしょうか。またできるのなら,どのように行うのでしょうか。

pls回帰分析を学ぶに当たり,岩田先生のpls回帰入門やAcreMaker様のホームページを見ましたが,有意かどうか扱っていませんでしたので,質問させて頂きました。
 
使用環境は、R2.11.1です。OSはWindowsVista です。
 
よろしくお願いします。

複数条件でベクトルを抽出する際に、条件が定数ではなくベクトルだったときについて

Saito (2010-08-27 (金) 08:52:15)

いつもお世話になっております。
似た様な質問はいくつかあったのですが、意外にも、条件抽出が複数で、かつ固定値ではなく、ベクトルで条件抽出をしている例が見当たらなかったので、質問させてください。

ある参照列(座標と価がセット)があるときに、それを参照して、新しい座標から、それにマッチングする価を参照したいと思っています。以下がサンプルプログラムです。

set.seed(1)
a <- seq(1, 10)
b <- seq(1, 10)
c <- rnorm(100)

###参照列###
d <- data.frame(expand.grid(a=a, b=b), c=c)

###当てはめたい座標###
e <- data.frame(a=sample(a, 50000, rep=T), b=sample(b, 50000, rep=T))

###これでは上手く動かない###
subset(d$c, d$a==e$a & d$b==e$b)

###これでも###
subset(d$c, d$a%in%e$a & d$b%in%e$b)


つまり、eの座標軸が与えられたときに、dの座標軸と対応させて、dの三列目(c列)を引っ張って来たいのです。

for文でやろうと思えばできるのですが、実際はもっとeが大きくて、for文が実行スピード的に使えません。%in%も、&が入っていなければ使えたのですが、この場合上手く動作しないようです。おそらく単純な問題だと思うのですが、思うようにいきません。

どなたか、分かる方がいらっしゃいましたら、ご教示いただけると幸いです。

> Ord <- data.frame(ord=seq(nrow(e)))   # 本来の順序を示す作業用変数
> A <- merge(merge(cbind(e,Ord),d),Ord) # 欲しい結果
> AA <- A[,-1]                          # 作業用変数列が目障りなら
> str(e)
'data.frame':	50000 obs. of  2 variables:
 $ a: int  3 3 6 3 2 6 6 2 3 8 ...
 $ b: int  6 8 4 8 7 8 4 2 6 7 ...
> str(A)
'data.frame':	50000 obs. of  4 variables:
 $ ord: int  1 2 3 4 5 6 7 8 9 10 ...
 $ a  : int  3 3 6 3 2 6 6 2 3 8 ...
 $ b  : int  6 8 4 8 7 8 4 2 6 7 ...
 $ c  : num  0.3411 0.6107 -0.415 0.6107 -0.0392 ...
> str(AA)
'data.frame':	50000 obs. of  3 variables:
 $ a: int  3 3 6 3 2 6 6 2 3 8 ...
 $ b: int  6 8 4 8 7 8 4 2 6 7 ...
 $ c: num  0.3411 0.6107 -0.415 0.6107 -0.0392 ...

もしくは(こちらの方が直感的でわかり易い?)

> ed <- merge(cbind(e,Ord),d)
> str(ed)
'data.frame':	50000 obs. of  4 variables:
 $ a  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ b  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ ord: int  32810 46779 23088 6441 17984 15437 12878 47025 38836 21870 ...
 $ c  : num  -0.626 -0.626 -0.626 -0.626 -0.626 ...
> ed[order(ed$ord),][,-3]   # 上のAAと同じもの

関数"pls"を見つけることができませんでした

学生 (2010-08-26 (木) 22:48:35)

RでPLSを行おうとして追います。
PLSパッケージをインストールし,読み込んでも,「関数"PLS"を見つけることができませんでした」と返ってきます。
Rになれるための練習として,インターネット「岩田先生のPLS回帰入門」の「wine」についてPLSで分析しようとしているのですが上手くいきません。
何卒よろしくお願い致します。

図の軸の数値が指数表示になる

(2010-08-26 (木) 16:50:23)

図を作成した際に,x軸やy軸に添えられる数字が指数表示(1e-01,1e+01など)に自動的になります。これを回避したいのですが,どなたか方法を教えて頂けませんでしょうか?
何卒よろしくお願い申し上げます。

DEAの線形計画問題について

atsuo (2010-08-23 (月) 11:52:39)

optimやconstrOptimで解決できないかと、あれこれやってみたのですが、
以下の線形計画問題を解決したいのです。

� 2*θ−2*λ[1]−2*λ[2]−1*λ[3] >= 0
� 1*θ−1*λ[1]−2*λ[2]−2*λ[3] >= 0
�  1 −1*λ[1]−1*λ[2]−1*λ[3] <= 0
� λ[1]+λ[2]+λ[3] >= 1
� λ[1]>=0
� λ[2]>=0
� λ[3]>=0
� �〜�の条件を満たし、θを最小にする問題です。

欲しいのはθの値とλ[1]、λ[2]、λ[3]の値です。
これはDEA分析なのですが、パッケージのFEARでも準備されていない解θが欲しいのです。
環境は、R2.11.1です。OSはWindows7 です。

GARCHの解釈

(2010-08-20 (金) 22:42:44)

Coefficient(s)や、Hessian Matrixに標示される、muやomegaの解釈方法はどこを参照すればよいのでしょうか?ご教示頂ければ幸いです。
例) Coefficient(s):
  mu ar1 ma1 omega alpha1

なるほど、貴重なアドバイスを有難うございます。出てきたGarch項を推計式の誤差項に代入する場合には、同時方程式を使う形になるのでしょうか? テキストを見比べて、自分で式を考えてみたのですが、なかなかうまくいきません。 アドバイス頂けると幸いです。宜しくお願い致します。 (fgarchはうまくInstallできなかったため、そちらのマニュアルを確認することは、失念していました。調べ方が甘く、すみませんでした。 ご丁寧にアドバイスをいただいて、有難うございました。)-- z 2010-08-31 (火) 19:05:42

Rのインストールで,メッセージ言語が文字化け

rの初心者 (2010-08-20 (金) 16:28:28)

とあることでRを用いて多変量解析を行う必要が生じまして、ダウンロードしてインストールしました。ところがメッセージ言語は文字化けしているようで日本語で表示されません。ただし「ファイル」「ヘルプ」といった表示は正常です。どうすればよろしいでしょうか? よろしくお願いします。

TARCH has not been available?

z (2010-08-20 (金) 14:29:41)

RでTARCH(Threshold GARCH)をLoopで100社*10年分行おうと調べてみたのですが、該当する機能が見つかりませんでした。
RでTARCHするのと同じ結果を得られる機能がございましたらご教示頂けると幸いです。
宜しくお願いします。

forループのリストにベクトルを使ってはいけないのでしょうか?

青葉ほととぎす (2010-08-20 (金) 14:18:14)

最近、Rでプログラムを始めたものです(Vistaでversion 2.10.1を使用)。
あるプログラムを作っているときに不可解な結果がでて、いろいろ試した結果、forループのところに問題がることが分かってきました。

以下は、問題の部分だけを抜き出したものです。

初め実験1を行っていたのですが問題があり、実験2のように書き変えたところうまくいきました。

sim.rate <- seq(0, 1, by=0.01)
result.matrix <- matrix(0, nrow=length(sim.rate), ncol=2)
# 実験1
for(i in sim.rate)
{
    result.matrix[i*100, 1] <- result.matrix[i*100, 1]+1
}
# 実験2
for(i in 1:length(sim.rate))
{
    result.matrix[i, 2] <- result.matrix[i, 2]+1
}

この2つは以下のような異なる結果を出します。

       [,1] [,2]
  [1,]    1    1
中略
 [27,]    1    1
 [28,]    2    1
 [29,]    0    1
 [30,]    1    1
中略
 [56,]    1    1
 [57,]    2    1
 [58,]    0    1
 [59,]    1    1
中略
 [95,]    1    1

そこで質問ですが、forループのリストしてベクトルを使ってはいけないのでしょうか?なぜこんな違いが生じるのかも教えていただけると助かります。

Windows7 Home Premiumで最新のRインストールできず

ななしのごんべ (2010-08-12 (木) 19:53:54)

Windows7 Home Premiumで最新のRを入れたら、パッケージがインストールできません。「コンピューターにiconv.dllがないため、プログラムを開始できません。この問題を解決するには、プログラムを再インストールしてみてください。」というエラーウィンドウが出ます。どなたか解決策をご教示いただけるとありがたいです。

積分した値が計算に使えない

しょーじん (2010-08-11 (水) 14:36:14)

積分した値を計算に使おうと思ったところ,「二項演算子の引数が数値ではありません」と出てしまい,計算に使えません.
その値自体を調べると「0.5 with absolute error < 5.6e-15」とでます.
計算に使うためにはどうしたらよいのでしょうか?

単回帰のデータが認識されない?

z (2010-08-07 (土) 23:56:43)

OLSを試みています。
推計式はComp1=α+βcomp0なので、以下の要領で、Syntaxを書いたら、Errorが出ました。

--Syntax-----
 Result1<-lm(comp1~comp0)
--error--------
Error in eval(expr, envir, enclos) : object 'comp1' not found

利用しているデータは下記の要領で、Comp0-90 まで、Dayは1-150まで。回帰式は、Comp<n>=α+βComp0 を Comp1-90まで順にComp0で回帰する予定です。

day        comp0        comp1        comp2        comp3        comp4    
1     -0.010873538 -0.017208413  0.009237875  0.000000000 -0.012048193
2      0.007344111 -0.001945525 -0.010297483  0.013313609  0.000000000
3      0.016795163  0.031189084 -0.019653179  0.002919708  0.034146341
4      0.006113627  0.003780718  0.004716981  0.010189229  0.007075472
5     -0.002080986 -0.007532957  0.011737089  0.010086455  0.145199063
6      0.005460344  0.001897533 -0.012761021 -0.021398003  0.022494888

Comp2以降はLoopでの処理を試みる予定なのですが、Comp1についての式は間違っていないようなのにエラーが出るので、調べ方にも窮しています。

gqt<-lapply(1:88, function(i){
  gqtest( lapply(1:88, function(i){
    lm(eval(parse(text=sprintf("comp%i", i))) ~ eval(parse(text=sprintf("comp%i", 0))), data=dataset)
  }), data= dataset,fraction=50, order.by=NULL) })
 Error in terms.default(formula) : no terms component 

どうアレンジしたら良いのでしょう? 調べ方等、アドバイスを頂けると幸いです。 宜しくお願い致します。 -- z 2010-08-31 (火) 19:00:42

+が勝手に出て、構文が使えない。

z (2010-08-07 (土) 21:47:40)

二進も三進も行かず、困っています。データを読み込む行はエラーが出ないのですが、二行目から+が出て、どんな構文を打っても、+が続きます。どこが間違っているのでしょうか?
例えば、下記のようになります。

dataset<-read.csv("C:/Users/(myname)/Documents/(foldername) /(sheet name).csv”, header=T)
+ attach (dataset)
+

行列データの要素計算方法について

森の熊五郎 (2010-08-05 (木) 16:21:03)

X <- c(0.1, 0.2, 0.3, 0.4)
Y <- c(0.3, 0.5, 0.8, 0.5)
Z <- c(0.6, 0.5, 0.7, 0.9)
sum(X)
sum(Y)
A <- rbind(X, Y)
sum(A)
sum(Z)
B <- rbind(A, Z)
sum(B)

この式を繰り返し文 forを使って実行する方法を知りたい。

BASICでは例えば、A(20,10)の各要素を足し算する方法として

S=0
for I=1 to 20
  for J=1 to 10
    S=S+a(i,j)
  next
next

なる文にて計算が出来ますが、Rの場合、内側のループはsum文で処理をするとして、外側のループをFor文で制御仕様とした場合、どのように文を組めば良いでしょうか?ご教示下さい。上記の例ではsum文とrbind文と使えばよいのですが、rbind文を使わず、for文を使って組む方法を知りたいと思います。

memory.limiitについて

yoshi (2010-08-02 (月) 14:56:40)

memory.limit()について質問です。
R上で

memory.limit()
[1] 1535

memory.limit(T)
[1] 12

memory.limit(F)
[1] 10
となります・ここで

memory.limit(4000)
[1] 4000
とした後は

memory.limit()
[1] 4000

memory.limit(T)
[1] 12

memory.limit(F)
[1] 10
となります。ここで質問なんですが「memory.limit()」「memory.limit(T)」「memory.limit(F)」の違いはなんなのでしょうか?今現在Rのメモリの上限のことで非常に困っています。また、皆さんはこのようなコマンドをどこから学んでいるのでしょうか?もし、良いWEBページ等がありましたら教えてくれると幸いです。

パッケージ「R.huge」の読み込み方法

yoshi (2010-08-02 (月) 10:49:35)

パッケージ「R.huge」のパッケージの読み込みが上手くいきません。
私がやったことは
�CRANミラーサイトから「R.huge」のzipファイルをダウンロード
�RのGUIで「ローカルにあるzipファイルからのパッケージのインストール」で「R.huge」のインストール
→RのGUI上では「パッケージ 'R.huge' は無事に開封され、MD5 サムもチェックされました 」と表示
�RのGUIで「パッケージの読み込み」で「R.huge」を選択し実行
→RのGUI上で「要求されたパッケージ R.oo をロード中です
 エラー: パッケージ 'R.oo' をロードできませんでした
 追加情報: 警告メッセージ:
 In library(pkg, character.only = TRUE, logical.return = TRUE, lib.loc = lib.loc) :
 'R.oo' という名前のパッケージはありません 」と表示
このように�までは良いと思うのですが、�でエラーが表示されてしまい、パッケージ「R.huge」が読み込まれていないと思います。

これを解決する方法を自分でも探したのですが分らなかったため今回投稿さしていただきました。よろしくお願いします。

2進数表記した時、右から数えて最初にゼロが出現する桁数

aMC (2010-08-01 (日) 06:07:54)

たとえば、2→10なので1桁目、5→101なので2桁目、7→0111なので4桁目、などのように、2進数表記した時、右から数えて最初にゼロが出現する桁数を求める効率のよいやり方はありますでしょうか?ベクトルを与えると、何桁目かをベクトルで返してくれるようなものが欲しいのですが、forループを使うと時間が掛かり過ぎます。

b進展開(10進数→2進数、2進数→10進数)関数

qMC (2010-07-31 (土) 22:56:24)

10進数→2進数、あるいは2進数→10進数へと変換する関数はありませんでしょうか?

Rscriptでコードを実行するとグラフに文字化けが発生

りりぽん (2010-07-31 (土) 11:56:57)

Rscriptで次のように、test.Rを実行させると、グラフタイトルの日本語文字化けが発生します。Rscriptを使わずに、直接、test.Rを実行すると問題なく日本語が表示されます。どなたか、解決方法をご教授いただけたら幸いです。よろしくお願いいたします。

###### test.bat ###################
C:\R-2.11.1-x64\bin\Rscript test.R
####### test.R ####################
postscript(file="test.eps", horizontal=F, family="Japan1Ryumin")
plot(1:10, main="日本語")
dev.off()
###################################

使用フォーマット・R環境:
ソースコード;UTF-8(BOMなし)エンコード
コマンドプロンプト;cp932;MSゴシック

R version 2.11.1 (2010-05-31) 
x86_64-pc-mingw32

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

forを使って変数を作成したい

basel_firb (2010-07-29 (木) 15:18:47)

過去に「行列にfor文中の変数を含んだ名前をつける 」という質問があり、内容を読み、構文 eval(parse(text=...)を使用するということらしいのですが、よく分かりませんでした。

データフレームに、変数var1 var2 var3・・・var68というのがあり、for文を使って,h1=var1/var68,h2=var2/var68・・・と分子のvarを1〜67まで繰り返し計算したいのですが、どのようにすればよいのでしょうか?

Rのメモリ使用量

sano (2010-07-29 (木) 12:39:07)

こんにちは、hclustで階層的クラスタリングを行おうとすると、

サイズ 3.4 Gb のメモリブロックを割り当てることができません。

とエラーメッセージが出てしまいます。
OSはwindows7,Rのversionは version 2.11.1(64bit), メモリは48GB積んでいます。しかしwindowsのリソースモニターを見ると、11GB程度でまだ余裕があります。

> memory.limit(T)

と実行すると、

[1] 7363

と表示されるので、Rが物理メモリを7GBしか、OSから割り当てられていないということではないかと思います。
30GB割り当ててもらおうと、

> memory.limit(30000)

と入力すると、

[1] 49149
警告メッセージ: 
In memory.size(size) :  メモリー限界を減らすことができません。

と表示され、うまくいきませんでした。
CPUがデュアルCPUでメモリを片側24GBずつ積んでいるのですが、このことも何か関係があるのでしょうか。
いろいろ調べてみたのですが、どうにもいかず、投稿させて頂きました。
お教えいただければ、幸いです。よろしくお願い致します。

T,Fを文字列として出力するには

sayaka (2010-07-26 (月) 11:59:13)

いつも参考にさせていただいております。自分でも本やGoogleなどで調べたのですがどうもわからず、質問させていただきました。
以下のようなデータ"test_data.csv"があるとき

"T","T"
"T","F"
"T","T"
"F","F"

data <- read.csv("test_data.csv",header=F)

というようにread.csvを使ってデータを読み込むと、以下のようにTやFが文字列ではなく論理値として出力されてしまいますが、どのようにすれば文字列として認識させることが出来るでしょうか?

> data
     V1    V2
1  TRUE  TRUE
2  TRUE FALSE
3  TRUE  TRUE
4 FALSE FALSE

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

式が上手く評価されない?

em (2010-07-21 (水) 01:30:51)

runif()関数を使って一様乱数を生成し、0〜100の範囲に変換してみました。

x <- round(runif(100), 2) * 100

ここで、7の出ている回数を調べようと

length(x[x==7])

とやると、実際にはベクトルxに7が含まれていても、0になってしまします。
また、x[x==7]のような式の結果もnumeric(0)となってしまいます。
どこがまずいのでしょうか?

デンドログラムのデータNo.の表示指定

hashiro (2010-07-18 (日) 16:50:49)

「Rによる統計解析」を読みながらクラスター分析を行っています。
エクセルで、194の変数、160行のデータです。
R内に読み込んでデンドログラムまでは出来ました。
この時、デンドログラム最下段には1〜160までのNo.が表示されますが、
この数字をエクセルのa列の文字で表示したいのです。
plotの中にlabels= で指定すればいいみたいですが、指定の仕方が悪くエラーが出てしまいます。
指定方法をどうか教えていただけませんか。

相関ルールについて

t_y (2010-07-16 (金) 10:57:30)

今現在相関ルール抽出のパッケージ「arules」を用いているのですが、分からないことがあります。しかし、言葉で端的に言うのは難しいので例を交えて説明します。
買い物バスケットで以下の例があるとします。

TID	アイテム集合
1	{パン、牛乳}
2	{ハム、牛乳}
3	{ビール、たばこ}

このデータを「data」に代入します。

> data <- list(c("パン", "牛乳"), c("ハム", "牛乳"), c("ビール", "たばこ”))

そして transactions 形式のデータを「data.tran」に代入します。

> data.tran <- as(data1, "transactions")

次に相関ルール「data.ap」に生成します

> data.ap <- apriori(data.tran)

生成された相関ルールは以下の 3 つです。

  lhs         rhs        support
1 {パン}   => {牛乳}   0.3333333
2 {ハム}   => {牛乳} 0.3333333
3 {ビール} => {たばこ} 0.3333333

相関ルールは lhs⇒rhs で「lhs が存在する場合に rhs が存在する」という意味です。
Support は「特定のルールの数÷全てのルール」です。例えば{パン}→{牛乳}の場合、support=1/3 です。

ここでもしアイテム「パン」と「牛乳」が新たなアイテム「食べ物」に属するとします。
属した場合の相関ルールの例を載せます

  lhs         rhs        support
1 {食べ物}   => {牛乳}   0.6666666
2 {パン}   => {牛乳}   0.3333333
3 {ハム}   => {牛乳} 0.3333333
4 {ビール} => {たばこ} 0.3333333

このように「食べ物」の中に「パン」と「ハム」が含まれているため suppot=2/3 になります。
しかし、data2 に以下のように代入し、

> data2 <- list(c(“食べ物”, ”牛乳”), c("パン", "牛乳"), c("ハム", "牛乳"),
+               c("ビール", "たばこ”))

同様に相関ルールを作成すると、

  lhs         rhs      support 
1 {パン}   => {牛乳}      0.25
2 {ハム}   => {牛乳}      0.25
3 {ビール} => {たばこ}    0.25
4 {食べ物} => {牛乳}      0.25

のようになり support の値が変化してしまいます。
これを解消する方法を教えて下さい。

SensomineRの分散分析でのエラー

Sensory (2010-07-15 (木) 19:17:26)

分散分析を行うと出現するエラーについてです.R,Rcommander,SensomineRを使っています.SensomineRよりpanel performance(分散分析)を起動させた場合,Product35×Panelist5×descriptor9ではエラーは出現しないのですが,Product31×Panelist5×descriptor9では,"置き換えるべき項目数が,置き換える数の倍数ではありませんでした”というエラーが出現します.

何か解決策があれば教えて頂けないでしょうか.よろしくお願い致します.

Rによる、SPSSファイルからSTATAファイルへの変換の仕方について

大学生 (2010-07-12 (月) 20:25:01)

SPSSファイルを、STATAファイルに変換したいのですが。

library(foreign)
read.spss("datafile",use.value.labels=FALSE)

で、SPSSファイルを読み込んだあと、そのまま、STATAファイルに変換して、保存したいのですが。この後のコマンドがわかりません。

Rのバージョンは2.8.1で、MACです。どうかよろしくお願いします。

3次元グラフにおける座標点の表示

メジロウ (2010-07-09 (金) 21:24:59)

マトリックスを3次元グラフで描画した時、x、y、zの座標まで、グラフ内に頂点として表示されてしまいます

x、y、zの座標点を取り除き、純粋に数値の座標点のみを頂点として描画するにはどうすれば良いでしょうか?

打ち込んだコマンドラインは以下の通りです

jhin <-matrix(c(
1.1479, 0.8587, 0.8719,
-1.099, 0.353, 0.972,
-1.3226, -0.4005, -0.9295,
-0.4804, -0.7637, -1.3671,
-0.1707, -1.3382, 1.1815,
0.2106, 0.7653, -1.7421,
0.3761, 1.8047, 0.0608,
0.4473, -1.3014, 0.7972,
1.625, 0.7415, 0.5675,
1.5616, -0.7194, -0.4122),
nrow = 10, ncol = 3)
jhin
rownames (jhin) <- paste("in", 1:10, sep="")
colnames (jhin) <- c("X", "Y", "Z")
jhin
rg <-jhin
library(sna)
gplot3d(rg, thresh = 1.1, displayisolates = TRUE, suppress.axes = FALSE,
        displaylabels = TRUE, xlab = "dimension 1",
        ylab = "dimension 2", zlab = "dimension 3")

Dataframeからある特定の因子を持つ行を抜き出したい

hello r (2010-07-09 (金) 15:46:59)

ググッたり、ここのDatagram tipsを見たのですが、解決しなかったので質問させていただけると助かります。

CSVファイルから読み込んだファイルに因子が含まれています。

      t   s
1   3755  sd
2   3840  sd
3   3856  si
4   3884  sd
5   4011  si
6   4031  sd
7   4033  sd
8   4115  sd
9   4147  si
10  4157  si


ここから因子ごとにdataframeを分離したいと考えています。
因子でなければ

v <- dataframe[,"s"] == "sd"
dataframeSD <- dataframe[v,]

などとして抽出できると思うのですが、最初の因子の比較がDataframeの該当列の因子の取得方法が分からないため、マッチングが行えません。
どのように対処すればよいでしょうか。
もしよろしければご教授願えないでしょうか。
よろしくお願いします。

メトロポリス法によるパラメータ推定プログラムを自作するにあたって

Saito (2010-07-06 (火) 00:54:49)

いつもお世話になっています。
ググっても、過去ログを漁っても、持っている本を読んでも、どうしてもわからなかったので、質問させてください。
都合により、メトロポリス法によるパラメータ推定プログラムを自作しています。が、メトロポリス法によるパラメータ推定がうまくいきません。以下に例を示します。

set.seed(1)
iter <- 20000 

###真の値###
a <- 10
b <- 2
x <- seq(0.01, 1, length=500)
y <- b*x + rnorm(500, a, 1) 
###尤度の計算(!!!対数にはしていない!!!)###
cost_func2 <- function(a3, b3) { 
sum(dnorm(y, a3 + b3*x, 1))
}
###乱数と初期値の設定###
u <- runif(iter, 0, 1)
u2 <- runif(iter, 0, 1)
a4 <- b4 <- numeric(iter+1)
a4[1] <- b4[1] <- 5
lag <- 4

###メトロポリス法(?)の実行###
for(i in 1 : iter) { 
a4_2 <- a4[i] + runif(1, -lag, lag)
b4_2 <- b4[i] + runif(1, -lag, lag) 
###一つ前のサンプルで計算されるcost_funcの値が近すぎると、ほとんど採択される???###
a4[i+1] <- ifelse(u[i] < cost_func2(a3=a4_2, b3=b4[i])/
                         cost_func2(a3=a4[i], b3=b4[i]),
                         a4_2, a4[i])
b4[i+1] <- ifelse(u2[i] < cost_func2(a3=a4[i+1], b3=b4_2)/
                          cost_func2(a3=a4[i+1], b3=b4[i]),
                         b4_2, b4[i])
}
hist(a4[(iter/2):iter])
hist(b4[(iter/2):iter])


疑問点は二つです。メトロポリス法の場合でも、尤度を計算すると思いますが、対数尤度にするとマイナスになる場合があります。メトロポリス法では、棄却するか採択するかのときに、一つ前のパラメータのサンプルと現在のパラメータのサンプルの比をとります。しかし、比を取ってしまうと、マイナス同士が消しあってしまうこともあります。例えば、パラメータ-10が前回のサンプル、-11が今回のサンプルだとします。-10と-11の比は、1.1ですよね。同じく、10と11の比も1.1です。しかし前者は前回のサンプルよりも尤度は小さくなっており、後者は尤度が大きくなっています。しかし、比をとるだけだと、両者とも採択されてしまう気がします。そのため上記のプログラムでは対数にしていません。
もう一つの疑問は、メトロポリス法では、上記のように比を取った後、一定確率(u)で、尤度が小さい方向にサンプルされても、採択する場合があります。しかし、前回のサンプルと今回のサンプルの比が、例えば、先ほどの11、10でサンプルされる順番が逆だったとすれば、10/11=0.91となります。この0.9という数字は、ほとんど1に近く、ずっとこの比でサンプルが続けられるとすると、どんどん間違った方向へサンプルが続いていきます。

このようなことを避けるにはどうすればよいのでしょうか。
上記のプログラムですと、a4とb4の平均値がそれぞれ10と2になればよいのですが、うまくいきません。仮にb4を消して、切片モデルでやるとうまくいくのですが、それでは解決にならないので・・・。

どなたかわかる方がいましたら、ご教授のほどよろしくお願いいたします。
なお、環境はWindows7, R-2.11.1です。

地理的加重回帰分析のgwrの予測値

酔鯨 (2010-07-03 (土) 10:33:29)

spgwrの地理的加重回帰分析のgwrの予測値は、XXX$SDFのテーブル出力を見れば、predの項目(厳密には、ヘッダがずれているので、右に1項目ずらす必要がある。)を見ればよいことは判りました。しかし、この関数を使い予測するためには、被説明変数の値が必要です。つまり、何らかの方法で予測しなければなりません。通常の重回帰で予測をするためには、回帰係数と説明変数の値だけで良いです。地理的重回帰分析では、通常の重回帰の予測と同じように、回帰係数と説明変数だけで予測値を得る関数はないのでしょうか?

Rによるアクセスログ解析

tadashi (2010-07-02 (金) 15:05:54)

"Rでアクセスログ" と検索してもでてきません。Rでアクセスログ(にかぎらず、ログデータ)の解析をすることはそれほどないのでしょうか?
もし、awstats 等でやっているようなことをRで代替している事例がありましたら、お教えください。

統計解析関数Tipsのrepeated mesure aonovaのコード

ちゃーぴー (2010-07-02 (金) 10:04:01)

RjpWiki内の統計解析Tipsにある反復測定分散分析 (Repeated measured ANOVA)(http://www.okada.jp.org/RWiki/index.php?R%A4%CE%C5%FD%B7%D7%B2%F2%C0%CF%B4%D8%BF%F4Tips#content_1_7)を使おうと思い,コードをコピペしようとしました.しかし,print.rep.anovaのコードが,ペーストしている途中でエラーとなります.当方の環境WindowsXP×2,Windows7-64bitでRのバージョンは2.11.0です.なお,64bitマシンでは32bit版,64bit版R両方を使用しています.3台のマシン全てで同様の結果となります.下記のような感じです.どなたか対策をご存知の方いらっしゃいますでしょうか?

+         rownames(ttx)[rownames(ttx)=="Residuals"]<-sub("Error: (.*\)","Error(\\1)",names(tx)[i])
 エラー:  "Error: (.*\)"で始まる文字列の中で '\)' は文字列で認識されないエスケープです


Shaffer法の多重比較

tasosi (2010-07-01 (木) 18:43:27)

Shaffer法の多重比較に関するパッケージを探しています.
色々調べてみたのですが,ANOVA君を活用するという手段はあったのですが,パッケージでは見つかりませんでした.
多重比較で,paired-t-testのP値を補正するのに用いたいと思っていますが,パッケージとしてはないのでしょうか?

行列におけるNA値を前後の行の値から補完したい

R初心者 (2010-07-01 (木) 14:27:54)

はじめまして。
ネットワークの研究をしている大学院生です。データのグラフ化をしたくて一週間ほど前からRの勉強を始めました。
googleとこのwikiでググってみたのですが、分からなかったので質問させていただきました。
環境はMacOSX 10.5.8
Rのバージョンは2.9.0です。
−−−−−−
やろうとしていることの概要
2つのプログラムからの出力をまとめた欠損値を持つCSVを、欠損値を補完したCSVにして、それをlatticeで線グラフ化する。
補完せずとも、線グラフが切れない手法があればそれでもOKです。
−−−−−−
ネットワークのスループットとRTTをTcpdumpを解析して出力するプログラムを書きましたが、都合上RTTとスループットが別に出力されて次のような形式になっています。

       time, seqnum(シーケンス番号), rtt, throughput
       0.0, 3136068389, 0.0469200000006822,  
       0.0808820000020205, 3136070885, 0.0611399999979767,  
       0.0999999999985448,  ,  , 24960.0~
       0.142041999999492, 3136074629, 0.0901450000019395,  
       0.142041999999492, 3136072133, 0.0699139999996987,  
       0.19999999999709,  ,  , 37440.0
       0.21197000000393, 3136077125, 0.0549669999963953,  
       0.232200999998895, 3136079621, 0.0598790000003646,  
       0.266950000004726, 3136084613, 0.0701439999975264,  
       .....

これをRにCSVとして食わせると、次のようになります。

         time     seqnum      rtt throughput
     1     0.000000 3136068389 0.046920         NA
     2     0.080882 3136070885 0.061140         NA
     3     0.100000         NA       NA      24960
     4     0.142042 3136074629 0.090145         NA
     5     0.142042 3136072133 0.069914         NA
     6     0.200000         NA       NA      37440
     7     0.211970 3136077125 0.054967         NA
     8     0.232201 3136079621 0.059879         NA

このNA値を前後の行の値の中間値となるように補完したいです。
なにかそのような関数はありますでしょうか。
補完のやりかたは
「1,2,NA,4,5」-> [1,2,(2+4)/2=3, 4, 5]
のようなもので、連続した値が抜けた際は
「1,2 NA, NA, NA, 6」-> [1,2,3,4,5,6]
となってもらえるのが理想です。

最初は作ろうかと思ったのですが、既にありそうな気がしたので質問させていただきました。
よろしくお願いします。
−−−−−−
自己解決しました。 Vectorを補完する関数を書いたので、行列を「列で分解、適用、合体」することで補完できると思います。

complementNaVector <- function(vector) {
	posVector <- 1:length(vector)
	naPosVector <- posVector[is.na(vector)]
	notNaPosVector <- posVector[!is.na(vector)]
	
	for (index in naPosVector) {
		lowerVector <- notNaPosVector[notNaPosVector < index]
		greaterVector <- notNaPosVector[notNaPosVector > index]
		
		
		if(length(lowerVector) != 0 && length(greaterVector) != 0) {
			low <- lowerVector[length(lowerVector)]
			high <- greaterVector[1]
			distance <- vector[high] - vector[low]
			stepl2h <- high - low
			stepl2i <- index - low
			vector[index] <- vector[low] + (distance * stepl2i / stepl2h)
		} else {
			vector[index] <- NA
		}
	}
	return(vector)
}

行列の補完

complementNaMat <- function(matrix) {
	iterator <- 2:ncol(matrix)
	mat <- complementNaVector(matrix[, 1])
	for(index in iterator){
		vector <- complementNaVector(matrix[, index])
		mat <- cbind(mat, vector)
	}
	return(mat)
}
## 以下でよさそう
complementNaMat <- function(matrix) {
	for (index in 1:ncol(matrix)) {
		matrix[, index] <- complementNaVector(matrix[, index])
	}
	return(matrix)
}

データの重なり

森の熊五郎 (2010-06-27 (日) 01:02:10)

コールセンタに日々大勢の方から電話がかかってきます。1日に何回も、また続けて何日も。そのようなコールセンタにおいて、日々かかってくる電話が何件あって、その電話が何人によってかかっているのか(重複をなくして)。またその結果を、1日単位、1週間単位、1ヶ月単位、一年単位で調べようとしております。本来は電話番号を使って、調べるのですが、個人情報ということもあり、ここでは、仮に下記のような人名が書かれたデータベースがあって、それをもとに調べることとします。
C201004,C201005,C201006
Sato,Sato,Nishihara
Tanaka,Yamada,Sato
Yamada,Kojima,Nishijima
Inoue,Kitano,Inoue
Sato,Shinagawa,Yokosuka
Kitamura,Handa,
Yamamoto,Kitamura,
Ohta,Handa,
Inoue,Nishi,
Kitamura,Yokota,
Kobayashi,Sato,

csvのファイルに上記のようにデータベースがあります。
行方向は日にち、列方向はその日に利用されたお客様のお名前が保存されています。

call=read.csv("call2010.csv",header=TRUE)
call201004=call$C201004
call201004=call201004[!is.na(call201004)]
call201005=call$C201005
call201005=call201005[!is.na(call201005)]
call201006=call$C201006
call201006=call201006[!is.na(call201006)]
length(call201004)
length(unique(call201004))
length(call201005)
length(unique(call201005))
length(call201006)
length(unique(call201006))
Tcall=c(call201004,call201005)
length(Tcall)
length(unique(Tcall))
Tcall=c(Tcall,call201006)
length(Tcall)
length(unique(Tcall))~

というプログラムを作りました。
�今回は3日分しかないので、簡単なのですが365日分だとさすがにベタ書きですと大変なので、BasicでいうところのFor next分みたいな文を作り簡単に処理したいのですが可能なのでしょうか?
�文字列のベクトルの結合はc(X,Y)ではうまくいきません。どうすればよいのでしょうか?
�欲を言えばピボットテーブルのように、お客様毎にいつ何件利用されたのかをしりたい。
コールセンタの利用件数は月に何万件もあって膨大なので、excelでは処理できず、Rに挑戦しているのですが、まだまだ素人なので。お手数をおかけしますがご教示下さい。

neuralパッケージのインストールについて

hiro (2010-06-24 (木) 02:54:33)

はじめまして、Ubuntu10.04(64bit)にて、R:2.10.1を使っている者です。
現在、豊田秀樹氏の「データマイニング入門」という書籍を読みながら、Rを学んでおります。

この書籍の第2章の「鉛筆の数え方」のところで、『neuralライブラリ』を用いたサンプルコードが載っているのですが、『neuralライブラリ』をインストールしようにも、install.packages("neural")としても、パッケージを見つけることが出来ないため、インストール出来ずに困っております。

この書籍はそもそもOS:windowsXP、R:2.7.0を前提にして書かれたものなので、neuralライブラリはLinuxには提供されていないライブラリなのでしょうか?

あるいは、neuralライブラリはLinuxにも提供されていたが、古くなったため、今では使われなくなってしまったのでしょうか?

この書籍でRを学びたいので、出来ればneuralライブラリを使いたいのですが、良い解決策をご存知の方がいらっしゃったら、対処方法をご教授下さい。

R + VIM (Windows)

尼河童 (2010-06-23 (水) 00:12:42)

みなさまこんにちは。
ウィンドウズマシーンで R 2.10.1 を使っています。はい、バージョンアップします。

エディタは昔から vim を使っているのですが、簡単に R コードを vim から R に飛ばす方法を模索しています。ウィンドウズというのがネックになりそうです。

素敵な方法をご存知の方はご教授ください。

png出力の背景透過について

としろう (2010-06-20 (日) 07:21:25)

plot関数などで作成した図を、png形式で出力する際に、枠やプロットした点以外の背景が透かしのままなのですが、背景を透かしの代わりに白で指定することはできるのでしょうか?
例えば下のプログラムでは、背景が透かしたままで出力されてしまいます。
使用環境はwindows,R-2.11.1です。

plot(1:20)
pp<-recordPlot()
png("test.png",bg="white")
replayPlot(pp)
dev.off()

Rで(ハフマン)符号の複合化

R初心者 (2010-06-19 (土) 21:29:59)

各アルファベットが0/1のコードに(ハフマン)符号化されているデータ「0,1,1,0,0,0,1,1,1,0,0,1,1,0,0,1,1,0,0,0,0,1,1,0,1,1,1,1,0,0,0,1」があるとします.例としまして,A〜Iのデータが以下の0/1のコードに符号化されている状況を考えます.

文字符号
A01
B111
C110
D101
E001
F000
G1001
H10001
I10000


上記の表でデータを複合化すると「0,1,1,0,0,0,1,1,1,0,0,1,1,0,0,1,1,0,0,0,0,1,1,0,1,1,1,1,0,0,0,1」は「AHCAGICBH」に複合化されます.

さて,データを複合化するための方法の1つが,符号を2分木化し,C言語等のポインタで2分探索法で複合化する方法です.9個の葉を持つ2分木は,9 個(種類)の異なる文字を含む文字列を符号化するために用いることができます.このとき,2分木の各葉を各文字に対応させます.各文字の符号は2分木の根から対応する葉へのパスによって決定します.このとき、左へ向かった場合 0,右へ向かった場合 1 とします.

 center

しかし,Rではポインタがなさそうなので,仕方なく以下のような方法で複合化してみました. データをx,上記のリスト(符号,文字)をy,zに格納するプログラムを以下で生成します.

> x <- c(0,1,1,0,0,0,1,1,1,0,0,1,1,0,0,1,     # データ
+        1,0,0,0,0,1,1,0,1,1,1,1,0,0,0,1)     # 結果:AHCAGICBH
> y <- c("01","111","110","101","001",
+        "000","1001","10001","10000")        # ハフマン符号
> z <- c("A","B","C","D","E","F","G","H","I") # 対応する文字

次に,データを1文字ずつ読み込み,リストyの要素のいずれに合致しているかをチェックし,合致していたらその符号に対応する文字をresultに格納する,というものです.

> result <- c()                       # 複合結果
> tmp  <- ""                          # バッファ
> for (i in 1:length(x)) {
+   tmp <- paste(tmp, x[i], sep="")   # データから1文字読み込み
+   j   <- 1                          # 複合化:yのどの値に該当するかチェック
+   for (j in 1:length(y)) {          # 
+     if (y[j] == tmp) {              # yのj番目の値に該当した場合
+       result <- c(result, z[j])     # zに複合化した結果(j)を格納
+       tmp    <- ""                  # バッファを初期化
+       break                         # 繰り返し文から抜ける
+     }
+   }
+ }
> result
[1] "A" "H" "C" "A" "G" "I" "C" "B" "H"

(1) 上記のプログラムでは,文字の種類が増えてくると計算速度が遅くなってしまいます.Rで,木構造を扱うようなことは出来るのでしょうか.
(2) Rで,C言語のポインタのようなことは出来るのでしょうか.
ご教示戴けますと幸いです.

grep("c",list(c("","")) の値が 1 になる問題

松田紀之 (2010-06-18 (金) 10:04:51)

Mac OS-X (10.5.8)でR2.11.1で文字列処理を試みているうちに,題名にあげた問題が見つかりました.以下,幾つか確認できたことです:
(1) patternの文字が"c"以外のアルファベットなら正常に動く.
(2) 参照されるlistの内容が c("") と単一の場合も正常に動く.
(3) listではなく,grep("c",c("","")) なら正常に動く.
プログラム全体の都合上,list() を対象にしています.何故この問題が起こるのか,またその対処法を教えてください.

batch処理でfile名を自動的に生成し、saveしたい

nakamura (2010-06-17 (木) 19:18:15)

単純な例として装置からの多数のdataをreadで取り込んで、XDR形式などで保存したいのですが、write.tableとsaveのfile名の指定のところで困ってます。" "で変数( fnとかfile.name[i] )を囲むと、文字列として変数が認識されているようで上手く行きません。ご教授お願いします。単純な例は以下です。”plotでのmain titleの指定”を参考にしました。

file.name <- c("map_1.csv", "map_2.csv") # ファイル名のベクトル
par(ask=T)
for (i in 1:2) {
    fn <- file.name[i] 	# ファイル名
    dt <- read.csv(fn , header=T, na.strings="NA")   
#   plot(dt, main=fn) 		# main タイトル付きの描画
    write.table(summary(dt), "fn_sum.txt", quote=F)
    save(dt, file="file.name[i].dat")
}

「0,1」のビットファイルを読み込む方法

R初心者 (2010-06-17 (木) 12:01:56)

「0100101000010101・・・」のように延々と続くファイルがあり,それを1文字(1ビット)ずつ読み込むことを考えております.

> ff <- tempfile()
> cat(file=ff,"01000101000100100101\n")
> read.fwf(ff, widths=rep(1,20))
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20
0  1  0  0  0  1  0  1  0   0   0   1   0   0   1   0   0   1   0   1

(1) ファイルを読み込む際文字数(ビット数)が分からなければ引数widthsを指定するときに困ります….読み込むファイルのサイズや文字数を計算する関数はありませんでしょうか.
(2) 上記では関数read.fwf()で無理やり読み込んでいますが,もっと良い関数があれば教えてもらえませんでしょうか.
どうかよろしくお願いいたします.

単語間の距離を求めデンドログラムの構築まで

AKB (2010-06-15 (火) 23:49:30)

基データ(csv)としてテキスト(301)×単語(1134)のマトリクスがあり,単語間の距離を算出し,デンドログラムの構築までしたいのですが”non-square matrix”と当然返されてしまいます.手順としては

data <-read.csv("C:/RTest/Real Game.csv", header=T, row.names=1)
data # 表示
old.op <- options(max.print=999999)
myCosine <- function(x) {
  ret <- matrix(0, ncol=ncol(x), nrow=ncol(x))
  for(i in 1:ncol(x)) { 
    for(j in 1:ncol(x)) {
      ret[i, j] <- (x[,i] %*% x[,j]) / (sqrt(sum(x[,i]^2)) * sqrt(sum(x[,j]^2)))
    }
  }
  ret
}
myCosine(data)
plot(hclust(as.dist(data)))

です.基のデータから単語間の距離に基づいたデンドログラムを構築する方法で何かうまい方法は無いでしょうか?

データフレームの取り出し

akira (2010-06-12 (土) 22:53:35)

以下のデータフレームがあったとして

   a b c d e f g
1  q p q 4 6 7 6
2  q p p 5 4 5 4
3  p p p 5 4 5 7
4  q q q 3 5 4 6

a~cそれぞれの値がpに該当する行を取り出す方法を考えていたのですが、どうしても思いつきませんでした。

[a]
  d e f g
3 5 4 5 7
[b]
  d e f g
1 4 6 7 6
2 5 4 5 4
3 5 4 5 7
[c]
  d e f g
2 5 4 5 4
3 5 4 5 7

どなたかお分かりの方、丸投げで申し訳ありませんがご教示頂けないでしょうか?

配列とGIS(緯度経度)情報の重ね合わせ図

宗二 (2010-06-11 (金) 02:01:30)

シェープファイルと2次元配列から作った、図の重ね合わせ方で悩んでます。
具体的には、日本の海岸線のシェープファイル(ライン)を

plot(シェープファイル, xlim=c(西経度, 東経度), ylim=c(南緯度, 北緯度))

というようにプロットした後、2 次元配列として用意した日本周辺の気温データを、同じ緯度経度の範囲で重ねて作図したいのですが、2次元配列をimage、あるいはcontour関数で作図するときに、x軸とy軸の範囲を緯度経度で指定できないものでしょうか?
配列をプロットするときの軸の範囲が、配列の要素番号と対応するので

image(配列, xlim=c(西経度, 東経度), ylim=c(南緯度, 北緯度))

のようにすると、範囲が大きくずれて何も表示されなくなってしまいます。
なにか良い解決方法ありましたらお願いします。

プログラムの中断・再開

shannon (2010-06-11 (金) 00:54:55)

プログラム実行を途中で中断・再開させる関数、あるいはコマンドはあるのでしょうか?
WindowsのTinn-R(ver:1.19.4.7)から、ショートカットキーでRコンソール(ver:2.11.0)にプログラム全体を送り実行しているのですが、プログラムの途中で、変数の値を確認するために実行を中断し、確認後に実行を再開させる、ということをやりたいのです。
手動で範囲指定してコンソールに送る、以外でいい方法はないでしょうか?

Rによる日別の時系列分析に関する文献を探しています。

T (2010-06-10 (木) 16:22:17)

只今卒論のテーマとしてある飲食店の一店舗の様々な日別データをRによって時系列分析しようと考えていまして、日別データを扱った時系列分析に関する本や文献を探しています。
ご存じの方がいらっしゃれば教えていただけると幸いです。

2次元配列の回転

do-san (2010-06-10 (木) 01:08:00)

2次元のバイナリ(もしくは文字列)配列を、90度、もしくは180度回転させたいのですが、そのような関数はあるのでしょうか?

filled.contourの 重ね描きについて

vp (2010-06-06 (日) 21:35:28)

宜しくお願いします。
以下のように filled.contour で描いた図に abline で x=0.5, y=0.5 の線を重ね描きしようと思いましたが、座標が対応しませんでした。

data(volcano)
filled.contour(volcano, color = terrain.colors, asp = 1) # simple

par(new = T)
abline(h = 0.5, xlim = c(0, 1), ylim = c(0, 1))
par(new = T)
abline(v = 0.5, xlim = c(0, 1), ylim = c(0, 1))

何かうまい方法は無いでしょうか?

GLMのfamily指定によるEstimateの正負の逆転について

mm (2010-06-06 (日) 01:59:41)

こんにちは。最近RでGLMに取り組み始めました。
まずは,説明変数が応答変数によって正,負のどちらに影響しているかを調べることが目的です。
この応答変数は,連続変数で正の値を取るため,familyをGammaにしました。
するとEstimateの値が,正規分布を仮定した回帰分析と逆の正,負の関係になってしまい,解釈に困っています。
下の例のように,irisのデータで,familyをgaussianとGammaでGLMを行っても同じ結果です。
指定するfamilyによって,正負の方向が逆転するということはあるのでしょうか?

test <- glm(Petal.Length ~ Petal.Width, data=iris)            
summary(test)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.08356    0.07297   14.85   <2e-16 ***
Petal.Width  2.22994    0.05140   43.39   <2e-16 ***
#################### 
testG <- glm(Petal.Length ~ Petal.Width, data=iris, family=Gamma)         
summary(testG)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.532340   0.017386   30.62   <2e-16 ***
Petal.Width -0.172682   0.008915  -19.37   <2e-16 ***

Bland Altman plotについて

sakura (2010-06-05 (土) 13:13:57)

最近、調べ物をしていたら、Bland Altman plot のことが目にとまりました。信頼性評価に使われるそうですが、馴染みがありません。Rのlibraryに搭載されたものがあるのでしょうか?差分の傾向を知るというのが、ポイントのようですが・・・。ご教示いただければ幸いです。

while のエラー

while (2010-06-03 (木) 03:44:10)

while を使用し、行列 kara 内に while 内処理を length(mg) 回実行したいのですがエラーが出ます。
(df kara 内の candy, drop という変数(行名)それぞれに対応する colMeans を求め、kara 内に格納する)
エラーは以下のように出ます。

以下にエラー kara[, i] <- g : 
  置き換えるべき項目数が,置き換える数の倍数ではありませんでした 
mg <- c("candy","drop") # データフレーム ame の中の colname を 2 つ選び vector に
kara <-matrix(, length(mg))
i <- 0
while(i <- i+1 <= length(mg)) {
    g <- colMeans(subset(ame, ame[i] == "eat", select = (範囲を指定)))	
    kara[,i] <- g	
}

Kaplan-Meier 曲線が 1 本だけ

?アロバ (2010-06-02 (水) 19:29:27)

たびたび皆様にご指導いただきたく、よろしくお願いします。

source("http://aoki2.si.gunma-u.ac.jp/R/src/km_surv.R", encoding="euc-jp")
# 1 は A 群,2 は B 群を表す
group <- c(1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,1,1,1,1,1,
           1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,1,1,2,2,1,2)
# 1 は死亡,2 は 生存(打ち切り)を表す
event <- c(1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,0,1,
           1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0)
# 生存期間
time <- c(2,84,318,198,198,197,192,306,96,90,88,66,48,264,
          24,5,243,216,42,5,20,0,12,144,48,192,192,228,
          176,180,84,84,123,117,115,267,98,96,86,63,44,41,
          33,20,18,252)
a.group <- group == 1	
km.surv(time[a.group], event[a.group])
library(survival) # survival ライブラリーを使う
dat <- Surv(time[a.group], event[a.group]) # survfit で使うオブジェクトを作る
res <- survfit(dat)

上記にて入力したのですが、Kaplan-Meier 曲線が A 群 1 本だけのグラフになってしまいます。A, B 両群の曲線を一つのグラフに入れたいのですが。。。 よろしくお願いします。

行列にfor文中の変数を含んだ名前をつける

lo (2010-05-31 (月) 00:05:13)

for (i in 1:10){
の中で
計算結果を行列に格納する工程を書いた後
i番目の結果を格納した行列(仮にkekkaという名前とします)に、新たに、
hogeとi(1から10までの値をとる)をくっつけた、hogeiという名前をつけたいです。
この場合、
hogei<-kekka
とやっても、hoge1とならずに、hogeiとなってしまいます。
paste("hoge",i)<-kekkaとやると、エラーになります。
このような、ループの回数に応じた連続した名称をつけるには、どのようにしたらよいのでしょうか?

val.prob{Design}のexampleコードに関する質問

vi (2010-05-27 (木) 19:22:41)

質問です。Designパッケージの説明ページ
http://bm2.genes.nig.ac.jp/RGM2/R_current/library/Design/man/val.prob.html
において、Examplesとして、実際に使用するときのコードが載っています。
この中の、後半部分$Survival analysis examples 中に

# Survival analysis examples
# Generate failure times from an exponential distribution
set.seed(123)              # so can reproduce results
n <- 2000
age <- 50 + 12*rnorm(n)
sex <- factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4)))
cens <- 15*runif(n)
h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
t <- -log(runif(n))/h
label(t) <- 'Time to Event'
ev <- ifelse(t <= cens, 1, 0)
t <- pmin(t, cens)
S <- Surv(t, ev)

# First validate true model used to generate data
w <- val.surv(est.surv=exp(-h*t), S=S)
plot(w)
plot(w, group=sex)  # stratify by sex

上記のような記載があるのですが、私は、est.survというのが、時間tにおける各ケースの計算上の生存率であると理解しています。その場合、

w <- val.surv(est.surv=exp(-h*t), S=S)

の部分の、exp(-h*t)は、比例ハザードを用いて、(時間tにおけるbaseline hazard)^exp(h)とするのが正しいのではないかと疑問に思っています。そうはいっても実際に試すと、きれいなキャリブレーションプロットがかけるのですが。私のest.survの解釈自体が間違っているのでしょうか?どなたか、引数est.survとSについての解釈を教えていただけると幸いです。

デンドログラムの葉の並び替え

jun (2010-05-25 (火) 22:57:47)

クラス "dendrogram" に対するメソッドで葉を並べ替える
reorder(x, wts, ...)
というものがありますが、使い方がよくわかりません。

葉をV1、V2、V3、V4・・・・と昇順に並べたいのですが。
昇順に並べても枝は絶対にクロスしないような距離行列とクラスタリング結果になっています。

wtsとは重みと書いてありましたが何の重みですか?
よろしくお願いします。
使用環境はR version 2.10.0 (2009-10-26)でWindowsVistaです。

1からV1への距離は存在しないので過度に大きな数として99.00を使っています。その他の99.00も同じ意味です。
クラスタリングを単連結で行うため結果に影響は及ぼさないと考えてそのようにしました。

> r
      V1     V2     V3     V4     V5     V6     V7     V8     V9
1 99.000  0.045 99.000 99.000 99.000 99.000 99.000 99.000 99.000
2  0.045 99.000  0.071 99.000 99.000 99.000 99.000 99.000 99.000
3 99.000  0.071 99.000  0.042 99.000 99.000 99.000 99.000 99.000
4 99.000 99.000  0.042 99.000  0.037 99.000 99.000 99.000 99.000
5 99.000 99.000 99.000  0.037 99.000  0.059 99.000 99.000 99.000
6 99.000 99.000 99.000 99.000  0.059 99.000  0.111 99.000 99.000
7 99.000 99.000 99.000 99.000 99.000  0.111 99.000  0.071 99.000
8 99.000 99.000 99.000 99.000 99.000 99.000  0.071 99.000  0.091
9 99.000 99.000 99.000 99.000 99.000 99.000 99.000  0.091 99.000
> plot(hclust(as.dist(r),method="single"))

既定値つきの引数

ms (2010-05-25 (火) 20:56:33)

RjpWikiの『Rの関数定義の基本』を参考にさせて頂き、関数定義で変数を既定値にしているのですが、結果を見ると、既定値になっていないようです。
x[2]+x[3]=1 としたいのですが、この書き方ではだめなのでしょうか?

ご意見頂きたくよろしくお願い致します。
使用環境は R2.10.1, XPです。

fr <- function(x, y=1) {
  
  LL <- 0 
  
  pp <- x[4]*(x[2]*(Data[, 19])^x[1] + x[3]*(Data[, 21])^x[1])^(1/x[1])
  dp <- x[4]*(x[2]*(Data[, 20])^x[1] + x[3]*(Data[, 22])^x[1])^(1/x[1])
  y <- x[2] + x[3]

  Ppp <- exp(pp) / (exp(pp) + exp(dp))
  Pdp <- exp(dp) / (exp(pp) + exp(dp))
  
  Ppp <- (Ppp != 0)*Ppp + (Ppp == 0)
  Pdp <- (Pdp != 0)*Pdp + (Pdp == 0) 
  
  Cpp <- Data[, 4] == 1
  Cdp <- Data[, 4] == 0

  LL <- sum(Cuchi*log(Puchi) + Csoto*log(Psoto))
  return(LL)
}
  • 何を計算しようとしている関数なのか理解できません。何を計算していようがかまわないのだけど,やっていることが何のためにやっているのか分からない。関数が返す LL は return の前にある式で計算され,その計算式で引用される変数はそれ以前のどこにも出てこない(それ以前の計算式は何のために何を計算しているの?)
    「x[2]+x[3]=1 としたい」というのも,よくわからない。x の要素は 4 個?で,x[2] と x[3] を足すと 1 になるような関係式がある?そもそも,x はベクトルで引き渡さないといけないようにも見えない。いちいち x[1] みたいに引用しないといけないし。4 つの要素を x1, x2, x3 x4 として引き渡せば,function(x1, x2, x3=1-x2, x4, y) とすれば,少なくとも x3 は 1-x2 という規定値を持つことにはなるでしょう。x[2]+x[3] が 1 になるような組み合わせは無限にあり,規定値になり得ない。(x2+x3=1 になるような解を求めるというような場合もありますが,そのような場合は「規定値」ではなく,「制約条件」という。概念がまるで違う)。
    ついでながら,引数で渡される y は規定値 1 を持つが,にもかかわらず他の計算に使われる前に y <- x[2] + x[3] と代入されてしまう。これはあなたが言っている 「x[2]+x[3]=1 としたい」というのとは,まるで違う。引数で渡された x[2], x[3] の値を足して,それを y という変数に代入しているだけ。それによって,x[2] や x[3] の値が別のものになるということではない。-- 河童の屁は,河童にあらず,屁である。 2010-05-25 (火) 21:30:01
  • 不勉強で申し訳ありません。
    LL <- sum(Cpp*log(Ppp)+Cdp*log(Pdp)) の間違いです。計算したかったことは対数尤度を計算したのち、optim(初期値,fn,method,…)で関数の最大化を行いたかったのです。そこで、ヘッセ行列を用いるため、x[1]としています。x[2]+x[3]=1というのは、制約条件として入れたいのですが、このような場合、どのようにすれば、この条件を満たす答えが出るか教えて頂ければ幸いです。お手数掛けます。 -- ms 2010-05-25 (火) 22:15:02
  • x[2]+x[3]=1 が制約条件なら,x[3] を求めなきゃ良い。計算の途中で x[3] のところを (1-x[2]) とすればよいだけでは?つまり,求めるパラメータは x[1], x[2], x[4] の 3 個だけということ(順序はつめればよいが)。 -- 河童の屁は,河童にあらず,屁である。 2010-05-25 (火) 22:26:14
  • もし x[1],x[2] 等にさらに正値等の条件が付くなら、線形不等式制約下での最適化関数 constrOptim がよいかもしれません。RjpWiki 中に解説がありますのでキーワード検索。 -- 2010-05-26 (水) 22:26:14
  • ありがとうございます。コメント頂いたとおり、x[2]x[3]は正の条件が付きます。optimで行う方法はありますでしょうか?constrOptimには、hessian はないですよね?質問ばかりですみませんが、よろしくお願いします。 -- ms 2010-05-28 (金) 18:32:13
  • ありがとうございます。コメント頂いたとおり、x[2]x[3]は正の条件が付きます。optimで行う方法はありますでしょうか?constrOptimには、hessian はないですよね?質問ばかりですみませんが、よろしくお願いします。 -- ms 2010-05-28 (金) 18:43:07
  • constrOptim のソースコードの optim を呼ぶところに hessian=TRUE を加える。 -- 河童の屁は,河童にあらず,屁である。 2010-05-29 (土) 09:24:05

lme関数出力のvariance componentの標準誤差

aor (2010-05-21 (金) 11:28:09)

いつも勉強させてもらっています。
nlmeパッケージ、lme関数を使用してマルチレベル分析の勉強をしています。
この出力のうち、 VarCorr関数を使うと分散が出力できますが、その標準誤差を求めたいと思っています。よい方法をご存知の方がいらしたら教えていただけないでしょうか。

小野寺先生らが訳された「基礎から学ぶマルチレベルモデル」 (ナカニシヤ出版) や、石田先生らが訳された「RとS-PLUSによる多変量解析」などで勉強していますが、標準誤差を求める方法はわかりませんでした。
上記の書やインターネットでの検索をすると、「標準誤差はあまり意味がないから信頼区間を使おう」という意見があり、標準誤差の出力方法を記載しているものは見当たりませんでした。
ですが、その意見はその意見として、勉強のため標準誤差の出力方法も知りたいと考えています。
標準誤差を出力するRの方法、あるいは計算式などについて記載されている書籍やサイト (日本語か英語) についてでもお教えいただければ助かります。
どうぞよろしくお願いします。
環境はR2.11.0, Windows Vistaです。

# サンプル (lmeのヘルプから) 
library(nlme)
fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
summary(fm2)
VarCorr(fm2)
intervals(fm2)

行列の要素の比較について

yosshiii (2010-05-20 (木) 16:28:42)

連続で質問です。質問ばかりですみません。
行列で

	01.txt	02.txt	03.txt
A	10	3	7
B	1	0	0
C	10	3	7

の行列があり、”A”という単語がテキスト01に10回, テキスト02に3回, テキスト03に7回出ているという意味の行列です。
そこで「A」と「C」の各テキスト毎の出現回数が等しいとき, 「A, C」をペアとして抽出するプログラムを作りたいんですけど、効率的なやり方が思い浮かびません。意見お願いします。

データフレームの作成

yosshiii (2010-05-19 (水) 12:44:00)

フォルダ内全てのテキストファイルからテキストファイル毎に固有表現の出現頻度をカウントしデータフレーム化を実装したいと思っています。以下にイメージしている実行例を載せます。
フォルダ「test」に01.txt, 02.txt, 03.txtのテキストファイルがありそれぞれ「”アリ” “スズメ”」「”アリ” “リス”」「”アリ” “ハチ”」という固有表現が書かれていたとします。これをデータフレーム化すると

	01.txt	02.txt	03.txt
アリ	1	1	1
スズメ	1	0	0
リス	0	1	0
ハチ	0	0	1

のようにしたいと思っています。
そして今現在「count」というリストにcount1=(1 1 0 0), count2=(1 0 1 0), count3=(1 0 0 1) の数値をいれる所まで出来ているのですが、このリストの「count」を上記のようにデータフレーム化したいです。よって、フォルダ内のテキストファイル数が幾つであっても良いように実装したいので、意見お願いします。

ループを用いたPSファイルの自動生成

だいもん (2010-05-19 (水) 12:05:32)

以下のコマンドをループでまわすことで、ポストスクリプトファイルを自動的に生成しようと試みています。(sheet, ddがループの変数に依存して替わります。)

postscript(sheet, horizontal = FALSE, paper = "special", height = 6, width = 6,colormodel="rgb")
levelplot(dd[, 4] ~ dd[,1]*dd[,2], , color=TRUE)
dev.off()

まずこのコマンドをループを遣わずに実行すると、X11()デバイスにはカラーで表示されるものが白黒のポストスクリプトファイルしか作成されません。

さらにループで複数回まわすと、全てのPSファイルが白紙(ファイルサイズ=5kb)になってしまいます。この問題については、recordPlot(), replayPlot()を用いたコードを試してみても同じ結果でした。(ただしこちらではカラーで出る。)

当方の環境はWindows XPで、R2.11.0を利用しています。お知恵を拝借できれば幸いです。

ラベルに文字列+式

N/A (2010-05-19 (水) 12:02:00)

ラベルで
coefficient β[i]
としたい場合どのようにすればよいでしょうか?

ylab="Regression function expression(beta)[i]"
ylab=parse( "Regression function",expression(beta),"[i]" )

でもできなくて案がつきました・・・
お手数をおかけしますがよろしくお願いします。

R version 2.11.0のsummary.survfit関数

nene (2010-05-17 (月) 16:54:23)

Design, Hmisc, survival,cmprskの各パッケージを読み込んだ、R 2.11.0で、とあるコード(R2.3.1で動作確認されているようです)を読み込んで実行すると、関数 "summary.survfit" を見つけることができませんでした、と表示されてしまい、うまく動きません。RSiteSearch("summary.survfit")を行ってみると、summary.survfit関数はきちんと存在するようなのですが、現在は仕様が変わって別の関数になっているのでしょうか?R2.3.1のころの仕様変更歴なども調べてみたのですが、特に記載がありません。お手数をおかけしますがよろしくお願いします。

クォート内の文字列処理

BURBUR (2010-05-15 (土) 11:55:27)

フォルダ内のファイルを検索したい関数を実装したいのですが

> test3 <- function(x) {
	file<-list.files(path="./x")
	return(file)
}
> test3(フォルダ名)

と実行するとpath名をクオートのでくくっているため

list.files: './x' は読めないディレクトリです
というエラーメッセージが出ます。~
これを防ぐ方法を教えて下さい。~

ディレクトリ内のファイル処理

ヨッシ (2010-05-14 (金) 16:17:19)

perlではディレクトリ内のファイルを取得する場合
opendir my $dir, "ディレクトリのパス" or die "$!";

	while(my $name = readdir $dir){
		if($name =~ /.*\.sgml/){   #ディレクトリ内のファイルの名前
			open F1, "ディレクトリのパス/$name" or die "$!";
			my @lines = <F1>;
			
				}
			}
		}
	}
	closedir $dir;

2つの箱ひげ図を1つのグラフに

アロバ (2010-05-12 (水) 09:45:56)

par(mfrow=c(1,2))
x<-c(1.12,0,6.82,23.9,17.66,2.67,6.09,0.13,5,0.36,2.07,3.1,10.79,57.29,5.14,
     107.91,80.05,0,21.66,45.33,0,57.51,103.7,40.62,110.27,77.03,203.22,85.15)
boxplot(x, names=c("x")) #箱ひげ図
y<-c(0,1.74,5.36,0,1.12,7.7,0.85,26.52,4.52,0.33,42.44,37.23,2.55,2.65,61.59)
boxplot(y, names=c("y")) #箱ひげ図
par(mfrow=c(1,1))

2種類のデータxとyを上記で入力したところ、1つのファイルの中に別々に箱ひげのグラフが出来てしまいます。
質問1;この2つのグラフを共通の縦軸のグラフにするにはどのようにすれば良いか?
質問2;それぞれのグラフの下に"x"、"y"と名前を入れる事は可能でしょうか?

ご多忙のところ大変恐縮です、よろしくお願いします。
使用環境は下記です。

R version 2.11.0 (2010-04-22) 
i386-apple-darwin9.8.0 

locale:
[1] ja_JP.UTF-8/ja_JP.UTF-8/C/C/ja_JP.UTF-8/ja_JP.UTF-8

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

Kaplan meier 生存曲線から任意の時間における生存率

jb3 (2010-05-11 (火) 21:58:13)

2つ質問があります。
datatableに

time  event group
55    0     1
23    1     1
61    1     1
12    1     1
90    0     1

以下同様に200行の生存時間timeと、イベントeventに関するデータがあるとします。
これを、読み込んで、<- datasetとして、Rの解説の本に書かれていたのを真似て、

> fitteddata <- Surv(dataset$time, dataset$event)
> sf <- survfit(fitteddata~group, data = dataset)
> plot(sf)

と実行すると、kaplan meier 生存曲線を描くことはできたのですが
質問1
もともと1グループしかないので、群分けはいらないと思い、上式の~groupを消すと

Survfit requires a formula or a coxph fit as the first argument

というエラーが出ます。単一群の生存曲線を描くのに、~groupは常に設定が必要でしょうか?
質問2
任意の時間における、kaplan meier法による生存率を計算する関数というのはあるのでしょうか?

Rを覚えたてで、2冊の解説書をみながら悪戦苦闘しています。
よろしくお願いします。

Rコマンダーでの層別化

mtk (2010-05-09 (日) 20:52:56)

MS Windows XP sp3上でRcmdrを使っています。
ver.2.8.0では変数をcharacterに変更するだけで層別化の因子と認識されていました。
ver.2.11.0でRcmdr 1.5-4を起動したところ,これでは認識されませんでした。そこで

データ → アクティブデータセット内の変数の管理 → 数値変数を因子に変換

で指定したところ認識されました。これは仕様の変更なのかトラブルなのか分かりません。
もし仕様の変更なのでしたら層別化の便利な使い方があれば教えてください。

R commderの適合性検定

sousou (2010-05-07 (金) 15:55:04)

Rcmdrで適合性検定を行いたいですが、どうやってもうまく行かなかったんですが、どなたが助けていただけませんか?
操作手順は分かりますが(統計量→要約→頻度分布→カイ2乗適合度検定)、例えば、AとBの選択者が同じかとうかを検討することとしますが、Aの選択者は127名で、Bの選択者は71名です。Rcmdrでどのように操作するのかを教えていただければありがたいですが、よろしくお願いします。

シェルスクリプト上でRによる解析

まさお (2010-05-07 (金) 15:51:27)

シェルスクリプト上でRによるさまざまな処理を行うには,どうしたらよいのでしょうか?

#!/bin/csh
R 
でRは起動するんですけれど,そこで止まってしまいます.R内のコマンドはシェルスクリプト内でどのように書けばよいのでしょうか.

linuxでR

さっこ (2010-05-07 (金) 12:07:39)

linux(GNU)にRを無事make installできたのですが(suで,/usr/local/r-2.10.1/に),RとタイプしてもRが起動しません(自分のusernameで).cshrcなどに何らかのパスを通す必要があるのでしょうか?

Rでのファイル操作と繰返し文

きりんさん (2010-05-04 (火) 11:17:02)
 Rで統計解析を行い,その結果をそれぞれの物質名ごとのファイルに格納するちょっとしたツールを作成しています.
これにループ処理を用いて,完成させたいのですがdir.create(“Ala”)などの「“”」に囲まれた部分を配列に置き換えるなどの操作が出来ずに困っています.
 他のプログラミング言語を用いずにRのみで行う方法はないでしょうか?
 拙い質問ではありますが,ご教授のほどを宜しくお願い致します.

データフレーム(ROC_table1.csv)
	Ala	Ethanolamine phosphate	Hypotaurine	Pro	Uridine	Glycolate
C	551.6935	29.13754	16.681	    173.8867	45.10912	50.05477
C	607.481	        0	13.91893	168.1533	42.37981	84.32064
C	506.6699	28.35182	16.64904	159.7178	43.60009	299.7703
C	530.9984	20.93959	11.42694	236.6067	42.33055	0
C	466.6074	19.3733	    14.7419	    176.2133	41.44249	180.9834
C	532.253	    18.66473	12.15251	172.6138	40.66632	181.5452
C	449.6792	9.951401	11.69104	130.8715	47.65735	0
C	403.2724	19.33708	0	    155.0835	47.37078	230.7595
C	423.5579	19.48771	11.05998	140.1291	38.48613	238.7741
C	463.7251	23.312	        0	126.3421	43.01137	248.8566
C	465.0218	28.58722	10.97927	137.6295	37.42445	177.0466
C	485.2845	35.14039	7.268277	143.0229	37.97251	223.5736
D	707.1228	25.1555	    11.25501	237.0575	37.70653	278.061
D	658.9739	31.52341	14.1557	    197.7347	33.66692	312.7555
D	475.4444	31.67994	11.48172	132.5917	41.70935	545.9441
D	613.6305	23.5282	    12.38427	164.597	    41.24212	63.77056
D	545.5581	33.57919	21.21189	174.7159	43.30637	419.7285
D	848.2421	32.92102	16.18306	248.8493	30.06797	402.9233
D	507.675	    25.12478	17.90848	185.7012	36.81457	315.3726
D	528.8479	20.10386	17.13785	196.1179	40.02536	126.3903
#データの読み込み
x <- read.csv("ROC_table1.csv", header=T)
data.frame(x)

#ROCソースコード
source("http://aoki2.si.gunma-u.ac.jp/R/src/ROC.R", encoding="euc-jp")

#患者群データを格納
Dise <- (x$X=="D")
disease.x <- x[Dise,"Ala"]

#健常者群データを格納
Cont <- (x$X=="C")
normal.x <- x[Cont,"Ala"]

#物質名のフォルダを作成
dir.create("Ala")

#新規フォルダへのディレクトリ変更
Direct <- getwd()
Direct1 <- paste(Direct,"Ala",sep = "/")
setwd(Direct1)

#フォルダに計算結果を格納
RR <- ROC0(disease.x, normal.x)
write.table(RR, file="Ala.data", sep="\t", row.names=FALSE, quote=FALSE)

#フォルダに図を格納
bmp(filename = "Ala.bmp", width = 480, height = 480)
ROC0(disease.x, normal.x)
dev.off()

#ディレクトリを元へ戻す
setwd(Direct)

新しい目盛りとかぶる

さつまいも (2010-05-03 (月) 20:23:20)

初めにヒストグラムを作ってから、X軸の目盛りを変えると最初の目盛りとかぶってしまいます。
どうしたら最初の目盛りを消すことができますか?

英語でテキストマイニングをする方法。語の頻度,形態素について

ピースケ (2010-04-29 (木) 17:51:37)

英語でテキストマイニングをする方法を探しています。。語の頻度,形態素について分析するパッケージはありますでしょうか。

R-2.11.0でtcltkウィジェットの日本語が文字化け

せーだ (2010-04-29 (木) 12:19:03)

Win7 Pro/WinXP Pro SP2を使っています。
R-2.11.0で両方とも再現しましたので投稿します。

以下のサンプルは、R-2.10.1までは問題なく動作していました。
しかし、2.11.0では、コンボボックス内の日本語だけが文字化けします。

require(tcltk)
tclRequire("BWidget")
tt <- tktoplevel()
tkgrid(tklabel(tt,text="好きな果物は?"))
fruits <- c("林檎","蜜柑","バナナ","梨")
comboBox <- tkwidget(tt,"ComboBox",editable=FALSE,values=fruits)
tkgrid(comboBox)

OnOK <- function()
{
    fruitChoice <- fruits[as.numeric(
          tclvalue(tcl(comboBox,"getvalue")))+1]
    tkdestroy(tt)
    msg <- paste(fruitChoice,"は美味しいよね",sep="")
    tkmessageBox(title="Fruit Choice",message=msg)
}
OK.but <-tkbutton(tt,text="   OK   ",command=OnOK)
tkgrid(OK.but)
tkfocus(tt)

2.11.0の新機能説明の項目の中に

Package tcltk now sends strings to Tcl in UTF-8: this means
that strings with a marked UTF-8 encoding are supported in
non-UTF-8 locales.

とあったので、これが影響しているのでしょうか。
また、このような文字化けは、BWidgetやTktable等、Tcl/Tkの拡張ライブラリに由来するウィジェット全般に発生しているのではないかと思います(Tktableでも、2.11.0ではセル中の日本語が化けました)。

もし対処法をご存知の方がいらっしゃれば、ご教示頂きたく存じます。
宜しくお願い致します。

エラー処理

josephine (2010-04-28 (水) 18:34:02)

大量のデータをS-W検定するスクリプトを作っています。
不適当な標本データが原因のS-W検定のエラーを上手く処理して結果をファイルに出力させたいのですが、try()では上手くゆきませんでした。

入力が適切な例

> a
[1] 0 0 0 1 0 0 0
> try(shapiro.test(a), silent=TRUE)

        Shapiro-Wilk normality test

data:  a 
W = 0.453, p-value = 4.136e-06

> s_shapiro_a <- try(shapiro.test(a), silent=TRUE)
> s_shapiro_a

        Shapiro-Wilk normality test

data:  a 
W = 0.453, p-value = 4.136e-06

入力が不適切な例

> b
[1] 0 0 0 0 0 0 0
> s_shapiro_b <- try(shapiro.test(b), silent=TRUE)
> s_shapiro_b
[1] "Error in shapiro.test(b) : all 'x' values are identical\n"
attr(,"class")
[1] "try-error"

try-errorの場合にはNAをs_shapiro_bに代入したいのですが、恐れ入りますが、アドバイスを頂けないでしょうか?よろしくお願いいたします。

contourplotを白黒で書きたい

orange (2010-04-26 (月) 16:32:03)

グラフィックス参考実例集にあるcontourplot("Cube Root Ozone (cube root ppb)"と同様な等高線図を,色付き(赤←→緑)ではなく,白黒で描きたいのですが,どのようにすればよいでしょうか.どなたかお教え頂けないでしょうか.よろしくお願いします.

計算結果の抽出方法

uribo (2010-04-25 (日) 22:43:10)

forループで計算を1000回繰り返したのですが、その計算結果の最後の値だけを1000個分抽出する式をつくりたい場合どうすればよいでしょうか?
恐らく意、はじめのforループ計算の値の最後の値のみを出す式を1000回繰り返すとはおもうのですがその最後の値の抽出方がわかりません。もしよろしければどなたかご教授お願いいたします。

多次元項目反応理論のパッケージ

totoro (2010-04-24 (土) 00:33:22)

Rで多次元項目反応理論(多次元IRT)または多次元カテゴリカル因子分析をしたいのですが、パッケージはありますでしょうか?

解析例などを示しているページがありましたら、教えていただけると幸いです。

latticeで縦軸目盛の位置を左にそろえたい

itok (2010-04-19 (月) 20:46:34)

VistaでR-2.8.0を使っております。最近latticeを勉強しております。
教科書の例題に以下のようなヒストグラムがありました。

data(Chem97, package = "mlmRev")
library(lattice)
histogram(~gcsescore | factor(score), data = Chem97)


出力をみると、縦軸目盛が、1行目では右に、2行目では左についています。
他の例題をみても、このように互い違いになるのが標準設定のようです。
希望としては、両方の行で、縦軸目盛を左に統一したいのですが、可能でしょうか。
ヘルプを読んでみましたが、どのパラメータを設定すればよいのかわかりませんでした。
どうかよろしくおねがいいたします。

主成分分析のスコアプロットでグループごとに色分けできますか?

にわか (2010-04-18 (日) 00:35:23)

主成分分析のスコアプロットで、データのグループごとに色分けして表示させたいのですが、可能でしょうか?
irisで
.PC <- princomp(~PL+PW+SL+SW, cor=TRUE, data=iris)
scoreplot(.PC,labels="names",cex=0.7,col=2)
このようにラベル表示はできるのですが、グルーピングが分かりやすいようにデータの1:50と51:100、101:150を色分けして表示できるとありがたいのです。
どなたかご教示頂けたら助かります。
よろしくお願いします。

GLMでLSMEANの計算

Kita (2010-04-14 (水) 15:33:17)

SASの場合、GLM解析を行った場合、Least Square Mean(LSMEAN)が計算されるとしていますが、RでLSMEANを求めるためにはどうすればいいのですか?

よろしくお願いします。

すでにCRANから消えたパッケージをインストールする方法について

Saito (2010-04-08 (木) 23:03:40)

下のほうで、adapt関数以外で重積分をする方法についてで質問した者です。
少し前の話ですので、新しく質問をさせていただきました。

あれから色々と勉強して、以下のようなサイトを見ながらやっていたのですが、どうにも巧くいきません。
http://www.okada.jp.org/RWiki/?RcmdrPlugin%C4%B6%C6%FE%CC%E7
http://www.okada.jp.org/RWiki/?%BB%E4%C5%AA%A5%D1%A5%C3%A5%B1%A1%BC%A5%B8%BA%EE%C0%AE%CB%A1
http://www.okada.jp.org/RWiki/?Windows%A4%C7R%A4%CE%A5%D1%A5%C3%A5%B1%A1%BC%A5%B8%A4%F2%BA%EE%C0%AE%A4%B9%A4%EB

まず今回はすでにあるtarファイルを使いますので、最初から作る必要はないと考え、adaptパッケージのtarファイルをこのサイトから直接C直下にダウンロードしました。そしてコマンドプロンプトに

cd c:\
R CMD INSTALL adapt_1.0-4.tar.gz

としましたが、以下のようなメッセージが出てしまい先に進めません。

c:\>R CMD INSTALL adapt_1.0-4.tar.gz
* installing to library 'C:\Users\Saito\Documents/R/win64-library /2.11'
* installing *source* package 'adapt' ...
** libs
  DLLフャ...
x86_64-w64-mingw32-gfortran     -O2  -c adapt.f -o adapt.o
make: x86_64-w64-mingw32-gfortran: Command not found
make: *** [adapt.o] Error 127
   ... ョケ
ERROR: compilation failed for package 'adapt'
* removing 'C:\Users\Saito\Documents/R/win64-library/2.11/adapt'
* restoring previous 'C:\Users\Saito\Documents/R/win64-library /2.11/adapt'


これはフォートランの命令が読めない、というメッセージなのでしょうか?文字化けも気になるのですが・・・。ちなみにPathは、

c:\>Path
PATH=c:\Rtools\bin;c:\Rtools\perl\bin;c:\Rtools\MinGW\bin;
C:\Program  Files (x86)\HTML Help Workshop;
C:\Program Files\R \R-2.11.0alpha-x64\bin;
C:\Program Files (x86)\MiKTeX 2.8\miktex \bin;
C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;
C:\Windows\System32\WindowsPowerShell\v1.0\

となっています。

ちなみにOSはWindows7、Rのバージョンは、R-2.11.0alpha-x64となっています。何か64bit版特有のことをしないといけないのでしょうか。fortranは持っていませんので、買えと言われても少し困ってしまいます。
自分では手詰まりになってしまって、解決できませんでした。もし何か解決策をご存知でしたら、教えていただけると幸いです。

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

もし何かお気づきの点がありましたら、ご助言いただけると幸いです。 -- Saito 2010-04-09 (金) 00:12:37

相関行列からクラスタ解析

KM (2010-04-07 (水) 16:28:46)

相関行列を類似度行列としてクラスタ解析したいと考えています。
hclust()にas.dist(x)という形で読み込みたいのですが、hclustに入力する距離行列は非類似度行列(値が小さい因子から統合される)とのことです。
−1x相関係数を非類似度として読み込んでみましたが、hang=-1にすると因子ラベルが0に揃ってしまいます。

以下の2点いずれかができれば解決すると思うのですが、いかがでしょうか?
・類似行列のクラスタ解析ができる方法を実行
・hclustでラベルを揃える位置を0以外の位置に指定する

よろしくお願いします。

optim関数について

あ〜る (2010-04-07 (水) 11:19:57)

Rのoptim関数と同じ関数をjavaで計算したく、javaのlibraryで探しています。
また、javaで完結しているコード(他の言語を呼び出したりしない)を探しています。

実際に計算したいのは、L-BFGS-B法ですが、javaのlibraryでは見つけられませんでした。

もしご存知の方がいらっしゃったらご教示ください。
よろしくお願いいたします。

maptoolsを使った地図グラフ作成について

kojiro_i619 (2010-03-31 (水) 13:50:28)

お世話になります。中澤先生の本から、以下

gnm <- readShapePoly("gunma.shp")
gnmdata <- gnm$att.data
gnmpoly <- map2SpatialPolygons(gnm, region.id=att(gnm, "region.id"))
## <−これがうまくいきません。
gunmadata <- gnmdata
gunmapoly <- gnmpoly
aged <- read.delim("agedprop.txt")
gunmadata <- merge(aged, gunmadata, sort=F, by="JCODE")
DD <- gunmadata$AP2006
classes <- cut(DD, seq(min(DD), max(DD), length=5), include.lowest=T)
cols <- topo.colors(4)
plot(gunmapoly, col=cols[ordered(classes)], xlab="", ylab="", axes=F)
legend(max(x)-0.3*(max(x)-min(x)), min(y)+0.1*(max(y)-min(y)),
       legend=names(table(classes)), cex=0.6,fill=cols)
title("群馬県市町村の65歳以上高齢者割合4区分(2006年)")
text(x, y, gunmadata$CITY1, cex=0.5, pos=1, offset=0)

などと、作ったのですが、途中の命令が、ないようです。
見本のコードでも、教えていただければ、幸いです。

forループの高速化

rbegginer (2010-03-30 (火) 21:27:26)

いつも勉強させていただいております。
以下のプログラムについて、高速化のコツ、ヒントなどをいただけないでしょうか。

以下のようなプログラムで、1千万の母集団から5サンプルの抽出を1000回行い、外れ値を除去して平均値を計算する、ということを行っております。
5サンプルが終わったら6,7,8...100とサンプルサイズを変えてまた平均値を計算...というようにしています。さらに基準や数値を変えて実行しようと思っています。

実はサンプル抽出を1000回ではなく10000回行いたいのですが、このプログラムが大変時間がかかり、5000回でやってみると50分程度かかります。
時間がかかっていると部分、高速化できるような書き方がありましたら教えていただけないでしょうか
こちらのTipsでforは時間がかかるという記述を読みましたが、sapplyにする方法も思いつかず悩んでおります。
どうぞよろしくお願いします。環境はWindowsVista, R2.10.1です。

pop <- round(rnorm(10000000, mean=300, sd=20)+(300*rexp(10000000)))
ssz <- c(5,6,7,8,9,10,15,20,25,35,50,100) # 抽出するサンプルの数
nrecv <- 0
nremv <- 0

for(j in 1:length(ssz)) {
   for(i in 1:1000){
      svct <- sample(pop, ssz[j])
      gmean <- mean(svct)
      gsd <- sd(svct)
      cfv.u <- gmean+(3*gsd) # 外れ値のカットオフポイント設定
      cfv.l <- gmean-(3*gsd) # 外れ値のカットオフポイント設定
      tmp <- svct
      tmp[tmp>cfv.u] <- NA # 外れ値除去
      tmp[tmp<cfv.l] <- NA # 外れ値除去
      svct.el <- tmp
      nrecv[i] <- mean(svct.el, na.rm=T) # 外れ値除去後サンプルの平均
   }
   nremv[j] <- mean(nrecv) # 1000個の平均値の平均
}

岡田昌史編 The R Book CD-ROM について

kojiro_i619 (2010-03-29 (月) 13:10:05)

CD−ROMの内容をロードするWEBサイトがありましたら、教えてください。

anova(lm(count~spray))でdim(X) must have a positive length エラー

にわか (2010-03-28 (日) 00:06:36)

winXPでR2.10.1を使っています。
線形モデルの結果にaov、anovaを用いた分散分析を行うと、anova行うとエラーがでます。aovではでません。
いろんな所(たとえばhttp://www1.doshisha.ac.jp/~mjin/R/15.html)に、anovaで行う事例が載っていますので、当方の何かに問題があると思うのですが、見当がつきません。
漠然とした質問で恐縮ですが、どのようなところに問題ありそうでしょうか。
よろしくご教示いただけたら幸いです。

mac osx 10.6.2 で RPostgreSQLパッケージをロードできません。

jiro (2010-03-25 (木) 23:11:55)

OS: mac osx 10.6.2( intel mac core 2 duo )
R ver: 2.10.1(32bit、64bit ともにインストール)
PostgreSQL ver: 8.3.1

上記の環境にて、R言語からPostgreSQLに接続しようとしています。
DBIのver 0.2-5 を install.package(DBI) コマンドにてインストールし、シェルからコマンドをたたき、RPostgreSQL をインストールしました。

$> sudo R CMD INSTALL RPostgreSQL_0.1-6.tar

途中、

rawToChar(magic[1:5])中で警告がありました:
文字列をnullに切り詰めました: 'Rpo\0\0'

という警告がでたもののインストールには成功したようです。
さて、Rのコンソールから、

> library( RPostgreSQL )

と実行すると下記のメッセージが表示され、パッケージのロードに失敗します。
64ビット版で失敗しているようなので、この問題の解決方法か、32bit 版の RPostgreSQL のインストール方法をお教え頂けないでしょうか。

よろしくお願いします。

Error in dyn.load(file, DLLpath = DLLpath, ...) : ~
  共有ライブラリ '/Library/Frameworks/R.framework/Resources/library/RPostgreSQL/
     libs/x86_64/RPostgreSQL.so' を読み込めません 
 dlopen(/Library/Frameworks/R.framework/Resources/library/
     RPostgreSQL/libs/x86_64/RPostgreSQL.so, 6):
     Symbol not found: _PQbackendPID
 Referenced from: /Library/Frameworks/R.framework/Resources/library/
     RPostgreSQL/libs/x86_64/
   RPostgreSQL.so
 Expected in: flat namespace
in /Library/Frameworks/R.framework/Resources/library/RPostgreSQL/libs/x86_64/RPostgreSQL.so 
エラー:  'RPostgreSQL' に対するパッケージもしくは名前空間のロードが失敗しました

PLS回帰用のNIRデータの構造(再)

にわか (2010-03-24 (水) 16:07:09)

標記の題名で11でお世話になった者です。
クリアに解決しましたので、投稿しておきます。
説明変数がマトリックスの場合には、protect関数(?) のI()を使えば、yarnのデータ構造のように、列毎に別々に変数を指定せずにまとめることができ、しかもデータフレームとしておけるようです。
例)

y <- c(2, 4, 6)
z <-matrix(1:12, 3, 4)
mydata <- data.frame(Y=y, Z=I(z))

以上、お世話になりました。

コンソールの内容を保存するコード

rbeginner (2010-03-20 (土) 14:07:10)

いつも勉強させていただいています。

コンソールに表示される命令と出力を全て保存する関数はあるのでしょうか?
メニューからファイル -> ファイルを保存 を選ぶとlastsave.txtというファイルとして保存されますが、これをコードで行いたいと考えています。
理想はsink(... append=T, split=T) のように逐一保存したいのですが、sink関数だと出力しか保存されないようです。

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

GLMによる予測値の予測区間の求め方

saka (2010-03-19 (金) 13:29:49)

GLMでガンマ分布とlogリンク関数を使って,将来の魚の漁獲量の予測値とその予測値の予測区間を求めたいのですが,予測区間をどのように求めたら良いのか分かりません。

推定値の信頼区間ではなく,それより広い範囲になることが多い予測区間の方です。

来年の予測漁獲量は○トンくらいだけど,95%の確率で○〜○トンの範囲になるということを知りたいのです。

GLMで解析するデータは下のようなものです。

これに,例えば新たに説明変数1,2,3がそれぞれ(12,5,1500)が得られたときに,予測漁獲量を「predict」を使って求めるところまではできるのですが,予測区間の求め方を教えてください。
よろしくお願いします。

目的変数 説明変数1 説明変数2  説明変数3
  7553       11.1       0.6        2240
 14897       13.0       1.6        1292
 30957       12.6       6.2         972
 35056        9.7      26.2        2655
 17436       11.4       9.4         547
  4656       14.0       0.2         820
  6879       12.3       0.0        1223
 41853       13.9       4.4        2043
 16308       12.1       4.8        1875
 11293       11.8       3.0        1363
  3964       10.2       9.8         236
  7993       13.1       0.5        1236
  6748       11.5       0.2        2144
  3522       12.2       0.0         758
 15395       10.5       7.6        1601
  7567       11.3       6.0        1175
  5823       10.7       4.4         880

以上です。

対応のある比率の検定?

にわか (2010-03-17 (水) 23:31:15)

ある植物を植えたときに、処理Aと処理Bで活着率に差があるかどうかを、時期を変えて調べた次のようなデータがあります。
各時期に処理A,Bとも30株ずつ植えたうち、表の株数が活着し、残りは枯死しました。

<各30個体を植えたときの活着した株数>
      時期1  時期2  時期3  時期4
処理A   14株   13株   10株   15株
処理B   15株   14株   13株   16株

この場合、処理A、Bのどちらが活着に有効かを検定する場合、各時期ごとには活着の成功か失敗かという問題だと思うので、
prop.test(c(10,19),c(30,30))
というように比率の検定をすればいいと思うのですが、
全体で見る場合、単に合計で 52/120 と 58/120 の比率の検定をするのか、時期を反復と見て対応のあるt検定を行ってよいのか悩んでいます。
あるいは、「対応のある比率の検定」のようながあるのでしょうか。
どなたかご教示下さい。
よろしくお願いします。

adapt関数以外で重積分をする方法について

Saito (2010-03-13 (土) 15:10:03)

あちこち調べ、自分でも試したのですが、わからなかったので質問させてください。

最近(といっても数ヶ月前)にCRANからadaptパッケージが消えたようです。
今まで重積分をするときにはadaptを使って積分していたので、消えた後も過去のRのバージョンのlibraryからadaptを持ってきて使っていました。

しかし、2.10.0以降になると以下のようにre-installを求められます。

> library(adapt)
エラー: package 'adapt' was built before R 2.10.0: please re-install it


re-installしようにもCRANにないのでできません。そこで仕方なく他の方法で重積分ができないか探したのですが見当たりませんでした。

どなたか、解決策をご存知でしたらご教授願えないでしょうか。
なお、環境はWindows 7、x64 R-2.11.1 Pre-releaseです。

prcompを用いた固有値・固有ベクトルの計算について

初心者 (2010-03-10 (水) 00:24:41)

prcompを使用して固有値・固有ベクトルを算出しているのですが、どうも結果が意図したものと異なってしまいます。
・Rの出力

> A <- matrix(c(0, -4, 4, 1, 4, -3, 1, 2, -1), nrow=3)
> A
     [,1] [,2] [,3]
[1,]    0    1    1
[2,]   -4    4    2
[3,]    4   -3   -1
> prcomp(A, scale=F, center=F)
Standard deviations:
[1] 5.579673e+00 9.312648e-01 1.184667e-16

Rotation:
            PC1        PC2        PC3
[1,] -0.7093545 -0.6210515 -0.3333333
[2,]  0.6431273 -0.3767530 -0.6666667
[3,]  0.2884500 -0.6872788  0.6666667

・参考にしたURL
http://www004.upp.so-net.ne.jp/s_honma/urawaza/eigenvector.htm
このページによると固有値は0,1,2となるはずなのですが、何か使い方等で間違っているのでしょうか?
Windows XPで R2.8.1 を使用しています。

Rの正規表現について

初心者 (2010-03-10 (水) 00:04:04)

Windows VistaのR-2.10.1です.(Rの他のバージョン,他のOSでもうまく動いて欲しいです)

例えば,「何章」や「何節」もしくは「番号のみ」(いずれも1〜2桁で,全角の数字の可能性もあります)から番号だけにするようなとき,

section <- c("1", "23", "4部", "56部", "7章", "89章")
section <- c(section, "1", "23", "4節", "56節")
sub("([[:digit:]]{1,2})[^[:digit:]]?", "\\1", section)

でしてみましたが、結果が

 [1] "1"    "23"   "4部"  "56"   "7章"  "89"   "1"   "23" "4節" "56"

のようになってしまいました.
どのようにしたら,うまく番号だけ抜き出せるのでしょうか? また何が原因なのでしょうか?? 宜しくお願いいたします.

listオブジェクトのgrepについて

R初心者 (2010-03-09 (火) 10:36:33)

R-2.10.1をWindows Vistaで使用しています。

mylist <- list(alph=letters, Alph=LETTERS, num=as.character(0:9))
grep("A", mylist, ignore.case=FALSE) # (1)
grep("C", mylist, ignore.case=FALSE) # (2)
grep("c", mylist, ignore.case=FALSE) # (3)

ここで、(1)と(2)などpatternが(3)の "c" 以外のgrepの挙動は私の望んだものとなっています。
しかし、(3)のgrepの挙動がおかしい(私が望んでいるものではない)です・・・。
本来ならば 1 が帰ってきて欲しいのですが、 c(1,2,3) が帰ってきてしまいます。
1 と帰ってくるような手立てはございますでしょうか? また、これはバグなのでしょうか?
grep以外の関数での解決方法でも結構です。よろしくお願いいたします。


添付ファイル: fileimage.png 1705件 [詳細] filehenokappa.png 1872件 [詳細] filer-cluster.png 1893件 [詳細] filescatterplot3d2.R 1388件 [詳細] fileomega-alpha.png 1704件 [詳細] fileexample999.png 1743件 [詳細] filesamplelengtherror.xls 1306件 [詳細] filetree.gif 1599件 [詳細] filejhin2.png 1900件 [詳細] fileplot3d.png 1836件 [詳細] filest20100614-2.png 1879件 [詳細] filemiss_filled_contour2.jpg 1861件 [詳細] fileglid.png 1831件 [詳細] filematplot.png 1820件 [詳細] filesampleperlerror.xls 1289件 [詳細] fileinstall-adapt.png 1843件 [詳細] filevennDiagram.png 1889件 [詳細] filescatterplot3d.png 1868件 [詳細] file3dplot.png 1757件 [詳細] filereorder.png 1780件 [詳細] filetest.txt 412件 [詳細] filevennDiagram2.png 1611件 [詳細] fileeval-parse.png 1794件 [詳細] fileapprox.png 1721件 [詳細] filejhin.png 1718件 [詳細] filemiss_filled_contour.jpg 1847件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2015-03-01 (日) 01:15:59