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

スレッド

Q&A (初級者コース) 過去の目次


●● 目次 ●●


svmでplotした結果からsupport vectorを取り除き等高線のみにしたい

ひろ (2013-02-17 (日) 23:00:54)

kernlabパッケージのksvmを使用した下記コードを実行してplotすると確立分布図が現れますが,その図からサポートベクター(SV)(○・●,△・▲のプロット)を取り除く方法はありますでしょうか?(e1071パッケージでもoKです)

library(kernlab)
# number of data points
n <- 150
# dimension
p <- 2
sigma <- 2  # variance of the distribution
meanpos <- 0 # centre of the distribution of positive examples
meanneg <- 3 # centre of the distribution of negative examples
npos <- round(n/2) # number of positive examples
nneg <- n-npos # number of negative examples
# Generate the positive and negative examples
xpos <- matrix(rnorm(npos*p,mean=meanpos,sd=sigma),npos,p)
xneg <- matrix(rnorm(nneg*p,mean=meanneg,sd=sigma),npos,p)
x <- rbind(xpos,xneg)
# Generate the labels
y <- matrix(c(rep(1,npos),rep(-1,nneg)))
# Train a nonlinear SVM with Gaussian RBF kernel and default parameters
svp <- ksvm(x,y,type="C-svc",kernel='rbf')
plot(svp,data=x)

scatter3dの軸の設定ができない

rkoma (2013-01-11 (金) 20:22:54)

R初心者です。よろしくお願い致します。
Rで作成する3次元散布図 (scatter3d)では、軸タイトルなどへのExpression関数への適用や、軸の最大・最小値の設定や、軸を対数目盛にすることは可能でしょうか?
Rcmdr の「3次元散布図」で作ったグラフをひな型にして、スクリプトウィンドウに追加して修正を行っています。具体的には、Rコマンダーから3次元散布図を作成して、

scatter3d(x軸の値, y軸の値, z軸の値, surface=FALSE,
          residuals=TRUE, bg="white",
          axis.scales=TRUE, grid=TRUE, ellipsoid=FALSE,
          xlab=" ", ylab=" ", zlab=" ")

となったところに log="x", log="z", xlim=c(0.1,100), zlim=c(0.1,100) およびExpression関数によるx,y,zlabのテキストの修正(上付き、下付き文字の追加)を加えています。特にエラーメッセージは出ないのですが、できませんでした。
2次元散布図では可能です。また、plot3d()では、最大値、最小値の設定はできたのですが、対数軸への変更はできませんでした。
 R初心者、プログラミング初心者なので、全く原因がわかりません。単純に機能がないというだけでしょうか?Help()やRjpWikiの過去の質問なども含め、できるだけ調べてみたのですが、よくわかりませんでした。
 また、scatter3d()やplot3d()に限らず、「最大値・最小値を設定できて軸を対数目盛りにできる3次元散布図」を作成する方法がありましたら、教えていただけると助かります。
OSはWindows2003、Rはver.2.15.2、Rcmdrはver.1.9-2を使用しています。
ちなみに、類似の質問をYahoo知恵袋にも投稿しています。そちらで解決しましたら、その旨ご報告します。
恐縮ですがご回答よろしくお願い申し上げます。

大学でパッケージがインストール出来ない

taipapa (2013-01-11 (金) 18:21:30)

某大学のパソコン(マック mountain lion)にR 2.15.2をインストールしました.その後にパッケージを色々入れようとすると,どのパッケージでも,以下の様なメッセージでうまくいきません.

URL 'http://essrc.hyogo-u.ac.jp/ 省略 /pls_2.3-0.tgz'
    を試しています 
Content type 'text/html' length 182 bytes
開かれた URL 
==================================================
downloaded 182 bytes
エラー:  ファイル ‘/var/folders/5j/省略/pls_2.3-0.tgz’
                は OS X バイナリパッケージではありません 
追加情報:   警告メッセージ: 
'tar' がゼロでない終了コード 1 を返しました  
tar: Unrecognized archive format
tar: Error exit delayed from previous errors.

なにをしても182バイトしかダウンロードしてくれません.

大学のネットワーク管理者に確認してもプロキシは使用していないとのことですので,プロキシの問題ではないと思います.また,幾つかのcranのミラーサーバーを試してみましたが,全てダメです.Bioconductorのパッケージをインストールしようとしても同様の症状です.マックのユーザーインターフェースの一つであるRパッケージインストーラでは,パッケージの一覧の取得は可能です.ですから,サーバーにはつながって中身も見えています.なのに,選択してインストールしようとすると上記のエラーでことごとくダメです.自宅ではなんの問題もなくR packageのインストールができているだけに,ほとほと困っております.うちの大学側のネットワーク絡みの問題なのでしょうが,かなり調べても同様の報告に行き当たることができず,こちらに相談する次第です.どなたかアドバイスいただければありがたいです.

コードがハイライトされた状態で印刷したい

たく (2013-01-11 (金) 11:21:32)

エディタでコードを色づけしてくれるものはいろいろありますが、ハイライトされた状態のカラーのまま印刷できるエディタや方法をご存じでしたらどなたかお教えいただけないでしょうか。環境はWin7を常用しています。

年次別に違うプロットマーカーを使用したい

Montecarlo (2013-01-07 (月) 03:15:01)

何度か試行錯誤を重ねてみましたが、うまくいかなかったため皆様のお知恵を拝借したく投稿します。
よろしくお願いいたします。

セッションインフォメーションは以下の通りです。

> sessionInfo() 
R version 2.15.2 (2012-10-26)
Platform: x86_64-w64-mingw32/x64 (64-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     

今したいことは題名にもあります通り、年次別に違うプロットマーカーを使用したいと思っています。
現在は以下のプログラムを走らせています。

year <- c(rep(2008, 10), rep(2009,10), rep(2010, 10), rep(2011,10))
x <- rnorm(40)
y <- rnorm(40)
data <- data.frame(year, x, y)

plot(data$x, data$y, type = "n")

temp <- subset(data, data$year == 2008); points(temp$x, temp$y, pch = 0)
temp <- subset(data, data$year == 2009); points(temp$x, temp$y, pch = 1)
temp <- subset(data, data$year == 2010); points(temp$x, temp$y, pch = 2)
temp <- subset(data, data$year == 2011); points(temp$x, temp$y, pch = 5)

しかし、最後の4行を関数化しようと思い以下の通り変更してみましたが、すべてのデータが全てのプロットマーカーで上書きされていくだけで思い通りにいきません。

point.add <- function (year) {
   temp <- subset(data, data$year == year)
   if (year == 2008) { pch = 0 }
   else if (year == 2009) { pch = 1 }
   else if (year == 2010) { pch = 2 }
   else if (year == 2011) { pch = 5 }
   points(temp$x, temp$y, pch =pch)
}

point.add(2008)
point.add(2009)
point.add(2010)
point.add(2011)

うまくいかない原因がわからないので、どなたかご教授いただけないでしょうか。
よろしくお願いいたします。

巨大なcsvファイルの一部を読み込みたい

HIROi (2013-01-04 (金) 02:48:26)

初めて質問させていただく初心者です。read.csvを使って、カラム数23,行数600万件程度のcsvファイルを読み込もうとしたら、memory sizeの限度らしく失敗しました。実行結果(一部省略)は以下の通りです。

EQ <- read.csv("VeryLargeData.csv")
エラー:  サイズ 45.4 Mb のベクトルを割り当てることができません 
追加情報:   14 件の警告がありました (警告を見るには warnings() を使って下さい) 
> warnings()
警告メッセージ: 
1: In scan(file, what, nmax, sep, dec, quote, skip, nlines,  ... :
  Reached total allocation of 1535Mb: see help(memory.size)
 # 中略・・・全く同じメッセージの繰り返し
9: In type.convert(data[[i]], as.is = as.is[i], dec = dec,  ... :
  Reached total allocation of 1535Mb: see help(memory.size)
 # 後略・・・全く同じメッセージの繰り返し、14.まで。

使用したいカラムは「キーとなるコード1カラム+倍精度小数2カラム=合計3カラム」のみなので、「列を間引いて」読み出せれば用は足ります。しかし、最も柔軟にファイルを読み込めると聞くscanのヘルプを見ても、「最初の数行を間引く」方法はあっても、「カラムのうちいくつかだけを読み出す」方法というのはどうも見当たらないように思われます。どなたかご助言をお願いいただけないでしょうか。何卒よろしくお願いします。
投稿における注意事項

R言語での多クラス(マルチクラス)分類 boostingについて

tomo (2013-01-02 (水) 13:00:50)

閲覧ありがとうございます。

私はR言語で学習をしている初心者です。
その中で、irisデータを用いてコマンドの学習をしているのですがエラーが出てしまい先に進めなくなってしまいました。
以下にコマンドを載せますものが私の用いたものです。

library(rpart)
library(adabag)
data <- iris
ndata <- nrow(data)
#乱数指定
set.seed(101)
#学習データ(data.learn)とテストデータ(data.test)に分ける
ridx <- sample(ndata, ndata * 0.5)
data.learn <- data[ridx,]
data.test <- data[-ridx,]
#3-fold crossvalidation、弱識別器の数を10とし、
#学習データに対しboostingを行ったものをdata.adaCvに代入
data.adaCv <- boosting.cv(Species ~ .,data = data.learn,v=3, mfinal = 10)
#テストデータに学習データを照らし合わせる。
resultPredict <- predict(data.adaCv, newdata = data.test, type="class")

というコマンドですが「predict」を行った際に

以下にエラー UseMethod("predict") : 
'predict' をクラス "character" のオブジェクトに適用できるようなメソッドがありません

というエラーが出てしまい、ネットや参考文献を見てもわからず投稿させていただきました。
このエラーの解決方法がお分かりのかたがいらっしゃいましたらよろしくお願いします。
お分かりの方いらっしゃいましたら例文等を添えてコメントいただけたら幸いです。よろしくお願いします。

タイムゾーン無指定による時間差計算の誤差

S (2012-12-28 (金) 11:29:56)

質問を書き込んでいるうちに自己解決しましたので、質問というより、同じ失敗をする人がでないように報告します。

> Sys.Date()
[1] "2012-12-28"

今日は12月28日です。1月27日まで何日あるのか計算しました。

> difftime("2013-01-27", Sys.Date())
Time difference of 29.625 days
> difftime("2013-01-27", "2012-12-28")
Time difference of 30 days

Sys.Date()を使った方はきっちり30日になりません。

> str(Sys.Date())
 Date[1:1], format: "2012-12-28"

Sys.Date()の出力はDateクラスで"2012-12-28"を持っているだけです。

> difftime("2013-01-27", as.Date("2012-12-28"))
Time difference of 29.625 days

2012-12-28をDateクラスにすると再現します。units = "days"を付けても変化はありません。

> difftime("2013-01-27", as.Date("2012-12-28"), tz = "JST")
Time difference of 30 days

タイムゾーンを指定すると30日に戻りました。バックエンドで時間差計算を行ってそのままモデルに放り込んでいるとなかなか気がつかないエラーだと思います(もしかすると、私だけ?)。

スクリプトファイル(拡張子R)の関連付け

yoro (2012-12-27 (木) 09:03:34)

お世話になります。
Windows7 (32或いは64bit)にて R 2.15.1 を使いっています。

お聞きしたいのは、hogehoge.R というようなスクリプトファイルをダブルクリックすることで、Rを起動させスクリプトを実行させる事は可能でしょうか?(バッチモードとしてではなく)

宜しくお願いします。

交互作用の自由度が統計ソフトによって異なるのはなぜでしょうか

こさく (2012-12-26 (水) 18:42:17)

いつもお世話になっていますが、初めて投稿します。統計も素人ですがご容赦ください。

2群(a、b群各30名)で、それぞれ直線回帰できるようなデータがあるとします。実際のデータでは平均値には有意な差はなく、a、b群での交互作用(傾きの差)が主な関心事です。

下記は作ったデータですが、yをxとgrpでモデル化しようとしています。下記のモデル1-4のうちモデル1を採用したいと思いましたが、残念ながらgrp:xの有意差が出ませんでした。

ただ結果を見るとgrp:xの自由度が2になっています。浅学ですが、下記モデル2に示されるようにx、grpの自由度がそれぞれ1で、grp:xの自由度も1*1=1でよいのではと思ってしまいます。

他の統計ソフト(JMP)で解析すると、モデル1におけるgrp:xの自由度は期待通り1になっていました。近くの統計の先生にもモデル1におけるgrp:xの自由度は1でいいのではないかと言われました。

都合のいい結果を採用してもよいのですが、もう少し掘り下げることができればと思い投稿させていただきました。

何かコメント等ご指摘いただければ幸いです。

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

set.seed(1234)
# grp "a" のデータ
x1 <- runif(30)
y1 <- runif(30)+x1*0.3

# grp "b" のデータ
x2 <- runif(30)
y2 <- runif(30)-x2*0.3

x <- c(x1, x2)
y <- c(y1, y2)

# grp (2群のラベル)
grp <- factor(c(rep('a', 30), rep('b', 30)))

library(car) # type3 を指定するのは、問題とは関係ないですが

# ここからが問題
# モデルの固定因子を変えると、交互作用項の自由度が変わる

# モデル 1、grp:x の自由度 2
> Anova(lm(y ~ grp+grp:x), type=3) 
Anova Table (Type III tests)

Response: y
            Sum Sq Df F value    Pr(>F)    
(Intercept) 1.5509  1 20.7241 2.904e-05 ***
grp         0.1585  1  2.1178   0.15118    
grp:x       0.3955  2  2.6423   0.08006 .  
Residuals   4.1908 56                      

# モデル 2、grp:x の自由度 1
> Anova(lm(y ~ x+grp+grp:x), type=3) 
Anova Table (Type III tests)

Response: y
            Sum Sq Df F value    Pr(>F)    
(Intercept) 1.5509  1 20.7241 2.904e-05 ***
x           0.3939  1  5.2640   0.02554 *  
grp         0.1585  1  2.1178   0.15118    
x:grp       0.2338  1  3.1245   0.08257 .  
Residuals   4.1908 56                      

# モデル 3、grp:x の自由度 1
> Anova(lm(y ~ x+grp:x), type=3) 
Anova Table (Type III tests)

Response: y
            Sum Sq Df F value    Pr(>F)    
(Intercept) 2.0111  1  26.357 3.573e-06 ***
x           1.2837  1  16.823 0.0001318 ***
x:grp       2.1557  1  28.251 1.843e-06 ***
Residuals   4.3493 57                      

# モデル 4、grp:x の自由度 2
> Anova(lm(y ~ grp:x), type=3) 
Anova Table (Type III tests)

Response: y
            Sum Sq Df F value    Pr(>F)    
(Intercept) 2.0111  1  26.357 3.573e-06 ***
grp:x       2.4071  2  15.773 3.533e-06 ***
Residuals   4.3493 57
  • それぞれ,Anova(lm(y ~ ほげほげ), type=3) というのを summary(lm(y ~ ほげほげ)) というようにやってみれば,どのようなパラメータ,幾つのパラメータが推定されているかがわかり,その数と Df がどういう関係になっているかがわかるでしょう。 -- 河童の屁は,河童にあらず,屁である。 2012-12-26 (水) 19:33:36
  • コメントありがとうございます。summaryで出てくるパラメータの数とDfが対応していることがわかりました。なぜモデルにより交互作用項のパラメータの数が変わるかが謎ですが、とっかかりができたのでもう少し勉強してみます。またわからなかったらすいませんがアドバイスお願いいたします。ありがとうございました。 -- こさく 2012-12-27 (木) 16:43:49

最尤法による共分散構造分析

りすりす (2012-12-26 (水) 12:17:20)

お世話になっております。

windowsでRを使用しています。
Rで、最尤法による共分散構造分析ができるパッケージはどれなのでしょうか?
semパッケージは、デフォルトで最尤法なのでしょうか?

どなたかわかるかたがいらっしゃれば、ご教授をお願いいたしす。

boxplotの図に斜線を入れたい

Saito (2012-12-21 (金) 21:29:04)

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

boxplot()の図に斜線を入れたいのですが、どうすれば良いでしょうか。boxplot(runif(100), runif(200), runif(300))のそれぞれのグラフで、斜線の種類も変えられると嬉しいですが、とりあえず斜線が引きたいです。

どなたかご存知の方がいらっしゃいましたら、ご教授をお願い致します。

barplotにおけるy軸の値

BCD (2012-12-19 (水) 23:02:56)

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

barplot関数を使ってグラフを作成した時にy軸に自動的に付加される値を取得したいです。
というのも小数点が多くつく場合、y軸が見づらいので10-3乗として表記をしたいからです。
barplotを使う前に元のデータに10の3乗をかけ算しても良いのですが、上記の方法は無いのか知りたいのです。

どうかご教示願います。

mean関数

sori (2012-12-19 (水) 14:24:56)

データ解析の勉強を始めたのですが、なぜ警告メッセージが出るのかがわからないので困っています。どなたかよろしくお願いします。

> x <- read.delim("clipboard")
> x
   学生名 性別 順位 国語 数学 英語
1    佐藤    1    5   64   48   78
2    鈴木    1    6   51   65   62
3    高橋    1    4   57   78   68
4    田中    1   14   38   62   42
5    渡辺    0   15   43   78   57
6    伊藤    0   13   52   73   53
7    山本    0    3   58   45   50
8    中村    0    1   72   36   71
9    小林    0   16   65   48   59
10   斉藤    1   17   53   58   35
11   加藤    0   12   30   55   53
12   吉田    0   18   56   72   49
13   山田    1    2   85   65   73
14 佐々木    1   19   69   62   66
15   山口    0   20   52   68   45
16   松本    0   11   55   52   63
17   井上    0    8   45   54   56
18   木村    1    7   61   32   88
19     林    1   10   43   57   43
20   清水    1    9   76   88   97
> attach(x)
> 性別
 [1] 1 1 1 1 0 0 0 0 0 1 0 0 1 1 0 0 0 1 1 1
> class(性別)
[1] "numeric"
> class <- factor(性別)
> class(性別)
[1] "numeric"
> 性別 <- factor(性別)
> class(性別)
[1] "factor"
> mean(x)
学生名   性別   順位   国語   数学   英語 
    NA   0.50  10.50  56.25  59.80  60.40 
 警告メッセージ: 
1: mean(<data.frame>) is deprecated.
 Use colMeans() or sapply(*, mean) instead. 
2: In mean.default(X[[1L]], ...) :
   引数は数値でも論理値でもありません。NA 値を返します 
> mean(国語)
[1] 56.25
> mean(x[,4:6])
 国語  数学  英語 
56.25 59.80 60.40 
 警告メッセージ: 
mean(<data.frame>) is deprecated.
 Use colMeans() or sapply(*, mean) instead. 

floor(2.3*100)

so (2012-12-19 (水) 13:50:56)

下記計算結果が導出される理由が分からず困っています。
ご教示いただけると幸いです。

> floor( c( 230, 2.3*100, 2.3*10*10 ) )
[1] 230 229 230
> floor( c( 230000, 2.3*100000, 2.3*100*1000, 2.3*1000*100 ) )
[1] 230000 229999 229999 230000

軸のラベルの向きを斜めに変更する方法

BCD (2012-12-17 (月) 12:44:46)

いつもお世話になっております。
グラフを描く時に、lasオプションを指定すればラベルの向きを変更することができますが、軸に対して斜め(45℃くらい角度をつける)にラベルを配置することはできないのでしょうか?

どうかご教示願います。

staxlab01.png

(多クラスのカテゴリデータでベイジアンネット)dealのautosearchでエラー

yamaguchi (2012-12-16 (日) 17:27:44)

はじめまして。いつもお世話になっています。
初投稿という事で至らない点もあるかと思いますが、よろしくお願いします。

dealパッケージで多クラスのカテゴリデータでベイジアンネットワークを作成しようとしたところ、autosearch関数で下記のエラーが発生しました。

>autosearch(dat.pre_net, dat.dummy, dat.prior)
以下にエラー apply(prior$jointalpha, didx, sum) :                                                     
 dim(X) must have a positive length

どのように解決したらよいか、ご助言をいただけますか。

以下実行手順です。

# サンプルとして成績データ作成
eng <- c("A", "B", "C", "C", "B", "D")
math <- c("B", "B", "C", "A", "B", "D")
sci <- c("A", "A", "B", "D", "B", "B")
res <- c("A", "B", "C", "C", "B", "C")
dat <- data.frame(eng, math, sci, res)
dat
  eng math sci res
1   A    B   A   A 
2   B    B   A   B
3   C    C   B   C 
4   C    A   D   C
5   B    B   B   B
6   D    D   B   C 

青木氏のサイトのqt3, make.dummyを参考にして、カテゴリデータをダミーデータに変換。
上記サンプルではダミーデータに変換しなくてもネットワークを描画できたのですが、自分のデータではjointprior関数を実行する際に "too large than due to tequnical limitations in R"というエラーが起き、オブジェクトのサイズを減らすためにダミーデータに変換することにしました。

dat.int <- dat
dat.int[,1:ncol(dat)] <- lapply(dat, as.integer) # カテゴリを数値に変換
cat <- sapply(dat.int, max) # 各列が持つ水準の個数を保存
cname <- unlist(sapply(dat.int, levels)) # 水準の名称を保存
dat.colnames <- colnames(dat.int)
dat.dummy <- make.dummy(dat.int)
colnames(dat.dummy) <- paste(rep(dat.colnames, cat), cname, sep=".")
           # 列名.水準 というカラム名を設定

library('deal')
dat.pre_net <- network(dat.dummy)
dat.prior <- jointprior(dat.pre_net, 20)
dat.dist <- learn(dat.pre_net, dat.dummy, dat.prior)
# 有効辺ノードの設定
res_list <- grep("res.*", colnames(dat.dummy))
           # このノードからは有効辺を引かない
[1] 12 13 14
# res → res以外への有効辺を引かないようにする
dat.col_num <- 1:length(colnames(dat.dummy))
         # 各ノードを数字のベクトルで取得
dat.no_res <- length(dat.col_num) - length(res_list)
         # resがつかないカラム名の個数を取得

a <- sort(rep(res_list,length(colnames(dat.dummy))))
    # res_listの各要素をカラム数分繰り返したベクトルを
    # ソート(有効辺が出てはいけないノード)
b <- 1:length(colnames(dat.dummy)) # 有効辺が入るノード
dat.ban <- cbind(a, b)
# 12,12などの重複を削除
dat.ban_df <- as.data.frame(dat.ban)
colnames(dat.ban_df) <- c("a", "b")
dat.ban_df[dat.ban_df$a == dat.ban_df$b,] <- c(NA, NA)
dat.ban_df <- na.omit(dat.ban_df)

# 有効辺を引かないノードリストの設定を適用
dat.pre_net <- getnetwork(dat.dist)
banlist(dat.pre_net) <- as.matrix(dat.ban_df)
# モデル探索
dat.posterior <- autosearch(dat.pre_net, dat.dummy,
                            dat.prior)
 以下にエラー apply(prior$jointalpha, didx, sum) :                                                     
 dim(X) must have a positive length

debug関数でautosearch()の処理を追って見ましたが、私には原因がわかりませんでした。
読みづらい点もあるかと思いますが、対処策・アドバイスのほど、よろしくお願いします。

また、実行環境は下記のとおりです。

>sessionInfo()                                                                                        
R version 2.14.1 (2011-12-22)                                                                          
Platform: x86_64-pc-linux-gnu (64-bit)                                                                 
locale:                                                                                                
 [1] LC_CTYPE=ja_JP.UTF-8       LC_NUMERIC=C                                                           
 [3] LC_TIME=ja_JP.UTF-8        LC_COLLATE=ja_JP.UTF-8                                                 
 [5] LC_MONETARY=ja_JP.UTF-8    LC_MESSAGES=ja_JP.UTF-8                                                
 [7] LC_PAPER=C                 LC_NAME=C                                                              
 [9] LC_ADDRESS=C               LC_TELEPHONE=C                                                         
 [11] LC_MEASUREMENT=ja_JP.UTF-8 LC_IDENTIFICATION=C                                                    
attached base packages:                                                                                
 [1] stats     graphics  grDevices utils     datasets  methods   base                                   
other attached packages:                                                                               
 [1] deal_1.2-35

package(e1071)におけるsvmについて

t.h (2012-12-13 (木) 17:28:14)

はじめまして。よろしくお願いします

題名の通りですが、package(e1071)におけるsvmについて質問があります。
私の最終的な目的は「1つのクラスに対して間違いを許さない2クラス判別を行うsvmをつくる」です。
(クラスAとBがあるとして、テストデータを入れたときにAがBと判定されることは許すが、BがAと判定することは許されないSVM、という意味です)

(1)このsvmに訓練データとしてデータを入れるとき、その順番は関係あるのでしょうか?
(e.g)「クラスA,クラスBという順番で入っている行列を訓練データと与えたとき」と、「クラスB,クラスAという順番で入っている行列を訓練データと与えたとき」に構成されるSVMは違うものになるのでしょうか?

(2)このsvmはデフォルトで実行すると、ソフトマージンsvm(最大マージンsvm)のようですが、ハードマージンsvmにするにはどのようにすればよいでしょうか?

(3)目的を達成するために、デフォルトのSVMからどのようにすればよいでしょうか?

(4*)デフォルトのsvmに3クラス以上のデータを訓練させた場合は1vs1式でのクラス分けなのか、1vsOtrhers式でのクラス分けなのかを教えていただきたいです

段々と抽象性の高い質問になっておりますが、((4)だけ別系統の質問です)最終的な目的を達成するにあたり、有用であればどの様な回答でも構いません。

見当違いの質問、読みづらい文章かもしれませんが、ご回答よろしくお願いします。

SEM の summary

しまりす (2012-12-07 (金) 19:14:06)

お世話になります。

SEMパッケージの練習のため、下記URLの「食物とガン」の共分散分析を自分でやってみました。
http://wwwhum.meijo-u.ac.jp/labs/hh002/r/content/sem.html

共分散の結果は出せたのですが、最後のsummary()による適合度の出力ができませんでした。
summary(データ名)と入力しても、何も出力されず、エラーメッセージも出ません。

以下がプログラムと出力結果です。

> # ガンデータのSEM #
> 
> setwd("C:/Users/a2122062/work")
> gandat <- read.table("gan.csv", sep=",", header = TRUE)
> 
> # 相関行列を求める #
> soukan <- cor(gandat)
> soukan
            calory      niku       nyu      sake    daicho  chokucho
calory   1.0000000 0.7731109 0.7155250 0.6341241 0.7394903 0.8100869
niku     0.7731109 1.0000000 0.7227997 0.4515933 0.8526245 0.6718607
nyu      0.7155250 0.7227997 1.0000000 0.2321334 0.5787253 0.4785180
sake     0.6341241 0.4515933 0.2321334 1.0000000 0.5881328 0.6095649
daicho   0.7394903 0.8526245 0.5787253 0.5881328 1.0000000 0.7523879
chokucho 0.8100869 0.6718607 0.4785180 0.6095649 0.7523879 1.0000000
> 
> # 共分散する#
> library(sem)
 要求されたパッケージ MASS をロード中です 
 要求されたパッケージ matrixcalc をロード中です 
 警告メッセージ: 
1:  パッケージ '‘sem’' はバージョン 2.15.2 の R の下で造られました  
2:  パッケージ '‘matrixcalc’' はバージョン 2.15.2 の R の下で造られました  
> modelA.gan <- specify.model()
1: 洋食傾向 -> calory ,b11 ,NA
2: 洋食傾向 -> niku ,b21 ,NA
3: 洋食傾向 -> nyu ,b31 ,NA
4: 洋食傾向 -> sake ,b41 ,NA
5: 消化管ガン -> daicho ,NA ,1
6: 消化管ガン -> chokucho ,b62 ,NA
7: 洋食傾向 -> 消化管ガン,g21,NA
8: calory <-> calory ,e1 ,NA
9: niku <-> niku ,e2 ,NA
10: nyu <-> nyu ,e3 ,NA
11: sake <-> sake ,e4 ,NA
12: daicho <-> daicho ,e5 ,NA
13: chokucho <-> chokucho ,e6 ,NA
14: 洋食傾向 <-> 洋食傾向 , NA , 1
15: 消化管ガン <-> 消化管ガン , delta2,NA
16: 
Read 15 records
 警
'specify.model' is deprecated.
Use 'specifyModel' instead.
See help("Deprecated") and help("sem-deprecated"). 
> 
> semgan <- sem(modelA.gan , soukan, N=47)
> std.coef(semgan, digit=3)
          Std. Estimate                           
1     b11    0.90353083       calory <--- 洋食傾向
2     b21    0.88860193         niku <--- 洋食傾向
3     b31    0.72052466          nyu <--- 洋食傾向
4     b41    0.61993727         sake <--- 洋食傾向
5            0.89590512     daicho <--- 消化管ガン
6     b62    0.83980749   chokucho <--- 消化管ガン
7     g21    0.98052591   消化管ガン <--- 洋食傾向
8      e1    0.18363205         calory <--> calory
9      e2    0.21038661             niku <--> niku
10     e3    0.48084422               nyu <--> nyu
11     e4    0.61567779             sake <--> sake
12     e5    0.19735401         daicho <--> daicho
13     e6    0.29472337     chokucho <--> chokucho
14           1.00000000     洋食傾向 <--> 洋食傾向
15 delta2    0.03856894 消化管ガン <--> 消化管ガン
 警告メッセージ: 
'std.coef' is deprecated.
Use 'stdCoef' instead.
See help("Deprecated") and help("sem-deprecated"). 
> summary(semgan)
> 

なぜなのでしょうか?
もしかしたらとても簡単なミスかもしれないのですが、どうぞよろしくお願いいたします。

lavaanパッケージによる多母集団分析のGFI,AGFIの出力

さとし (2012-12-04 (火) 21:46:32)

お世話になります。Windows7のR 2.15.2を使っております。

lavaanパッケージで多母集団分析を行いました。

lav.model.1 <- lavaan(model.1, data=d1, model.type="sem",
group="group", int.ov.free=TRUE, auto.var=TRUE, auto.fix.first=TRUE,
fixed.x=FALSE, auto.cov.lv.x=FALSE,  auto.cov.y=FALSE)
summary(lav.model.1, fit.measure=TRUE, standardized=TRUE)

このようなオプションでsummaryを出力したところ,色々な適合度が表示されるのですが,GFIやAGFIが出力されません。

そこで質問なのですが,

1.GFIやAGFIの出力を可能にするオプションはあるのでしょうか?
2.GFIやAGFIの算出に必要な共分散行列を出力することはできるのでしょうか?
3.あるいは,出力された他の適合度指標からGFIやAGFIを式展開などで導出することは可能でしょうか?

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

x軸ラベルの回転について

れこ (2012-12-04 (火) 15:45:58)

お世話になります。WindowsのR 2.15.1を使っています。

y <- rnorm(20)
x <- factor(rep(c("AAAAA","BBBBB"),10))
stripchart(y ~ x, vertical=TRUE)

このとき、x軸のラベルを反時計回りに90度回転させたグラフにしたいのですが、どのようにすればよいでしょうか?

宜しくお願いします。

一定時間が過ぎると新規個体群が入ってくる生物の生活史シミュレーションの高速化をしたい

Saito (2012-12-04 (火) 11:34:01)

いつもお世話になっております。
ある生物が、時間が経つに連れて大きくなり、一定時間が過ぎると子供を産んで(子の数は固定)次の新しい年級群が加入し、自然に死んだりたまに捕獲される、という生活史をシミュレーションしたいと思っています。ただ、プログラムに問題があり、大変時間がかかります。下記は私が作った例なのですが、新規に年級群が入るというイベントをif文で書いてしまっています。 実際には、時間も年級群の数ももっとべらぼうに多いですし、パラメタの値を色々変えて検討したいので、素早く計算出来るよう今よりも高速化をしたいと考えております。しかし、私の文法がへたくそなせいで、これ以上は高速化が出来ません。どなたかご教授頂けますと幸いです。

Size <- seq(0.1, 30, by=0.5) #個体のサイズ分布
year <- 40 #時間
class <- 3 #生物の年級群の数

M <- catch <- matrix(0, ncol=year*class, nrow=length(Size))

###Catchは既知のデータ###
catch[, 5] <- rpois(length(Size), lambda=2)
catch[, 10] <- rpois(length(Size), lambda=3)
catch[, 56] <- rpois(length(Size), lambda=4)
catch[, 77] <- rpois(length(Size), lambda=5)
catch[, 105] <- rpois(length(Size), lambda=6)

SS1 <- matrix(NA, nrow=length(Size), year*class)
SS2 <- matrix(NA, nrow=length(Size), year*class)
SS3 <- matrix(NA, nrow=length(Size), year*class)

m1 <- runif(year*class, 0, 1) #平均成長の値(時間変化する)
m2 <- runif(year*class, 5, 10)
m3 <- runif(year*class, 15, 20)

sd1 <- runif(year*class, 1, 5) #平均成長の分散(時間変化する)
sd2 <- runif(year*class, 1, 3)
sd3 <- runif(year*class, 1, 2)

surv1 <- runif(year*class, 0.9, 0.99) #生存率(時間変化する)
surv2 <- runif(year*class, 0.9, 0.99)
surv3 <- runif(year*class, 0.9, 0.99)

###シミューション開始###
M[1, 1] <- S2 <- S3 <- 40000
for(i in 2 : (year*class)){
if (i <= year) {
A1 <- dnorm(Size, mean=m1[i], sd=sd1[i])
G1 <- sample(Size, sum(M[, i-1]), prob=A1/sum(A1), rep=T)
M[, i] <- 
SS1[, i] <- rbinom(rep(1, length(Size)),
                   table(cut(G1, breaks=c(Size, max(Size) + 0.5))), 
                   prob=surv1[i])

M[, i] <- M[, i] - catch[, i]
M[which(M[, i] < 0), i] <- 0
S1 <- M[, i]

} else { #次の年級群が加わってくる
if (i > year & i <=(year*2+13)) {
A1 <- dnorm(Size, mean=m1[i], sd=sd1[i])
A2 <- dnorm(Size, mean=m2[i], sd=sd2[i])

G1 <- sample(Size, sum(S1), prob=A1/sum(A1), rep=T)
G2 <- sample(Size, sum(S2), prob=A2/sum(A2), rep=T)

S1 <- SS1[, i] <- rbinom(rep(1, length(Size)),
                         table(cut(G1, breaks=c(Size, max(Size) + 0.5))), 
                         prob=surv1[i])
S2 <- SS2[, i] <- rbinom(rep(1, length(Size)),
                         table(cut(G2, breaks=c(Size, max(Size) + 0.5))), 
                         prob=surv2[i])

S1 <- S1 - round(catch[, i]*(SS1[, i]/(SS1[, i]+SS2[, i])))
       #それぞれの年級群の数に比例した値で捕獲
S2 <- S2 - round(catch[, i]*(SS2[, i]/(SS1[, i]+SS2[, i])))
S1[which(is.na(S1))] <- 0
S2[which(is.na(S2))] <- 0

M[, i] <- S1 + S2
M[which(M[, i] < 0), i] <- 0
} else { #また次の年級群が加わってくるが、加わるタイミングが違う
if (i > (year*2+13)) {
A1 <- dnorm(Size, mean=m1[i], sd=sd1[i])
A2 <- dnorm(Size, mean=m2[i], sd=sd2[i])
A3 <- dnorm(Size, mean=m3[i], sd=sd3[i])

G1 <- sample(Size, sum(S1), prob=A1/sum(A1), rep=T)
G2 <- sample(Size, sum(S2), prob=A2/sum(A2), rep=T)
G3 <- sample(Size, sum(S3), prob=A3/sum(A3), rep=T)

S1 <- SS1[, i] <- rbinom(rep(1, length(Size)),
                         table(cut(G1, breaks=c(Size, max(Size) + 0.5))), 
                         prob=surv1[i])
S2 <- SS2[, i] <- rbinom(rep(1, length(Size)),
                         table(cut(G2, breaks=c(Size, max(Size) + 0.5))), 
                         prob=surv2[i])
S3 <- SS3[, i] <- rbinom(rep(1, length(Size)),
                         table(cut(G3, breaks=c(Size, max(Size) + 0.5))), 
                         prob=surv3[i])
S1 <- S1 - round(catch[, i]*(SS1[, i]/(SS1[, i]+SS2[, i]+SS3[, i])))
S2 <- S2 - round(catch[, i]*(SS2[, i]/(SS1[, i]+SS2[, i]+SS3[, i])))
S3 <- S3 - round(catch[, i]*(SS3[, i]/(SS1[, i]+SS2[, i]+SS3[, i])))

S1[which(is.na(S1))] <- 0
S2[which(is.na(S2))] <- 0
S3[which(is.na(S3))] <- 0

M[, i] <- S1 + S2 + S3
M[which(M[, i] < 0), i] <- 0
}}}}
###シミューション終わり###
###得たいもの1###
SS1[is.na(SS1)] <- 0
SS2[is.na(SS2)] <- 0
SS3[is.na(SS3)] <- 0
SS <- SS1 + SS2 + SS3 
SS[, c(5, 10, 56, 77, 105)]
   #(既知データである)捕獲があった時点での
   #捕獲前の推定個体数の総和

###得たいもの2###
M #捕獲を考慮して計算を続けたときの推定個体数 

対応した2つのデータフレームに一括でペアT検定とWilcoxon符号順位検定を行い、各列の因子とp値一覧を出力したい。

のぶ (2012-12-03 (月) 23:04:10)

お世話になります。

Rは始めたばかりの初心者です。
ご教示頂けると幸いです。

下記のようなデータフレームX、Yの、各A-A',B-B',C-C',D-D'間の、ペアT検定とWilcoxon符号順位検定を行い、各列の因子とp値一覧を出力したいと考えております。

X <- data.frame(A = c(1,5,6,4,2), B = c(2,1,7,4,3) ,
     C = c(1,NA,NA,NA,NA) , D = c(3,8,NA,4,1))
Y <- data.frame(A' = c(NA,3,2,8,1), B' = c(3,2,5,4,9) ,
     C' = c(NA,2,8,3,10) , D' = c(1,3,NA,2,NA))

検討したい因子がかなり多く、p値の一覧として出力した上で、BH法にてFDRを計算する予定です。
なので、

	検討できた数	p-value(paired-T.test)	p-value(wilcoxon)
A-A'	4	0.6892	 0.7127
B-B'	5	0.4144	0.5807
C-C'	0	NA	NA
D-D'	3	0.09547	 0.1736

のように、要約された結果一覧を得たいです。

for (i in 1:4) print(t.test(X[,i],Y[,i],paired=TRUE))
for (i in 1:4) print(wilcox.test(X[,i],Y[,i],paired=TRUE))

とするところまではたどり着きましたが、各検定結果において、何と何を比較したか分からない上にまとまっておらず、さらに全てNAの列があるとそこでスクリプトがとまってしまいます。
何か良い方法はありますでしょうか。
お手数をおかけいたしますが、よろしくお願いいたします。

mimRが手ごわい?!

悩めるポー (2012-12-01 (土) 18:40:47)

R version 2.13.0(2011-04-13)
mimR_2.6.2
gRbase_1.4.4の環境下で、mimRとMIM3.2.0.7のコネクションがうまくいかないのは何故でしょうか?
The mimR Package
January 31, 2006をなぞっても、、、。

webページの更新日時取得

asap (2012-11-30 (金) 21:30:28)

webページの内容は

 src <- readLines("http://www.jma.go.jp/jp/")

のようにして取得できますが、そのページの更新日時を取得できる方法はありますか?
例えばfirefox(Win版)であれば、右クリック>ページの情報を表示 を選択すると
そのページの更新日時が表示されますが、それをRで取得したいと思っています
よろしくお願いします

stripchartで対応のある点を線でつなぐ

COM (2012-11-26 (月) 10:02:04)

お世話になります。
stripchart() で因子別に一次元散布図を描き、対応するデータを線でつなぎたいと思っています。

y <- rnorm(20)
x <- factor(c(rep(1, 10), rep(2, 10)))
stripchart(y~x, vertical=T)

このとき、因子1の1番目と因子2の1番目、因子1の2番目と因子2の2番目・・・
というように点を線でつなぐにはどのようにすれ- d ばよいでしょうか?
宜しくお願いします。

factor データを strptime で"POSIXt" に変換したときにデータ数が変わる??

やまだ (2012-11-24 (土) 06:24:36)

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

x
[1] 2010/02/02 00:15:00 2010/02/02 01:15:00 2010/02/02 02:15:00 2010/02/02 03:15:00
[5] 2010/02/02 04:15:00 2010/02/02 05:15:00 2010/02/02 06:15:00 2010/02/02 07:15:00
[9] 2010/02/02 08:15:00 2010/02/02 09:15:00 2010/02/02 10:15:00 2010/02/02 11:15:00
[13] 2010/02/02 12:15:00 2010/02/02 13:15:00 2010/02/02 14:15:00 2010/02/02 15:15:00
[17] 2010/02/02 16:15:00 2010/02/02 17:15:00 2010/02/02 18:15:00 2010/02/02 19:15:00
20 Levels: 2010/02/02 00:15:00 2010/02/02 01:15:00 ... 2010/02/02 19:15:00

(y <- strptime(x, "%Y/%m/%d %H:%M:%S"))
[1] "2010-02-02 00:15:00" "2010-02-02 01:15:00" "2010-02-02 02:15:00" "2010-02-02 03:15:00"
[5] "2010-02-02 04:15:00" "2010-02-02 05:15:00" "2010-02-02 06:15:00" "2010-02-02 07:15:00"
[9] "2010-02-02 08:15:00" "2010-02-02 09:15:00" "2010-02-02 10:15:00" "2010-02-02 11:15:00"
[13] "2010-02-02 12:15:00" "2010-02-02 13:15:00" "2010-02-02 14:15:00" "2010-02-02 15:15:00"
[17] "2010-02-02 16:15:00" "2010-02-02 17:15:00" "2010-02-02 18:15:00" "2010-02-02 19:15:00"

length(x); length(y)
[1] 20
[1] 9

上記のようなコマンドを実行し、length(y)が9となるのはなぜでしょうか。どうぞご教授ください。

barplotやbwplotにおける因子(項目)の並べ替え

atr (2012-11-20 (火) 19:25:14)

お世話になります。

lab <- sample(rep(c("亜","い","u","エ","O"), length.out=10), 10) #サンプル
sp <- runif(10, min=0, max=7)
DF1 <- data.frame(lab, sp)
MEAN1 <- tapply(DF1$sp, DF1$lab, mean)
barplot(MEAN1)

にて描画される棒グラフのバーの並び順を例えば左から「"亜","い","u","エ","O"」とする方法はないのでしょうか。
低い検索スキルのせいかRjpwiki内では見つかりませんでした。ご教授頂けると幸いです。

SVMの結果がsummaryとtableで違う

(2012-11-13 (火) 16:00:14)

はじめまして。よろしくお願いします。
今、問題になっているコードは以下の部分です。

library(e1071)
(eeg.svm <- svm(NAME~.,data=zz,cross=7))#svm
summary(eeg.svm)
pred2 <- predict(eeg.svm, zz)
table(pred2,zz[,1])

zzは、1列目NAMEに人物名(イニシャル)が入っており、2〜1201列目に各個人の波形が入っている、50×1201のデータフレームです。

このデータをsvmにかけて、個人の判別を行おうとしています。
そこで、summary(eeg.svm)を実行し結果を見ると

Call:
svm(formula = NAME ~ ., data = zz, cross = 7)

Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.0008333333  

Number of Support Vectors:  50

 ( 10 10 10 10 10 )
 
Number of Classes:  5 

Levels: 
 m N o u y

7-fold cross-validation on training data:

Total Accuracy: 62 
Single Accuracies:
 42.85714 85.71429 42.85714 57.14286 85.71429 42.85714 75 

と出てきて、Accuracyが62%だとわかります。
では、その内訳がどうなっているのか、と予想からクロス表を作るため

pred2 <- predict(eeg.svm, zz)
table(pred2,zz[,1])

を実行してみると、

pred2  m  N  o  u  y
    m  8  0  0  0  0
    N  2 10  1  0  0
    o  0  0  9  0  0
    u  0  0  0 10  0
    y  0  0  0  0 10

と出てきます。このクロス表によればAccuracyは94%です。
なぜ、このように違った結果が出てしまうのでしょうか?

見辛い文章になっていたら申し訳ありません。(←編集し直しましたよ)
ご回答よろしくお願いします。

「NULL に対して属性を設定しようとしました」のエラーメッセージの原因を教えてください

佐藤 (2012-11-07 (水) 22:38:03)

初めまして。今回、初めて投稿させて頂きます。
現在とあるシステムのエラー原因が分からず困っております。
ご指南を宜しくお願いします。

<状況>
4年前からあるシステムを運用していて、Rのソースを一切変えず、これまで運用をしてきたのですが、先日下記エラーが発生しました。想定しないデータが連携された為、以下のエラーが発生するようになったと思われるのですが、エラー内容をみたところ、何が原因か分かりません。

下記のエラー内容を見て、原因が分かれば教えて頂けないでしょうか?

<ログファイルからのエラー直後の箇所を抜粋>

if (nrow(tDatStdModelPros) > 0) {

    res <- by(tDatStdModelPros, factorList,
              function(x) analyticalPerspectiveByStdModelFunc( x,
              stdPredictPros[(stdPredictPros[,"COLUMN_3"] == x[1,"COLUMN_3"]),],
              indexInfo[x[1,"COLUMN_3"],"COLUMN_11"],
              stdModelProsList[x[1,"COLUMN_2"],], startYM, 0))
    # 結果が List なので matrix に変換
    colName <- c("COLUMN_1",
                 "COLUMN_2",
                 "COLUMN_3",
                 "COLUMN_4",
                 "COLUMN_5",
                 "COLUMN_6",
                 "COLUMN_7",
                 "COLUMN_8",
                 "COLUMN_9",
                 "COLUMN_10") 
    predictValProsStd <- matrix(unlist(res), ncol=10,
            byrow = TRUE, dimnames = list(NULL, colName))

    predictValProsStd[,"COLUMN_9"] <- 3
    predictValProsStd[,"COLUMN_6"] <- 2

    rm( res)
}
以下にエラーmatrix(unlist(res), ncol = 10,
  byrow = TRUE, dimnames = list(NULL,  : 
        NULL に対して属性を設定しようとしました
実行が停止されました~

Jaccard係数(類似度)計測後のデンドログラム描画方法を教えてください

mura (2012-11-07 (水) 14:18:02)

はじめまして。R初心者の大学院生で、修士論文のデータ整理で困っております。ご指南を宜しくお願いします。

実は、Jaccard(ジャッカード)係数で約10個のデータの相互の結びつきの強さを計算し、横書きのデンドログラムで描画したいと考えています。
スケールは、ジャッカード係数なのでデータは最大値は1.0で最小値は0となり0.XXXで表現され、数値が大きい(1に近い)程結びつきが強いということで樹形図は結節されていきます。

今、10個のデータ間のジャッカード係数を計算してある縦横の表(対照行列)がJaccard.CSVという名称のCSVファイルがあるとき、Rでどのようにすれば横書きの樹形図を書くことができるでしょうか?
Rについては超初心者で、本を見ながら練習をはじめたばかりです。
大変お手数ですが関数とその他引数も含めてどのようにすれば良いかご指南願えれば幸いです。
どうぞ宜しくお願い致します。

日頃、デンドログラムを描くことがないので、今までplot.hclust()ってhorizオプションがないことを知らなかった。何か理由があるのだろうか?

インストール済みのパッケージのアーカイブ化(zip形式)

わnがねぇー (2012-10-27 (土) 11:23:59)

 最近、CRAN より、インストールしたパッケージが、CRANサイトより消えてしまいました。

 新しいPCに、Rをインストールしているのですが、このパッケージも利用する必要があります。

 ンストール済みのパッケージをアーカイブ化(zip形式)して、新しいRの環境で利用する方法はあるのでしょうか?

長さの異なるベクトルを渡した時の"関数の引数"と"配列の要素番号"で展開される順序

satto (2012-10-26 (金) 02:26:34)

初めて投稿させて頂きます
以下のコードを実行した時に表示される要素番号とその要素が指す値が実際のものと一致しません。

xn <- 2
yn <- 3
zn <- 4
a <- c(1:(xn*yn*zn))
dim(a) <- c(xn, yn, zn)
f1 <- function(x, y, z) {
	cat("a[x,y,z] : ")
	cat(
		sprintf("[%d,%d,%d]=%2d,",
			x, y, z, a[x,y,z])
	)
	cat("\n")
}
f1(1:xn, 1:yn, 1:zn)

おそらく、sprintf()の書式指定子に渡される順序と、配列の[]内でベクトルが展開される順序が一致していないのだと思いますが、このような表示を行うために何か良い方法はないでしょうか。

回帰分析の繰り返し処理のやり方について

悩むぽー (2012-10-24 (水) 20:15:37)

お世話になります。
処理を繰り返すには、for ループを使えばいいんだろうな?と思うのです。
例えば、回帰分析の結果を変数の組み合わせで、数十回繰り返す場合、どんな処理をするのが、より効率的かと悩んでいます。文法をしっかりと理解しないまま使うことに課題はあるのですが、どなたかご教示ください。よろしくお願いします。

変数ABとvariable1の回帰分析をするとして、変数ABとn個の異なる変数の回帰分析を繰り返す処理です。下記のように、commandを繰り返し入力すれば出来ることですが、、、。

RegModel.1 <- lm(AB~variable1, data=dataset)
summary(RegModel.1)
RegModel.2 <- lm(AB~variable2, data=dataset)
summary(RegModel.2)
RegModel.3 <- lm(AB~variable3, data=dataset)
summary(RegModel.3)

年次別データの欠損年の自動補間

atr (2012-10-23 (火) 19:20:06)

お世話になります。
下記のように、'00〜'02年、'04年は測定実施、'03年測定せず、というデータに対し、年毎の平均値を棒グラフ化(または折れ線グラフ化)したい場合に

yr <- sample(rep(c(2000:2002, 2004), length.out=10), 10) #サンプル
sp <- runif(10, min=0, max=7)
DF1 <- data.frame(yr, sp)
MEAN1 <- tapply(DF1$sp, DF1$yr, mean)
barplot(MEAN1)

上記barplot(MEAN1)のように'02年の次が'04年とならず、下記barplot(MEAN2)のように、測定値のない'03年がゼロ等となるようにデータ自動補間またはプロット時に'03年が自動補間する方法はないでしょうか?

DF2 <- DF1
DF2[(nrow(DF2)+1), ] <- c(2003, 0) # '03年平均ゼロ化(便宜上)
MEAN2 <- tapply(DF2$sp, DF2$yr, mean)
barplot(MEAN2) #狙いの出力

測定ナシの年のデータを手動追加で作成、はデータが膨大であり出来れば避けたいです。以上宜しくお願いします。

パッケージvcdのモザイクプロットに関するいくつかの疑問

vcd初心者 (2012-10-23 (火) 17:43:10)

お世話になります。
Windows 7上でR i386 2.15.1を使って、vcd_1.2-13を試しています。
パッケージvcdのモザイクプロットに関して、いくつか疑問があります。

library(vcd)
mosaic(Titanic)

とすると、描画領域全体に広がるプロットが作成されますが、次のように仮想データを入力して描画すると、

tab <- matrix(c(10, 30, 40, 60, 70, 30), nrow=3, byrow=T)
dimnames(tab) <- list(Size=c("S", "M", "L"), Gender=c("Male", "Female"))
mosaic(tab, split_vertical=T)

正方形の領域にしかプロットされません。
これを描画領域全体に広げるにはどうしたらよいのでしょうか。

また、度数を各セルに表示しようとして、次のように入力すると、

mosaic(tab, split_vertical=T, labeling=labeling_cells(text=tab))

度数は表示されるのですが、変数名と水準名のラベルが消えてしまいます。
最低、水準名は表示したいのですが、どうしたらよいのでしょうか。

さらに、ラベルを日本語で表示しようとして、次のように入力すると、

> dimnames(tab) <- list(サイズ=c("S", "M", "L"), 性別=c("男", "女"))
> mosaic(tab, split_vertical=T)
 以下にエラー grid.Call.graphics(L_downviewport, name$name, strict) : 
  Viewport 'cell:サイズ=S' was not found

とエラーになってしまいます。日本語が通らないようなのですが、何か対処法はありますでしょうか。

以上、長々と申し上げましたが、どなたか解決策を御教示いただければ幸いです。

関数の引数の値の有効性の確認方法の教授願い

ひろ (2012-10-23 (火) 12:02:25)

関数の引数(パラメータ)aの値がvalidか否かで処理をif分岐することを考えています.しかし,aの値の存在,validが正しく判断できずに困っています.

というのもexists("a")とするとTRUEが返りますが,他is.na,is.null,in.nan等で値の有無を調べると「wrapup 中にエラーが起こりました 'a'が見つかりません」とエラーが返ってきます.aは存在するのに,アクセス出来ない原因がわかりません.お知恵をお貸し頂けると助かります.

■ 追加ソースコード

test<-function(a,b){
	browser()
	missing(a)		# TRUE
	missing(b)		# FALSE
	print(exists("a"))	# TRUE
 	print(exists("b"))	# TRUE
 	print(ls(a))		# [1] "a" "b"
	print(ls(b))		# 正常にオブジェクト一覧表示
	print(search(a))	# wrapup 中にエラーが起こりました
                               #  使われていない引数 (a)
	print(search(b))	# wrapup 中にエラーが起こりました
                               #  使われていない引数 (b)
	is.na(a)		# wrapup 中にエラーが起こりました
                               #  オブジェクト 'a' がありません
	is.na(b)		# FALSE
}
test(b=1)

添え字が許される範囲か否かを確かめる方法

ひろ (2012-10-23 (火) 10:46:51)

添え字が許される範囲か否かを確かめる方法はありますでしょうか?

■ 試行内容:
SVMといった機械学習を用い,弁別予測・正解結果より精度を算出しようとしています.

ところが,結果の行列の構成が結果で異なります. 下記の様に,"正常","異常"の2つの行列が通常現れる

pred   異常 正常
  異常 100     0
  正常    0  100

ところが予測結果が悪いと次のように"異常"のみ現れ,"正常"行がないことがあります.結果res["正常",]とすると"wrapup 中にエラーが起こりました 添え字が許される範囲外です"とエラー.

pred   異常 正常
  異常 100   100

hclustの各labelsの枝の長さ

kin (2012-10-23 (火) 07:59:59)

お世話になります。

hclustで樹形図を作成しています。

今回の質問は、各ラベルから延びている枝の高さを得たいのです。
→ 各ラベルに直接接続している枝の高さです。

dat.hc <- hclust(dat, method="ward")

とした場合、dat.hc$heightで各クラスターごとの枝の高さが得られ、dat.hc$merge でクラスター形成履歴?が判るので、そのラベル(マイナス値)と対応させるのかなと考えていますが、いまいちスマートなスクリプトが思い浮かびません。

アドバイスなど頂ければ幸いです。
宜しくお願いします。

K-meansの結果出力

Micheal (2012-10-23 (火) 05:09:56)

dataset <- 1:100
x <- matrix(dataset, 50, 2)
w <- km <- kmeans(x, 2)

の結果得られるClustering vectorを出力したくて

write(w, file="result.txt", ncolumns=1)

としたのですが、うまく出力する事が出来ません。ご示唆を頂ければと思っています。よろしくお願い致します。

mvpartによる決定木描画時の凡例について

SKYLINE (2012-10-20 (土) 21:54:00)

http://mjin.doshisha.ac.jp/R/19.pdfを参考に mvpart 関数を用いて決定木を描画しました。

library(mvpart)
iris.rp <- mvpart(Species~., data=iris)
plot(iris.rp, uniform=TRUE, branch=0.6, margin=0.05)
text(iris.rp, all=T, use.n=TRUE)

ところが、上記のサイトのような凡例 (色と品種の対応) が表示されません。
どうすれば凡例が表示されるのでしょうか?
ご教授お願い致します。

環境
MacOSX 10.8.2
R 2.15.1
RStudio 0.96.331
mvpart 1.6

数字変換

アトム (2012-10-19 (金) 07:41:20)

100x250のマトリックス(dim(a)=100, 250)の中に1から9までの数字が入っています。この中から3だけを選んで10に変換したいと思い、

b <- gsub(3, 10, a)

を実行したのですが上手くいきません。ご示唆を頂ければと思っています。

0から始まるデータの読み込み

ひろ (2012-10-18 (木) 13:01:50)

何時も参考にさせて頂いています。R初心者です。x<-c(034, 0223, 098, 076, 054)というデータを読み込もうと思ったら0が表示されずに、>x 34, 223, 98, 76と読み込まれてしまいます。数字は数値ではありません。0も一緒に読み込む事は出来るのでしょうか?

半角と全角混在文字列のバイト数

kin (2012-10-18 (木) 10:44:24)

お世話になります。
日本語のフレーズを指定文字数で改行したいと考えています。
しかし、ncharでは半角全角も1文字と判断するので半角と全角混在文字列では、指定文字数で改行を入れても1行の見た目の長さが異なってしまいます。
バイト単位で文字数を得るような関数はありませんか?
宜しくお願いします。

出現回数の計算法

Micheal (2012-10-18 (木) 09:04:32)

R初心者です。ランダムな組み合わせの数字が入ったデータフレームを作成しました。

y <- c(897, 34565, 32, 1,345, 324, 98765, 808, 897, 808, 808, 456, 3, 454, 3, …)
length(y)=30000

この中から同じ組み合わせの数字が幾つ出現しているかを調べ表示したいのですが、解決方法はありますでしょうか?得たい出力は、

(組み合わせの数字,出現回数)

808, 2
1345, 35
32, 3
…

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

異なる列数のデータフレームの連結

ひろ (2012-10-15 (月) 15:12:54)

いくつかの列は同じ列名で,幾らか異なる列名のdata.frameを2つ行連結したいのですが,列数が一致しないとエラーが返ってきて実現できません.解決方法がないでしょうか?

例えば,data.frame1とdata.frame2があるとします.2つの行連結すると,共通の列名が有るところは通常通り連結して,無い箇所は"NA"といった値で補完され,結果data.frame.Answerに成ることが希望です.

data.frame1:

    A    C    D
1   1    2    3
2   4    5    6

data.frame2:

    A    B    C    D
1   7    8    9   10
2  11   12   13   14

data.frame.Answer(<- data.frame1+data.frame2)

    A    B    C    D
1   1    NA   2    3
2   4    NA   5    6
3   7    8    9   10
4  11   12   13   14

コンソール画面の更新

kin (2012-10-14 (日) 07:25:10)

時間がかかる処理の進行状況がわかるようにcatにてコンソールに出力しています。
しかし、コンソール画面の行数が26行とすると、そこまでの出力が溜まらないと画面に表示されないようです。
RGUI設定のbuffer console by deault?のチェックは外しています。

現在は、事前に改行(\n)を26個くらい出力することで対応していますが、他に良い方法はありますでしょうか?

geoRでのバリオグラム作成について

cet (2012-10-14 (日) 01:46:07)

はじめまして。いつもこちらを参考にさせていただいています。
R初心者なのですが、どうしても研究でKrigingによる空間補間が必要になり、現在、R version 2.15.1(Win7/64bit)でgeoRのパッケージを利用しています。
緯度、経度、値の情報が入ったcsvを読み込みchl200708というデータフレームを作成しました。(緯度範囲:30-55N,経度範囲:140-180E データ行数:57690)

このような形式のデータを使用しております。

> head(ttt200708)
      lon    lat     ttt
1 140.042 54.958 0.73995
2 140.083 54.958 0.74817
3 140.125 54.958 0.74817
4 140.167 54.958 0.77717
5 140.208 54.958 0.77717
6 140.250 54.958 0.78524

(7~57690行目は省略いたします) このデータをパッケージgeoRコマンドを用いて

testgeo <- as.geodata(ttt200708)

でgeodataの形式にすることはできました。
しかし、距離行列を計算しようとすると”負の長さのベクトルは許されません”というエラーが出てしまいます。
同様に、バリオグラムを作成しようとすると

> vario1 <- variog(testgeo,max.dist=0.1)
variog: computing omnidirectional variogram
 以下にエラー dist(as.matrix(coords)) :  負の長さのベクトルは許されません 

となってしまいます。緯度・経度の並び順が悪いのかと思い、昇順・降順を入れ替えてみたりしましたが、それでも同様のエラーになってしまいます。
オンラインでもこの間違いの解決方法を探したのですが、みつかりませんでした。

勉強不足なのはわかっているのですが、どこかにこの解決方法が載っていますでしょうか?
もし、解決方法をご存知でしたら、教えてください。
どうぞよろしくお願いいたします。

tickの長さ

MY (2012-10-12 (金) 02:50:04)

いつもお世話になっております。
plot関数などの作図関数を使ってグラフのx軸、y軸の目盛りの線の長さを調節したいのですがどのようにすれば良いでしょうか?目盛りの太さを調節するにはlwd.ticksオプションをいじればよいことは分かったのですが...
どうぞよろしくお願い致します。

合計値

Micheal (2012-10-11 (木) 07:40:43)

R初心者の者です。

dataset=c(1,9,2,7,72,3,7,99,21,20,11,12,122,45,198,334,43,39,34,45…)
(要素数は100個)というデータセットにおいて、前から順に5個づつ合計していきたいのですが、計算方法が分かりません。イメージは(1+9+2+7+72, 3+7+99+21+20, 11+12+122+45+198, ….)です。どうぞ、宜しくお願い致します。

各行のカウントの仕方について

sheet (2012-10-09 (火) 23:11:19)

24行x500列の時系列データ(=data)で、1〜24行の各行が時間、1〜500列が人物の位置データになっています。1=家、2=オフィス、3=買い物…10=それ以外、で表しています。 (data[,1]=(10,10,10,10,10,10,10,10,4,4,2,2,2,2,2,2,2,2,2,2,3,4,4,10)、の様な)

これらのデータから各要素(1から10)の各時間毎の総計(rowSums(data==1), rowSums(data==2)….)を取り出し、更にその各時間毎(各行毎)に「新たにその行動を起こした人の数だけを数えたい」と思っています。例えば

data[,1]=(1,10,10,10,10,10,10,4,4,4,2,2,2,2,2,2,2,2,2,2,3,4,4,10)
data[,2]=(10,10,10,10,10,10,10,10,4,4,4,4,2,2,2,2,2,2,2,2,3,4,4,10)
data[,3]=(10,10,10,10,3,3,4,4,4,4,2,2,2,2,2,2,2,2,2,2,3,4,4,10)

というデータセットにおいては、第一行のrowSums(data==10)=2, rowSums(data==1)=1になると思います。

ここでdata[,1]では第一行に1が初めて現れているので、3つのデータセット全てにおける第一行の1という要素の合計数は1(rowSums(data==1)=1)。反対に、data[,2]とdata[,3]においては10という要素が現れているので、第一行の10という要素の合計は2となります(rowSums(data==10)=1)。

しかし第二行を見てみると、data[,1]では10という要素が「新たに」現れているので、第二行の10という要素に関してはカウントするのですが、data[,2]とdata[,3]の第二行は10となっていて、第一行から連続しているので「新しく行動を起こした人」とは考えず、カウントしない事とします。この場合、第二行のrowSums(data==1)は1、第二行のrowSums(data==10)は2となる様に計算したいのですが、どの様にすれば良いのか検討が付きません。

何かご示唆を頂ければと考えております。よろしくおねがいします。

時系列データ変換

Pen (2012-10-09 (火) 08:18:53)

時系列データを作っています。各列に個人のIDを振り、各行には毎分毎に、その個人から得られたデータを保存してあります。データは100人から得られたので、列の総計は100。行は1日を分毎に保存してあるので、24(h) x 60 (min) = 1440行となっていて、それを行列(=a)としてRに読み込みました。ここから時系列データにする為に、ts関数を使って

a <- ts(a, start=0, frequency=1440)

としたら、

第一行:0,000000000
第二行:0,0006944444
第三行:0,0013888889
第四行:0,002083333
....

となってしまいました。
この数字は一体何なのでしょうか?
普通の時刻(00:00:00, 00:01:00, 00:02:00)で表示する事は可能なのでしょうか?
何かアドバイスを頂ければ幸いです。

計算結果の出力について

SP (2012-10-09 (火) 07:56:54)

PP.test の繰り返し計算で得られた結果を出力したいのですが、どうにも上手くいきません。

for (i in 1:3) {
  b <- PP.test(a[, i])
  print(b)
}

という計算式を作り、得られた結果は

Phillips-Perron Unit Root Test
data:  a[, i] 
Dickey-Fuller = -3.5585, Truncation lag parameter = 7, p-value = 0.03666

	Phillips-Perron Unit Root Test
data:  a[, i]
Dickey-Fuller = -3.8536, Truncation lag parameter = 7, p-value = 0.01631

	Phillips-Perron Unit Root Test
data:  a[, i] 
Dickey-Fuller = -7.0615, Truncation lag parameter = 7, p-value = 0.01

となりました。この P-value だけをまとめて出力したいと思っています。ファイル形式は何でも構いません。宜しくお願いします。

散布図のプロットを分散化する方法

yy (2012-10-04 (木) 13:17:40)

対応分析や多次元尺度構成法などで散布図を作成しています。

データや解析手法により座標面の一部分に集中してプロットされることがしばしば発生しています。

この際に、作成されたプロットデータをスプリングモデルのような感じでノードの近いものは斥力が発生し、遠いものは引力が発生するみたいな感じで、ノードを分散させてみたいと考えていますが具体的な関数や手法の見当が全くついていません。

なにか、アドバイスをいただければ幸いです。

引数名の表示

taked (2012-10-03 (水) 18:58:04)

関数の引数名を表示させたいのですが、可能でしょうか?
下はうまくいっていない例ですが、このような関数を実行したときに引数名(オブジェクトの変数名)である"y"をcat変数などで表示させたいです
yはベクトルです。よろしくお願いします。

 y <- c(1.9, 0.8, 1.1, 0.1,-0.1,4.4,5.5,1.6,4.6,3.4)
 Data <- function(x){
 #  cat(names(x))
 }
 Data(y)

線形判別分析をしたいのですが(エラー:variables are collinear)

(2012-10-03 (水) 15:42:52)

はじめまして。
線形判別分析を用いて、個人認証システムを実装しようと考えています。
元になる波形はcsv形式で保存されており、
30000行×10列で、1〜5列目が人物yの波形データ、残りが人物mの波形データです。

eegMY <- matrix() #空の行列を用意
eegMY <- matrix(scan("csvファイルのパス",sep=',',n=30000),nrow=30000,ncol=10,byrow=T) #行列にデータを読み込む
eegMYt <- t(eegMY) #転置する
name <- matrix(, nrow=10, ncol=1) #個人を判別するための情報をy,mとして、行列を用意する
name[1:5,1] <- "y"
name[6:10,1] <- "m"
x <- cbind(name,eegMYt) #データと個人識別行列を結合する
z <- x[1:10,1:101] #※1
y <- as.data.frame(z) #ldaは行列形式は対応していないのでデータフレームに型変換する
library(MASS)
ldadata <- lda(V1~.,data = y) #実行できない(エラー:variables are collinear)

※1 データが大きすぎるため実行できないようなので、とりあえず範囲を区切っています。
現在このようなコードになっています。
V1に用意した個人認証用の"y"か"m"が入っているので、それをもとにldaを実行したいのですが、

In lda.default(x, grouping, ...) : variables are collinear

という警告メッセージが表示されます。
どのようにすれば、このプログラムが完成に向かうでしょうか?

この様な、本格的なコードを書くのは初めてなので、見当違いの質問かもしれませんが、ご教授いただけたら幸いです。よろしくお願いします。

1つのグラフを複数ページPDFにしたいのですが

yds (2012-09-27 (木) 13:50:39)

hclustでデンドログラムを作成しpdfに出力しています。
pdf(family="Japan1GothicBBB", file="dendorogramTI.pdf", width=11.69*3, height=8.27)
td.hc <- hclust(td.d, method="ward")
plot(td.hc, cex=0.7, hang=-1, xaxs="i", yaxs="i")
dev.off()

良し悪しは別として、デンドログラムの項目(要素)が多くなり、例えばA4を横3枚に並べたようなサイズになってしまいます。
これをA4で3ページの1つのPDFとして出力することはできないでしょうか?

setwd() 後の、Rコードのソースを読み込み… での問題

yy (2012-09-26 (水) 08:57:07)

windows7 32bit R i386 2.15.0 を使用しています。

Rのショートカットの作業フォルダは以下のようにしてデスクトップにRフォルダを作成しています。

"C:/Users/hogehoge/Desktop/R"

その配下にcsvフォルダと、scriptフォルダを作成して、各々csvファイルや拡張子rのスクリプトファイルを格納しています。
ユーザーは、ファイル → Rコードのソースを読み込み… でscriptフォルダに移動後に任意のrスクリプトファイルを読み込みます。
rのスクリプトファイルの最初には、下記のように記載されていてcsvファイル選択後、解析処理をさせています。

library(svDialogs)
dirBase <- getwd()
setwd(paste(dirBase, "/csv" , sep=""))
csvfile <- dlgOpen(title="対象のCSVファイルを指定してください。", filters=dlgFilters[c("All"), ])$res
# 別のスクリプトを連続して実行する場合に備え、scriptフォルダに移動
setwd(paste(dirBase, "/script" , sep=""))

その後、Rのコンソール画面で getwd() すると、

> getwd()~
[1] "C:/Users/hogehoge/Desktop/R/script"

のように正しく変更されています。
◆ 今困っていることは、ファイル → Rコードのソースを読み込み…をクリックしても、最初の時はダイアログのフォルダが C:/Users/hogehoge/Desktop/R/csv となっています。
その後、ファイルを選択せずにダイアログを閉じ、もう一度、同じようにファイル → Rコードのソースを読み込み…を行うと、C:/Users/hogehoge/Desktop/R/script になっておりrスクリプトファイルを選択できるようになります。
◆ これを最初のRコードのソースを読み込み…で C:/Users/hogehoge/Desktop/R/script が選択できるようにしたいのですが、何か手立てはありますでしょうか?

Rの起動時の文言で気になるところ

Montecarlo (2012-09-25 (火) 23:10:06)

いつもお世話になっております。
新しいパソコンへRをインストールしたところ、以前のパソコンでは表示されなかった文言が表示されるようになってしまいました。

 [以前にセーブされたワークスペースを復帰します] 

この文言の意味するところが分からず、何か影響があるのでしょうか。 「ワークスペースを○○したい」のページは見つけましたが、あまりヒントになりませんでした。
ご存知の方がいらっしゃいましたら、なにとぞお教えいただけますようお願いいたします。 なお「sessionInfo()」は以下の通りです。

> sessionInfo()
R version 2.15.1 (2012-06-22)
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   
  • 前回 R を終了したときの状態を回復したと言うことで,気にする必要は全くありません。R を終了するときに定義した変数や関数など(要するにワークスペース)を保存し,次に R を立ち上げたときにそこから読み込んで,前の変数や関数を引き続き使えるようにすると言うことです。もし,終了時にそれ以前に定義した変数や関数などをチャラにしたいということなら,終了時に「ワークスペースのイメージファイルを保存しますか?(OS により文言は異なる)」に対して「保存しない」を選べばよいのです。 -- 河童の屁は,河童にあらず,屁である。 2012-09-26 (水) 06:44:50
  • 河童の屁は、河童にあらず、屁である様 ご回答ありがとうございました。気にしなくていいということがわかり安心いたしました。 -- Montecarlo 2012-09-26 (水) 23:34:21

正規表現にマッチした個数

yy (2012-09-20 (木) 21:15:57)

お世話になります。
たとえば

x <- "最新のコメントが付いた記事"

とした場合、正規表現を "(最新|記事)" とした場合に「マッチした個数:この場合2個」を得たいのですが、regexprや grep ではできなさそうですが何か適当な関数はありますでしょうか?
よろしくお願いします。

3次元分布図のプロットごとに日本語ラベルをつけるには

a.i (2012-09-19 (水) 04:32:46)

plot(x, y)
text(x, y, label)

での2次元分布図ではきちんと日本語表示されますが、

plot3d(x, y, z)
text3d(x, y, z, label)

でのrglパッケージで3次元分布図を作成すると日本語の部分だけが表示されないか文字化けしています。
きちんと日本語部分を表示する方法はあるでしょうか。
よろしくお願いします。

表作成

MY (2012-09-18 (火) 11:58:43)

Rを使って表を作成しpngなどのファイルに出力してグラフ化したいのですが、そういったことができる関数はありますか?
よろしくお願いします。

表示される有効桁数(2.0や4.0といった表示方法)

taked (2012-09-18 (火) 11:12:50)

小数点以下の表示をそろえるため、2.0や4.0を得たいのですが、調べても分からなかったので教えてください
よろしくお願い致します
下の例ではround関数を使っていますが、round関数でなくても、4.0(数値)の入力に対して、4.0を得られれば良いです

 > round(4.01,digit=0)
 [1] 4
 > round(4.01,digit=1)
 [1] 4
 > round(4.01,digit=2)
 [1] 4.01
 > 4.0
 [1] 4

biplotの色

超初心者 (2012-09-14 (金) 15:52:22)

いつも掲示板を参考にしています。

x<-matrix(c(2,2,3,6,2,0,6,7,9,1,4,1,4,2,1,3,5,6,3,2,6,4,5,7,5,1,1,7,3,2,0,5,0,8,1,8),9,4)
rownames(x)<-c("みかん","いよかん","はっさく","ぽんかん","桃","すもも","白桃","20世紀梨","新高梨")
colnames(x)<-c("あ町","い町","う町","え町")

library(MASS)
x.coa<-corresp(x,nf=4)
biplot(x.coa)


上記のようなサンプルで出力されるbiplot図に対して、柑橘類(みかん、いよかん、はっさく、ぽんかん)はオレンジ色、(桃、すもも、白桃)はピンク色、(20世紀梨、新高梨)にはグリーンをあてたいのですが、どうやったらできるのでしょうか?

よろしくお願いします。

Rstudio for Mac で,HRML に含まれる PNG 画像の日本語文字化け

(2012-09-11 (火) 15:22:23)

Mac OS X 10.8.1, R 2.15.1,Rstudio 0.96.330 を使っています。
R コンソールから pdf("file.png"...) のようにして書き出したファイル中の日本語は文字化けしませんが,Rstudio で,*.Rmd から *.html を生成すると,figure フォルダに作られる *.png ファイルの日本語は文字化け(トーフ)します。当然,できあがった html ファイル中の *.png 画像も文字化けしています。
Mac で作った *.Rmd を Windows に持って行って Rstudio で *.html にすると文字化けはしません。
Rstudio for Mac の場合には,フォントファミリーが日本語になっていないのだと思いますが,フォントファミリーを設定するオプションがわかりません。
R コンソールから knit2html から始めて debug でトレースしてみた限りでは,knit 関数の中にある process_file の中(あるいは更にその奧か)で *.png が作られるようですが,やはり,フォントファミリーの設定のタイミングが全くわかりませんでした。対処法はあるでしょうか。

画像に使われるフォントはpar(family="フォント名")で指定します.Macでヒラギノ角ゴシックを使う場合,

 par(family="HiraKakuPro-W3")@

とします.

ldaの結果からクロス表を作りたいのですが要素数が合いません

ウエモト (2012-09-10 (月) 16:15:33)

初めて質問させて頂きます。宜しくお願いします。
Rを使って線形判別法(lda)によってデータの分類・予測をしたいと思っております。

library(MASS)
group <- c("y","y","y","y","y","m","m","m","m","m")#人物の属性
wave <- c(y1,y2,y3,y4,y5,m1,m2,m3,m4,m5)#各要素は1列30000行の行列
traindata <- data.frame(GROUP=new{2012-09-17 (月) 21:48:32};
  • group,WAVE=wave)#データフレームにする
    (ldadata<-lda(GROUP~.,data = traindata))#ldaにかける
    plot(ldadata)#学習データの第一判別関数得点の分布
    apply(ldadata$means%*%ldadata$scaling,2,mean)#判別得点
    table(group,predict(ldadata)$class)#判別結果をクロス表にする←失敗
    このような状態です。実行すると最終行で、要素が一致せず、失敗してしまいます。
    原因は、groupが要素数10なのに対し、predict(ldadata)$classが要素数30000だということに起因しています。
    自分としては、データフレームtraindata中のwaveにおける、各要素である各行列がどちらのグループに属するのかを調べたいのですが、どうすればいいかで困っています。
    google検索などを使ってみましたが、自分の力不足でどのようにすればいいのかわかりません。
    Rの知識が少なく、低レベルな質問かとも思いますが、ご教授願います。
    見辛い文面でしたら申し訳ありません。

censored logit/probit model

taka (2012-08-20 (月) 21:31:47)

お世話になります。
Rでcensored logit modelあるいはcensored probit modelを利用したのですが、ただのlogit model、probit model、またCensored regression modelのTobit modelについては、こちらのwikiの「Rでエコノメトリクス」のページなどを参照させていただき、それぞれ

glm(y ~ x1 + x2 + x3, family=binomial(link="logit"))
glm(y ~ x1 + x2 + x3, family=binomial(linknew{2012-10-04 (木) 13:35:44}; ~

="probit"))

vglm(y ~ x, tobit(Lower=0), trace=TRUE)

などの関数を使うということまではわかったのですが、censored logit/probit modelではどのような関数を使えば良いのかがgoogleなどで調べられませんでした。
もしご教授いただけましたら幸いです。
当方、Rの基礎どころか統計の基礎も危うく、見当違いな質問でありましたら誠に申し訳ございません。

xtableライブラリで小数点の桁数を指定

meet (2012-08-07 (火) 22:35:28)

xtableライブラリで、列ごとに小数点の桁数を指定する(2列目は小数点第一位まで、3列目は第二位まで)ことはできますか?また1列目に行番号をつけないようにできますか?ご教授お願いします。

library("xtable")
sex    <- c("F", "F", "M", "M", "M")
height <- c(158.0, 162.5, 177.3, 173.3, 166.3)
weight <- c(51.22, 55.33, 72.44, 57.44, 64.33)
( x    <- data.frame(SEX=sex, HEIGHT=height, WEIGHT=weight) )
mytable <- xtable(x) 
print(mytable, type="html")

行列からデータフレームを作りたい

ランゲル・ハンス (2012-08-04 (土) 11:00:04)

いつも掲示板を参考にしています。
次のような行列Mがあります。

m <- c(16, 15, 7, 14, 9, 11, 8, 4, 1, 16, 12, 17, 10, 9, 19, 6)
M <- matrix(m, 4, 4)
M

行列Mの行をX、列をY、値をVとし、行列Mの左下が(X, Y)=(1, 1)になるように次のようなデータフレームを作りたいと思います。

X 	Y	V
1	1	14
2	1	4
3	1	17
4	1	6
1	2	7
2	2	8
3	2	12
4	2	19
1	3	15
2	3	11
3	3	16
4	3	9
1	4	16
2	4	9
3	4	1
4	4	10

方法をご教示いただけないでしょうか?
上の例は4行4列ですが、本当はn行m列でデータフレームを作りたいと思います。どうぞよろしくお願いします。

画像の特定色のピクセル数を求めたい

しまだ (2012-08-03 (金) 23:43:47)

お世話になります。
画像を与えた時、特定の色のピクセル数を得たいのです。そこで:

library(biOps)
data <- readJpeg("c:/hoge_20120101.jpg")
x <- c(1:352)
y <- c(1:288)
blue <- c(0)
grey <- c(0)
white <- c(0)
black <- c(0)
yellow <- c(0)
for ( i in x ){
for ( j in y ){
if (( data[i,j,1]==0 ) && ( data[i,j,2]==0 ) && ( data[i,j,3]==254 ))
{ blue <- ( blue + 1 ) }
if (( data[i,j,1]==192 ) && ( data[i,j,2]==192 ) && ( data[i,j,3]==192 ))
{ grey <- ( grey + 1 ) }
if (( data[i,j,1]==255 ) && ( data[i,j,2]==255 ) && ( data[i,j,3]==255 ))
{ white <- ( white + 1 ) }
if (( data[i,j,1]==0 ) && ( data[i,j,2]==0 ) && ( data[i,j,3]==0 ))
{ black <- ( black + 1 ) }
if (( data[i,j,1]==255 ) && ( data[i,j,2]==255 ) && ( data[i,j,3]==0 ))
{ yellow <- ( yellow + 1 ) }
}}

(1)処理を速くしたい
スクリプトを実行したところ、動いたのですが、処理が遅いように思います。 コーディングが素人で下手くそだと思うのですが、どこを修正すれば良いのでしょうか。
(2)png画像を扱いたい
biOpsは、png画像を読み込めないようです。オリジナルはpngのため、無理矢理jpegに変換しているのですが、pngが読み込めて同様の処理ができるパッケージはありませんでしょうか。

○○したいができない

IonaZUNNNNNNN (2012-08-03 (金) 23:28:40)

マニュアルをちゃんと嫁
そもそも、そんなことをしたいというのは変態(そんなことする必要はない)
どうしてもしたければ、関数(プログラム)を書け
関数(プログラム)を書けないなら、あきらめよう

プログラム中で、S1 - Snまでのベクトルのデータフレームの作成法

tt (2012-08-02 (木) 19:01:06)

n は、プラスの整数です。nが決まっていれば、以下でOKなのですが、異なった数字の時(実際にはせいぜい、5-10程度)には、プログラム中では対応できないですよね。

mS <- data.frame(S1,S2, ---,Sn)

一つの方法として、S1のみのデータフレームを作り、そこに1個ずつベクトルを追加していく方法はできるのですが、いかにも非効率ですよね。もう少しすっきりした方法は無いでしょうか? 
よろしくお願いします。

Rstudio serverでinstall.packages()

taked (2012-08-01 (水) 17:41:01)

Rstudio serverで、インストールできるパッケージとできないパッケージがあります。解決法があればご教授下さい。
CentOS 6.2、 rstudio-server-0.96.228-x86_64 です。
例えばggplot2はインストールできましたが、XMLなどいくつかのパッケージがうまくインストールできません。
パッケージがR自体のバージョン もしくは64bitに対応してないのでしょうか?

 > install.packages("XML")
 Installing package(s) into ‘/home/●○●○/R/library’
 (as ‘lib’ is unspecified)
 trying URL 'http://cran.cnr.Berkeley.edu/src/contrib/XML_3.9-4.tar.gz'
 Content type 'application/x-gzip' length 923501 bytes (901 Kb)
 opened URL
 ==================================================
 downloaded 901 Kb
 
 * installing *source* package ‘XML’ ...
 ** package ‘XML’ successfully unpacked and MD5 sums checked
 checking for gcc... gcc
 checking for C compiler default output file name... a.out
 checking whether the C compiler works... yes
 checking whether we are cross compiling... no
 checking for suffix of executables... 
 checking for suffix of object files... o
 checking whether we are using the GNU C compiler... yes
 checking whether gcc accepts -g... yes
 checking for gcc option to accept ISO C89... none needed
 checking how to run the C preprocessor... gcc -E
 No ability to remove finalizers on externalptr objects in this verison of R
 checking for sed... /bin/sed
 checking for pkg-config... /usr/bin/pkg-config
 checking for xml2-config... no
 Cannot find xml2-config
 ERROR: configuration failed for package ‘XML’
 * removing ‘/home/●○●○/R/library/XML’
 Warning in install.packages :
 installation of package ‘XML’ had non-zero exit status
 
 The downloaded source packages are in
 	‘/tmp/Rtmp5ZwZF3/downloaded_packages’

subset の仕様について

tadashi (2012-07-25 (水) 01:10:16)

こんな現象があります。
下記のようなsubset の使い方は間違いのようです。どこかにこの件について書いてあるでしょうか

> subset(test,test[,"flag"]==flag)


> test<-data.frame(flag=sample(c("a","b","c"),rep=T,10))
> x<-"a"
> subset(test,test[,"flag"]==x)
  flag
1    a
6    a
9    a
> flag<-"a"
> subset(test,test[,"flag"]==flag)
   flag
1     a
2     c
3     b
4     b
5     c
6     a
7     c
8     c
9     a
10    b

二元配置分散分析での residuals

分散分析 (2012-07-20 (金) 11:26:29)

Rで二元配置分散分析をしたのですが、Residuals が表示されないことがあります。その理由を教えていただけると助かります。

> a <- c(15.65, 14.40, 17.50, 15.90, 20.30, 19.45, 20.45, 19.05)
> 要因1 <- factor(c("T1", "T1", "T2", "T2", "T3", "T3", "T4", "T4"))
> 要因2 <- factor(c("Y1", "Y2", "Y1", "Y2", "Y1", "Y2", "Y1", "Y2"))

> summary(aov(a ~ 要因1 * 要因2))
            Df Sum Sq Mean Sq
要因1        3 34.026 11.3421
要因2        1  3.251  3.2513
要因1:要因2  3  0.151  0.0504

よろしくお願いします。

hoge.Rdataをダブルクリックで起動したときのファイル名取得

ふしぎ (2012-07-19 (木) 17:53:11)

例えば、hoge.RdataをダブルクリックしてRを起動した後で、このファイル名"hoge.Rdata"を取得する関数や方法はありますでしょうか?
list.files() で作業フォルダ内のファイル名は取得できますが、他にも.Rdataの拡張子を持つファイルがある場合、hoge.Rdataを一意に決定できずに困っています。解決策やヒントなど、ご教示いただければ幸いです。
動作環境は、Win XP、R version 2.14.2です。

shapeファイルのポリゴンの座標値の取り出し方

蓮井 (2012-07-15 (日) 06:45:20)

win-xp、R2.15.1、maptools0.8-16、sp0.9-99を使っています。
shapeファイルの各ポリゴンについて、
 ・面積を求める
 ・各ポリゴンがある座標値を含むかどうかを確認する
の2つの作業を実施したく、各ポリゴンの座標値を取り出し方についてご教示いただきたく、よろしくお願いします。
readShapeP最新 cat(sprintf(br;特別なことをするわけですから,なにもしないでいつでも好きなときに利用できるというわけにはいきません。でしょうね。 -- 河童の屁は,河童にあらず,屁である。 olyでポリゴンを読み、スロットpolygonsを表示させると、

> map <- readShapePoly("sample.shp")
> map1[1,] @ polygons
> 
> An object of class "Polygons"
> Slot "Polygons":
> [[1]]
> An object of class "Polygon"
> Slot "labpt":
> [1] -32959.94 -64662.22
>  <中略>
> Slot "coords":
>            [,1]      [,2]
>  [1,] -32919.99 -64656.18
>  [2,] -32920.01 -64659.99
>   <中略>
>  [13,] -32919.99 -64656.18
>   <後略>

となりますが、このcoordsの値を取り出すためにはどうすればよろしいでしょうか。
「R の S4 クラス、メソッド入門」は読んでみたのですが、よくわからず…。
よろしくお願い申し上げます。

sendmailRでメール送信したい

kbys (2012-07-13 (金) 21:23:22)

sendmailRというパッケージでgmailを使ってメール送信する方法について教えてください。マニュアル(http://cran.r-project.org/web/packages/sendmailR/sendmailR.pdf)を参考にしたのが下のコードですがうまくいきません。よろしくお願いいたします。"oooooooooo@gmail.com"は仮のアドレスです。

library(sendmailR)
sendmail(from="oooooooooo@gmail.com", to="oooooooooo@gmail.com", subject="件名",
         msg="本文", control=list(smtpServer="smtp.gmail.com"))
以下にエラー socketConnection(host = server, port = port, blocking = TRUE) : 
コネクションを開くことができません 
追加情報:   警告メッセージ: 
In socketConnection(host = server, port = port, blocking = TRUE) :
smtp.gmail.com:25 cannot be opened

コマンドプロンプトによるrscript.exeの実行について

taked (2012-07-10 (火) 15:15:07)
お世話になります.rscriptでRスクリプトを実行して、グラフを保存したいのですがうまくできません.
以下のプログラムを、task.Rとしています.

 library(ggplot2)
 ratings <- qplot(rating, data=movies, geom="histogram")
 ggsave(ratings, file="fig1.png")


Rコンソールにコピペして実行するとグラフが保存できるのですが、コマンドプロンプトで
  カレントディレクトリ>rscript task.R
として実行するとグラフが保存ができていません.コマンドプロンプトにはRコンソール実行時とほぼ同様のメッセージがでます.実行できているようですがグラフは保存されません.
  Loading required package : methods
  Saving 7 x 7 in image
  stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
誤りなどあったらアドバイス頂けないでしょうか、よろしくお願いします.Win VISTA, R 2.15.1です

多次元の同時分布を求め、ある組み合わせを取る確率を求めたい

Saito (2012-07-03 (火) 18:50:19)

いつもお世話になっております。
連続値を取る確率変数が複数種類あった場合に、その経験的な同時分布を求め、確率変数がある組み合わせをとる時の同時確率を求めたいです。
例えば、

> library(np)
> x <- rbeta(100, shape1=10, shape2=1)
> y <- rgamma(100, shape=1)
> z <- rnorm(100)
> fit<-npudens(~x+y+z)
> plot(fit)

とすると、3変量の場合の同時確率分布が得られるようですが、その後に例えば、x=0.1, y=0.5, z=-0.1となるような確率を計算をしたいと考えています(出来ればベクトルの形で与えられるようにして、x=c(0.1, 0.2, 0.3), y=c(0.1, 0.2, 0.5), z=c(0.2, 0.4, -0.1)などと与えて、3つの同時確率が出てくると嬉しいです)。

また、上記の計算はおそらく連続値を丸める必要があると思われるのですが、そうではなく、x=0.9-1.1、y=0.4-0.6、z=-0.2-0.2といったように、特定の区間に入る確率も計算したいです。

ネットを調べていると、同じような質問はあったのですが、
https://stat.ethz.ch/pipermail/r-help/2010-November/260994.html
それに対する返答が私には理解できず、手づまりとなっています。

その他、色々と調べたのですが、二変量の場合や、正規分布に従う場合などはいくつか例があるようなのですが、特に分布形を仮定しない多変量となると資料が少なく途方に暮れております。

どなたか分かる方がいらっしゃいましたらご教授頂けると幸いです。
どうぞよろしくお願い致します。

最頻出データの抽出方法

KANA (2012-07-02 (月) 20:15:29)

Rにて下記のイメージように、ひとつのテーブルにある最頻出の読み方だけを抽出したいのですが、可能でしょうか。

AAAは、あああが3件・えーえーえーが2件のため「あああ」を抽出
UUUは、うううが2件・ゆーゆーゆーが1件のため「ううう」を抽出

■データ内容イメージ
テーブル名:table1

アルファベット,読み
AAA,あああ
UUU,ううう
AAA,あああ
EEE,えええ
AAA,えーえーえー
UUU,ゆーゆーゆー
EEE,いーいーいー
AAA,えーえーえー
EEE,いーいーいー
UUU,ううう
AAA,あああ

■抽出結果イメージ
AAA,あああ
UUU,ううう
EEE,いーいーいー

以上、よろしくお願いいたします。

軸の上付き文字による指数表記について

ymbr (2012-06-30 (土) 16:35:57)

お世話になっております。
2e+05などの大きな値を実数表記する際にoptions(scipen=○○)を使って『200000』と表示のではなく、上付き文字の指数を用いることで『2×10^5(実際は上付き)』のような書式で表示することはできますでしょうか?
よろしくお願い致します。

メモリのエラーについて

taked (2012-06-29 (金) 17:46:33)

お世話になります.lm()で 線形モデルによる回帰を行いたいのですが、データフレームが大きいためかエラーとなります.

 > memory.size(max = FALSE)
 [1] 14.01
 > df <- read.csv("data.csv", as.is=T)
 > memory.size(max = FALSE)
 [1] 471.71
 # (中略)
 > memory.size(max = FALSE)
 [1] 747.44
 > LM <- lm(df$y ~ . , data = df )
  エラー:  サイズ 268.0 Mb のベクトルを割り当てることができません 
 > memory.size(max = FALSE)
 [1] 1210.4
 >memory.limit(size = NA)
 [1] 2047
 > object.size(df)
 229790744 bytes

(質問1) もともとの"data.csv"は約100MB(250列×20万行です)です.これを読込すると、object.sizeは約220MBとなりRが約450MBのメモリを使うのは正常でしょうか?
(質問2) lm()の実行前のメモリ使用量は約747MBですが、lm()などを実行するにはやはりメモリ(memory.limit=2047MB)が不足なのでしょうか?何か他の原因が考えられるでしょうか?
解決方法や回避策など(変数等を減らす以外で)アドバイスを頂けると幸いです.よろしくお願い致します.
R 2.14.1 , --max-mem-size=2047MB, Windows vista 32bit, メモリ 3.00 GB, 起動させているソフトはRとテキストエディタのみです.

2変量ガウスコピュラ関数について

esta (2012-06-29 (金) 05:07:05)

いつもお世話になっております。
2変量ガウスコピュラ関数について、確率密度関数をグラフで表示したいのですが、2変量ガウスコピュラ密度関数の中に標準正規分布の逆関数があります。この逆関数を表現できず困っています。標準正規分布の逆関数はどのように表現すればよいのでしょうか。

総当り法(ロジスティック回帰分析)でのエラー

hiro (2012-06-20 (水) 11:56:54)

Rを始めたばかりです。
現在、ロジスティック回帰分析にかけるために、総当り法でAICからモデル選択を行なっています。
変数が20個と多いので、all.logisticの結果はすべて表示されず、デフォルトだとdeviance順に20個しか表示されません。なので、print.all.logisticを使ってAICの小さい順にソートしたいのですが、エラーが出てしまいます。

> print.all.logistic(obj, sort.by=c("AIC"))
 以下にエラー data.frame(obj$deviance, obj$AIC, obj$formula) : 
   オブジェクト 'obj' がありません 

どのようにすればよいのでしょうか。
よろしくお願いいたします。

Rだけでデータクリーニング

tadashi (2012-06-17 (日) 18:15:17)

実際のデータを読み込むときには、データクリーニングをする必要があります。
そういった作業をRだけできるだけやる方法を書いてあるページはないでしょうか?
Rクックブックには、そういったのはRは不得意だから、やらないほうがいいとかいてありますが、もしかしたらあるのではないかと思っての質問です。
とりあえず、scan で全部読み込んで、それからいろいろ処理しましょうといような内容になるかと予想しています。

関数の最大値の求め方

esta (2012-06-15 (金) 16:30:16)

はじめまして。R歴1週間の初心者です。

f<-function(x)exp(-3.0346+0.0558*x-0.0003*x*x)/{1+exp(-3.0346+0.0558*x-0.0003*x*x)}
plot(f,0,250)

これを実行するとfのグラフが表示されます。mathematicaで最大値とそのときのxの値を求めたところ、x=93のときf(93)=0.391765が最大値として得られました。これをRで求めるにはどのようにしたらよいのでしょうか?

The R Tipsを一通り読んでみましたが分かりませんでした。
どうかよろしくお願いいたします。

doMC 並列処理

BCD (2012-06-14 (木) 11:06:29)

いつもお世話になっております 。
Linux、R2.1.4環境下でlibrary(doMC)のregisterDoMC()を使ってforeach~%dopar%{で並列処理を行っているのですが、Rがクラッシュしてしまい処理を終えることができません。解決策をご存知の方がいらっしゃいましたら御教授願います。

RData 形式 日本語を windows,mac,linux で共通してつかえるのか?

tadashi (2012-06-12 (火) 10:08:39)

前も聞いているのですが、できるんでしょうか? UTF-8で、CSVで保存すれば、なんとでもなるだろうというのはわかっているのですが、バイナリで共通してもてると便利だと思っています。それほどめんどうでなく、いい方法はないのでしょうか?

X軸ラベルの変更方法について

TK (2012-06-09 (土) 21:27:13)

Rコマンダーで折れ線グラフを作成した結果思うようなグラフができたのですが、X軸の目盛が0,50,100....と50刻みで表示されました。これを0,30,60,90,....360まで30刻みにしたいのですが、作成方法を詳しくお願いします。
R初心者ですのでプログラムの書き込み場所手順等、詳しくお願いします。
今は、Rコマンダーでデーターセットしてグラフ作成する事しかわかってませんので宜しくお願いします。

移動ブロック・ブートストラップについて

[[E&T]] (2012-06-06 (水) 21:34:40)

 Efron and Tibshiraniの"An Introduction to the Bootstrap"でブートストラップの勉強をしています。同書のp.99-102の「8.6 移動ブロック・ブートストラップ」に関して、質問があります。
 時系列データからランダムに3つの観測値からなるブロックを復元抽出して、新しい時系列を200系列作成します。それに対して1次の自己回帰モデルを当てはめ、回帰係数βのブートストラップ推計値を求めようというものです。
 これを追試しようとして、library(boot)のtsbootを適用すべき問題であろうと推測し、以下のようなコードを試してみました。(対象とする時系列は、library(bootstrap)のlutenhorm$V4です。)

library(bootstrap)
library(boot)
block.beta <- function(tsb){
 ar.fit <- ar(tsb, order.max=1)
 as.vector(ar.fit$ar)
}
lut.block3 <- tsboot(lutenhorm$V4, statistic=block.beta, R=200, l=3, sim="fixed")

すると、以下のようなエラーメッセージが出て、うまくいきません。

以下にエラー tsboot(lutenhorm$V4, statistic = block.beta, R = 200, l = 3,  : 
  replacement(置き換え)の長さが0です 

何か、根本的なところで私自身に誤解があるのかもしれませんが、原因がよく分かりません。どなたか解決策を御教示いただければ幸いです。

ARで伝達関数のゲインやコヒーレンス・位相を出したい

Rbeginer (2012-06-05 (火) 11:10:03)

Rで時系列解析を考えている初心者です。

AR(自己回帰モデル)を用いた時系列解析に関する質問です。
RでARを用いた伝達関数のゲインやコヒーレンス、位相を計算することを検討しています。
調べていてもどうもたどり着けません。
Rでそもそも可能でしょうか、またパッケージを導入すれば可能なのでしょうか?
ご教授いただければ幸いです。

グラフの軸の目盛を内向きにしたい

sato (2012-06-05 (火) 10:39:04)

R初心者です.
Rで折れ線グラフや棒グラフなどの単純な作図を行なっています.
軸の目盛りですが,デフォルトでは外向きとなっていると思います.
これを内向きにしたいと思って,いろいろ調べたのですが,よくわかりません.
ご教授くださいますようお願いいたいします.

グラフ内に線を書き込みたい

BCD (2012-06-02 (土) 21:07:42)

いつもお世話になっております。
グラフ内における座標について質問があります。

ベクトルがいくつかあり、それらをfor文で繰り返しを用いてそれぞれをplot関数で折れ線グラフを作成しています。各々のグラフにおいて作図領域の中の同じ位置にlines関数を用いて線を表示したいのです。それぞれのグラフでは座標位置が異なるのでlines関数に渡すxとy座標を統一することができません。for文を使って繰り返しの中で自動的に線を加えたいので各グラフにおける座標を取得することができれば良いと思うのですが方法が分からないのです。

何かご存知の方がいらっしゃれば御教授願います。

動画像のコマ送りとショートカットキー

Rstudent (2012-05-31 (木) 09:29:52)

Rのグラフィックデバイスについて質問があります。

Mpeg等の動画像を手動で一定のフレーム数進め、その都度、動画像上のある点の座標を取得したいと考えております。
座標の取得はlocator()を利用すれば可能だと思いますが、以下について方法がありましたら、ご教示願います。

1.動画像を手動で進める方法
animationパッケージのani.start()を用いてHTML上で再生すると、コマ送りのボタンが表示されますが、このようなことをRのグラフィックデバイス上でできないものでしょうか。

2.ショートカットキーを用いてグラフィックデバイスにコマンドを与える方法
動画像を進める際にショートカットキーでコマ送りを行ったり、ショートカットキーで変数に値を入れたりすることができると非常に効率的になるのですが、そのような方法はあるのでしょうか。

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

コンマ秒以下の経過時間の計算について

timesaver (2012-05-21 (月) 18:48:18)

コンマ秒以下の経過時間を計算したい場合、例えば
as.difftime(c("0:3:20", "0:3:20.5"),units="sec")
とすると結果は
Time differences in secs
[1] 200 200
attr(,"tzone")
となり秒以下は無視されます。これを解決する簡単な方法はどなたかご存じ
ありませんか?過去記事やパッケージはないかと探したのですが見当たりませんでした。

ランダム効果の推計に関して

optimal (2012-05-20 (日) 01:03:58)

期間10年、9カ国、説明変数10から成るパネルデータを分析しています。分析対象は、経済学でいうところのグラビティモデルです。

プーリング推計と固定効果推計については出力できたのですが、ランダム効果推計を行おうとすると、下記のようなエラーがでてしまいます。
パネルデータには欠損値はありません。どなたか、対処法をご教示くださいませんでしょうか?なお、期間10年対象国21カ国(上記9カ国含む)のパネルデータ、期間10年対象国12カ国のパネルデータについては推計が出来ました。

> north.pl <- plm(form, data=north, model="pooling")
series    are constants and have been removed
> north.wi <- plm(form, data=north, model="within")
series    are constants and have been removed
> north.ra <- plm(form, data=north, model="random")
series    are constants and have been removed
 以下にエラー if (sigma2$id < 0) stop(paste("the estimated variance of the",  : 
   TRUE/FALSE が必要new{2012-07-14 (土) 04:53:48};
  • tadashiさん、Ionaさんありがとうございます。pythonのsmtplibをrJythonから利用する方法を調べてこのサイト(http://r.789695.n4.nabble.com/Email-out-of-R-code-td3530671.html)で、特にファイル添付できるnutterb 氏の投稿を参考にしています。送信はできたのですが本文に日本語を送信すると???となって受信されます。エンコードを明示的に指定しないといけないことが分かり修正を試みているのですがpythonも初心者なものでうまくいっておりません。大変お手数おかけしますが、nutterb 氏の投稿の中で修正箇所・コードについて教えて頂けないでしょうか? -- kbys なところが欠損値です 使用環境は以下の通りとなります。
    R version 2.15.0 (2012-03-30)
    Platform: x86_64-pc-mingw32/x64 (64-bit)

    以上につき、ご教示お願い致します。

テキストとグラフをPDF出力する方法について

taked (2012-05-18 (金) 17:40:50)

Rを使って、簡単なレポートを作成したいと思ってます。
pdf()関数で、グラフはpdfに出力できます。
例えば、任意のテキスト、データフレーム、グラフを並べてpdfに出力することは可能でしょうか?ご教授お願いします。

df <- data.frame(x=c(1:5), y=sample(100,5))
pdf("c:/test.pdf")
  plot(df)
dev.off()

アクセスログの解析は、Rでは面倒なのでしょうか?

tadashi (2012-05-17 (木) 14:28:27)

Apacheのcombined形式のログをRで読み込んで処理しようと思っていますが、読み込んだあとの処理(正規表現)で結構苦労しています。rubyでやれば超簡単なのですが、Rでは非常に難しいのでしょうか? やりたいことは、IPアドレスやエージェントなどの項目への切り分けです。一行のアクセスログを綺麗に分解したいのです。非常に面倒かどうかを知りたいのです。

epsファイルの作成でエラー

Toy (2012-05-16 (水) 22:03:33)

Windows XPとWindows 7でR version 2.15.0を使用しています。
プロットのタイトルやラベルに日本語を使用して、これをepsファイルに保存すると、

警告メッセージ: 
1:  'mbcsToSbcs' 中の '繧ウ繝ウ繝励Λ繧、繧「繝ウ繧ケ' で変換に失敗: <e3>をドットで置き換えました
(以下略)

のようなワーニングが大量に出て、作成されたepsファイルを見ると、確かに日本語の部分がドットに置き換えられています。
これは、既知のバグなのでしょうか?
それとも、何か対処法があるのでしょうか?

地域R勉強会

tadashi (2012-05-12 (土) 14:18:50)

いろいろあるみたいです。たとえば、札幌R勉強会とかあるみたいですが、情報ご存知のかた、教えてください。

ファイル出力の方法を教えて下さい。

のざわ (2012-05-03 (木) 15:18:29)

soundと言うライブラリで、色々な音が作れる事が判りました。(だだしWindowsマシンだけ)。ここで質問なのですが、折角作ったサウンドファイルをファイルとして出力し保存したい場合のやりかたが判りません。ファイル出力の方法を教えて下さい。宜しくお願いします。

mimRのエラー

タカ (2012-05-03 (木) 06:37:33)

いつもお世話になっております。
グラフィカルモデリングを行うために、mimRをインストールしたのですが、説明書通りに打ち込んでもエラーになってしまいます。

> data(carcass)
> gmd.carc <- as.gmData(carcass)
> m.main <- fit(mim(".", data=gmd.carc))
MIM returned NULL or NA, can not continue 以下にエラー   
strsplit(mimFormula, "/") :  文字列でない引数です 

原因を御存じの方がいらっしゃいましたら、ご指導いただけると助かります。
よろしくお願いいたします。

使用環境は、
R version 2.15.0 (2012-03-30)
Platform: i386-pc-mingw32/i386 (32-bit)

dataframeの変数名を変数で指定したい

タカ (2012-04-29 (日) 21:24:21)

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

データフレームを作成する際、変数名を他の変数で指定したいのですが、

> x #xはデータフレーム
  var
1   5
2   4
3   3
4   2
5   1
> y="var2" #変数名をvar2にしたい
> z=c(1,1,1,1,1)
> x=cbind(x,y=z)
> x
  var y
1   5 1
2   4 1
3   3 1
4   2 1
5   1 1


上のように変数名をvar2として追加したいのですが、yになってしまいます。
変数で指定することはできないのでしょうか?
よろしくお願いいたします。

テストデータ作成

tadashi (2012-04-27 (金) 16:33:04)

POSデータの詳細行のようなテストデータをつくろうと思っています。簡単にRでできることなのでしょうか?ある範囲の文字列(orange,lemon,apple) があってそれに何個というのをランダムにつけます。下記のようなデータです。

orange 2
lemon 3
apple 5
lemon 2
apple 4

のようなサンプルデータです。for などつかわずにできそうな気がするのですが、できるのでしょうか?
やりたいことは、文字列一覧からランダムに文字列を選んで、それにランダムに自然数をわりあてます。Rだとなんかさくっと出来る方法があるような気がする気がしています。

複数グループの散布図の回帰直線を描くには?

yamda (2012-04-22 (日) 11:51:20)

例えば,1月〜12月に対し,変数A(90,86,88,82,75,67,55,46,36,21,17,10)と変数B(72,61,52,49,39,36,30,27,23,20,16,12)が与えられているときに,AとBの回帰直線を同時に描くにはどのようにすればよいでしょうか?ご教示いただけますと,幸甚です。

文字列とデータフレーム操作

hayashi (2012-04-19 (木) 20:51:21)

変数名y1,y2,y3,y4,y5,.....y100 があり,"y1 + y4","y1 + y2 + y4 +y5"のような文字列が与えられたときに,文字列に出現した変数名をカウントしたいのですが,どのようにすればよいか教えていただけないでしょうか

最終的に欲しい結果は,次のような結果です

  y1 y2 y3 y4 y5....y100
1  1  0  0  1  0....0
2  1  1  0  1  1....0

df <- - data.frame(y1=0, y2=0, y3=0, y4=0, y5=0)
r1 <- "y1 + y4"
r2 <- "y1 + y2 + y4 +y5"
df
x1 <- c(1, 0, 0, 1, 0)
x2 <- c(1, 1, 0, 1, 1)
df2 <- rbind(df, x1)
df2 <- rbind(df2, x2)
df2

数値の最大表示桁数

BCD (2012-04-19 (木) 18:42:45)

非常に桁数の大きな計算をしています。例えば以下のような計算式です。

> (1e+307)*(1e+307)
[1] Inf

しかしながら、計算の桁数が大きいと計算結果が全てInfになってしまい実際の数値を知ることが出来ません。どうすれば計算結果をInfでなく数値として表示させることができるのでしょうか?よろしければ御教授願います。

パッケージのインストールの不都合

さと (2012-04-18 (水) 09:58:44)

4月から新しい環境で仕事を始めたのですが、インターネットに接続してRに新しいパッケージをインストールしようとしたところ、

> utils:::menuInstallPkgs()~
 --- このセッションで使うために、CRANのミラーサイトを選んでください --- 
 以下にエラー m[, 1L] :  次元数が正しくありません 

なるエラーメッセージが出てインストールができませんでした。以前の職場ではそのような不都合はなかったのですが、原因がわからず困っております。
新しい職場のシステム管理者はR使いではないので原因がわからなかったのですが、何か解決のヒントをご教授いただければ幸いです。

実行デモ

教師 (2012-04-15 (日) 23:05:22)

講義で一連のR命令を実行していく様子を,プロジェクターでスクリーンに映したいと考えています.命令を書いたソースファイルを読み込み,一行ずつ実行し,その都度解説をしたいと考えています.
何か適当な工夫,もしくはそれ用のパッケージをご存知なら教えてください.
例えば,par(ask=TRUE) でひき続く画像出力を段階的に出力するイメージです.

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

ノンパラメトリックとパラメトリック的なモデルの比較

北原 (2012-04-13 (金) 18:09:05)

Rとより一般的な統計の質問かもしれません。
AICはパラメトリックモデルの比較に使う指数だと知っていますが、ノンパラメトリックとパラメトリックモデルの比較にはどのようなものが使えるんでしょうか(Rで計算できる指数の中で)?

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

Windows版Rのインストールガイドのおすすめ

教師 (2012-04-09 (月) 19:40:44)

今年度からRを使った講義を受け持つことになった理工系大学の教員です.もっぱら Ubuntu Linux でしかRを使った経験がないのですが,おそらく学生は Windows を使っているものが多いのではと考えています.Windows 版 (日本語化)Rのインストールは今ではほぼ自動化されているのでしょうか.それとも多少の手動による設定が必要でしょうか.Windows のバージョンでかなり違うのでしょうか.

Windows 版 R のインストールについて現在最も参考になるサイト,もしくは本を教えてください.この Wiki で見つかるインストール法が最良でしょうか.

よろしくお願いします.

95% CI の算出について

Hiro? - 2 (2012-04-09 (月) 12:22:04)

cmprsk_2.2-1 を利用して、競合リスクを考慮したevent のCumulative incidence を下記のようなプログラムで計算しました。しかし、これですと 95% CI が計算されません。どのようなコマンドを加えるとよいでしょうか。ご指摘いただいたように実際のデータをつけて再投稿いたします。
どうかお願いします。

install.packages("cmprsk")
library(cmprsk)
MYDATA
(症例は25例で、days:観察期間、censor:eventあり=1、eventなし=0、競合リスク=2)
(下記データをexcel上に作成)
   days censor group
1     3      1     1
2     4      1     1
3     4      1     1
4     5      1     1
5     6      1     1
6     5      1     1
7    11      1     1
8     6      1     1
9     2      1     1
10   19      2     1
11   21      2     1
12    3      1     1
13    2      2     1
14    3      1     1
15    3      1     1
16    7      1     1
17   28      0     1
18    5      1     1
19    5      1     1
20    3      1     1
21    3      1     1
22    4      1     1
23    5      1     1
24    4      1     1
25    3      1     1

MYDATA <- read.delim(pipe("pbpaste"), header=TRUE) #(ここで上記excelデータをコピー)
Event <- MYDATA
Event$censor <- ifelse(MYDATA$censor == 1, 1, 0)
Comp <- MYDATA
Comp$censor <- ifelse(MYDATA$censor == 2, 1, 0)
Mixed <- MYDATA
Mixed$censor <- ifelse(MYDATA$censor == 1, 1, ifelse(MYDATA$censor == 2, 2, 0))
result <- cuminc(MYDATA$days, MYDATA$censor, MYDATA$group, cencode=0)
result2 <- timepoints(result, Event$days)
result2$"est"
                  2      3      4      5      6      7     11   19    21    28 (days)
group 1 0.04 0.32 0.48 0.68 0.76 0.80 0.84 0.84 0.84 0.84(目的とする事象の累積%)
group 2 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.08 0.12 0.12(競合リスク)

国勢調査をR化したけど、

tadashi (2012-04-06 (金) 00:58:14)

windows 用のまずつくってみました。
http://www.e-stat.go.jp/SG1/estat/List.do?bid=000001034991&cycode=0
の1番目のです。
ダウンロードは下記からできます。解凍してから読み込んでください。
http://it.isogaya.co.jp/koku001.zip
文字コードがshift-jis ぽいので、Linux ,Mac では読めないと思います。共通して読めるようにするのはどうしたらいいのでしょうか?

> load("koku001.RData")
> names(koku001)
 [1] "地域コード"     "地域識別コード" "境域年次.2010." "境域年次.2000."
 [5] "地域名称"       "人口22"         "人口17"         "人口増減数"    
 [9] "人口増減率"     "面積"           "人口密度"      

下記は、全国の人口を求めています。全国の集計を除いて、都道府県で人口を集計しています。当たり前ですが、全国の集計と同じ結果になります。

t <- subset(koku001, 地域識別コード == "a")
t[1, 6]
t <-t[2:48,]
sum(as.numeric(t$人口22))

read.table での文中ダブルクオーテーションの処理

tadashi (2012-04-05 (木) 17:57:36)

アクセスログをread.table で読み込んでいますが、結構トラブルが発生して取り込めなくなっています。簡素化すると下記のようになります。

"a"     "b"     "c"
"aa"a"  "bbb"   "ccc"

read.table(file = "test.tsv", sep = "\t", header = TRUE) 

で読み込むと、

In read.table(file = "test.tsv", sep = "\t", header = TRUE) :
  'test.tsv' の readTableHeader で不完全な最終行が見つかりました

というメッセージがでます。データ作成をRuby でやっているので、そこで処理すればいいともいえるのですが、R の中ではダブルクオーテーションがはいっていて読み込みエラーがでるのは、どうやって処理するのが、よく行われているのでしょうか? 

衆議院小選挙区と国勢調査のdataset

tadashi (2012-03-26 (月) 00:17:16)

公開されているデータから、いろいろ面白い分析ができると思っています。衆議院小選挙区と国勢調査を結びつけて分析するような動きは日本であるのでしょうか? 現在、Rで小選挙区と国勢調査の地域コードをリンクしたRのdataをつくりました。こういったことに興味ある方はいらっしゃるのでしょうか?

パワースペクトルの両対数図

北原 (2012-03-21 (水) 21:39:27)

Rの初心者のものです。
Rでパワースペクトル解析(spectrum)を行い、plotで図を描いています。ところで、横軸(x軸)も対数をとった図にしたいんですが、どうすればいいんですか?
教えてください。

三次元座標の移動について

順天太郎 (2012-03-21 (水) 09:30:32)

R 及び数学には初心者です。
質問させていただきます。
幾つかの点がxyz座標上にあります。
その点を P (1,2,3), Q (3,-2,4), R (-,2.-5, 8), 1:5, sep=, S(-2,4, 8), T (-2,-4, -3)とします。
あるxyz軸に対してか平面”(平面Aとする)を設定し、その平面がA1x + B1y + C1z + D1=0とした時、A1 = -0.01、B1 = 0.05, C1= -0.03, D1 =1 すなわち  -0.01x + 0.05y – 0.03z + 1 = 0で表されるとします。
2. 座標(?)全体を回転し、平面Aが平面xyと重なり、かつP, Q, Rの三点の中心を通る様にしたいです。(平面A’)。
3. この回転運動を行ったとき、上記P,Q,R,S, Tの新しい位置はどの様になるでしょうか?またその式はどの様に求められるのでしょうか。
 宜しくお願いたします。
 順天太郎

igraphを用いたグラフ化

パン屋 (2012-03-20 (火) 20:57:18)

こんばんは。
igraphを用いてグラフを書きました。

c1 <- c("a", "a", "a", "b", "b", "c")
c2 <- c("b"," c", "d", "c", "d", "d")
c3 <- c(5, 10, 1, 0, 100, 20)
tmp1 <- data.frame(c1 = c1, c2 = c2, c3 = c3)
library(igraph)
tmp2 <- graph.data.frame(tmp1[1:6,], directed = F)
E(tmp2)$weight <- c3
print(tmp2, e = TRUE, v = TRUE)
plot(tmp2, vertex.label = V(tmp2)$name,
     layout = layout.fruchterman.reingold,
     edge.label = E(tmp2)$weight)

としたのですが、重み(c3)が結果に反映されてないようです。どのようにしたらよいかご教授いただけないでしょうか。よろしくお願いします。

R version 2.14.2 (2012-02-29)
Platform: i386-pc-mingw32/i386 (32-bit)
other attached packages:
[1] igraph_0.5.5-4

mtext で空白を挿入

サワボノリ (2012-03-16 (金) 23:49:52)

こんにちは。

plot(0, xlab="")
mtext(side=1, at=0.8, text="panda score=0", line=-4, adj=1)
mtext(side=1, at=0.8, text="cat score=1+", line=-5, adj=1)

とした時に 0 と 1 が揃うように表示したいのですが、やり方がわかりません。family="mono" であれば "panda score=0 " とするだけで良いのですが、 family="sans" の時はそうシンプルにはできません。
"+" と同じ幅の空白を挿入する方法がありましたら教えてください。
latex の \phantom{+} と同じようなものはないのでしょうか?

system()のスペース処理について

Cent (2012-03-13 (火) 13:43:33)

こんにちは。Ceut OS 5.1, R 2.13.2を使用しています。

名前にスペースが含まれているディレクトリにsystem()でunixコマンドを実行したいのですが、スペースが邪魔でエラーになってしまいます。簡素化した例を挙げると

#"aaa bbb"というディレクトリにUNIXコマンドを実行する
dir_name <- "aaa bbb"
dir.create(dir_name)
path <- paste(getwd(),"/",dir_name, sep="")
system(paste("ls ",path," > ",path,"/ls_res.txt",sep=""))

この場合、以下のエラーメッセージが表示されます。

ls: bbb: そのようなファイルやディレクトリはありません
ls: bbb/ls_res.txt: そのようなファイルやディレクトリはありません

gsub()でディレクトリ名の"aaa bbb"の半角スペースをエスケープしてからsystem()に渡せば良いと思ったのですが、

gsub("aaa bbb","aaa\ bbb",path)

gsub("aaa bbb","aaa\\\\ bbb",path)

としてみましたが、 "/root/aaa\ bbb"という文字列をつくり出すことができません。どのようにすればsystem()で処理できるでしょうか?

宜しくお願いします。

行単位の随時処理

SAS ユーザ (2012-03-10 (土) 14:10:23)

Rで行単位にテキストファイル※を読み込み随時処理を行うことは可能でしょうか?
目的は一度にメモリに乗り切らないような大量データの処理です。
AWKやSASのデータステップをイメージしています。
readlines、及びscanは一度に全ファイルを読み込みに行くので目的に沿いません。

※可能であればデータフレームに対しても同様の処理を行いたいです。

データフレームの繰り返し処理

tadashi (2012-03-09 (金) 17:36:51)

アクセスログをRで分析するために、データフレームの繰り返し処理をしたいのですが、ある程度、詳しく書いたページはないでしょうか?
for (適当な変数名 in データフレーム) { データフレーム名$変数 etc}
で使えるようなのですが、普通は使わないのでしょうか?

集計表の作成

のす (2012-03-09 (金) 06:21:30)

よろしくお願いします

   A  B  C  D
1 赤 青 黄 白
2 赤 黄 白 白
3 青 白 白 黒
4 白 黄 赤 黒
5 黒 白 青 青
6 赤 赤 赤 赤

このような表から

   A B C D
赤 3 1 2 1
青 1 1 1 1
黄 0 2 1 0
白 1 2 2 2
黒 1 0 0 2

このようなクロス集計表を得たいのですが、どのようにしたら出来るのでしょうか。
いままでtableを使って列ごとに処理した後、エクセルで成型していました。

環境は Win7 64, R2.14.1です。

Rでグラフなどの出力をMathematicaのノートブック風にできないでしょうか?

kejuyan (2012-03-08 (木) 00:06:40)

Rでグラフなどの出力をMathematicaのノートブック風にできないでしょうか?
私の環境は、Windows7,R-2.12です。
色々調べたところTeXmacs + Rという組み合わせで実現できるような記述も見つけたのですが、設定の仕方が分かりませんでした。
TeXmacs + Rに限らず、Windows7でRの出力をMathematica風に出来る方法があれば、具体的な設定方法も含めてご教授宜しくお願いします。

データフレームを正規表現を条件にして抽出する

tadashi (2012-03-05 (月) 15:26:37)

apache のアクセスログのデータフレームで、user-agent のところで、bot の文字がはいっているデータ(列)を抽出したいのですが、どうするのがいいのでしょうか、正規表現で、簡単にできると思うのですが、書き方がわかりません。
subset(データフレーム、条件(正規表現))のように書くと思われるのですが、可能ではないでしょうか?

ファイルが読み込めない

ユーザー (2012-03-05 (月) 11:11:10)

read.tableコマンドでCSVファイルを読み込もうとしましたがう下記のような警告がでてうまくいきません。どうすればよいでしょうか?
以下、コマンドと警告です。

> Ar1TestData <- read.table("C:\Users\g-fujimoto\Dropbox\Ar1TestData.csv",header=TRUE)
エラー:  "C:\U"で始まる文字列の中で8進文字なしに'\U'が使われています

WindowsからのイベントID(E) 63のエラー警告について

saito (2012-03-04 (日) 10:20:33)

OS Windows7 x64を使用しておりますがイベントID(E)63のエラー報告がされます。メッセージは次の様です。"C:program files\R\r-2.14.1\Tc\bin64\tk85.dll"のアクティブ化コンテキストの生成に失敗しました。マニフェストまたポリシーファイル"C:program files\R\r-2.14.1\Tc\bin64\tk85.dll"行9のエラーです。要素"assemblyIdentity"の属性"processerArchitecture"に無効な値"x64"が指定されています。
上記のようなエラーメッセージですが気に病むような問題なのか無視しても構わないのか判断に迷っております。何方か解決策をご教授願えないでしょうか?お忙しい中申し訳ありませんが宜しくお願い致します。

パッケージのソースコードについて

みのむし (2012-03-02 (金) 15:24:49)

Rを使って計量経済学の分析をしています。

最近ある研究者が自分の開発した手法をRで使うためのパッケージを開発して、学術目的に限り無料で配布しています。

私は彼の手法を少し改良したものを実際の分析で使いたいと考えています。
そこで、彼の作成したパッケージに手を加えたいのです。

ところが、ダウンロードしたパッケージを見ると、拡張子がrdbやrdxとなっており、もととなったコードを読むことができません。

Rのパッケージのソースコードを読むにはどうすれば良いのでしょうか?
どなたかご教示いただければ大変幸いです。

データの抽出

tsujimo (2012-03-01 (木) 00:41:25)

最近Rを触り始めました

#用語がおかしかったらお知らせください,なるべく改めて行きたいと思います

データフレームからデータを抽出するのに以下を試みましたが,結果がどうも解せません
c(1,3)で抽出しようとしたものは何が抽出されたのでしょうか?
どういうデータが取りこぼされたのでしょうか?

#結果のデータを見比べて見ても特に規則性は見えませんでした

> nrow(raw[which(raw["Q1"]==c(1), useName=TRUE, arr.ind=TRUE),])
[1] 8526
> nrow(raw[which(raw["Q1"]==c(3), useName=TRUE, arr.ind=TRUE),])
[1] 14750
> nrow(raw[which(raw["Q1"]==c(1,3), useName=TRUE, arr.ind=TRUE),])
[1] 11546

ある行"Q1"を選択肢に複数の不連続な選択番号で抽出したかったのですがどうすれば良かったのでしょうか?
複数の不連続な選択番号(ここではc(1,3))を関数の引数として受け取りたいとも考えています
ループ版は作れたのですが悲しい位遅いです
以上よろしくお願いします

セル平均と切片

そとの (2012-02-28 (火) 10:15:35)

ネットでみつけたデータを使って回帰分析の勉強をしています。

d <- read.table("http://personality-project.org/r/datasets/R.appendix2.data", header=T)
lm1 <- lm(Alertness ~ Gender * Dosage, data=d)

上の回帰分析で求められる切片と下の計算で得られる平均が同じになることはわかったのですが、両者のsdの違いが説明できません。

mean(d$Alertness[d$Gender=='f' & d$Dosage=='a'] )
summary(lm1)$coef[1,2]
sd(d$Alertness[d$Gender=='f' & d$Dosage=='a'] )

このふたつは同じになるものではないのですか?直接的なRの質問でありませんが、よろしくおねがしいます。

pisa2009データ 150メガバイトに必要はメモリは?

tadashi (2012-02-25 (土) 10:38:36)

http://pisa2009.acer.edu.au/downloads.php の生徒の回答(テキストで、非圧縮で1G程度)をRに読み込むと、150メガくらいになります。64bit のwindows で、メモリは、何Gぐらいあると快適に動くのでしょうか?

lmtest, zooのロード

akira uegaki (2012-02-22 (水) 23:10:18)

Mac OS X (10.7.2)ですが、パッケージlmtestをインストール後、これをロードしようとすると、「エラー: パッケージ '‘zoo’' をロードできませんでした」となって、ロードできません。したがって、DWテストができません。Windowsで同じことをしたときは簡単に成功したのですが。いろいろネット上で探ってみたのですが、わかりません。どなたか、ご教示ください。

dataset パッケージ

tadashi (2012-02-21 (火) 18:16:56)

dataset パッケージ library/datasets/html/00Index.html にいろいろデータがあるのですが、少し詳しく解説しているページがあったら教えてください。

Rで開発したソフトの配布

ひろ (2012-02-20 (月) 23:09:37)

Rで開発したソフトを配布したいのですが,ソースコードが丸見えなため,どうにかできないかと悩んでおります.
1.Rの難読化か機械語にコンパイル方法ありますでしょうか?私が調べてみたところ存在せずでした.(RCC(R2C++)は入手できず,Pコンパイラやらはありましたが...)
2.compilerパッケージでコンパイルしてみましたが,ソースコードが添付されてしまいます.ソースコードやコメントを除去する方法はございますでしょうか?英語のマニュアル読んでみましたが,それらしきは見当たりませんでした.もしかしたら,英語苦手なので見落としている可能性はありますが.

windowsでUTF8の対応について質問

ひろ (2012-02-20 (月) 15:10:24)

UTF8に依存した文字コードを含む文字列データをR言語を使って処理しようとしています.Linuxと違ってWinではUTF8の読み込める物の,エラーメッセージが文字化けしたりと困っておりまして,そこでお聞きしたいのが,

1. WindowsでもUTF8はShift-JISの同様に簡単に処理できますか?UTF8の文字列を扱う度iconvを使って変換が必要でしょうか?
2. Mac/Linuxへの移行を検討した方がいいでしょうか?
3. Windows版Rは近い未来にUTF8に対応するか否か?

モデルに基づくクラスター分析(mclust)"hc for EEE model is not currently supported"

おちあい (2012-02-20 (月) 04:47:08)

モデルに基づくクラスター分析の勉強のために、以下のスクリプトをあるウエブサイト(http://mjin.doshisha.ac.jp/R/29/29.html)から使用しています。

> install.packages("mclust")
> iris2<-iris[51:150, 1:4]
> library(mclust)
> plot(EMclust(iris2))
> mhc <- hc(modelName="EEE", data=iris2) # もとはEEE

ここまで来ると

 以下にエラー hcEEE(data = iris2) : 
 hc for EEE model is not currently supported

が表示されます。
他のモデル、EII、VII、VVVではこのエラーは表示されません。
EEEのモデルを採用し処理を続けるにはどうしたらよいのでしょうか。アドバイスを頂けたらと思います。
使用環境
R version 2.14.0 (2011-10-31)
Platform: x86_64-pc-mingw32/x64 (64-bit)
Windows 7

パッケージが読み込めない

ひらやま (2012-02-18 (土) 10:40:35)

分布モデル作成に有用なdismoパッケージを使っているのですが、先日から突然このパッケージが読み込めなくなってしまいました。読み込もうとすると「応答なし」となって固まってしまいます。その他のパッケージは問題なく読み込むことができ、dismoパッケージのみこの症状が出ます。
以下のことは試してみましたが解消されません。
・パッケージの更新
・パッケージの削除、再インストール
・旧バージョンのパッケージの使用
・Rをアンインストールし、再インストール 
また、昨日まで問題なくdismoパッケージが読み込めていた別PCでも、今朝になると全く同じように読み込めなくなっていました。昨夜から今朝にかけて行っていた作業内容はパッケージdismo内のmaxent()というコマンドを使って、架空データで10000回の分布モデル作成でした。
どのように対処していいものか悩んでおります。ご助言いただけますと助かります。よろしくお願いします。
使用環境
OS: Windows 7
R version: 2.14.1
対象package: dismo

barplot

ただばやし (2012-02-04 (土) 14:51:06)

積み上げグラフを書いていて、バーを塗りつぶした上に模様を描きたいと考えています。(例えば群集組成を示す図で生物の属ごとに色を与え、種ごとに模様を変えたい)
barplotで図を書いているのですが、angleとcolで模様を変えたり模様に色を付けることはできるのですが、バー自体を塗りわけることができず困っております。
ご助言頂けますと助かります。
よろしくお願いします。

ARFIMAのdについて

のざわ (2012-02-04 (土) 08:36:34)

宜しくお願いします。

以下のようなプログラムで和分を実数と整数(1)にしてsummaryを見ているのですが相関係数もAICも全く変わりません。
何かバグでもあるのでしょうか?

library(fracdiff)
AP.d <- fdGPH(AirPassengers)$d
AP.fra <- fracdiff(AirPassengers, nar = 3, dtol = AP.d, nma = 1)
summary(AP.fra)
AP.fra <- fracdiff(AirPassengers, nar = 3, dtol = 1.0, nma = 1)
summary(AP.fra)

arima関数エラーをtry()でキャッチする方法について

Yuki (2012-02-01 (水) 16:10:38)

多数の時系列データをARIMAモデルで予測したいと考えています。
try関数によりエラーとなった場合は次のデータに行くようなプログラムを作成しようとしていますが、try関数を使ってもうまく戻り値を捉えることができず困っています。対処方法をご存じの方がいらっしゃいましたらお教え頂けますか?
私が引っかかっている箇所は以下になります。
利用データ:(エラーとなった1データのみ記載しています。実際はループして処理をしています。)

myts <-ts(c(92,85,77,178,174,161,93,81,120,136,84,117,92,108,87,97,82,99,
            142,205,142,75,60,92,113,93,128,64,72,71,89,98,103,127,97,101,
            86,81,85,57,79,86,83,99,89,132,136,89,85,107,50,22,97,148,102,
            79,79,80,40,99,178,77,97,89,81,73,126,129,160,98,23,86,105,113,
            120,106,105),freq=7, start=c(7,13))

arimaで処理が正常に終わる場合、以下のように、

> (class(try(fit <- arima(myts,order=c(1,1,1),method="CSS-ML"),silent=true)))
[1] "Arima"

が返ってきますので、下記のようにifelse関数を用いてmyresに戻り値を格納し、myresの値を利用して、FALSEの場合、例外処理に行くようにしたいと考えました。

> (myres <-ifelse(class(try(fit <- arima(myts, order = c(1, 1, 1),
   method = "CSS-ML"), silent = true)) == "Arima", TRUE, FALSE) )
[1] TRUE

正常にARIMAで予測できる場合は、期待通り「TRUE」が返ってきます。 しかし、エラーとなるmethod="ML"で実行すると、期待した「FALSE」が返ってきません。

> (myres <-ifelse(class(try(fit <- arima(myts, order = c(1, 1, 1),
   method = "ML"), silent = true)) == "Arima", TRUE, FALSE) )
 以下にエラー value[[3L]](cond) :  オブジェクト 'true' がありません 
 追加情報:   警告メッセージ: 
In arima(myts, order = c(1, 1, 1), method = "ML") :
  possible convergence problem: optim gave code=1

のようにエラーとなりますが、myresの値は更新されず、以下のように前の値のままです。

> myres
[1] TRUE

念のため、myresをrmしてから確認すると、以下のようにmyresに値が帰っていないことがわかりました。

> rm(myres)
> (myres <-ifelse(class(try(fit <- arima(myts, order = c(1, 1, 1),
   method = "ML"), silent = true)) == "Arima", TRUE, FALSE) )
 以下にエラー value[[3L]](cond) :  オブジェクト 'true' がありません 
 追加情報:   警告メッセージ: 
In arima(myts, order = c(1, 1, 1), method = "ML") :
  possible convergence problem: optim gave code=1
> myres
 エラー:  オブジェクト 'myres' がありません 

利用環境は以下のとおりです。

Rのバージョン:2.14.1
OS:WindowsXP SP3
利用パッケージ:
stats     graphics  grDevices utils     datasets  methods   base  

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

expressionについて

まあくん (2012-01-26 (木) 10:32:57)

expression()の使い方ですが、

expression(paste(bgroup("(", atop(n, x), ")"), p^x, q^{n-x}))

temp <- "q^{n-x}"
expression(paste(bgroup("(", atop(n, x), ")"), p^x, temp))

のように書き換えたいのですが、うまくいきません。どのようにすればうまくいくでしょうか?

コレスポンデンス分析の散布図表示

のす (2012-01-23 (月) 18:02:07)

投稿ルールも読まず、汚い投稿をしてしまい重ね重ね申し訳ありません。
お詫びします。
先ほどの、質問ですが、以下の様な処理を行った場合、「小学生…」のグラフと、「イチゴ…」のグラフが重ねて出力されます。
ところが、「Excelで学ぶコレスポンデンス分析」(高橋信 著)の123ページには「大学生、高校生…」のグラフと、「なし、いちご…」のグラフとを別々に描いてもかまわないと記述されています。
これをRで行い、2つのグラフが得られるようにしたいのですが、よくわかりません。
よろしくお願いします。
環境はwin7(64)
R ver 2.14.1です。

library(MASS)
y <- data.frame(小学生 = c(4, 3, 4, 7, 3), 中学生 = c(5, 9, 6, 7, 4),
                高校生 = c(6, 2, 3, 7, 6), 大学生 = c(7, 1, 4, 6, 7))
rownames(y) <- c("イチゴ", "みかん", "ぶどう", "りんご", "なし")
y.coa <- corresp(y, nf=2)
biplot(y.coa, main="好きな果物調べ")

陰関数のplotについて

るち (2012-01-23 (月) 15:11:05)

はじめまして。
今T^2管理図の管理限界楕円を図に出したいのですが、その前に陰関数を図に表せるかが知りたいと思っています。

円であれば半径の指定さえすれば表示できることはわかるのですが、楕円になる陰関数での図式化の方法があれば教えていただけたらと思います。

npパッケージのnpcmstestでbwsで具体的にbandwidthを指定すると「仮引数 ”bws” が複数の実引数にマッチしました」というエラーがでてしまいます

ななみ (2012-01-21 (土) 21:55:03)

はじめまして。npパッケージでモデルの検定をしています。

model <- lm(Y ~ X, y = TRUE, x = TRUE)
X <- data.frame(X)
npcmstest(model = model, xdat = X, ydat = Y, bws = 1)

のようにバンド幅を指定すると

「仮引数 ”bws” が複数の実引数にマッチしました」

というエラーがでてしまいます。
どなたか、これに対する解決方法をご存じでしたらお教え頂きますようお願い致します。

関数内で、引数を文字列として使用したい

masa (2012-01-19 (木) 14:13:00)

お世話になります。
ヘルプやこのサイトを検索して見つけられず、非常に基本的な点ではと思うのですが、ご教授ください。
下記のような関数で階層化クラスタリングを行おうとしています。dataset の名称を文字列として、ところどころ(# here の行)で使いたいのですが、下記のままだとオブジェクトの内容そのものが入ってしまいます。

myhc <- function(varset, dataset) {
	hc <- hclust(dist(model.matrix(varset, dataset)), method= "ward")
	plot(hc, main= paste("Cluster Dendrogram for Solution", dataset, sep=""),  # here
	     xlab= paste("Observation Number in", datasetname), 
	     sub="Method = Ward; Distance = Euclidian", 
	     hang=-1)
	rect.hclust(hc, k = 5, border = "red")
	dev.print(jpeg, filename = paste(dataset, ".jpg"), width = 800, height = 600) # here
}
sessionInfo()~ [#w973c77e]
R version 2.14.1 (2011-12-22)
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] splines   graphics  stats     utils     grDevices tcltk     base

other attached packages:
[1] RODBC_1.3-4      Rcmdr_1.8-1      relimp_1.0-3     car_2.0-11
[5] survival_2.36-10 nnet_7.3-1       MASS_7.3-16

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

merge()でbyに複数列指定する方法

zou (2012-01-16 (月) 16:33:26)

merge()で、by= に複数列指定する方法があれば教えて下さい。

(x <- data.frame(
   name1 = c("tanaka", "tanaka", "tanaka", "suzuki", "suzuki", "suzuki"),
   year1 = c("2000", "2001", "2005", "2000", "2005", "2008"),
   test1 = 1:6))
(y <- data.frame(
   name2 = c("tanaka", "tanaka", "suzuki", "tanaka", "suzuki", "suzuki"),
   year2 = c("2000", "2001", "2000", "2005", "2005", "2008"),
   test2 = c(11, 12, 14, 13, 15, 16)))
 # m1 <- merge(x, y, by.x = "name1", by.y = "name2", all = TRUE))

ほしい結果は次のものです。よろしくお願いします。
なお、ソートして結合するやり方ではなく、merge()を使いたいです。

tanaka 2000 1 11
tanaka 2001 2 12
tanaka 2005 3 13
suzuki 2000 4 14
suzuki 2005 5 15
suzuki 2008 6 16

Rcmdrへの日本語を含むファイルパスのデータのインポート

mamiy (2012-01-16 (月) 11:52:48)

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

R2.13.2を使用していましたが、最近、R2.14.1に変更しました。
R2.13をアンインストールしてから、R2.14.1をインストールしたのですが、Rcmdrにて、日本語を含むパスのファイルのデータインポートでエラーが出ます。

[3] エラー:  
  <86><e3><82>ケ繝育畑繝<87>繝シ繧ソ/test.csv", header=TRUE, sep=",",
  na.strings="NA", dec=".", strip.white=TRUE)に不正なマルチバイト文字があります

R2.13.1でもR2.13.2でもエラーはなく、日本語を含むファイルパスのものも、データインポートできていたのです。何か設定などしなければいけないのでしょうか。
データのインポートを実行した際にRコマンダーのスクリプトウィンドウに表示されるスクリプトを選択して、右クリックで「実行」を選択すると、エラーは出ません。アクティブデータセットには、セットされませんが。
ご存じの方、おられましたら、お教えください。よろしくお願いいたします。
sessionInfo()の結果は

R version 2.14.1 (2011-12-22)
Platform: i386-pc-mingw32/i386 (32-bit)

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

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

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

forループとplotの組み合わせを高速化したい

宗二 (2012-01-15 (日) 16:02:34)

よろしくお願い致します。下のプログラムは、2列×10行の配列bを作り、要素が999の行で区切られたそれぞれの部分を線プロットしています。

a1 <- seq(1, 10, by = 1)
a2 <- seq(1, 5.5, by = 0.5)
b <- cbind(a1, a2)
b[1, ] <- b[4, ] <- b[10, ] <- 999
separator <- which(b[, 1] == 999)
for (i in 1 : (length(separator) - 1)) {
    data <- b[(separator[i] + 1) : (separator[i + 1] - 1), ] 
    if (i == 1) {
        plot(data, xlim = c(1, 10), ylim = c(1, 10),
             type = "l")
    } else {
        par(new = TRUE)
        plot(data, xlim = c(1, 10), ylim = c(1, 10),
             type = "l", axes = F, ann = F)
    }
}

このサンプルのように区切りの数が少なく、forループの数が少なければ問題ないのですが、区切りの数が莫大になると、作図するのにかなり時間がかかってしまいました。そこで、forループを使わないで高速化できないか考えたのですが、具体的な方法が浮かばず悩んでいます。何か高速化するいい方法、もしくはアイディアがあればご教授お願いします。

ロジスティック回帰?

竹下 (2012-01-15 (日) 08:03:21)

こんにちは。
response variable が割合 (p) なのですが、 lm( log(p / (1-p)) ~ x ) とするのは統計学的に良いやり方なのでしょうか?

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

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)

添付ファイル: fileoutput.png 926件 [詳細] fileRnw.png 1398件 [詳細] fileerror.png 1234件 [詳細] fileexample.png 1417件 [詳細] filestaxlab01.png 1678件 [詳細] filechclust01.png 1427件 [詳細] fileloglog2.png 1491件 [詳細] filehtml.png 1187件 [詳細] fileerror2.png 1221件 [詳細] filenorth.csv 813件 [詳細] filecmdscale.png 1247件 [詳細] filebarplot2.png 712件 [詳細] fileexample-cp932.Rnw 679件 [詳細] filebakaexcel.png 749件 [詳細] fileスクリーンショット 2012-06-09 22.08.55.png 570件 [詳細] fileloglog.png 1441件 [詳細] filealpha.png 1289件 [詳細] fileスクリーンショット 2012-11-20 20.03.00.png 884件 [詳細] filebiplot.png 754件 [詳細] filejibundesettei.png.png 784件 [詳細] filejibundesettei.png 1320件 [詳細] filefactor.png 1709件 [詳細] filehist.png 1688件 [詳細] fileboxplot2.png 1477件 [詳細] filerotated.png 1702件 [詳細] filebarplot1.png 1316件 [詳細] fileboxplot.png 1594件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2015-03-01 (日) 01:15:59