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

新規投稿はできません。




半角数字データを因子型にできない

お恥ずかしい? (2007-01-31 (水) 02:23:49)

> 一度、自分で検証して、検索して、頭を整理して、みんなにわかりやすく質問した方が良いと思いますよ。投稿法も。 -- 2007-01-31 (水) 05:47:26
ご迷惑おかけしました。
質問を整理します。
(1)数値型から因子型への変更は,下記のご教示の通りです。幼稚な質問をしてしまいお詫びします。色々ありがとうございました。
> factor(x)としたって,それを変数に付値しないと結果は保存されないでしょう。z$X <- factor(z$X) とやってみてくださいな。
> データといえど,全角数字を使うのは好ましくない。2番目のsumary の結果の X と Y の数字のフォントの違いに注意。 -- 2007-01-31 (水) 08:43:01

> x <- c(1,2,3,4,5)
> y <- c("1","2","3","4","5")
> z <- data.frame(X=x,Y=y)
> summary(z)
       X      Y    
 Min.   :1   1:1  
 1st Qu.:2   2:1  
 Median :3   3:1  
 Mean   :3   4:1  
 3rd Qu.:4   5:1  
 Max.   :5         
> z$X <- factor(z$X)
> summary(z)
 X      Y    
 1:1   1:1  
 2:1   2:1  
 3:1   3:1  
 4:1   4:1  
 5:1   5:1  

(2)

> x <- c(1,2,3,4,5)
> y <- c("1","2","3","4","5")
> z <- data.frame(X=x,Y=y)

の段階で, edit(z) としてデータエディタを出し,xの列の属性をnumericからcharacterに変えても反映しないのですが,操作を間違えているでしょうか;
(3)
> 外部データのハンドリングが良くないということはないでしょう。あなたが,よく理解していないからではないかと思います。具体的理由なしに,そのような判断を提示されるのはいかがなものかと思いますよ。
 すみません。CSVファイルから読み込むときにたとえデータが半角英数字であっても,「この列はカテゴリーで読む」とかの指定があったらいいのになぁ,と思ったのです。

  • > edit(z)としてデータエディタを出し,xの列の属性をnumericからcharacterに変えても反映しないのですが
    z <- edit(z) のようにしましたか。単に edit(z) とやったって,変更の結果が z に付値されないでしょう?エクセルや何かと違って,セルを変更したって,それを変数に付値しないと何にも残りません。
    > 「この列はカテゴリーで読む」とかの指定があったらいいのに
    あってもいいけど,読みこんだ後に,「この列はfactorに変換する」ということでもよいですよね。factor だったら,1,2,3 を hi, lo, med のようにわかりやすく変換することもできるし,順序を入れ替えることだってできるわけです。実に細かい操作ができますよね。
    投稿方法はまだまだ改善が必要(^_^;) 変に編集すると,表示の体裁がおかしくなる。必要な場合には,他の記事がどういう風に表現されているか見比べてから。よくわからないときにはへたに手を出さない。 -- 2007-01-31 (水) 11:57:35
  • help(read.csv) を良く読む、特にオプション colClasses の意味を。良く読まないで的外れな(質問はともかく)批判はしない。 -- 2007-01-31 (水) 12:48:34
    ファイルtest.table (CSV ファイルということだから)
       height, weight, sex
    No1,    168,     63,   1
    No2,    172,     75,   1
    No3,    170,     60,   2
    
    > x <- read.csv("test.table",
             colClasses=c("character","numeric","numeric","factor"), header=TRUE)
    > x
        height weight sex
    No1    168     63   1
    No2    172     75   1
    No3    170     60   2
    > str(x)
    `data.frame':   3 obs. of  3 variables:
    $ height: num  168 172 170
    $ weight: num  63 75 60
    $ sex   : Factor w/ 2 levels " 1"," 2": 1 1 2
  • クラスがバラバラでなく,factor にすべきものが少数なら,読んだ後に選択的に factor 関数を使う方が私は好きだ。 -- 2007-01-31 (水) 12:54:42
  • エクセルと違って,データが格納されているメモリー中のデータを直接いじっているのではなくて,複製をいじっているんだと言うことをお忘れなく。factor(x) だけしても x に変化がないのも同じ理由。変更を加えたコピーを元のオブジェクトに付値して初めて変更の成果が見られるのだ。 -- 2007-01-31 (水) 12:58:20
  • ご迷惑をお掛けいたしましたこと,深くお詫び申し上げます。やっと扱いの初期段階がクリアできそうです。懇切なご教示,皆様有難う御座いました。 -- お恥ずかしい? 2007-01-31 (水) 17:10:03
  • 迷惑ではないですよ。エクセルからRという世界に踏み込んだ一歩は素敵なことです。「属性は前でも後でも自由に、好きなだけ変更できる」Rには他にも素敵な機能がいっぱいです。これからもっとRにはまってください。そして情報を提供して下さい。 -- 2007-02-01 (木) 12:39:31

判別分析:関数qda()について

mu? (2007-01-26 (金) 01:44:56)

はじめまして.
判別分析の勉強をしているものです.
関数lda()を用いれば線形判別係数を求められますが, 関数qda()を用いたときも係数は求められるのでしょうか. ご存知の方いらっしゃいましたらご教授お願いいたします.

  • 実際にやってみれば分かるでしょう。 -- 2007-01-26 (金) 08:38:18
    > library(MASS)
    > ? qda
    > ? predict.qda
    > example(qda)
  • qdaでldaの係数が求まるか,という質問でしょうか。やってみれば,としか言いようがない。 -- 伊太利屋次郎? 2007-01-26 (金) 12:14:29
  • まったく同じ問題に直面しているが,答えがわかった? -- lene? 2007-01-31 (水) 21:16:06
  • ということは,あなたもやっては見なかったと言うこと?(プログラムソースも見るとよいでしょうけど) -- 2007-01-31 (水) 21:54:40

PDFファイルへの出力時に日本語の部分が文字化けする

R初心者X? (2007-01-24 (水) 15:39:40)

こんにちは.初歩的な質問で失礼致します.

R2.4.1を新たにインストールして,次のデータセットを実行しました.

d <- c(0:100)
v1 <- 10000 * exp(-0.05*d)

pdf(file="c:/test.pdf") 

plot (d,v1,xlab="時間", ylab="価値", 
main="", sub="指数関数", pch=22, 
ylim=c(0, 10000), type="l", lwd=2)

dev.off()

しかし,PDFファイルに出力しようとすると,

Warning messages:
1: font width unknown for character 0x90 
2: font width unknown for character 0x90 
3: font width unknown for character 0x90 
4: font width unknown for character 0x90 

という警告が出て,日本語で表記した部分が文字化けしたものがPDFファイルに出力されました.もし,お心当たりの方がいらっしゃいましたら,対処法について,ご教示頂ければ幸いです.

なお,R2.4.1はこのWikiを通じてダウンロードしたものです.
起動時の初期画面の文字化けは,Rconsole, Rdevga, Rprofile.siteの3つのファイルを上書きすることによって解消されましたが,PDFファイルの日本語部分は文字化けするようです.お忙しいところ,恐縮ですが,以上,宜しくお願い申し上げます.

  • 日本語化掲示板は読んだかな? -- 2007-01-24 (水) 15:46:02
  • 早速のコメント有難うございます.日本語化掲示板で,図形出力デバイス pdf( ) の文字化の記事を確認し,舟尾先生の2006-12-05 (月)の3つのファイルを再度上書き保存してみたのですが,同じエラーが生じました. -- R初心者X? 2007-01-24 (水) 16:17:57
  • Rconsole等を上書きしなくても、手元の環境では以下のコードで文字化けしません。 -- takahashi? 2007-01-24 (水) 17:25:21
    pdf(file="c:/test.pdf",family="Japan1") 
    plot (d,v1,xlab="時間", ylab="価値",  main="", sub="指数関数", 
      pch=22,  ylim=c(0, 10000), type="l", lwd=2)
    dev.off()
    familyの引数については、?postscriptFontsで。
  • takahashi様.ご教示,大変有難うございました.familyというオプションに対して引数にJapan1を指定すると,文字化けしないことを確認致しました.また,ご指摘いただきましたヘルプを元に,Japan1HeiMin?やJapan1GothicBBBなどの他の引数も試してみましたが,こちらも正常に描画されました.重ねて御礼申し上げます. -- R初心者X? 2007-01-24 (水) 18:05:27
  • pdfの日本語が文字化けするので日本語を使わないようにしていましたがようやく解決策(family="Japan1")を発見しました。意外とこの事に気付いていない人って多いのではないでしょうか。私が検索下手ってだけかもしませんが… -- ishi? 2009-03-26 (木) 18:17:05

変数名に変数?

R初心者? (2007-01-24 (水) 06:57:46)

大変、初歩的な質問ですみません。
データを別のデータに保存する際、データ名の変更を簡単にするために数字さえ変更すればよい、という感じにしたいです。
イメージですが、下の例ですと、コンソール画面で"ver_2"と打てば"1 2 3 4 5"と表示されるようにしたいです。どのようにすればよいでしょうか?
重ね重ね、初歩的な質問ですみません。

data<-c(1,2,3,4,5)
ver_num<-2
ver_ver_num<-data
  • こういう時のためにリストという便利なデータ型があるんですよ。 -- 2007-01-24 (水) 07:44:28
    > ver <- list(NULL)
    > ver_no <- 0
    > ver[[ver_no <- ver_no + 1]] <- 1:5
    > ver[[ver_no <- ver_no + 1]] <- matrix(1:4, 2,2)
    > ver[[ver_no <- ver_no + 1]] <- letters[1:5]
    > ver
    [[1]]
    [1] 1 2 3 4 5
    
    [[2]]
         [,1] [,2]
    [1,]    1    3
    [2,]    2    4
    
    [[3]]
    [1] "a" "b" "c" "d" "e"
    
    > ver[[1]]
    [1] 1 2 3 4 5
    > ver[[2]]
         [,1] [,2]
    [1,]    1    3
    [2,]    2    4
    > ver[[3]]
    [1] "a" "b" "c" "d" "e"
    もしどうしても ver_2 などという命名法がお気に入りなら、次のような関数を用意してつかう。変数は大局的環境中の変数になります。
    > foo <- function(var_name, ver_no, data) {
                      nam <- paste(var_name, "_", ver_no, sep="")
                      assign(nam,data,env=.GlobalEnv) }
    > foo("ver", 1, 1:5); ver_1
    [1] 1 2 3 4 5
    > foo("ver", 2, matrix(1:4,2,2)); ver_2
         [,1] [,2]
    [1,]    1    3
    [2,]    2    4
    > foo("ver", 3, letters[1:5]); ver_3
    [1] "a" "b" "c" "d" "e"
    > foo("another.ver", 1, 1:5); another.ver_1
    [1] 1 2 3 4 5
  • お節介ですが x<-data などとせず x <- data と空白を入れる癖を付けておいた方が良いですよ。x<-1 は論理値が要求される場面では比較と解釈されます。なにより他人に(ということは自分にも)読みづらい。 -- 2007-01-24 (水) 09:33:49
  • 脇道だが,x<-1が比較演算と解釈されるのはどんなときだろうか? -- 2007-01-24 (水) 11:37:00
    > x<-9
    > x
    [1] 9
    > if(x<-1) print("x is less than -1 ??")
    [1] "x is less than -1 ??"
    > x
    [1] 1
  • 失礼、逆です.論理判断 x < -1 をしているつもりなのに、代入 x <- 1 とされてしまう。 r-help でも何度も見掛けるパターンです. -- 2007-01-24 (水) 12:48:46
  • 分かりました!どうもありがとうございます。また、<- と <(-a) の危険も学ばさせていただきどうもありがとうございます。 -- 初心者? 2007-01-24 (水) 16:38:23

binary relationから行列の作成

学生? (2007-01-22 (月) 14:29:30)

こんちには。行列の解析をしています。外部で作成した行列ファイルをRに読み込ませる時に これまでは、事前に行列の形にしてRに読み込ませていました。 ただ、行列情報を普段はbinary relationの形で持っているので、 Rにこのまま読ませて、R内部で行列に変換するとういうことが したいのですが、どのようにすればよいでしょうか?

>binary relation 形式
A<tab>B<tab>0 
A<tab>C<tab>1 
B<tab>C<tab>2

>matrix 形式
  A B C
A 1 0 1
B 0 1 2
C 1 2 1
  • もし記号の種類と数が予めわかっていて、ファイルがタブ区切り、ということなら(もっとエレガントな一般解がありそうですが)。 -- 2007-01-22 (月) 15:22:17
    > x <- scan("/tmp/test.data", what=list(character(0), character(0), double(0)), sep="?t")
    Read 3 records
    > x
    [[1]]
    [1] "A" "A" "B"
    
    [[2]]
    [1] "B" "C" "C"
    
    [[3]]
    [1] 0 1 2
    > mat <- matrix(1, 3,3, dimnames=list(c("A","B","C"), c("A","B","C")))
    # 永続付値演算子 <<- の使用に注意
    > lapply(1:3, function(i) mat[x[[1]][i], x[[2]][i]] <<- mat[x[[2]][i], x[[1]][i]] <<- x[[3]][i])
    > mat
      A B C
    A 1 0 1
    B 0 1 2
    C 1 2 1
  • ひとつの方法 -- takahashi? 2007-01-22 (月) 15:34:52
    > d<-read.table("filename.tsv",as.is=T)
    > n<-sort(unique(c(d[,1:2],recursive=T)))
    > r<-matrix(1,length(n),length(n),dim=list(n,n))
    > for(i in 1:dim(d)[1])r[d[i,1],d[i,2]]<-r[d[i,2],d[i,1]]<-d[i,3]
    > r
      A B C
    A 1 0 1
    B 0 1 2
    C 1 2 1
  • 記号がA,B,C...ということだとして,以下のように(長いけど変なことはしていない) -- 2007-01-22 (月) 17:54:21
    > df <- read.table("test.dat", header=FALSE, as.is=TRUE)
    > df
      V1 V2 V3
    1  A  B  0
    2  A  C  1
    3  B  C  2
    > # A, B, C,... を 1, 2, 3, ... に変換する関数
    > sufix <- function(x) sapply(x, function(L) which(LETTERS==L))
    > 
    > # 行列のサイズ n を自動で計算
    > n <- which(LETTERS==rev(names(table(c(df[,1], df[,2]))))[1])
    > x <-diag(n) # 対角行列を作る
    > suf <- cbind(sufix(df[,1]), sufix(df[,2])) # 添え字行列を作る
    > x[suf] <- df[,3] # 該当する要素位置に数値を格納
    > x <- x+t(x) # 下三角行列にコピーして
    > diag(x) <- 1 # 対角成分は1にする
    > colnames(x) <- rownames(x) <- LETTERS[1:n] # 名前を付けて完成
    > x
      A B C
    A 1 0 1
    B 0 1 2
    C 1 2 1

適応度の高低を等高線図で描く方法

moriuchi? (2007-01-21 (日) 01:12:29)

はじめまして。
花形質が適応度に与える影響を調べています。
各個体ごとに、花弁長(X)および、雌しべ長(Y)を計測し、
適応度としては種子数(Z)を用います。
図としては、X軸に花弁長、Y軸に雌しべ長をとり、
適応度の高低を等高線図で示せないかと考えています。
contour()関数を使用すればよいのかと思うのですが、
どのように記述すれば、このような図を描けるのでしょうか?
初歩的質問で申し訳ないのですが、ご教授いただけないでしょうか。

  • example(contour) -- 2007-01-21 (日) 10:46:38
  • ? contour -- 2007-01-21 (日) 14:11:48
  • contour 関数は昇順の x,y に対して、各格子点 c(x[i], y[i]) ごとに z[i,j] が与えられているような状況を想定していますから、今の場合は使えない(?) なぜ単純に z を x,y に回帰しないのですか。その上で各格子点上の値を予測すれば、contour 関数が使えるでしょう(example(nlm) の中の例が参考になるでしょう。)。もしくは、少し想定外になりますが、kriging 法を使うということもありえます。パッケージ gstat, geoR 参照。 -- 2007-01-21 (日) 14:25:42
  • 早速ご教授いただきましてありがとうございます。回帰を実施し、Z=a1X+a2Yの関数をcontour(Z)として当てはめればよいでしょうか? -- moriuchi? 2007-01-21 (日) 17:13:44
  • X,Y のどのような関数にあてはまるかが分かっていれば,回帰分析によりパラメータ推定を行い,X, Y のグリッド点における推定値 Z を計算すると contour(X, Y, Z) で等高線を描くことができるでしょう。Z=a1X+a2Y では単なる複数の直線になるでしょうけど。
    以下に例を挙げてみましょう。z = exp((x-a)^2+(y-b)^2) という関数に従うとして,データから a, b を非線形最小二乗法により推定して,その推定値に基づいて x, y グリッドの z を計算して等高線を描きます。 -- 青木繁伸 2007-01-21 (日) 17:31:25
    > df <- structure(list(x = c(0.3, 0.4, 0.5, 0.7, 1, 1.4, 1.5, 1.6, 2, 
    + 2.1, 2.5, 2.7, 3, 3.2, 3.3), y = c(1.3, 1.4, 1.5, 1.8, 2, 2.3, 
    + 2.5, 2.7, 3, 3.1, 3.5, 3.6, 4, 4.2, 4.5), z = c(324, 167, 90, 
    + 23, 7, 2, 2, 1, 1, 1, 2, 2, 7, 18, 51)), .Names = c("x", "y", 
    + "z"), row.names = 1:15, class = "data.frame")
    > ans <- nls(z~exp((x-a)^2+(y-b)^2), df, start=list(a=2, b=3))
    > a <- ans$m$getPars()[1]
    > b <- ans$m$getPars()[2]
    > cat("a=", a, "  b=", b, "?n")
    a= 1.982297   b= 3.01758 
    > x <- seq(0,4,length=100)
    > y <- seq(1, 5, length=100)
    > z <- outer(x, y, function(xi, yj) exp((xi-a)^2+(yj-b)^2))
    > contour(x, y, z, levels=c(1.1, 2, 5, 10, 100, 1000))
    contour-example.png
  • おそらく一番手軽なのはパッケージ akima を使い、不規則空間データを空間スプライン補間することでしょう。以下の例が参考になる? -- 2007-01-21 (日) 19:11:02
    > library(akima)
    > xdata <- runif(30); ydata <- runif(30); zdata <- rnorm(30)+0.2*xdata^2-0.5*xdata*ydata+0.9*ydata^2 # 例示用人工データ
    > data <- data.frame(x=xdata, y=ydata, z=zdata)
    > with(data, contour(interp(x,y,z)))
    akima.interp.png
  • 上にある df データフレームのデータを akima を使って等高線を描こうとしてやってみましたが,うまくいかないようですね。
    data データフレームの例は,zdata <- rnorm(30)+0.2*x^2-0.5*x*y+0.9*y^2 というようなモデルがあるなら,添付画像のようになるが,akima で当てはめられた等高線図とは似てもにつかない。 -- 2007-01-21 (日) 20:11:35
    akima-example2.png
  • 失礼。最初の例の rnorm 以降の x,y は xdata,ydata でなければなりません。x=x,y=y,z=z という指定が初心者には混乱するかと思い、最後に名前を変えたのが中途半端でした。直しておきました。なお、当然ながらスプライン補間はモデル抜きの現象的な当てはめですから、結果をどう解釈するかは解釈者次第。 -- 2007-01-21 (日) 20:26:02
  • 皆様ありがとうございました。大変参考になりました。じっくりと読みチャレンジしたいと思います。 -- moriuchi? 2007-01-21 (日) 21:00:20

C.FORTRAN言語をR言語の変換する方法

初心者? (2007-01-20 (土) 12:38:21)

FORTRAN言語をR言語に取り込む機能はあるのでしょうか。
もしご存知の方いましたらご教授お願いしたします。

  • CRAN にある公式マニュアル「Writing R Extensions」を熟読。 -- 2007-01-20 (土) 14:30:14
  • FORTAN で書いたプログラム(サブルーチン)を呼び出して利用したいということ?
    Rから他言語利用を見るべし- 2007-01-20 (土) 22:01:56

64ビット版R for Windows Vista

Maerklin? (2007-01-17 (水) 01:44:31)

Windows Vistaがそろそろリリースされますが、Windowsの64ビット用のR言語バイナリは出るのでしょうか?

線形判別分析における定数変数の扱い

matuo? (2007-01-16 (火) 17:40:17)

3次元レーザスキャナで取得した30種の物体の形状データを、線形判別分析で判別させようとしています。
形状データを一定の方向からの距離画像に変換し、その画像のそれぞれの画素の値を変数としていて、 変数が1176個、データ数が1種30枚の計900あるデータフレームを学習データとしています。

data.lda <- lda(train.data, label) #labelはtrain.dataに対応するラベル

と実行すると、

以下にエラーsprintf(ngettext(length(const), "variable %d appears to be constant within groups",  : 
文字オブジェクトには %s 書式を使ってください

とエラーが出てしまいます。
レーザスキャナで物体データの取得がされなかった点は距離画像での値が0になり、 画像の端にあたるデータが定数変数になっているのですが、 ldaに渡すデータに定数変数が含まれているとエラーが出力されるようです。
手動でそういった変数を取り除くと上手く判別されるのですが、 常に画像の全く同じ場所に物体データがあるとは限らず応用性に欠けます。
関数内部で定数変数の処理を行うことは出来ないのでしょうか。

  • ああ,エラーメッセージは,エラーメッセージを出力しようとしたプログラム部分にエラーがあったという,とんでもない状態を示していますね。
    それはさておき,lda の内部にはそのような処理を行う部分は含まれていないのでは?
    ldaに渡す前に,前もってデータをチェックして,全て同じ値を持つような変数を削除してからldaに渡すようにする関数を噛ませてやればよいだけだと思います。
    しかしまあ,変数が1176で,ケースが900だと,colinear になって,解は求まらないんじゃないでしょうか? -- 2007-01-16 (火) 18:05:33
    variables are collinear in: lda.default(x, grouping, ...) 
    実に簡単なことです。以下のようなプログラムを実行すると,全く同じ値が入っている第3列が削除されているのが確認できるでしょう。
    x <- matrix(rnorm(900*1176), 900)
    colnames(x) <- paste("V", 1:1176, sep="")
    x[,3] <- 777
    x <- x[,apply(x, 2, var) != 0] # 必要なのはこれだけ!
    x[1:10, 1:5]
  • エラーメッセージが気になりますが、R のバージョンは最新でしょうね。 -- 2007-01-16 (火) 18:52:44
  • 質問者の環境はわかりませんが,R2.4.1 でも,あの通りのようですよ。 -- 2007-01-16 (火) 18:57:41
  • 例も示していただいてわかりやすい回答をありがとうございます。
    エラー文がそのような意味だったとは思いもしませんでした。
    やはり事前にデータを抜き出す必要があるのですね。そうすると実際に物体の判定を行うときに学習データと同じ列を削除する必要がありますね。
    もっと分析方法を検討しようと思います。RはR-2.4.1を使っています。 -- matuo? 2007-01-16 (火) 18:54:42
  • 気になることが。エラーメッセージの和訳は各パッケージ毎に必要ということでしょうか。質問者の例で一部が日本語なのは、その部分が標準パッケージが受け持っている部分だから、ということですかね。 -- 2007-01-16 (火) 19:56:31
  • 昨日報告があったそうです.次のリリースで直るとの事. -- なかま 2007-01-16 (火) 21:30:07
  • 翻訳はパッケージ毎に必要になります. 今回のはVRの中のメッセージで,英語の他に仏語の翻訳があります. 日本語で出たのは, ldaがメッセージを吐こうとしたらそんなフォーマットでは吐けないとR本体側のメッセージ(日本語)が出てきたと言う事です.(エラーがエラーになって出ています) -- なかま 2007-01-16 (火) 21:49:54

使う事の出来ない関数

サンダース? (2007-01-16 (火) 10:41:40)

library(KernSmooth?)の中に、bkde2Dという関数があり、その関数の中身を
print(bkde2D)で見てみました。
すると、その関数の中でlinbin2Dという関数が使われていたのですが、この関数の中身を見ようとヘルプやらネットで検索してみましたが、どうしても見ることが出来ません
また、このlinbin2Dという関数を使おうとしても、関数 "linbin2D" を見つけることができませんでしたとなって、使う事が出来ません
どうすればこの関数を利用できるのでしょうか?誰か分かる方がいれば教えてください、お願いします

  • KernSmooth:::linbin2D
    • takahashi? 2007-01-16 (火) 10:45:51
  • すばやい回答ありがとうございます!解決しました -- サンダース? 2007-01-16 (火) 10:48:37
  • Rgonzuiも -- 2007-01-16 (火) 13:36:27

sample

sarita? (2007-01-16 (火) 07:14:20)

超初心者です。sampleという関数の使い方がわかりません。sample=runif(75,min=-100,max=0)とすると、ランダムに75個の数字が出てくるわけなんですが、データファイルの特定の列を指定して、そこからある条件にあったデータをサンプリングする場合、そのデータファイルの指定をどこでしていいのかわかりません。質問の意図がわかっていただけるでしょうか。。。

  • 質問の意味はとてもわかりにくいですが、もし本当に sample() 関数の使い方がわからないのだとして、一例をあげます。-- 2007-01-16 (火) 07:45:24
    > sample(1:10, 5)  # 1:10 から5個を非復元無作為抽出
    [1]  9 10  5  8  3
    > sample(1:10, 20, rep=TRUE) # 1:10 から20個を復元無作為抽出
    [1]  1  2  6  4  3  7  2  3  7  2 10  1 10  9 10 10 10  5  2  5
  • たぶん出てきた「ランダムな75個の数値」をインデックス指定に使う方法がわからないということでしょうが...データファイルを読み込んでできたデータフレームがd, その列名がcol だとすると、d$col[sample]でいいんですよ。入門書とかRによるデータ解析とかを少しは時間を作ってきちんと読んだほうがいいでしょう。 -- 2007-01-16 (火) 08:03:52
  • 質問は,わかりにくいと言うより,意味不明です。質問にあるsampleは,関数じゃなくて,あなたが作ったオブジェクト(75個の一様乱数)に付けた名前ですね。また,後半は前半と全く関係ないのでは?「データファイルの特定の列」ではなくて「データフレームの特定の列」だろうし,「データファイルの指定」というのが,データファイルからデータフレームに読み込むときにデータファイルの指定をどこでするのかということなら,read.table 関数を調べてください,ということになるでしょう。= 「ランダムな75個の数値」は,-100〜0 という負の数値ですからね,データフレームから取り除く要素を指定したいんでしょうかね?
    タイトルの付け方も不適切だし。 -- 2007-01-16 (火) 08:08:29
  • ありがとうございます! -- sarita? 2007-01-16 (火) 08:22:43
  • 意味不明の質問ですみませんでした。 -- sarita? 2007-01-16 (火) 08:33:49
  • で,本当は(わかりやすく書けば)何が知りたかったんですか?記事のタイトルくらいは編集できますよね。 -- 2007-01-16 (火) 08:47:27

下付け文字

初心者を早く卒業したい? (2007-01-14 (日) 18:53:11)

下付け文字を含むグラフのタイトルを作りたいのですが、うまくいきません。

環境はMac os 10.48、Rのバージョンは2.31です。

標準のエディタまでは下付け文字の入力ができる(他のエディタからコピペして)のですが、コンソールでは表示がうまくいかずに普通の位置に文字がきてしまいます(グラフも同様です)。

解決法をご存じの方がいらしたらお願いします。

  • > 標準のエディタまでは下付け文字の入力ができる(他のエディタからコピペして)のですが
    それは,エディタが下付・上付文字を表示する機能を持っているからでしょう。グラフィックにおける文字表示はそのような機能を持っていないと言うことでしょう。R のコンソールだって,そんな表示はサポートしていない。添付の画像参照(右はエディタ。ちなみに,他のエディタからコピペしなくてもできる。それを実行したコンソールが左)。どこかに数式を表示する関数について書いてあったと思いますので,それを代用するということで我慢する。 -- 2007-01-14 (日) 19:19:34
    supsub.png
  • たとえばplot(ほにゃほにゃ,main=expression(AA[1]))とか...demo(plotmath)とかするといろいろわかって幸せかもしれません。 -- 2007-01-14 (日) 19:54:44
  • その方法も試してみたのですが、私の環境ではAA[1]のように表示されてしまいうまくいきませんでした。 -- 初心者を早く卒業したい? 2007-01-14 (日) 20:03:44
  • できないわけはないと思いますけど。
    関係あるかどうかは分からないが,なぜ,R2.4.1 にしないのだろうか?
    まずは,環境を整えるのが先決。添付は,以下の実行結果。 -- 2007-01-14 (日) 20:52:21
    hist(rnorm(1000),main=expression(AA[1]))
    demo6.png
  • ありがとうございました。title( " " )でタイトルを書く方法にこだわりすぎていました。 -- 初心者を早く卒業したい? 2007-01-14 (日) 21:00:59
  • 当然ながら,title でも描けますよ。ただ,title で描いたっていいけど,冗長でしょう。 -- 2007-01-14 (日) 21:07:34
    > hist(rnorm(1000), main="")
    > title(main=expression(AA[1]))
  • "The R Tips" にはこの件についての解説がありましたから、R1.9.1の段階ですでにplotmathの実装はあるはずです。まさかtitle("AA[1]")とかやってませんよね..."初心者を早く卒業したい"のならば、場当たり的な対応はせずに適切な解説をしっかり読んだほうが今後のためにもいいかもしれません。 -- 2007-01-15 (月) 02:25:06
  • 「The」のついていない方の R-Tips です。-- 2007-01-15 (月) 16:25:06

ARIMAモデルの定数項

rwiki? (2007-01-09 (火) 18:47:05)

arima で推定しようとしたら intercept が出てこなかった。なぜ?

デンドログラムの線間隔および線太さ指定方法

leau? (2007-01-05 (金) 11:10:03)

昨日はどうもありがとうございました。‘R’入門者ですのでおてやわらかにお願いいたします。
1.デンドログラム最下段の線間隔を広げたいのですが、コマンドの例をご教示ください。
2.デンドログラムの線の太さを小さくしたいのですが、コマンドの例をご教示ください。

  • 1. 図に赤で示した間隔ですか?
    clust.png
    そうだとすれば,昨日貰った回答を試してみましたか?
    >-ウィンドウサイズを広げてみましたか? -- 2007-01-04 (木) 10:47:07
    >-(補足)単純にマウスで、たとえばWebブラウザのウィンドウのサイズを変更するときのようにウィンドウの端をつかんで広げればよいのです. -- 2007-01-04 (木) 10:49:21 2. 昨日の回答の「...」を調べましたか。ヒントは par 関数の説明のところに出てくる lwd ですが,読んでください。細くできればいいですが。
    入門者だからって,説明はわかるようにやらないと。 -- 2007-01-05 (金) 13:09:51
  • 上図の間隔を広げたいのです。それをコマンドのスクリプトで行いたいのですが。昨日いただきました回答はすべて試してみましたが、データ数が500程度あるために重なってしまうためです。フォントも0.05まで試しましたが、それでも重なります。 -- leau? 2007-01-05 (金) 13:33:30
  • 500個のデータでデンドログラムを描くことが適切かどうかも考えてみましたか?大局的な傾向を探るのが目的だとしたら、cutreeで枝落としすべきでしょう。昨日から見ていて、leauさんの質問は、最初の質問と本当にやりたいことの間にずれがありそうなので…。cutreeの使い方はオンラインヘルプで自助努力してください。ここに書き込んで回答を待つ間に、自分で調べて実際に試してみたほうが早い(ことの方が多い)場合があります。 -- 2007-01-05 (金) 14:09:36
  • 図は単に画面上にウィンドウとして表示するだけではなく、PDF形式やPNG形式の画像として出力することもできる、ということを思い出せば、重なる問題はすぐに解決できると思います。線の太さはすべての線を変更するなら簡単ですが、一部だけとなると大変そうですね。 -- 2007-01-05 (金) 14:16:02
  • いろいろありがとうございました。参考になりました。しばらく自助努力いたします。 -- leau? 2007-01-05 (金) 14:35:22

クラスター分析におけるフォントサイズ変更方法について

leau? (2007-01-04 (木) 10:34:38)

クラスター分析で出力されるデンドログラムで、最下段にデータNo.が表示されますが、このときのデータNo.のフォントを小さくするにはどのようにしたらよいのでしょうか。
現状では、データの数が多いと、データNo.が重なり合って判読不能です。

  • ウィンドウサイズを広げてみましたか? -- 2007-01-04 (木) 10:47:07
  • (補足)単純にマウスで、たとえばWebブラウザのウィンドウのサイズを変更するときのようにウィンドウの端をつかんで広げればよいのです. -- 2007-01-04 (木) 10:49:21
  • 本当に,文字サイズを小さくしたいのなら,以下のように,plot 関数で cex を指定してやればよいです。? hclust で表示される plot の引数に ... があるのが,くせ者です。描画色を変えることだってできるんですよ。-- 2007-01-04 (木) 11:52:14
    > hc <- hclust(dist(USArrests), "ave")
    > plot(hc, cex=0.5, col="red")
  • ありがとうございました。上記コマンドで、フォントが小さくなる様子はわかりましたが、取り込んだデータ(Dataset)を用いて、ウォード法、ユークリッド距離で解析し、フォントを小さくする指定方法を教えていただきたいのですが。 -- leau? 2007-01-04 (木) 14:38:25
  • > 取り込んだデータ(Dataset)を用いて、ウォード法、ユークリッド距離で解析し、フォントを小さくする指定方法を教えていただきたい
    はて?何をどうしようと希望されているのか,さっぱりわかりませんが? -- 2007-01-04 (木) 14:41:25
  • 説明不足ですみません。コマンドの例> hc <- hclust(dist(USArrests), "ave") > plot(hc, cex=0.5, col="red")において、上記条件に合う引数の指定方法が知りたいのですが。 -- leau? 2007-01-04 (木) 15:08:30
  • distはデフォルトでユークリッド距離行列を返すし、ウォード法が使いたいのなら"ave"の部分をmethod="ward"にするだけのことです("ave"はmethod="average"の省略形)。少しはhelp(hclust)やhelp(dist)を見ましたか? -- 2007-01-04 (木) 15:35:53
  • フォントを小さくする方法を聞いているのだから,dist や hclust はちゃんと使えていると思っていましたが,買いかぶってましたようで。残念ですね。
    上の方もおっしゃっているように,あなたがまずすべきことは,dist と hclust のオンラインヘルプを読むことでしょう。そこまで,手取り足取り教えてもらえると思わない方がいいです。 -- 2007-01-04 (木) 16:02:09
  • お陰様でうまくできました。早速オンラインヘルプを見てみます。 -- leau? 2007-01-04 (木) 17:03:57

独立標本データの分析

yoshizaki? (2007-01-03 (水) 20:26:11)

独立標本データの分析を行いたいと思っております。
http://aoki2.si.gunma-u.ac.jp/R/two_sample.html
青木先生のページにあるのですが、

男女 血液型  データ数  平均値  標準偏差
男 A   
男 B
男 AB
男 O
女 A
女 B
女 AB
女 O
全体

ともう一歩踏み込んで出力をさせたいと思っているのですがうまくいきません。将来的には、群の数をもっと増やして分析したいと思っておりますので、どうかご教授下さい。

  • 複数の群変数で多元分類をして,各群ごとのデータ数と平均値と標準偏差を計算したいということでしょうか。同所にあります breakdown.R
    http://aoki2.si.gunma-u.ac.jp/R/breakdown.html
    ではいかがでしょうか。breakdown という名前は SPSS 由来なので,ちょっと連想できないかもしれませんね。件の場所にはたくさんありすぎまして,http://aoki2.si.gunma-u.ac.jp/R/Rstat.pdf の第8章あたりをご覧いただくとよろしいかもしれません。 -- 青木繁伸 2007-01-03 (水) 20:30:48
  • 有難うございます。一瞬で解決しました。まだまだ調査不足だったみたいですね。これからも勉強させていただきます。 -- yoshizaki? 2007-01-03 (水) 20:58:21

2行2列のchisq.testの結果が?

リュウ? (2007-01-01 (月) 21:18:35)

chisq.test(matrix(c(435,165,265,135),2,2,byrow=T))の結果です

        Pearson's Chi-squared test with Yates' continuity correction

data:  matrix(c(435, 165, 265, 135), 2, 2, byrow = T) 
X-squared = 4.1716, df = 1, p-value = 0.04111

手計算でするとカイ2乗は4.464となります
2,2行列では結果が定義に基づく手計算とずれてしまいます
何故ずれるのかご教授ください よろしくお願いします

Rでロジスティック回帰分析

R初心者? (2006-12-29 (金) 06:32:09)

The R Bookのp.308で解説されているCalibrationで、表13.23のようなQuintileのテーブルの作り方をご存知のお方はいらっしゃるでしょうか?
ネットでググってみたんですが、ロジスティック回帰モデルに関してこのようなテーブルの作成法を解説しているサイトはありませんでした。
是非ともご教授お願いします。

  • この質問に答えようとする人は,少なくともその本を持っていいる人に限られるということですね。上手な質問を試みてみましょうよ。 -- 2006-12-30 (土) 09:34:15
  • 他の本ならともかく、このWikiから生まれたThe R Bookの所持をこのWikiで前提と置いても無問題と思いますが。 -- 2006-12-30 (土) 18:01:55
  • ですよね。でもね,年末で,今,その本に,アクセスできないんですよ。さすがに二冊用意しておくと言うこともできなかったのです。また,いかに高名な本であっても,あることを前提にしての質問というのは,非難はできないでしょうが,回答確率は低くなっても(どれくらい低くなるかは別として)やむを得ないのでは? -- 2006-12-30 (土) 19:11:34
  • 普通に推測確率をfittedとかpredict(x,type="response")とかで算出して、その値に応じてあてはめた症例群を層別して、EstimatedとObservedを書いていけばよいのではないでしょうか。自分で関数を書くとよい練習になると思います。 -- 2007-01-04 (木) 11:25:17
  • 色々なご回答ありがとうございました。当方で試してみたところ、table関数とcut関数、plogis関数、predict関数を組み合わせたものを使えば、全く同じとは言わなくても、まあまあのCalibrationのテーブルを作れるようです。色々な助言ありがとうございました。 -- R初心者? 2007-01-05 (金) 22:17:00
  • せっかくなので、出来上がった関数を掲載していただけると後から勉強する人に役立つと思います. -- 2007-01-06 (土) 16:45:42
  • 関数と言うほどのものではありませんが、例えばglm関数を使ってあるデータに対して
    logisticmodel <- glm(y ~ x, family=binomial, data)
    と計算した後、その結果(logisticmodel)を利用して、
    table(cut(plogis(predict(logisticmodel)), br = c(0,.2,.4,.6,.8,1)), data$頻度)
    等と入力すると、The R Bookのp.308で解説されているCalibrationに近いテーブルが出力されます。(もっとも本の方は度数ですが)
    初心者同士で敢えて解説しますと、predict関数と言うのは、原則、線形回帰の為の関数で、従って、
    predict(logisticmodel)
    の計算で出力される数値はロジスティック曲線の定義域の部分を出力するようです。すなわち、これのままでは確率として用を成さないので、定義に従ってロジスティック累積分布関数に変換してやる必要性が生じます。
    以上より、ロジスティック累積分布の関数plogisを使って、
    plogis(predict(logisticmodel))
    と出力します。これで確率となるわけです。
    次にcut関数ですが、これは任意のデータをbr=で指定したとおり区分分けをするようです。基本的な書式は以下のようです。
    cut(データ, br=c(・・・))
    今回の場合は、『データ』の部分にロジスティック回帰での確率が入り、また、五分位の為、区分けをbr = c(0,.2,.4,.6,.8,1)とします。当然br=以下は別に任意で構いません。
    ここで区分けがされるので、あとはtable関数を用いて、今までの作業をした部分と、ロジスティック回帰の目的変数にしたかったデータ上の頻度を2つ利用してクロス集計すれば完成、となります。-- R初心者? 2007-01-08 (月) 14:46:23
  • ありがとうございました! 参考になります. -- 2007-01-10 (水) 23:31:25

コックス回帰のパワーアナリシス

サンプルサイズ? (2006-12-28 (木) 01:00:34)

コックス回帰のパワーアナリシスを行ってくれる関数、パッケージはあるのでしょうか?

  • パワーアナリシスというのは,検定力についてあれこれ検討するわけで,コックス回帰の検定力というのは何を指すのでしょうか? -- 2006-12-30 (土) 09:35:08
  • NCSS PASSなどにはCox proportional hazard model のサンプルサイズ推定計算機能もあって、そのヘルプには参考文献も載っているはずですが、それに相当するもののRでの実装は私が以前調べたときにはなかったと思います。式自体は簡単だったので自分で計算した覚えがあります。 -- 2007-01-04 (木) 11:08:18

時系列処理の高速化の仕方

? (2006-12-27 (水) 15:52:43)

R初心者です。
時系列データの処理をする上で、なかなか高速化がうまく行きません(R最適化のTIPSは読みましたが)。
行いたい処理の簡略化したものは以下のようなものです。

目的:ある指定された期間の最大値および最小値をすべての期間で算出

サンプル1
10000個のデータに対して前後15個のデータを含む31個のデータのうち最大値と最小値の計算をすべてのデータに対して(初めと終わりの15個は除く)適応する。

#初期値
N=10000
x=diffinv(rnorm(N))
span=15
#関数定義例
test1 <- function(x,span){
n=length(x)
span_max=NULL
span_min=NULL
for(i in (span+1):(n-15)){
 temp_data=x[(i-15):(i+15)]
 span_max[i]=max(temp_data)
 span_min[i]=min(temp_data)
}
outdata=cbind(span_max,span_min)
}
#実行速度のチェック
system.time(test1(x,span))
[1] 0.67 0.03 0.72   NA   NA

といった感じです。
行列を使って少しだけ高速化すると、
サンプル2

test2 <- function(x,span){
n=length(x)
temp_data=matrix(NA,n,span*2+1)
for(i in (span+1):(n-span)){
  temp_data[i,]=x[(i-span):(i+span)]
}
span_max=apply(temp_data,1,max)
span_min=apply(temp_data,1,min)
outdata=cbind(span_max,span_min)
}
system.time(test2(x,span))
[1] 0.53 0.03 0.58   NA   NA

となります。
しかし、より高速化する必要があり、このプログラムでforループで作成しているような行列を高速に作成する方法を教えて頂けませんでしょうか?

  • 行列必要なくないですか? -- takahashi? 2006-12-27 (水) 17:28:38
    test3 <- function(x,span){
      n=length(x)
      span_max<-span_min<-NULL
      sapply((span+1):(n-span),
    function(i){span_max<<-max(x[(i-span):(i+span)]);span_min<<-min(x[(i-span):(i+span)])})
      span_max<-c(rep(NA,span),span_max,rep(NA,span))
      span_min<-c(rep(NA,span),span_min,rep(NA,span))
      outdata=cbind(span_max,span_min)
    }
    
    test4 <- function(x,span){
      n=length(x)
      span_max<-c(rep(NA,span),sapply((span+1):(n-span),function(i)max(x[(i-span):(i+span)])),rep(NA,span))
      span_min<-c(rep(NA,span),sapply((span+1):(n-span),function(i)min(x[(i-span):(i+span)])),rep(NA,span))
      outdata=cbind(span_max,span_min)
    }
    
    test5 <- function(x,span){
      n=length(x)
      span_max<-span_min<-NULL
      sapply((span+1):(n-span),function(i){d<-x[(i-span):(i+span)];span_max<<-max(d);span_min<<-min(d)})
      span_max<-c(rep(NA,span),span_max,rep(NA,span))
      span_min<-c(rep(NA,span),span_min,rep(NA,span))
      outdata=cbind(span_max,span_min)
    }
    
    t2<-system.time(test2(x,span))
    t3<-system.time(test3(x,span))
    t4<-system.time(test4(x,span))
    t5<-system.time(test5(x,span))
    print(t2)
    print(t3)
    print(t4)
    print(t5)
    手元のマシンは遅いので
    [1] 2.92 0.02 2.97   NA   NA
    [1] 1.46 0.00 1.47   NA   NA
    [1] 1.45 0.00 1.45   NA   NA
    [1] 0.95 0.00 0.97   NA   NA
    という感じでした。
  • 少し外しているかも知れませんが、この計算速度で不満ということは、同じような計算を何度もしたいか、それとも実際はもっと長い時系列を扱いたいか、のどちらかでしょう。前者で、しかも時系列の長さ、スパンの長さがいつも同じ場合は次のような補助行列(これはこれで工夫と時間がかかりそう)を一度作っておくと、その後の毎回の作業が早くできます。つまり、xx <- matrix(x[X], nn,2*span+1) と置くと xx[i,] が i 番目の max,min を求めるべき時系列の連続する31個のデータになるような添字行列 X を予め作っておくというアイデアです。xx の各行の最大最小値を求めるには関数 max.col を使います。min.col 関数はないので、max.col(-xx) としています。以下の例の X は参考のための例で本当の X ではありません。しかし計算時間は基本的に同じはずです。X の作り方は、行列Tips大全の行列をベクトルとしてアクセスする方法辺りをみるとヒントになるかも知れません。 -- 2006-12-27 (水) 21:34:30
    N=10000
    x=diffinv(rnorm(N))
    span=15
    n <- length(x)
    nn <- (n-15)-(span+1)+1
    X <- matrix(sample(1:n, (2*span+1)*nn, rep=TRUE), nn, 2*span+1)
    testA <- function(x, span){
      xx <- matrix(x[X], nn, 2*span+1)
      resmax <- xx[cbind(1:nn, max.col(xx))]
      resmin <- xx[cbind(1:nn, max.col(-xx))]
      cbind(resmax, resmin)
    }
    > system.time(testA(x, span))
    [1] 0.074 0.014 0.089 0.000 0.000
  • max.colがあったんですね. 私はmatrixのトリックを使ってみました. -- なかま 2006-12-28 (木) 11:35:32
    testI <- function(x,span){
        suppressWarnings(XX<-matrix(c(x,rep(NA,span*2+1)),nrow=(length(x)+span*2),
                                    ncol=span*2+1)[(span*2+1):(length(x)),])
        cbind(XX[cbind(1:dim(XX)[1],max.col(XX))],XX[cbind(1:dim(XX)[1],max.col(-XX))])
    }
    > system.time(testI(x, span))
    [1] 0.052 0.008 0.063 0.000 0.000
  • takahashiさんありがとうございます。 -- ? 2006-12-28 (木) 11:44:40
    試してみたのですが、test2とtest3〜5では出力が違ってきます。
    間違っていたと思われる部分を自分で修正してみました test5の訂正
    test5 <- function(x,span){
      n=length(x)
      span_max<-span_min<-NULL
      sapply((span+1):(n-span),
           function(i){d<-x[(i-span):(i+span)];span_max<<-max(d);span_min<<-min(d)})
      span_max<-c(rep(NA,span),span_max,rep(NA,span))
      span_min<-c(rep(NA,span),span_min,rep(NA,span))
      outdata=cbind(span_max,span_min)
    }
    test5_ver2 <- function(x,span){
      n=length(x)
      span_max<-span_min<-NULL
      sapply((span+1):(n-span),
       function(i){d<-x[(i-span):(i+span)];span_max[i]<<- max(d);span_min[i]<<-min(d)})
      span_max<-c(span_max,rep(NA,span))
      span_min<-c(span_min,rep(NA,span))
      outdata=cbind(span_max,span_min)
    }
    という感じでよいでしょうか?
    しかし、
    > system.time(test5(x,span))
    [1] 0.27 0.01 0.28   NA   NA
    > system.time(test5_ver2(x,span))
    [1] 0.74 0.00 0.74   NA   NA
    となってしまいます。
    どのように修正すればよいでしょうか?
  • 申し訳ないです。間違ってました。test5_ver2のようにするとベクトルのサイズ変更に時間がかかると思うので、最初に確保しておくか、以下のようにすればいいと思います。 -- takahashi? 2006-12-28 (木) 12:24:57
    test6 <- function(x,span){
      ## 新しいコード
      n=length(x)
      ret<-sapply((span+1):(n-span),function(i){d<-x[(i-span):(i+span)];c(max(d),min(d))})
      outdata<-cbind(c(rep(NA,span),ret[1,],rep(NA,span)),c(rep(NA,span),ret[2,],rep(NA,span)))
    }
    
    print(system.time(test5(x,span)))
    print(system.time(test5_ver2(x,span)))
    print(system.time(test6(x,span)))
    
    [1] 0.28 0.01 0.30   NA   NA
    [1] 1.76 0.06 1.84   NA   NA
    [1] 0.18 0.00 0.18   NA   NA
    こんな感じです。
  • まとめ -- ? 2006-12-28 (木) 15:53:19
    みなさんありがとうございます。
    かなり速くなりますね。

    とりあえずtestAの答え??
    testA_ver2 <- function(x, span){
      n <- length(x)
      nn <- n-2*span
      X <- suppressWarnings(matrix(1:n,n+1,2*span+1)[1:n,])
      xx <- matrix(x[X],n,2*span+1)
      resmax <- xx[cbind(1:n, max.col(xx))][1:nn]
      resmin <- xx[cbind(1:n, max.col(-xx))][1:nn]
      cbind(resmax, resmin)
    }
    って感じですかね?
    最終的な結果は以下の通りです。(計算量を多くしました)
    N=50000
    x=diffinv(rnorm(N))
    span=15
    system.time(test1(x, span))
    [1] 12.42  3.71 16.16    NA    NA
    system.time(test6(x,span))
    [1] 1.09 0.10 1.22   NA   NA
    system.time(testA_ver2(x, span))
    [1] 0.29 0.08 0.38   NA   NA
    system.time(testI(x, span))
    [1] 0.24 0.03 0.26   NA   NA
    なかまさんのモノはシンプルで非常に早いですね。
    みなさんすばらしいですね。
    もっと勉強します。。。
  • 中間さんの超絶技巧に私同様ついていけない人も多い(ばかり?)と思うので、簡略化した例で解説。 -- 2006-12-28 (木) 22:27:53
    n <- 10; span <- 3 とします。要点は
    > XX
         [,1] [,2] [,3] [,4] [,5] [,6] [,7]
    [1,]    7    6    5    4    3    2    1
    [2,]    8    7    6    5    4    3    2
    [3,]    9    8    7    6    5    4    3
    [4,]   10    9    8    7    6    5    4
    という行列を作ることです。すると例えば x[XX[1,]] は最初に最大最小を求める副時系列 x[7:1] になります。XX を作る中間マジックは、まずベクトル
    c(x,rep(NA,2*span+1) = c(1,2,3,4,5,6,7,8,9,10,NA,NA,NA,NA,NA,NA,NA)
    を作り、R にこれを dim=c(16,7) の行列にせよ、と命令します。すると、R は困って、文句(suppressWarnings 関数はこれを無視)をいいながらもベクトルをリサイクルしながら、せっせと行列に変えていきます。結果は
    > X<-matrix(c(x,rep(NA,span*2+1)),nrow=(length(x)+span*2),ncol=span*2+1)
    Warning message: 
    行列のデータ長 [17] が行数  [16] を整数で割った、もしくは掛けた値ではありません
    > X
          [,1] [,2] [,3] [,4] [,5] [,6] [,7]
     [1,]    1   NA   NA   NA   NA   NA   NA
     [2,]    2    1   NA   NA   NA   NA   NA
     [3,]    3    2    1   NA   NA   NA   NA
     [4,]    4    3    2    1   NA   NA   NA
     [5,]    5    4    3    2    1   NA   NA
     [6,]    6    5    4    3    2    1   NA
     [7,]    7    6    5    4    3    2    1
     [8,]    8    7    6    5    4    3    2
     [9,]    9    8    7    6    5    4    3
    [10,]   10    9    8    7    6    5    4
    [11,]   NA   10    9    8    7    6    5
    [12,]   NA   NA   10    9    8    7    6
    [13,]   NA   NA   NA   10    9    8    7
    [14,]   NA   NA   NA   NA   10    9    8
    [15,]   NA   NA   NA   NA   NA   10    9
    [16,]   NA   NA   NA   NA   NA   NA   10
    となります。ここで不要な行 1:6, 11:16 を捨てれば目的の行列が得られます。お見事。教訓:何とかとリサイクル規則も使いよう。NA 値を使った利点を生かし、まず X を求めてから
    XX <- X[complete.cases(X),]
    としても良かったわけです。

AIC以外での次数決定

蟹玉? (2006-12-19 (火) 14:56:13)

現在、時系列モデルのひとつであるARモデルを用いて予測を行っていますが、デフォルトであるAIC以外の情報量基準(BICなど)を利用するにはどのように記述すればよいのでしょうか?
ヘルプ[help(ar)]にはAICを使わない方法は書いてありましたが、他の情報量基準の使い方は書かれていませんでした。

  • どこかのパッケージに用意されている可能性はありますが、基本的には自分でプログラムすることになるでしょう。 -- 2006-12-19 (火) 15:59:30
  • AIC=ー2(モデルの最大尤度)+2k (k はモデルのパラメータ数)、BIC=-2(モデルの最大尤度)+log(n)k (n はデータ数)ですから、BIC=AIC-2k+log(n)k の関係があります。ですからAIC値がわかればBIC値はすぐに計算できます。 -- 2006-12-19 (火) 18:11:46
  • BICを求めるならlibrary(fSeries)でもよいのでは? -- 2006-12-19 (火) 18:30:32
  • 質問者は BIC を如何にして求めるかを聞いているのではなく,ar 関数において AIC の代わりに AIC を使う方法を聞いているのかと思いますが? -- 2006-12-19 (火) 21:39:30
  • Rのソースコードのsrc/library/stats/R/ar*.Rの中でAICを計算してるところをBICを計算するように書き換えて関数定義を上書きすれば可能です。 -- takahashi? 2006-12-19 (火) 22:05:42
  • ありがとうございます。仰られる通りでして、ar関数でのBICの使い方を質問したかったのです。 -- 蟹玉? 2006-12-20 (水) 15:17:55
  • 上記のar*.RファイルはRの解凍後ディレクトリに発見しましたが、関数定義を書き換えた後はどのように反映すればよいのでしょうか? -- 蟹玉? 2006-12-20 (水) 15:21:20
  • パッケージングするのが王道ですが、書き換えた関数群のRソースコードのファイルを source("...")で読み込むのが手っ取り早いでしょう。この際、名前空間の関係でexportされてないシンボルが幾つかあって実行時エラーになります(R_eurekaとか)。そいつらはstats:::R_eurekaとしてあげます。 -- takahashi? 2006-12-20 (水) 15:50:28
  • 要するに、以下のコードを走らせればよいです。BICの計算は間違ってるかも。
    BICでモデル選択してることがわかるかと思います -- takahashi? 2006-12-20 (水) 16:07:36
    ar.yw.default <-
        function (x, aic = TRUE, order.max = NULL, na.action = na.fail,
                  demean = TRUE, series = NULL, ...)
    {
        if(is.null(series)) series <- deparse(substitute(x))
        ists <- is.ts(x)
        x <- na.action(as.ts(x))
        if(ists)  xtsp <- tsp(x)
        xfreq <- frequency(x)
        x <- as.matrix(x)
        if(!is.numeric(x))
            stop("'x' must be numeric")
        if(any(is.na(x))) stop("NAs in 'x'")
        nser <- ncol(x)
        if (demean) {
            xm <- colMeans(x)
            x <- sweep(x, 2, xm)
        } else xm <- rep(0, nser)
        n.used <- nrow(x)
        order.max <- if (is.null(order.max)) floor(10 * log10(n.used))
                     else round(order.max)
        if (order.max < 1) stop("'order.max' must be >= 1")
        xacf <- acf(x, type = "covariance", lag.max = order.max, plot = FALSE,
                    demean = demean)$acf
        if(nser > 1) {
            ## multivariate case
            snames <- colnames(x)
            A <- B <- array(0, dim = c(order.max + 1, nser, nser))
            A[1, , ] <- B[1, , ] <- diag(nser)
            EA <- EB <- xacf[1, , , drop = TRUE]
            partialacf <- array(dim = c(order.max, nser, nser))
            xaic <- numeric(order.max + 1)
            solve.yw <- function(m) {
                # Solve Yule-Walker equations with Whittle's
                # generalization of the Levinson(-Durbin) algorithm
                betaA <- betaB <- 0
                for (i in 0:m) {
                    betaA <- betaA + A[i + 1, , ] %*% xacf[m + 2 - i, , ]
                    betaB <- betaB + B[i + 1, , ] %*% t(xacf[m + 2 - i, , ])
                }
                KA <- -t(qr.solve(t(EB), t(betaA)))
                KB <- -t(qr.solve(t(EA), t(betaB)))
                EB <<- (diag(nser) - KB %*% KA) %*% EB
                EA <<- (diag(nser) - KA %*% KB) %*% EA
                Aold <- A
                Bold <- B
                for (i in 1:(m + 1)) {
                    A[i + 1, , ] <<- Aold[i + 1, , ] + KA %*% Bold[m + 2 - i, , ]
                    B[i + 1, , ] <<- Bold[i + 1, , ] + KB %*% Aold[m + 2 - i, , ]
                }
            }
            cal.aic <- function() { # omits mean params, that is constant adj
                det <- abs(prod(diag(qr(EA)$qr)))
                return(n.used * log(det) + 2 * m * nser * nser)
            }
            cal.resid <- function() {
                resid <- array(0, dim = c(n.used - order, nser))
                for (i in 0:order) {
                    resid <- resid + x[(order - i + 1):(n.used - i),
                                       , drop = FALSE] %*% t(ar[i + 1, , ])
                }
                return(rbind(matrix(NA, order, nser), resid))
            }
            order <- 0
            for (m in 0:order.max) {
                xaic[m + 1] <- cal.aic()
                if (!aic || xaic[m + 1] == min(xaic[1:(m + 1)])) {
                    ar <- A
                    order <- m
                    var.pred <- EA * n.used/(n.used - nser * (m + 1))
                }
                if (m < order.max) {
                    solve.yw(m)
                    partialacf[m + 1, , ] <- -A[m + 2, , ]
                }
            }
            xaic <- xaic - min(xaic)
            names(xaic) <- 0:order.max
            resid <- cal.resid()
            if(order > 0 ) {
                ar <- -ar[2:(order + 1), , , drop = FALSE]
                dimnames(ar) <- list(1:order, snames, snames)
            } else ar <- array(0, dim=c(0, nser, nser),
                               dimnames=list(NULL, snames, snames))
            dimnames(var.pred) <- list(snames, snames)
            dimnames(partialacf) <- list(1:order.max, snames, snames)
            colnames(resid) <- colnames(x)
        } else {
            ## univariate case
            r <- as.double(drop(xacf))
            z <- .Fortran(stats:::R_eureka, ## ここも変更
                          as.integer(order.max),
                          r, r,
                          coefs=double(order.max^2),
                          vars=double(order.max),
                          double(order.max))
            coefs <- matrix(z$coefs, order.max, order.max)
            partialacf <- array(diag(coefs), dim=c(order.max, 1, 1))
            var.pred <- c(r[1], z$vars)
            xaic <- n.used * log(var.pred) + 2 * (0:order.max) + 2 * demean
    
            ## ここから
            tmp<-data.frame(order=0:order.max,AIC=xaic)
            xaic <- n.used * log(var.pred) + log(n.used) * (0:order.max) + 2 * demean
            tmp$BIC<-xaic
            print(tmp)
            ## ここまで
            
            xaic <- xaic - min(xaic)
            names(xaic) <- 0:order.max
            order <- if (aic) (0:order.max)[xaic == 0] else order.max
            ar <- if (order > 0) coefs[order, 1:order] else numeric(0)
            var.pred <- var.pred[order+1]
            ## Splus compatibility fix
            var.pred <- var.pred * n.used/(n.used - (order + 1))
            if(order > 0)
                resid <- c(rep(NA, order), embed(x, order+1) %*% c(1, -ar))
            else resid <- as.vector(x)
            if(ists) {
                attr(resid, "tsp") <- xtsp
                attr(resid, "class") <- "ts"
            }
        }
        res <- list(order=order, ar=ar, var.pred=var.pred, x.mean = drop(xm),
                    aic = xaic, n.used=n.used, order.max=order.max,
                    partialacf=partialacf, resid=resid, method = "Yule-Walker",
                    series=series, frequency=xfreq, call=match.call())
        if(nser == 1 && order > 0)
            res$asy.var.coef <-
                solve(toeplitz(drop(xacf)[seq_len(order)]))*var.pred/n.used
        class(res) <- "ar"
        res
    }
    
    print(ar(lh))
    
    実行結果
    > source('tmp.r')
       order       AIC       BIC
    1      0 -56.12519 -56.12519
    2      1 -73.43620 -71.56500
    3      2 -73.89383 -70.15143
    4      3 -74.43186 -68.81825
    5      4 -72.94150 -65.45669
    6      5 -71.21907 -61.86306
    7      6 -69.43864 -58.21144
    8      7 -67.96236 -54.86395
    9      8 -65.96929 -50.99968
    10     9 -65.69066 -48.84985
    11    10 -63.69097 -44.97896
    12    11 -61.89799 -41.31478
    13    12 -59.94707 -37.49266
    14    13 -57.97006 -33.64445
    15    14 -56.38814 -30.19133
    16    15 -56.99202 -28.92400
    17    16 -55.08691 -25.14769
    
    Call:
    ar(x = lh)
    
    Coefficients:
         1  
    0.5755  
    
    Order selected 1  sigma^2 estimated as  0.2079 

rep.aovは何処にある?

高井? (2006-12-16 (土) 09:55:35)

反復測定データの分散分析を考えており,ここを検索したら「Rの統計解析関数Tips」でrep.aovという関数の記述が見つかりました。そこで,CRANのパッケージ,R Site,仲間さんのCRAN Code searchを試してみましたがいずれも検索できませんでした。rep.aovは何処にあるのでしょうか?

  • Rの統計解析関数Tips」にあります。 -- 2006-12-16 (土) 10:35:18
  • rep.aov で,ググったら,一発目に出てきます。そのあと,ページ内で検索というブラウザの機能を使ったら,目的のところまで,瞬間移動できますが。 -- 2006-12-16 (土) 17:39:28
  • ああ〜。あのコードが関数なんですね。どなたかの自作関数ということでしたか。恥ずかしい。 -- 高井? 2006-12-16 (土) 20:24:41

yatex-modeとnoweb-minor-modeの両立

S? (2006-12-13 (水) 17:42:12)

emacs上でRnwファイルを編集しているのですが、noweb-minor-modeを実行すると、勝手にtex-modeがロードされます。使い慣れたyatex-modeを使いつつ、noweb-modeで楽に編集したいのですが、yatexとnowebを両立した設定例など助言をいただければ幸いです。VineLinux3.2 + emacs-21.3-0vl7 + emacs-ess-5.3.2です。
よろしくお願いします。このWikiのSweave関連ページには目を通しましたが、解決には至りませんでした。

  • わたしのは, 5.3.0ですが,
    (defun Rnw-mode ()
      "Major mode for editing Sweave(R) source.
    See `noweb-mode' and `R-mode' for more help."
      (interactive)
      (require 'ess-noweb)
      (noweb-mode 1)
      (noweb-set-doc-mode 'yatex-mode)
      (noweb-set-code-mode 'R-mode))
    とか書いてあります.-- なかま 2006-12-13 (水) 19:52:11
  • なかまさま、お陰様で私の環境でもうまくいきました。感謝申し上げます。 -- S? 2006-12-14 (木) 11:10:54

制約条件つき関数の描画または同時描画

figure? (2006-12-13 (水) 15:10:10)

はじめまして。関数の描画についての質問です
(1)H(x) >0で1 x<=0で0の描き方。(line文以外)
(2)円をその半径を超えた領域内に描くとき、consoleにNANを出させないで描く方法?
(3)大きさの違う二つ円を円のなかに円が含まれる場合にmatplotで描く方法?
大変初歩的な質問ですみませんが、よろしくお願いします.
xの条件付きで関数定義をしても、plot()させることができませんでした。

  • (1)curve(ifelse(x>0,1,ifelse(x<=0,0,NA)),-1,2) -- takahashi? 2006-12-13 (水) 16:18:31
  • (3)matplot(data.frame(sin(2*pi*0:180/180),0.5*sin(2*pi*0:180/180)),data.frame(cos(2*pi*0:180/180),0.5*cos(2*pi*0:180/180)),type="l") -- takahashi? 2006-12-13 (水) 16:18:56
  • (2)plot(sin(2*pi*0:180/180),cos(2*pi*0:180/180),xlim=c(-3,3),ylim=c(-3,3),type="l") いまひとつ質問の意味がわからないけど・・・ -- takahashi? 2006-12-13 (水) 16:24:24
  • (1)は,plot.stepfun(0, 0:1) あと,好きなだけお化粧を。(2), (3) は,いみふめ -- 青木繁伸 2006-12-13 (水) 19:17:22
  • どうも有難うございます。うっかりしていたのは y=sqrt(r^2-x^2)を使うとばかり思い込んでいたのでif 文がベクトルに使えないので困っていました。ifelseがあること、sin,cosで描けば良いことをうっかりしました。ところで、ついでにお伺いしたいのですが、このように絵がいた場合、円の中を好みの色で塗りつぶせるのでしょうか?円の中に含まれる小円の外を塗りつぶしたりRGB指定で。あまりごてごて命令を書かないで塗りつぶしたいものですがいかがでしょうか -- figure? 2006-12-14 (木) 20:57:52
  •  -- takahashi? 2006-12-14 (木) 23:23:37
    plot.new()
    polygon(0.5+0.5*sin(2*pi*0:180/180),0.5+0.5*cos(2*pi*0:180/180),col="red",border=F)
    polygon(0.5+0.25*sin(2*pi*0:180/180),0.5+0.25*cos(2*pi*0:180/180),col="white",border=F)
  • 以下のようなものでも,「ごてごて」してるとは思いませんし,かえって簡潔では? -- 2006-12-15 (金) 00:13:41
    circle <- function(cx, cy, f, col)
     {
    	theta <- seq(0, 2*pi, length=400)
    	polygon(cx+f*sin(theta), cy+f*cos(theta), col=col, border=FALSE)
    }
    plot.new()
    circle(0.5, 0.5, 0.5, "red")
    circle(0.5, 0.5, 0.25, "white")

ARモデルのAICについて

おやじ? (2006-12-12 (火) 17:39:53)

AR()で予測式をつくり、そのAICを$AICで返した時に、最小値が0となっているのですが、これは最小となったAICを基準として差分をとっているんでしょうか?教えてください、お願いします。

  • ARという関数は私は見つけられなかったのですが,arという関数でしたら,そのhelp に,
    Order selection is done by AIC if aic is true. This is problematic, as of
    the  methods here only ar.mle performs true maximum likelihood estimation.
    The AIC is computed as if the variance estimate were the MLE, omitting
    the determinant term from the likelihood. Note that this is not the same
    as the Gaussian likelihood evaluated at the estimated parameter values.
    In ar.yw the variance matrix of the innovations is computed from the
    fitted coefficients and the autocovariance of x.
    と書いてあるのですが,これがあなたの質問の回答なんでしょうかね。 -- 2006-12-12 (火) 17:47:45
  • その通りです。AIC 値は差だけが意味を持ちますから、わかりやすく最小値がゼロになるように調整してるわけです。例えば R のソースの ar.mle.R のコードを見ると、各当てはめ次数毎に AIC 値をベクトル xaic に入れた後、xaic <- xaic - min(xaic) と差分をとっていますから、最小値は当然ゼロになります。 -- 2006-12-12 (火) 20:27:17
  • ありがとうございます!!やはり差をとっているんですね!すっきりしました! -- おやじ? 2006-12-13 (水) 16:15:37

閾値を設定したクラスタリング解析について

学生? (2006-12-12 (火) 11:44:01)

関数hclust()を使って階層的クラスター分析をする際にcutree()関数を使って事前にクラスタ数を指定して個体がどこのクラスターに属するかに関する情報を返してくれますが、これとは別に事前に閾値を設定してその閾値で、クラスター情報を返すにはどうすればよいのでしょうか?

  • このページの下の方に(目次参照)「クラスター図を横向きで描画したい」という質問があります。hclust はデンドログラムを書くために必要な情報を全て返してくれるので,それを使って,どんなことでもできると思いますよ。 -- 2006-12-12 (火) 12:19:06
  • なぁんだ。? cutree しなかったんですか?ちゃんと,h という引数があるじゃありませんか。てっきり,サポートされていないかと思いましたよ。 -- 2006-12-12 (火) 12:34:09

repeated-measures ANCOVA

s.y.? (2006-12-11 (月) 15:52:15)

皆様、
repeated-measures ANCOVAをRで行うための簡単な方法はありますか。今、具体的に問題に遭遇しているわけではないのですが、自分はRだけでやっていけるのか、それともSPSSを買った方がよいのか、判定するための一助としたいと思っています。よろしくお願い申し上げます。

  • Repeated measured ANOVA については,Rの統計解析関数Tipsに記述があるようですね。 -- 2006-12-11 (月) 16:10:25
  • Repeated-measures ANOVAではなく、repeated-measures ANCOVAに関する情報を探しています。 -- s.y.? 2006-12-11 (月) 16:25:07
  • わかっています。ですから,「Repeated measured ANOVA については」と書きました。あなたは,自分で google 等を検索したのでしょうか? -- 2006-12-11 (月) 16:41:15
  • はい、検索はしたのですが、わかりませんでした。Rに関する理解が足りないという以前に、ANCOVAに関する理論的な理解が足りないのかも知れませんが、何かしらご教示いただけましたら幸いです。 -- s.y.? 2006-12-11 (月) 16:54:31
  • 困った時の「R site search」( http://finzi.psych.upenn.edu/search.html )。 -- 2006-12-11 (月) 19:30:31
  • R site searchを使用してみましたが、やはりわかりませんでした。Rだけでやって行くには力不足なのかも知れません。 -- s.y.? 2006-12-12 (火) 00:40:31

loess平滑化について

MU? (2006-12-09 (土) 12:38:20)

皆様のお教えをいただけるようにお願い申し上げます.
体重成長の結果を平滑化しようとしています.
使用環境は次の通りです.
PCはWindows XPで、R version 2.2.1. パッケージは methods, stats, graphics, grDevices, utils, datasets, base が入っています.
エクセルから読み込んだマトリックスBWを元に、図を作成しました

uc.png

そして、RjpWiki アーカイブスの〔Rの基本パッケージ〕中の〔平滑化関数一覧〕の中から「2.2散布図とLoess平滑化の同時プロット」のまねを試みたのですが、predictionのところで「存在しません」とエラーが出てしまいます.
具体的には次の通りです.

> BW<- read.delim("clipboard")
> x<- BW[,1]
> y<- BW[,2]
> scatter.smooth(x,y,span=2/3,degree=1,family=c ("symmetric", "gaussian"),
      xlab=deparse(substitute(x)), ylab=deparse(substitute(y)),
      ylim=range(y,prediction$y), evaluation=50)
以下にエラーrange(y, prediction$y) : オブジェクト "prediction" は存在しません

この場合の、解決策を教えていただけないでしょうか.
勉強が足りない状態で、お伺いすることをお許しください.よろしくお願いします.

  • 2.1 の lowess 関数の戻り値を prediction に付置しているのでしょう。prediction がないといわれているのだから,それ以上の理由はないですね。 -- 2006-12-09 (土) 13:40:12
  • たぶんfamilyの引数でもエラーが出そうですね。 -- 2006-12-09 (土) 15:31:52
  • family の指定法は,尋常ではないですが,エラーにはなりません。それよりは,lowess の方の xy とか o の定義はどこにも書かれていないのでそちらで苦労するでしょう(適切に設定するだけでいいので,書かれていなくてもいいのですが)。 -- 2006-12-09 (土) 18:16:52
  • 取り急ぎ、お礼を申し上げます。まだ理解が進まず、解決できていないのですが、頑張ってみます。 -- MU? 2006-12-09 (土) 18:21:40

実行ファイルが文字化け?です。

-0.889? (2006-12-03 (日) 00:13:33)

実行ファイルをダウンロードし、出来たショートカットをクリックすると、変な文字ばかりです。アンインストールをし、再度ダウンロードしたのですが、結果は同じでした。Windowsがいけのいのですか?
それから、マクネマー検定のことですが、不一致の数だけが問題で、一致している数はいくらであろうと問題ないのですか?

  • 変な文字と言われても確認のしようがないと思います。実例を出してみてください。 -- まさひと? 2006-12-03 (日) 09:13:00
  • 実行ファイルが文字化け,ではなくて,できたショートカットをクリック(Rの起動)した後,出てくるウインドウの中の表示が文字化けということでしょう?アンインストールしたりする前に,インストールについての説明をもう一度読む方がよろしいかと思います。この質問,どこか別のところ(あるいはここ)にしてましたよね。違うコメントも付いていたと思うんだけど,回答者自ら削除したのかな? -- 2006-12-03 (日) 14:38:26
  • この質問の下3つ目ですね。やたら,スレッドを増やさないようにしましょう。『トップページの「Rのインストール」は読んだのか』についての,お返事もなかったようですが。 -- 2006-12-03 (日) 14:43:35
  • 一度に関連の無い質問を複数するのも行儀が悪いですね。 -- 2006-12-03 (日) 17:46:54
  • ご迷惑をおかけして申し訳ありません。 フリーのソフトウエアをダウンロードして、使わせていただくのも初めてですし、このような書き込みをするのも始めての経験です。書き込みのルールがあるなんて知りませんでした。申しあけありませんでした。
    トップページの「Rのインストール」からR-2.3.1をダウンロードして、ショーカットから起動させてみたのですが、文字化け?です。「THe R Project」には「exe」ファイルは見つけることは出来ませんでした。そして、教えていただいた「tukuba」のサイトからR-2.4.0をダウンロードして、「R]から起動させたウインドウの中の表示が文字化けでした。
    それをお見せしようとコピーして貼り付けましたら「Rはフリーソフトウェアであり、「完全に無保証」です。」日本語に変わってしまいました。訳が分かりません。起動させたウインドウの表示は文字化けのままです。
    長くなりすみません。 -- マキヤン? &new{2006-12-03 (日) 22:08:40}
  • 勝手ながら整形させていただきました.ご質問の件ですが,以下のサイトの Windows版 (10) の操作をしてください.おそらく症状が治ります.

    http://cwoweb2.bai.ne.jp/~jgb11101/files/cart/cart.html

    この分だと,R を使い始めた後もいろいろと分からないことが出てくるでしょうから,データ解析環境R とか R の基礎とプログラミング技法 あたりで R の基礎を学んで下さい.それでも分からないことがあれば質問してください.教えて君 にならないように. -- 通りすがり? 2006-12-03 (日) 22:36:52
  • いろいろなご教示をありがとうございました。がんばってみます。-- マキヤン? &new{2006-12-05 (火) 00:05:12}

Linux版Rにパッケージをなんでもかんでも追加する方法

Fedore Core5使い? (2006-12-02 (土) 16:44:50)

Linux(Fedore Core5)に、「追加パッケージをなんでもかんでも追加する」を参考にして、パッケージをなんでもかんでも追加しようとしたのですが、

# 1. CRANから一覧を入手
packs<-CRAN.packages(contriburl=contrib.url("http://cran.md.tsukuba.ac.jp/"))
# 2. 一覧をもとにダウンロードとインストール
install.packages(packs[,1],contriburl=contrib.url("http://cran.md.tsukuba.ac.jp/"))

エラーが出て上手く行きません。。。
(「contriburlは廃止予定です」というエラーは無視してよいと思うのですが)
対処法を教えていただければうれしいのですが。。。

  • 例えば、どんなエラーになるのでしょう。一般には、一つのパッケージをインストールしようとしてもエラーになることはあります(コンパイル環境やら、依存性やらで)。 -- 2006-12-02 (土) 19:46:49
  • 御指導通り、「一つのパッケージをインストールしようとしてもエラー」になりました。。。
    > install.packages("Rcmdr")
    --- Please select a CRAN mirror for use in this session ---
    Loading Tcl/Tk interface ... 以下にエラーdyn.load(x, as.logical(local),  
    as.logical(now)) :
           共有ライブラリ '/usr/lib/R/library/tcltk/libs/tcltk.so' を読み込めません
    libtk8.4.so: 共有オブジェクトファイルを開けません: そのようなファイルやディレ クトリはありません
    エラー:.onLoad は 'tcltk' のための 'loadNamespace' に失敗しました
    ひとつずつインストールするのが良いようですね。御指導ありがとうございました。-- 2006-12-02 (土) 19:46:49
  • tkのライブラリが無いといってますよ. yum install tk をすれば解決しませんか.? -- なかま 2006-12-04 (月) 10:28:58

RODBCとEXCELファイル・ACCESSファイルの読み込みについて

-0.889をハンドルネームにするのは・・・? (2006-12-02 (土) 10:33:42)

EXCELファイルやACCESSファイルを読み込む際の話ですが、Windows版のRですと、パッケージ「RODBC」のodbcConnectExcel?()関数などで読み込むことが出来ますが、MacOSX版Rですと、パッケージ「RODBC」の中にはodbcConnect()関数しかなく、EXCELファイルやACCESSファイルを読み込むことが出来ません・・・。(ちなみに、当方のMacにOffice Macは入っております)CSVファイルなどに落としてread.csv()関数で、という手もあるのでしょうが、これらのファイル形式のデータをodbcConnect()関数などで読み込む方法をどなたかご存知でしたらご教授いただけますでしょうか??
あと、これは興味本位の質問なのですが、Linux版RでEXCELファイルやACCESSファイルを読み込む場合はどうすればよいのでしょう??(Open Office に ODBC ドライバが入っていたりするのでしょうか??)
ご指導よろしくお願い致します。

  • レスがついていないようなので,まず最近のMacOSXに日本語が通るODBCマネージャ(windowsでいうところののデータソース(ODBC))はあるのでしょうか?あるのなら、ドライバはiODBCあたりにありませんか?LinuxについてもiODBCあたりを探してみては? -- okinawa 2006-12-12 (火) 10:28:08
  • Windows版で読み込めるんですか?僕もRcmdr経由でRODBC試してみたんですが、Excelとのコネクトは失敗しています。コードをエディタにコピペして実行してはみたんですが、やはり上手く動きませんでした。相変わらずcsvでExcelとはやり取りせざるを得ません。 -- R初心者? 2006-12-29 (金) 06:35:14
  • Mac,Linuxは実機が手元にないのでどこまで可能か不明ですが、Windowsでは設定が正しければ可能です。SQL自体もクエリー対象のアプリによっては、テーブルの指定が.だったり_だったり、セパレータも"だったり'だったりしますので、コピペですぐ動くことはまれです。私の知りうる限りMacでは、ODBCドライバは別途購入する必要があったり、日本語のSQL抽出がうまく通らなかったりしていたと思いますが・・・。 -- okinawa 2006-12-30 (土) 11:48:43
  • R2.4.1&RODBCver1.1-7,WinXPsp2でc:?test.xls作成し、コンソールから下記のように入力したらちゃんと読み込めましたよ。 -- okinawa 2006-12-30 (土) 17:32:43
    > library(RODBC)
    > channel <- odbcConnectExcel("C:/test.xls")
    > channel
    RODB Connection 7
    Details:
     case=nochange
     DBQ=C:?test.xls
     DefaultDir=C:?
     Driver={Microsoft Excel Driver (*.xls)}
     DriverId=790
     MaxBufferSize=2048
     PageTimeout=5

実行ファイルのダウンロードについて(ダウンロードファイルが何か変です)

-0.889? (2006-11-30 (木) 21:57:15)

Rのダウンロードを教えてください。The R projectに行っても実行ファイルがみつかりません。

  • http://cran.md.tsukuba.ac.jp/bin/windows/base/ へ行きましょう。たぶん Windows ユーザですね? -- 2006-11-30 (木) 23:57:15
  • 実行ファイルのダウンロードについて教えていただきありがとうございました。ダウンロードできたのですが何か変です。初心者用のQ&Aに書きますので、また教えてください。何か、相当大変みたいです。 --マキヤン? (2006-12-02 (土) 01:02:02)
  • 真の初心者です。R-2.4.0-win32をダウンロードして、ショートカット「R」から起動させたのですが、初期画面が得体の知れない文字ばかりで、何か変です。インストールする際にコンポーネントを確認したのですが、「version for east asian languages」は入っていないようでした。Rが使えるようになるまでに、大変な苦労が待っていそうです。手取り足取り教えてください。-- -0.889? 2006-12-02 (土) 01:13:26
  • トップページの「Rのインストール」は読んだんですか? 手取り足取りなんて... 揚げ足なら取ってくれるでしょうね。  -- 偽の初心者? 2006-12-02 (土) 05:52:44
  • Macintoshだともう少し簡単で,少なくともクリックして文字化け画面が出るなんてことはない。今すぐ電気屋さんへいって,Macintoshを買いましょう。 -- 2006-12-02 (土) 08:36:36

同じ大きさの複数のプロットを重ねて描画

初心者です? (2006-11-22 (水) 23:26:49)

いつもお世話になっています.
共通の横軸を持つプロットを複数積み重ねて描画したいのですが,

theta<-seq(-pi,pi,by=0.1)
par(mfrow=c(3,1),mar=c(0,4,0,0),oma=c(5,4,4,1))
plot(theta,sin(theta),type="l",xlab="",ylab=expression(sin(theta)),xaxt="n")
axis(1,tick=T,labels=F)
plot(theta,cos(theta),type="l",xlab="",ylab=expression(cos(theta)),xaxt="n")
axis(1,tick=T,labels=F)
plot(theta,tan(theta),type="l",xlab=expression(theta),ylab=expression(tan(theta)))
mplot1.png

とすると,一番下のプロットに付けたxlabが消えてしまいます.そこで

theta<-seq(-pi,pi,by=0.1)
par(mfrow=c(3,1),mar=c(0,4,0,0),oma=c(5,4,4,1))
plot(theta,sin(theta),type="l",xlab="",ylab=expression(sin(theta)),xaxt="n")
axis(1,tick=T,labels=F)
plot(theta,cos(theta),type="l",xlab="",ylab=expression(cos(theta)),xaxt="n")
axis(1,tick=T,labels=F)
par(mar=c(5,4,0,0) plot(theta,tan(theta),type="l",xlab=expression(theta),ylab=expression(tan(theta)))
mplot2.png

としてみましたが,これでは最後にプロットしたグラフの大きさが上の二つとは変わってしまいます.
どなたかうまい方法をご存知であれば教えていただけないでしょうか
環境 : R version 2.5.0 Under development (unstable) (2006-11-20 r39944)
i386-pc-mingw32

  • 最初のプログラムの方で,2行目のmarの設定を変えればいいだけではないでしょうか?
    mar=c(0,4,0,0)じゃ,xlab を描くだけの余白がないんですよ。3,4番目の0も必要とあればもう少し大きくしないといけない場合もありますよ。 -- 2006-11-22 (水) 23:33:53
    par(mfrow=c(3,1),mar=c(4,4,0,0),oma=c(5,4,4,1))
  • 回答ありがとうございます.説明不足だったのですが,三つのプロットの横軸は共通なので,上二つのプロットの横軸にはtickだけ入れ,なおかつそれらの下にくるプロットにぴったりと重ねたいんです.そうするとmar=c(0,4,0,0)としか書けないような気がするんですが・・・
    二つ目に書いたプログラムでは三つ目のプロットのxlabの入るスペースをpar(mar=c(5,4,0,0))として確保したんですが,その分プロットの高さが低くなってしまったという次第です -- 初心者です? 2006-11-22 (水) 23:42:10
  • あ,そうか。それじゃ,3つのグラフの間に隙間ができるというわけなんだね? -- 2006-11-22 (水) 23:47:35
  • そうなんです.残された方法はmtextで下の余白にxlabを書き込む,とかになるんでしょうか -- 初心者です? 2006-11-23 (木) 00:04:45
  • でもね。縦軸の目盛りの数値が重なってみっともないとか言うこともあるわけで。
    文字の大きさのバランスも変だったので,私の環境で調整してみると以下のようになった。このような感じかな? -- 2006-11-23 (木) 00:02:30
    theta<-seq(-pi,pi,by=0.1)
    par(mfrow=c(3,1),mar=c(0,4,0,0),oma=c(4,1,1,1),xpd=NA,cex=3)
    plot(theta,sin(theta),type="l",xlab="",ylab=expression(sin(theta)),xaxt="n")
    axis(1,tick=T,labels=F)
    plot(theta,cos(theta),type="l",xlab="",ylab=expression(cos(theta)),xaxt="n")
    axis(1,tick=T,labels=F)
    plot(theta,tan(theta),type="l",xlab="",ylab=expression(tan(theta)))
    text(0, -110, expression(theta))
    multi.png
  • 縦軸(野ラベル,目盛り)は左右交互に描くとよいかな。(上の例だと,cos の目盛りは右に描く)。((それも変か)) -- 2006-11-23 (木) 00:13:52
  • ありがとうございました.参考にしてやってみました.
    theta<-seq(-pi,pi,by=0.1)
    par(mfrow=c(3,1),mar=c(0,4,0,0),oma=c(5,4,4,1),las=1)
    plot(theta,sin(theta),ylim=c(min(sin(theta))*1.5,max(sin(theta))*1.5),type="l",xlab="",
        ylab=expression(sin(theta)),xaxt="n",yaxt="n")
    axis(1,tick=T,labels=F)
    axis(2,tick=T,labels=T,at=pretty(sin(theta)))
    plot(theta,cos(theta),ylim=c(min(cos(theta))*1.5,max(cos(theta))*1.5),type="l",xlab="",
        ylab=expression(cos(theta)),xaxt="n",yaxt="n")
    axis(1,tick=T,labels=F)
    axis(2,tick=T,labels=T,at=pretty(cos(theta)))
    plot(theta,tan(theta),ylim=c(min(tan(theta))*1.5,max(tan(theta))*1.5),type="l",xlab="",
        ylab=expression(tan(theta)),xaxt="n",yaxt="n")
    axis(1,tick=T,labels=T)
    axis(2,tick=T,labels=T,at=pretty(tan(theta)))
    mtext(expression(theta),side=1,line=3)
    mplot3.png
    あとは仰られるように文字サイズの微調整も必要そうですね. -- 2006-11-23 (木) 00:48:17
  • 水色の背景で表示される部分は自動改行されず、ブラウザーで表示するとスクロールが必要になります。適宜強制改行(複数行に分ける)をお願いします。(直しておきました。) -- 2006-12-02 (土) 06:14:52

Rmapのインストールについて

まめ? (2006-11-22 (水) 14:54:07)

Rmapを用い、地図を作成したいのですが、手順どおり作業を行ってもうまくいきません。
ご存知の方教えてください。

使用環境は以下です。
PC:Windows XP

> sessionInfo()
Version 2.3.1 (2006-06-01)
i386-pc-mingw32

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

手順は以下の通り行いました。
1、以下のファイルを入手
・shapefile.dill
・proj.dill
・Rmap_1.1.0
2、shapefill.dillおよびproj.dillをRインストール先のbinファイルに保存
3、RGuiを起動し、『パッケージ』から『ローカルにあるzipファイルからのパッケージのインストール...』を選択し、Rmap_1.1.0.zipを読み込む

以上を行うと以下のように示されます。

> utils:::menuInstallLocal()
updating HTML package descriptions

そしてRmapパッケージをロードすると、エラーが返されます。

> library(rmap)
以下にエラーlibrary(rmap) : 'rmap' という名前のパッケージはありません

そこで、再び『パッケージ』に入り、『パッケージの読み込み...』のパッケージ一覧からRmapを選択すると以下のようなメッセージが返されます。

> local({pkg <- select.list(sort(.packages(all.available = TRUE)))
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
以下にエラーlibrary(pkg, character.only = TRUE) : 
       'Rmap' は有効なパッケージではありません。バージョン < 2.0.0 でインストールされたかも知れません

以上が私の行った作業とその結果です。 初心者のため、何が問題点なのか分かりません。お手数ですが、どなたかお教えいただけませんでしょうか?

  • RmapってR1.9.Xまでしか使えなかったような記憶があります。maptoolsを使うか、R1.9.Xでやってみてはいかがでしょうか? -- okinawa 2006-11-23 (木) 11:43:49
  • ありがとうございます。試してみます。 -- まめ? 2006-11-29 (水) 17:56:14

Rの予約語

ちょもらんま? (2006-11-14 (火) 13:11:51)

Rで使われているRのシステム自体が用いる変数名=予約語(data等々)の全リストが知りたいのですがどこかにありますでしょうか
またもし自分で作ったRスクリプト中でRの予約語と重複するような変数名を用いた場合、該当スクリプトの動作がおかしくなるというようなことはありますでしょうか?

  • http://cran.r-project.org/doc/manuals/R-lang.html#Reserved-words
    FALSE, TRUE, break, else, for, function, if, in, next, repeat, while, NULL, NA, Inf, NaN, ..., ..1, ..2 など
    これらは,変数名などとして使うことはできません(for <- 5 などとすると,文法エラーになる) -- 2006-11-14 (火) 16:20:10
  • これらはエラーになるから実害は無いでしょうが、たとえば sin <- function(x) cos(x) などとするとエラーにはならず困った事になる。 -- 2006-11-14 (火) 21:34:36
  • それは予約語でないので仕方ない。そういうことをする人が悪い。よく,Sと比較されるのは,SではT/FはTRUE/FALSEと同じく予約語だが,Rでは再定義可能なので,T<-0とかF<-345等とされていると困ったことになる。
    ちなみに,本当に sin を cos にしたいなら,sin <- cos だけで十分。よく使う関数を1文字で表したいなんていう向きには,r <- read.table とか w <- write.table とか好き勝手に再定義できます。 -- 2006-11-14 (火) 21:38:20
  • みなさま、ご教授ありがとうございました。疑問がすべて解決いたしました。大変感謝しておりますm(_ _)m -- ちょもらんま? 2006-11-15 (水) 10:38:21
  • おぉ、仲間を発見!私もTという名前の変数をうっかり作ってしまい、そのままイメージを保存してしまったので、以降のコマンドオプション指定で、はまりました。 -- 2006-11-15 (水) 12:20:09
  • 予約語使用即エラーにならない例も(将来は知らない)。-- 2007-01-04 (木) 13:41:29
    > ..1 <- 3
    > ..1
    [1] 3 

マトリクスからNAを除去する方法

ふみふみ? (2006-11-14 (火) 10:53:52)

数値(整数4桁、小数3桁)からなるデータフレーム、またはマトリックス中にデータが無いことを示す「NA」が行中、または列中の特定箇所に含まれてしまっている様な場合、その「NA」を含む、行、列を除き、「NA」を含まない、純粋な数値からなるデータフレーム、またはマトリックスを抽出したいと考えております。そのような事を可能にするような関数等存在しますでしょうか。本Wiki,色々なRのサイト、マニュアル等見ましたが適した関数が無く、質問させていただきました。
どうぞよろしくお願いいたします。

  • 捜すより作る方が早い -- 2006-11-14 (火) 11:16:04
    > # 例を作る
    > x <- matrix(1:30, 5, 6)
    > x[1,2] <- x[3,4] <- x[4, 6] <- NA
    > x
         [,1] [,2] [,3] [,4] [,5] [,6]
    [1,]    1   NA   11   16   21   26
    [2,]    2    7   12   17   22   27
    [3,]    3    8   13   NA   23   28
    [4,]    4    9   14   19   24   NA
    [5,]    5   10   15   20   25   30
    > # 適用   わずかにこれだけ
    > y <- x[!is.na(rowSums(x)), !is.na(colSums(x))]
    > y
         [,1] [,2] [,3]
    [1,]    2   12   22
    [2,]    5   15   25
  • ありがとうございます。本当に迅速な回答痛み入ります。自分のとてつもない勉強不足さと馬鹿さ加減にほとほと嫌気がさしそうですが、これからもがんばって勉強していきます。ご迷惑をおかけすることもあるとおもいますが、どうぞよろしくお願いいたしますm(_ _)m -- ふみふみ? 2006-11-14 (火) 11:51:48
  • x[complete.cases(x),complete.cases(t(x))] というのもあり -- 2006-11-14 (火) 21:58:02
    > x[complete.cases(x),complete.cases(t(x))]
         [,1] [,2] [,3]
    [1,]    2   12   22
    [2,]    5   15   25
    それでもって,びっくりしたことは,後者の方が効率がよいということ。
    > x <- matrix(1:30, 5, 6)
    > x[1,2] <- x[3,4] <- x[4, 6] <- NA
    > system.time(for(i in 1:100000) y <- x[!is.na(rowSums(x)), !is.na(colSums(x))])
               user          system           total   user.children system.children 
             24.156           0.480          27.397           0.000           0.000 
    > system.time(for(i in 1:100000) z <- x[complete.cases(x),complete.cases(t(x))])
               user          system           total   user.children system.children 
              4.914           0.087           5.540           0.000           0.000 
    ていうか,行列のサイズにも依存するかもね。こういうのは,判断が難しい。
    > x <- matrix(rnorm(1000000), 1000, 1000)
    > x[sample(1000,50),sample(1000,50)] <- NA
    > system.time(for(i in 1:100) y <- x[!is.na(rowSums(x)), !is.na(colSums(x))])
               user          system           total   user.children system.children 
             27.518           2.444          32.602           0.000           0.000 
    > system.time(for(i in 1:100) z <- x[complete.cases(x),complete.cases(t(x))])
               user          system           total   user.children system.children 
             25.149           4.482          32.726           0.000           0.000 

table()で集計した時にカウントがゼロの項目も0と出したい

<ふ>? (2006-11-09 (木) 22:24:40)

みなさま、

データの単純な集計にtable()をつかってますが、以下のような場合に困っています。
たとえば、xを、

> x <- c(1,1,1,2,2,4,4,5)

として、

> table(x)

を実行すると、

x
1 2 4 5 
3 2 2 1 

となりまして、データがある、1,2,4,5は、でてきますが、データがない、3は、出力にあらわれません。
しかし、グラフに表示するにあたって、項目3が、ゼロであることを表現したいのです。
つまり、

x
1 2 3 4 5
3 2 0 2 1

というような状態をつくりたいのですが、なにか方法はありませんでしょうか。
たとえば、項目の要素を、ベクトルで与えて(ここでは、c(1:5))、それに対して「集計する」というようなことができればいいのですが。
アドバイスよろしくお願いいたします。

  • table(factor(x, levels=1:5)) -- コラロド? 2006-11-09 (木) 23:50:35
  • ありがとうございます。既出なんですね、やはり...。いろいろさがしたのですが、見つからなくて、お聞きしてしまいました。たすかりました。それから、プログラム部分では「半角スペース」ですね。修正、お手数おかけしました。次回から、やれると思います。 -- <ふ>? 2006-11-10 (金) 00:07:30
  • 知っているといつか役に立つ(?)関数達 の (74) を参照。 -- 2006-11-10 (金) 06:41:20
  • 「知っていると..」(74)も拝見いたしました。ありがとうございます。 -- <ふ>? 2006-11-10 (金) 09:02:46

ファイルサーバーのファイルを取得する方法

Akira? (2006-11-08 (水) 17:59:29)

WindowsXPからLinux(ubuntu)に乗り換えようとしています。R-2.4.0patchです。
ファイルサーバにあるデータを読み込みたいのですが、サーバにはUserName?とPasswordが設定してあります。
WinXPではGUIでサーバにアクセスしておけば、R上でも読み込むことができます。
Linuxではうまくできません。Linux上でsetwd("//xxx.xxx.xx.xxx/abc")のようにアクセスする方法はありますでしょうか?
WinXPでは次のようにしています。

  1. //xxx.xxx.xx.xxx/abcにデータがあるとします
  2. GUIでサーバにアクセス(UserName?とPasswordを入力)
  3. Rを起動
  4. setwd("//xxx.xxx.xx.xxx/abc")で移動可能
  • mountしとけばよろしいのでは? -- 2006-11-08 (水) 18:21:12
  • 解決しました。ありがとうございます。ubuntuのDesktopにmountしてもダメで、Rを起動するterminalでmountしなければならなかったのですね。 -- Akira? 2006-11-08 (水) 21:07:19
    >sudo smbmount //xxx.xxx.xx.xxx/abc /mnt/abc -o username=abc

jpeg()、bmp()などのファイル番号の取得方法は?

<ふ>? (2006-11-07 (火) 17:11:15)

みなさま、

グラフを出力するにあたり
jpeg()
のように、ファイル名を指定しないでおくと、Rplot%03d.jpg のような書式で順番にファイル名をインクリメントして書き出してくれます。

で、この%03dにわたされている数値を参照する方法はありませんでしょうか。

現在やっているのは、
1)xtable()である項目をHTMLの表にする。
2)その項目のグラフを描画しファイルに書き出す。

これを、一本の report.html というようなファイルで束ねることです。
1)は、そのまま、sink(repot.html)で吐き出してくれるのですが、グラフのファイルを、<IMG SRC="....">で指定したいと思っております。次に書き出すファイル名(%03dに渡す番号)がわかれば、ファイル名を書き出せます。

アドバイスよろしくお願いいたします。

  • ファイル名に含まれる連番を自分で管理すればよろしいのではないでしょうか?
    以下のようなプログラムの断片では,ファイル名はFig201.jpeg, Fig202.jpeg, Fig203.jpeg のようになりますね。
    蛇足ですが,グラフに jpeg はあまりお勧めではないかも知れません。 -- 2006-11-07 (火) 18:39:48
    for (i in 1:3) {
    	file.name <- sprintf("Fig%03i.jpeg", i+200)
    	jpeg(file.name, width=400, height=300)
    	hist(rnorm(1000))
    	dev.off()
    }
  • 実は、一連のデータを集計しながらグラフを描いたり、削除したりします。その意味では、何番目かをスクリプトの側ではコントロールしたくないのです。で、実際に書き出しているデバイスから番号が得られれば、確実に追えると考えているのです。 -- <ふ>? 2006-11-07 (火) 19:30:37
  • あなたがその情報にアクセスするタイミングが今ひとつ明確でないのですが,jpeg ファイルに書き出した直後にそのファイル名を知りたいなら,以下のような関数を作っておく。知りたいときに関数を呼ぶと,ファイル名が返される。。。。 -- 2006-11-07 (火) 21:13:23
    fname <- function()
    {
    	a <- list.files()
    	sprintf("Rplot%03i.jpg", max(as.numeric(substr(a[grep("Rplot", a)], 6, 8))))
    }
    利用例は以下のごとし。ただ,私は,pdf 関数でやってみた。
    > fname <- function()
    + {
    + 	a <- list.files()
    + 	sprintf("Rplot%03i.pdf", max(as.numeric(substr(a[grep("Rplot", a)], 6, 8))))
    + }
    > pdf(onefile=FALSE)
    > hist(rnorm(1000))
    > hist(runif(1000))
    > fname()
    [1] "Rplot002.pdf"
    > plot(1:100, (1:100)^2)
    > plot(1:100, (1:100)^3)
    > fname()
    [1] "Rplot004.pdf"
    > plot(rnorm(100), rnorm(100))
    > fname()
    [1] "Rplot005.pdf"
  • なるほど、作業ディレクトリのファイル名から引き出してくるわけですね。関数までつくっていただきありがとうございます。先ほどスクリプトに組み込み、目的を達することができました。ありがとうございます。 -- <ふ>? 2006-11-07 (火) 22:49:45
  • 前述の方法よりかなりわかりにくいが、次の一行コマンドでも最新のタイムスタンプをもつファイル名を取り出せます (?)。参考まで。解説すると、glob2rx で Unix 風のファイル名のワイルドカード指定を正規表現に変換し、list.files でそれに合致するファイル名の情報データフレームを得、その生成時刻列を "mtime" で取り出し、一番最新の行番号を which.max で求め、さらに対応する行ラベル名を得ています、ややこしい(笑) -- 2006-11-07 (火) 22:53:03
    > row.names(a <-  file.info(list.files(pattern=glob2rx("Rplot[0-9]*.jpg")))["mtime"])[which.max(a[,1])]
    [1] "Rplot002.jpg"
    ちなみに得られる正規表現は以下の通りですから、これを直接 pattern に指定しても可。
    > glob2rx("Rplot[0-9]*.jpg")
    [1] "^Rplot[0-9].*??.jpg$"

read.fortran, read.fwf とマルチバイト

青木繁伸 (2006-11-06 (月) 18:26:00)

自分ではやらないことも,他人がやってしまうとその後始末が大変です。

私なら数値にコード化してデータ入力するのだが,日本語のままデータを入力してくれた学生がおりまして,read.fortran か read.fwf で読ませればよいやとタカをくくっていたのだが,これらの関数は,バイト数で読むのではなく文字数で読むということが判明。AWKかperlでスクリプトを組んで,CSVファイルにはき出すと何のことはないのだが,その学生はプログラムできないので,エクセルの「データ区切り」を使ってCSVファイルを作らせることになってしまった。
入力欄が5バイト分という半端な状態だったのであるが,これを偶数バイト(6バイト)にして全角3文字に満たないときは残りを全角空白で入力)なるようにエディタでちょいちょいといじって,読めるようにしてしまった。
read.fwfのソースをいじるのも気が進まない。

  • 一次的にC localeで読み込んでしまえばいいかと
    lc<-Sys.getlocale("LC_CTYPE")
    Sys.setlocale("LC_CTYPE","C")
    df<-read.fwf("hoge.dat",width=c(バイト,バイト,バイト))
    Sys.setlocale("LC_CTYPE",lc)
    . -- なかま 2006-11-06 (月) 20:23:28
  • なるほど。そんな手があるんですか。勉強になります。 -- 青木繁伸 2006-11-06 (月) 22:52:48
  • やってみました。残念ながら,読みこむことはできませんでした。最終行に関するwarningについては,cr/lf が付いているので問題ないはずです。 -- 青木繁伸 2006-11-07 (火) 10:16:54
    > lc<-Sys.getlocale("LC_CTYPE")
    > Sys.setlocale("LC_CTYPE","C")
    [1] "C"
    > con <- file(description="ID1001-2.txt", open="r", encoding="shift-jis")
    > df <- read.fwf(con, width=c(4,rep(1,13),6,2,1,1,6,2,rep(1,59), 2, 2, 1, 4), header=TRUE)
    Warning messages:
    1: 入力コネクション 'ID1001-2.txt'  に不正な入力がありました 
    2: 'ID1001-2.txt' に関する readLines で不完全な最終行が見つかりました 
    >  Sys.setlocale("LC_CTYPE",lc)
    [1] "ja_JP.UTF-8"
  • CP932環境でだけです. charwidth==byteでは無い環境ではread.fwfは救済不能です. -- なかま 2006-11-07 (火) 12:02:22
  • マックな人はあきらめてと言われたのですが,必要な人もいるかも知れないので,以下のようなものを提示しておきます。
    Mac OS X で,gawk (awk でもよいけど)がインストールされていることを仮定しています。
    perl を使う方がいいのだろうけど,私は perl をしゃべれない。 -- 青木繁伸 2006-11-08 (水) 13:46:53
    # fn.in 読みこみたいデータファイル名
    # width 各変数に読みこむ欄の桁数(バイト数!)
    # header 1行目に変数名があるかないか
    get.fwf <- function(fn.in, width, header=FALSE)
    {
    	nf <- length(width)
    	pos <- cumsum(width)+1
    	pos <- c(1, pos[-nf])
    	fn.out <- tempfile()
    	fn.awk <- tempfile()
    	con <- file(fn.awk, open="w", encoding="euc-jp")
    	cat("BEGIN {?n", file=con)
    	for (i in 1:nf) {
    		cat(sprintf("pos[%i]=%i; width[%i]=%i?n", i, pos[i], i,width[i]), file=con)
    	}
    	cat("}?n", file=con)
    	if (header) {
    		cat("FNR == 1?n", file=con)
    		cat("FNR > 1", file=con)
    	}
    	cat(sprintf("{ for (i = 1; i < %i; i++) {?n", nf), file=con)
    	cat("str = allblank(substr($0, pos[i], width[i]))?n", file=con)
    	cat("printf ?"%s??t?", str}?n", file=con)
    	cat("str = allblank(substr($0, pos[i], width[i]))?n", file=con)
    	cat("print str }?n", file=con)
    	cat("function allblank(str) {?n", file=con)
    	cat("for (j=1; j <= length(str); j++)?n", file=con)
    	cat("if (substr(str, j, 1) != ?" ?") return str?n", file=con)
    	cat("return ?"NA?"", file=con)
    	cat("}?n", file=con)
    	close(con)
    	# awk しかない人は ln -s awk gawk するか,下の gawk を awk に書き換え
    	cmd <- paste("gawk -f", fn.awk, fn.in, ">", fn.out)
    	print(cmd)
    	system(cmd)
    	df <- read.table(file(fn.out, open="r", encoding="euc-jp"), header=header)
    	unlink(fn.awk)
    	unlink(fn.out)
    	return(df)
    }
    # 利用例
    # df <- get.fwf("test.dat", width=c(4, rep(1, 13), 5, 2, 1, 4), header=TRUE)

pieをclockwise=TRUEでつかいたい

<ふ>? (2006-11-04 (土) 23:17:00)

みなさま、
pieを使って円グラフを書くときに、今までの習慣で、0時から時計回りに表示させています。そのために、clockwise=TRUEと引数を指定して記述しているのですが、だんだん面倒になりまして、(Rってオブジェクト指向言語なんだから)
pie2() = pie(closkwise=TRUE)

というように、新しいインスタンスをつくれないだろうかと思っています。
ですが、上記の形ではNGでした。

なにかやりようがありそうですが、どなたかご教示いただけませんでしょうか。

  • 例えば次のようにすればよさそうですが、もっと簡単な方法があってもよさそうな。 -- 2006-11-05 (日) 00:30:10
    > pie2 <- function(x, labels = names(x), edges = 200, radius = 0.8,
                     init.angle = 90, density = NULL, angle = 45,
                     col = NULL, border = NULL, lty = NULL, main = NULL, ...)
            {
              pie(x, labels, edges, radius, clockwise = TRUE, init.angle,
              density, angle, col, border, lty, main, ...)
            }
  • pie2 <- function(x, ...) pie(x, clockwise=TRUE, ...) だけでいいんじゃ? -- 2006-11-05 (日) 07:47:25
  • pie のソースを書き換えて,pie に定義し直すというのも。 -- 2006-11-05 (日) 07:48:42
  • ありがとうござます。個人的にはソース書き換えでやっていたのですが、学生達が使う環境全部に手を加えるわけにはいかない、という事情がありました。で、07:4:25のやり方で、できました。ありがとうございます。パラメータを非defaultの固定値で使いたい時は、この方法で定義すればいいわけですね。 -- <ふ>? 2006-11-05 (日) 10:08:34
  • なるほど。pie2 <- function(...) pie(...,clockwise=TRUE) だけでもよさそうですね。 しかし、いずれも pie2(x,clockwise=FALSE) とするとエラーになる。最初の例を次のように変更すれば大丈夫。 -- 2006-11-05 (日) 11:02:26
    > pie2 <- function(x, labels = names(x), edges = 200, radius = 0.8,
                       clockwise = TRUE, init.angle = if(clockwise) 90 else 0,
                       density = NULL, angle = 45, col = NULL, border = NULL,
                       lty = NULL, main = NULL, ...)
              { pie(x, labels, edges, radius, clockwise, init.angle,
                    density, angle, col, border, lty, main, ...)
              }
  • ああだ,こうだあるようですが,このような応用はオブジェクト指向とか言うことではなく,シュガーコートのレベルではないかと思います?(念のため,最初のコメンテーター 2006-11-05 (日) 07:47:25 です) -- 2006-11-05 (日) 21:06:46
    ああだ,こうだいいたいときには
    pie2 <- function(x, clockwise=TRUE, ...) pie(x, clockwise=clockwise, ...)
    とすれば完璧なんじゃないでしょうか。こうしておけば,clockwise=FALSE が指定されても,問題ないでしょう(最初に提示しなかったのは冗長だという理由からだけです)
    それと,pie2 の引数名にclockwiseを使いましたが,これは原理をわかりやすくするという意味ではcw=TRUE 等と書いて,関数定義部分ではclockwise=cw などとするのがおぎょうぎのいい定義の仕方でしょう。
    そのようにしなかったのは,利用者の便宜のためです。
    てうか,そんなような別の関数を定義するとか,そのような使い方を覚えておくとか言う面倒なことをするより,pie(なんだらかんだら, clockwize=TRUE) とやってね。。。というのが,簡単明瞭じゃないですか。そのような引数が存在するんんだなあ。ということも,わかるし。使い分けも。分かるし。シュガーコートというのは,個人的なもので,みんなに使ってもらうために作るのは善し悪しですよ。(誰にとってのシュガーコートかということ。人によってはよけいなお世話だし)

PLSRについて

kazu? (2006-11-02 (木) 15:21:46)

Rを用いて、PLSRを行いたいと思っています。
Rとpls.pcrのversionはそれぞれ、2.4.0、1.2-1です。 下記のページを参考にしています。
http://cse.naro.affrc.go.jp/iwatah/toukei/pls/2004/pls_intro04.htm
pls.pcrパッケージに含まれたサンプルデータではpls解析を行えることを確認しましたが、「手持ちのデータでPLS回帰分析」の部分がうまくいきません。
wine.txtを

> wine <- read.table("wine.txt")

で読み込み

> wine

で読み込めたことを確認しました。
次に

> Y <- scale(wine[,1:3])
> X <- scale(wine[,4:7])

でデータを切り分けました。
plsrを行うために、
http://cran.r-project.org/doc/packages/pls.pdf
を参考に

> wine.pls <- plsr(Y ~ X, 6, validation="CV")

と入力しましたが

以下にエラーmvr(Y ~ X, 6, validation = "CV", method = "kernelpls") : 
       Invalid number of components, ncomp

と表示されうまくいきません。ncompの値をいじってはみたのですが同様のエラーが表示されます。
どのような原因が考えられるでしょうか?

  • ncomp をいくつにしましたか?wineデータがどんなのかわからないので,正規乱数で適当にこしらえた場合には,ncomp が 1 〜 4 なら,エラーは表示されず解らしきものが出ますが。。。(1のときには別のエラーが出るがぁ)-- 2006-11-02 (木) 16:11:15
  • コメントありがとうございます。具体的にはこのようなデータです。
    	Hedonic	Goes_with_meat	Goes_with_dessert	Price	Suger	Alcohol	Acidity
    1	14		7		8	7	7	13	7
    2	10		7		6	4	3	14	7
    3	8		5		5	10	5	12	5
    4	2		4		7	16	7	11	3
    5	6		2		4	13	3	10	3 
    ncompの値は1〜10位の間で振ってみました。-- kazu? 2006-11-02 (木) 16:29:18
  • 元のデータが,5行のデータであるとは思いません。(示してくれたページにあるのは,最初の5行でしょう?)そうであるとしても,ncomp を1にしたら,少なくともエラーは表示されず,解?らしきものは表示されるが? -- 2006-11-02 (木) 16:57:30
  • 元データは下記アドレスからダウンロードしたものを用いましたが、5行のデータのようです。
    http://cse.naro.affrc.go.jp/iwatah/toukei/pls/
    ncompを1にすると解らしきものは確かにでますね。確認できました。ありがとうございます。
    validationまでの結果は上記のページのものと一致するようです。ですが、trainingはうまくいかないようですね。
    TRAINING: % variance explained
    以下にエラーdimnames(tbl) <- list(c("X", yvarnames), paste(1:object$ncomp,  : 
           'dimnames' の長さ [1] が配列の大きさと違っています
    また、上記のページにあるようにncompを3にすると
    > wine.pls <- plsr(Y ~ X, 3, validation="CV")
    以下にエラーLa.svd(XtY) : 'x' のに無限値か欠測値があります
    と表示されますね。。。??初心者過ぎてお恥ずかしいです。-- kazu? 2006-11-02 (木) 17:04:03

データ読み込みをしたあとの分析ができません

koike? (2006-10-31 (火) 23:28:31)

統計初心者です。因子分析を実際にやってみようと思いやってみたところ、データは読み込むことができるのですが

result <- pdf(dat)

と入力すると

Warning message:  
Communality >= 1. in: pdf(dat) 

とでてしまい、それ以上なにもできなくなってしまいます。
変数が多かったからいけないのかと何度も変えてみてもできません。
どうしたら、因子分析できるようになるのでしょうか?

pdf.datの中身は、適当にうちこのような数値です。

   バ ス ミ マ ク チ メ コ ア ラ か き キ ウ カ プ オ
1   7  7  8  3  9  9  6  8  3  2  1  7  9  3  3  9  7
2   7  8  9  9  9  9  2  7  7  9  5  4  7  8  9  9  4
3   7  4  3  3  6  4  7  7  6  3  3  4  6  3  3  6  7
4   9  6  6  5  8  6  8  9  6  5  4  4  9  5  4  8  9
5   9  5  7  5  6  8  4  4  4  6  9  6  6  6  5  9  6
6   5  7  5  5  5  7  5  8  5  8  9  5  9  5  8  8  9
7   9  7  6  3  7  9  4  6  5  2  2  6  7  5  5  9  9
8   7  7  6  8  7  6  5  8  7  5  6  5  6  7  9  5  6
9   7  7  4  8  7  7  4  7  6  3  5  5  7  7  6  8  6
10  5  5  9  5  8  8  5  6  3  3  3  8  7  2  8  8  6
11  8  1  8  9  7  9  3  9  4  7  4  4  9  7  6  6  7
12  9  4  7  5  8  4  5  7  5  6  9  4  4  5  8  9  5
13  7  7  6  8  8  9  6  7  7  6  6  5  9  7  5  9  7
14  9  6  4  6  9  9  3  6  8  9  4  6  7  8  4  9  9
15  8  7  5  5  5  9  4  5  4  3  9  4  4  6  5  8  6
16  7  4  7  4  7  3  3  7  5  3  6  6  6  6  6  4  6
17  6  4  2  2  7  6  4  5  5  8  3  4  6  6  5  7  5
18  8  7  5  5  6  8  3  5  6  3  7  3  5  4  5  7  9
19  5  9  8  5  7  5  5  5  9  4  7  6  6  7  7  8  7
20  6  6  6  8  8  6  1  5  7  8  4  5  7  7  7  7  6

もう一つですが、因子分析はやはり変量がサンプル数より多いと実行できないのでしょうか?

  • コメントありがとうございます。
    やろうとしていたことは、データを読み込んでの因子分析です。青木先生が書かれている因子分析を参考にしています。
    始めに読み込むデータファイル名を pdf.dataとしていました。しかし、
    dat <- read.table("pdf.data", header=TRUE)
    と打ち読もうとすると、ファイルがないとエラーがでてしまいました。どうして読み込んでくれないのかわからなく試行錯誤しpdf.datとdatファイルに直したところ読み込めたので、pdf.datとしました。読み込めた中身を確認するためdatと打つと、ちゃんと確認することができました。
    そして、因子分析を実行するため
    result <- pdf(dat)
    と入力したのです。その結果エラーが
    Warning message:
    Communality >= 1. in: pdf(dat)
    とでてしまったのです。
    「因子数が大きすぎるが大きすぎる、デフォルトで決まる数より少なくすれば」というのは、どういうことなのでしょうか?因子数が大きいとは因子の数が多いということでしょうか?またデフォルトで決まる数よりも少なくすることは、どうしたらできるのでしょうか?
    初心者な回答ですみません。また至らぬ点があるとは思いますが、よろしくお願いします。 -- koike? 2006-11-01 (水) 11:48:52
  • 貴方のやったことを整理しましょう。データファイル名は pdf.data だというんですね。実際にワーキングディレクトリにあるファイルは pdf.data ですか?たぶん違うのではないでしょうか?それと,因子分析のプログラムとか書いてありますが,第三者にはどこにあるプログラムかすぐにはわからないのではないでしょうかね。http://aoki2.si.gunma-u.ac.jp/R/pfa.html のことだと思いますが,そこに書かれている因子分析を行う関数名は pfa なんですけど,それを貴方は pdf という名前に書き換えたのですか?そうではないということなら,result <- pdf(dat) としたのがまず変ですよね。それは result <- pfa(dat) だったのではないでしょうか?さて,実際はその関数を適用したものであるとして,貴方は因子分析についてどのくらい理解しているのでしょうか。抽出すべき因子数とか,因子数が不適切だと共通性が1を超えることがあるとかについては,ご存じなんでしょうか?たぶん知らないのだとは思いますが,pfa 関数の説明を見ると,引数に factors というのがあり,それは「求める因子の数(指定しない場合には,固有値が 1 以上の数になる)」ということですよね。この引数で,因子数を決めればよいのでしょう。
    実際には以下のようになるでしょうね。 -- 2006-11-01 (水) 12:02:39
    > x <- read.table("pdf.data", header=TRUE)
    > result <- pfa(x, factors=3)
    > result
    $rotation
    [1] "Varimax"
    
    $correlation.matrix
                バ          ス          ミ          マ          ク          チ          メ
    バ  1.00000000 -0.21985732 -0.08560197 -0.01360948  0.09620525  0.21050276  0.02883936
    以下略。。。。

maptoolsで緯度経度の範囲をカッチリと指定・プロットするには?

浜田? (2006-10-31 (火) 09:56:56)

以下のスクリプトで、指定した範囲の地図を描画させたいのですが、余分な範囲まで、含まれてしまいます。指定の範囲(xlim,ylim)だけをカッチリとプロットさせるには、どうすれば良いのでしょうか?

library(maptools)
jpn <- read.shape("japan_ver60.shp")
png(file="test.png", width=800, height=800)
plot(jpn, xlim=c(135.0,136.0), ylim=c(34.0,35.0))
dev.off()

実行結果は、

test.png

となってしまいます。、(図が大きくて、すみません・・・・)

  • 軸目盛を計算させねばいいわけだから, par(mar=c(0,0,0,0));plot(jpn, xlim=c(135.0,136.0), ylim=c(34.0,35.0),axis=F) でどうでしょう. -- なかま 2006-10-31 (火) 10:51:42
  • png とplotの縦横比が異なるので、plotが地図の表示はアスペクト比1:1であるので表示範囲を調整したのではないでしょうか? -- la-la-la-la? 2006-10-31 (火) 10:57:10
  • ああ, やってみるとダメですね. んじゃ
    png(かっちりサイズ)
    plot.new()
    par(mar=c(0,0,0,0),usr=c(135.0,136.0,34.0,35.0)) # usrはかっちり指定してください
    plot(jpn, xlim=c(135.0,136.0), ylim=c(34.0,35.0),add=T,axis=F)
    # ワーニング等は気にしてはいけない.
    axisが怪しい...ですが-- なかま 2006-10-31 (火) 11:21:01
  • おお! 出来ました! 有難うございました。 別途用意していたクリッカブルマップ用のデータに重ねるので、カッチリしたものが必要なわけでした。 助かりました。 -- 浜田? 2006-10-31 (火) 12:08:53

本の内容教えてください

g? (2006-10-30 (月) 20:41:20)

RでGISの"関連情報"に「Springer より、R を使った本が出る予定」との記事があったのですが,これは,「Rの基礎とプログラミング技法」のことを指しているのだと思いますが,この本は,GISについての記述があるのでしょうか? amazonの目次紹介では,それらしき項目はなさそうなのですが,購入された方,教えてください.

  • 言葉足りずで申し訳ありません。「Rの基礎とプログラミング技法」とは別物です。今年来日されたRoger Bivand 氏等(他2名?)が執筆しているもので、空間統計についてのもので、Springerから出版予定です。現在、Springer本家には情報がないようです。 -- 海の見える放送局? 2006-10-31 (火) 10:07:34
  • 詳しい情報をありがとうございました. -- g? 2006-10-31 (火) 10:36:35
  • Springer は R を使った本のシリーズを計画中みたいですね。 -- 2006-10-31 (火) 20:54:04

WindowsでRを使わせる

青木繁伸 (2006-10-30 (月) 18:46:38)

私はWindowsがほとんどわからないのに(Macならある程度わかるんですけど),Windowsもあまり知らない初心者にRを使わせようとして四苦八苦しています。
「こんなときにWindowsではどうやるのだっけ?」例えば,ファイルの拡張子を表示するように設定するにはどうやれば良いのか。まるっきり知らないわけではないので,最終的にはキーワードなんかでググってみてやり方を思い出したり。
というグチはこれくらいにして,今日聞きたいのは,プロジェクト(課題)ごとにディレクトリを作り,あるプロジェクトのデータ解析はいつもそこで行うというように決めるとき。ただ,一つの課題が終わるまで,ずーーーと同じディレクトリばかりで仕事をしているわけではない。あちこちのディレクトリを渡り歩いて,少しずつ仕事をするような場合。
Rを起動して,どのプロジェクトの仕事をしようかと言うとき,まずやらなくてはならないのは「作業ディレクトリ」の変更。毎度やるのも面倒だという場合には,Rのディレクトリに .Rprofile を作って,その中に setwd("/foo/bar") と書いておけばよいが,違うプロジェクトの仕事をするときには,コンソールから setwd("/foo/baz") とすれば,移動はできる。しかし,R を終了して再度 /foo/baz で仕事を継続しようとしても,R は /foo/bar を作業ディレクトリにしてくれる(怒)。しばらくは /foo/baz で仕事をするとすれば,.Rprofile を書き換えないといけない。面倒だ。
こんなとき,皆さんはどうやっているのでしょうか。
なお,いろいろやっていて,どのディレクトリで仕事をしようかというのが決まれば,そのディレクトリを開いて,そこにある .RData をクリックすれば,そのディレクトリを作業ディレクトリとして R が立ち上がるということはわかりました。(もしかして,これが一番良いやり方なんでしょうか?)
もっと良い方法はないんでしょうか。

  • ショートカットのプロパティの作業ディレクトリをそもそも変えてしまうのが一番簡単だと思います. Rguiの引数に環境変数 R_PROFILE,LANGUAGE等を指定するのも手だと思います. あとはそれらの自動化と言うことが一番嬉しいかと思いますので, VBスクリプトなど組むのがいいんじゃないかと.
    // hoge.js
    var Shell, FS, Desktop,Desktopex, path, Link, Rpath;
    Rpath = "c:??Program Files??R??R-2.4.0??bin??Rgui.exe";
    Shell = WScript.CreateObject("Wscript.Shell");
    FS = WScript.CreateObject("Scripting.FileSystemObject");
    Desktop = FS.GetFolder(Shell.SpecialFolders("Desktop"));
    Desktopex = FS.CreateFolder(Desktop + "??演習");
    path = Desktopex.Path + "??R-2.4.0_ex01.lnk";
    Link = Shell.CreateShortcut(path);
    Link.TargetPath = Rpath ;
    Link.Arguments = "LANGUAGE=en";
    Link.WorkingDirectory = "C:??";
    Link.Save();
    ダブルクリックなり,バッチ(CMD.EXE)の中ではファイル名で実行してくれます. -- なかま 2006-10-30 (月) 20:25:54
  • ありがとうございます。検討して貰うことにします。 -- 青木繁伸 2006-10-31 (火) 12:09:35

論理値のベクトル同士の演算

げのまー? (2006-10-30 (月) 11:58:10)

同じ長さの論理値のベクトル同士から,各要素ごとにandをとってベクトルが得られるような方法をご教示ください。
例えば,以下の例でaとbからc(TRUE,FALSE,FALSE,FALSE)というベクトルが得たいのです。
ループで各要素ごとに処理する方法は知っておりますので,ベクトル全体で演算する方法を知りたいです。
以下の方法では,1個目の要素についての結果しか得られておりません。

> a<-c(TRUE,FALSE,TRUE,FALSE)
> b<-c(TRUE,TRUE,FALSE,FALSE)
> a && b
[1] TRUE
  • 下の方(下5つめ)にも同じような質問がありましたね。& と && の違いを調べてみましょう。 -- 2006-10-30 (月) 12:26:17
  • 解決しました。Tips大全の論理判断のところで&と&&の違いが分かるようにしておきました。 -- げのまー? 2006-10-30 (月) 12:50:56

UTF-8 の読み込み

93795号? (2006-10-30 (月) 10:40:41)

中国語を含むデータを使用する必要があるのですが,データをうまく読み込むことができずに苦戦しています。

「ファイル Tips 大全」の「特定の文字コードのファイルを read」を参考に encoding を UTF-8,改行コードを CR/LF として,以下の3行を含むテキストファイル(utf8.txt)を作成しました。

  石巻市
岩出山町
  岩沼市

そして,ファイルを読み込もうとしたのですが,エラーが出てしまいました。

> read.table(file(description="utf8.txt",open="r", encoding="utf-8"))
以下にエラーread.table(file(description = "utf8.txt",  : 
        ファイルの開始部分が空です
追加情報: Warning messages:
1: 入力コネクション 'utf8.txt'  に不正な入力がありました 
2: 'utf8.txt' の readTableHeader で不完全な最終行が見つかりました 

同じ内容のファイルを Shift-JIS で作成した場合は問題なく読み込めました。
この先,何らかの対処をするにしても,どこから手をつけたらよいのか見当がつきません。
皆さんのお知恵を拝借できれば幸いです。
使用環境等は下の通りです。よろしくお願いいたします。

> sessionInfo()
R version 2.4.0 (2006-10-03) 
i386-pc-mingw32 

locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;
   LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932

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

> iconvlist()
....
[363] "UTF-16"
[364] "UTF-16BE"
[365] "UTF-16LE"
[366] "UTF-32"
[367] "UTF-32BE"
[368] "UTF-32LE"
[369] "UTF-7"
[370] "UTF-8"
....
  • そのデータファイルにはヘッダーがないわけですから,
    read.table(file(description="utf8.txt",open="r", encoding="utf-8"), header=FALSE)
    としなければなりませんね。 -- 2006-10-30 (月) 11:46:00
  • ご回答,ありがとうございます。でも,header=F を指定してもやはり読み込めませんでした。 -- 93795号? 2006-10-30 (月) 12:21:04
  • 実際にやってみてからのコメントだったので,私のところではできた理由がわかりませんねぇ。 -- 2006-10-30 (月) 12:25:23
  • そうですか... read.table は header=F がデフォルトなので,header=F の指定の有無が問題ではないと思うのです。しかし,何が問題なのか,どこから手をつけるべきなのか,手がかりがわからなくて困っています。 -- 93795号? 2006-10-30 (月) 12:33:19
  • 2: 'utf8.txt' の readTableHeader? で不完全な最終行が見つかりました 。。。というのは,ファイルの最終行に改行コードが付いていないという理由からです。 -- 2006-10-30 (月) 12:35:50
  • ありがとうございます。最終行にも改行を入れて保存したつもりなのですが,ファイルがきちんとできていないのでしょうか?保存したファイルを16進で表示すると
    00000000: efbb bf20 2020 e79f b3e5 b7bb e5b8 820d
    00000010: 0a20 e5b2 a9e5 87ba e5b1 b1e7 94ba 0d0a
    00000020: 2020 20e5 b2a9 e6b2 bce5 b882 0d0a     
    のようになっています。しかし,情けないことに自分ではこれで問題ないのかどうかがよくわかりません。 -- 93795号? 2006-10-30 (月)
  • efbbbf(最初の3オクテット)はFFFEのばいとおーだーまーくをutf8に変換したら出てくる物で, utf-8上では無意味なコードです(オーダーいらないでしょ?). iconvはutf-8のBOMを処理しないので変なコードだと文句を垂れます. nkfやiconvが使えるなら除去可能ですが, よくわからないならエディターの設定項目とか,他のエディターで処理できないかやってみて下さい. -- なかま 2006-10-30 (月) 15:08:12
  • ”efbbbf”はそういう意味で,それがあるからダメなんですね。勉強になりました。私が作ってみたUTF-8ファイルと違うので,ワケワカメでした。UTF-8では使わない方がよいとされているBOMだソウですが,変わり者のエディタが付けちゃうんでしょうね。 -- 2006-10-30 (月) 17:24:12
  • みなさん,ありがとうございました。UTF-8のファイルは Shift-JIS のファイルを xyzzy というエディターでコード変換して作成したものです。読み込めない原因と対処の方針がわかったので,他のエディターなどでもう一度試してみます。ありがとうございました。 -- 93795号? 2006-10-30 (月) 17:34:14
  • ご報告。
    nkf を使って文字コードの変換を行ったところ,無事に読み込むことができました。ありがとうございました。
    しかし,肝心の中国語を含むデータではエラーが出てしまいました。
    > read.table(file(description="test.txt",open="r", encoding="utf-8"), header=T, fill=T)
      A B    C
    1 1 2 中国
    2 2 1 日本
    3 3 1 ?国
    4 4 2     
    Warning messages:
    1: 入力コネクション 'test.txt'  に不正な入力がありました 
    2: 'test.txt' の readTableHeader で不完全な最終行が見つかりました 
    4行目のC には「两国」と入力されているのですが,読み込んでくれませんでした。なかなか難しいです。 -- 93795号? 2006-10-30 (月) 20:13:05
  • あら〜中国語もあるんですか...nkfは日本語しか処理しません. マルチリンガルならutf-8環境を持ったMACなりUnixで処理するか, どうしてもWinなら台湾(CP950)に切替えてのBig5の文字数に期待をかけるか...Window環境ではUtf8は支援しませんので... -- なかま 2006-10-30 (月) 20:36:42
  • R に加えて文字コードや言語環境(?)についての勉強もしなくてはならないようで,かなり遠い道のりになりそうです。当面は中国語の部分を R 以外で処理することにして,ごまかしごまかし やっていこうかと思います。なかま先生をはじめ皆さん,本当にありがとうございました。 -- 93795号? 2006-10-30 (月) 20:53:03

R 2-4-0 へのXMLパッケージのロード

2-4-0?? (2006-10-29 (日) 15:46:45)

はじめまして。お世話になります。

Windows binaryのR2-4-0をWindows XP上で使っています。XMLパッケージ(XML_0.99-93)をロードしようとすると以下のエラーが出ます。(XEmacs-21.4.19、ess-5.2.8を使った結果ですが、Rguiでも同じでした)

> library(XML)
Error in dyn.load(x, as.logical(local), as.logical(now)) : 
	unable to load shared library  'c:/PROGRA~1/R/R-24~1.0/library/XML/libs/XML.dll':
  LoadLibrary failure:  The specified procedure could not be found.

Error: .onLoad failed in 'loadNamespace' for 'XML'
Error: package/namespace load failed for 'XML'

何か解決法はありますでしょうか? 2-4-0の変更点のページか何かで、”methodsを使っている関数は書き換えが必要”というような記述を見たような記憶があるのですが、それでしょうか?

また、ヘルプを呼び出そうとしても似たような以下のエラーがでます。

> ?c
Error in  print.help_files_with_topic("c:/PROGRA~1/R/R-24~1.0/library/base/chm/c") : 
	CHM file could not be displayed

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

  • 10/26に出た最新版(1.1-1)に更新すればいいだけの話ではないですか?こちら(R-2.4.0 for Win)ではそれで何の問題もなくロードできますよ。 -- 2006-10-29 (日) 19:06:48
  • ありがとうございます。どこにあるのですか? -- 2-4-0?? 2006-10-29 (日) 20:46:13
  • パッケージのアップデートとかいうメニューがあるんじゃないですか

作業スペースの表示について

永野? (2006-10-24 (火) 12:06:26)

作業スペースの読み込み後LOADされたコンソール画面を表示するにはどうしたらよいですか。

  • 「LOADされたコンソール画面」というのは、コマンド履歴のような物のことですか?だとしたら、例えば、Rをインストールしたディレクトリ(フォルダ)中の、.Rhistoryというファイルに入っています。作業スペースというのは(命名がよくないと思いますが)、実際には、作業履歴・環境ではなく、Rオブジェクトをシリアライズしたもののことなので、「LOADされたコンソール画面」というようなものは存在しないと思いますよ。LOADされるのはRのオブジェクトです。 -- T. I.? 2006-10-24 (火) 13:45:15
  • history() の事でしょうか. -- なかま 2006-10-24 (火) 14:29:08
  •    例えば新規コンソール画面で<3+3=6としてこの画面をC:?にBBBというファイル名で作業スペースとして保存するとします。それから作業スペースを読み込むでBBBを読み込むと画面に load C:¥BBBと表示されますがこのBBB即ち3+3=6を画面に表示するにはどうすればよいでしょうかということです。 -- 永野? 2006-10-25 (水) 16:56:49
  • プロンプトの向きが逆で,かつ,<3+3=6はエラーになりますが。なかまさんがコメントしたhistory()とは,違うのですか? -- 2006-10-25 (水) 18:16:43
  • はて、見たままの画面をセーブすることは(R の機能としては)できないはずですが。sink 又は capture.output は R のコマンドの(テキスト)出力をファイルに落としてくれますが、コマンドまで含めて記録してくれませんし、ロードできるような形式ではありません。一方、R のコマンド履歴は .Rhistory ファイルに記録され、history() 関数で一覧したり、上下矢印で順に過去の命令入力を再現できますが、出力そのものを記録するわけではありません。単に全てを記録したかったら、Unix-like OS では R | tee foo.txt と起動すれば、コマンド及びその出力全てがファイルに記録されますが、次のセッションでそのまま(画面に)再現できるわけではありません。 -- 2006-10-25 (水) 18:59:21
  • やっぱり、「作業スペース」というものに誤解があるんじゃないかなあと思います(私も同じような勘違いをしていたので)。保存できるのは、Rオブジェクトなので、3+3=6が保存したければ、x <- function() { 3+3 } のように、なんらかのオブジェクトに包んでやらないといけません。このxは保存できるので、3+3だとか、結果の6を後で得ることができます。 -- T. I.? 2006-10-25 (水) 23:23:08
  • いろいろとご教示ありがとうございました。 -- 永野? 2006-10-26 (木) 12:48:37
  • テキストエディタ経由なら「見たままの画面」をコマンド付きでそのまま保存できます。保存したファイル(テキストファイル)からコマンドを再実行することもできます。例えば、ESS。 -- 谷村 2006-10-27 (金) 10:45:45
  • そういう意味なら,コンソールは,入力・出力の両方とも,見たままをファイルに「保存」できますね。そのファイルをRのエディタで読んで,範囲を選択してcommand+returnで実行できますし。~作業スペースの保存・読み込みの対象は「.RData」ですね。 -- 2006-10-27 (金) 11:23:35

Rで扱える行列のサイズの上限

pahud? (2006-10-23 (月) 21:24:34)

Rで扱える行列のサイズの上限は,Rがインストールされているプラットホームによって異なるのでしょうか.それとも,Rが仕様として上限を設定しているのでしょうか.

  • 一種の興味からしりたいのでしょうか?それとも,他にもっとたくさんのメモリーを使えるプラットフォームがあるならそちらを使いたいということがあるのでしょうかねえ
    あなたが使っているのは,どのプラットフォームでしょうか。そして,行列のサイズの上限をチェックしたのは,どのような手段で行ったのでしょうか?それが分かれば,あなたと違ったプラットフォームでRを利用している方のコメントも得られるかもしれませんね。もっとも,Rのソースの,ここに,こういう風に書かれているから,分かるでしょう?というコメントも(なかまさんから?)得られるかもしれませんが。 -- 2006-10-23 (月) 21:52:26
  • 経験的にいえば,そんなぎりぎりの制限に挑戦しなくても,問題は解決可能なことが多いと思います。はっきりいえば,エレガントな解法が必要とするメモリ量は,シンプルな解法が必要するメモリ量の数分の1,極端な場合には数百分の1にもなる。実行時間もおおむねそれに比例する。。。
    どのような場合であろうと,あなたが何を知りたいのか,はっきり・具体的に書くのが吉。 -- 2006-10-23 (月) 22:38:51
  • 32bitなら例えば2GB(ベクトルデータサイズ)までなら作成可能です.しかし,計算結果を格納したりその他のデータも考えなければいけないので, それは無意味な上限でしょう. ちなみに64bitなら46340^2が(メモリーがあれば,現在の)上限になります. -- なかま 2006-10-24 (火) 00:22:08
  • ソースコードを見ると(array.c)、allocMatrix()で、行x列 <= INT_MAX をチェックして、それ以下になるようにしているので、INT_MAX=32767 (16bit)な環境なら、200x200行列は仕様上無理、ということになります。32bitでは、例えば、50000x50000行列は、仕様上無理、ということになります。 (行列が内部では配列として扱われること考えれば、あたりまえのことですが。) 実際には、メモリのほうが、先にたりなくなると思います。-- T.I. 2006-10-24 (火) 13:03:14
  • INT_MAXは2^31-1です(INT_MAXが16bitと言うのは...お年が...以下自粛). 16bit環境と言う物が仮りにあってもRは動かないでしょう. 仰せのとおり行列は内部ではただの一次元配列+ヘッダーのサイズが2G以下なら作れますが, そもそも使える全体のサイズが2Gだったり3Gだったりするので M<-A%*%B なら最大でも1G以下になります. ですので最適解としては32bitなら実メモリ次第ですし,64bitならメモリが16G以上あるなら,46340^2が上限と言うことです.決して200x200等(小さすぎ!)と小さな行列ではありませんです.はい. -- なかま 2006-10-24 (火) 13:40:37
  • なかまさん、すみません。どうも、更新するまえのバージョンに、コメントされているみたいです。(INT_MAX=32767は例です。質問者が、プラットフォームによって、上限が異なるのかどうか効いているので。16bitのコンピューターは私もほとんど使ったことありませんです(そんなに年寄りじゃない(^^;;; -- T.I.? 2006-10-24 (火) 13:53:37
  • がーん, 使った事ある...σ(T_T) -- なかま 2006-10-24 (火) 14:27:12
  • みなさま回答ありがとうございました.INT_MAXは知りませんでした.勉強になりました.なお,今回の質問は,冒頭で明示したように,乗り換えを意図した上でのものです. -- pahud? 2006-10-25 (水) 07:52:15

ifelseを使った文字列の置換について

matsu? (2006-10-23 (月) 11:56:05)

はじめまして。こんにちわ。
オブジェクトの指定について質問させてください。

> data8
    V1 V2 V3 V4 V5 V6 V7 V8
1   A  B  B  B  B  B  B  B
2   B  B  B  B  B  B  B  B
3   A  A  A  B  A  B  A  A
4   B  A  B  A  B  A  B  B
5   B  B  B  B  B  B  B  B
6   B  B  B  B  B  A  B  B
7   B  B  B  B  B  B  B  B

以上のようなデータがあって、例えば1列目と2列目の同じ行の要素を見て、A,AであればA、A,BであればB、B,AであればC、B,BであればDというように置き換えて新しいベクトルを作りたいので、まずは試しに2区分するために以下のように書きました

ifelse(((data8[,1] == "A") && (data8[,2] == "A")),"A","B")

しかし結果は

[1] "B"

と表示されベクトルを得ることができません。4区分にしても同じようになります。これを解決する方法はあれば教えて下さい。
よろしくお願いします。

  • & と && の違いを調べてみましょう -- 2006-10-23 (月) 12:32:50
  • できました。tipsをよく読めばできるような愚問でした。申し訳ありません。 -- matsu? 2006-10-23 (月) 13:54:20
  • 例解:factor(c("A","C","B","D"))[as.integer(data8[,2])*2-2+as.integer(data8[,1])] -- 2006-10-23 (月) 14:16:50
  • 別解:factor(diag(matrix(LETTERS[1:4],byrow=T,2)[data8[,1],data8[,2]])) -- 2006-10-23 (月) 14:22:19
  • 別解2:factor(paste(data8[,1],data8[,2]),labels=LETTERS[1:4]) -- 2006-10-23 (月) 14:58:28
  • 別解3:factor(LETTERS[factor(paste(data8[,1],data8[,2]))]) -- 2006-10-23 (月) 16:08:09
  • 色々な回答ありがとうございます。こんなにブラッシュアップできるものだとは・・・。奥が深い。 -- matsu? 2006-10-23 (月) 21:14:44

カテゴリごとにまとめたデータの整形(データベースとテーブルの相互変換)

take? (2006-10-17 (火) 19:51:40)

下記のように、カテゴリー値を含めてデータが一列になったものを、 データのカテゴリごとにまとめた形に整形したいのですが、 このような作業をしてくれる関数はありますでしょうか。
subsetなどを駆使すれば無理やりできますが、頻繁に使われますよね。 非常に初歩的な質問かと思いますが、適切なキーワードが見つかりませんでしたので どなたかお助けくださると助かります。

変換前(データベース形式)

yeartypevalue
2001A5
2002A10
2003A15
2001B1
2002B3
2003B5

変換後(テーブル形式)

AB
200151
2002103
2003155
  • そんなにしょっちゅう使うものでもないし一般的でもないので,処理を行う関数を定義して,せいぜい簡単に呼び出せるようにしておくだけで十分なんではないでしょうか?関数の仕様も汎用化できるほど一般的ではないと思いますが。あれこれ探したり,聞いたりしている間に関数を書いた方が早いかも? -- 2006-10-17 (火) 20:52:57
    > # 引数は3列からなるデータフレームとする
    > convert <- function(df)
    + {
    + 	df[,1] <- as.factor(df[,1])
    + 	df[,2] <- as.factor(df[,2])
    + 	m <- matrix(0, length(table(df[,1])), length(table(df[,2])))
    + 	for (i in 1:nrow(df)) m[df[i,1], df[i,2]] <- df[i,3]
    + 	rownames(m) <- levels(df[,1])
    + 	colnames(m) <- levels(df[,2])
    + 	return(m)
    + }
    > df <- data.frame(x = c(2001, 2002, 2003, 2001, 2002, 2003),
    + y = c("A", "A", "A", "B", "B", "B"),
    + z = c(5, 10, 15, 1, 3, 5))
    > convert(df)
          A B
    2001  5 1
    2002 10 3
    2003 15 5
  • 例えば year=2001 で type=A のものが二つ有れば例の様にはまとめられません。最初のデータフレームを x とすると y <- table(x); y[,,"15"] は value=15 の値に対する year,type の組合せごとの重複度を込めた度数分布テーブルを計算してくれますが、これはご希望のものではないし... -- 2006-10-17 (火) 21:37:58
  • 21:37:58 のコメントは,原質問者ではないのですね。そうなんですよ,仕様が曖昧というか全く不明確なので,実際の関数を書くのも容易ではない。
    たぶん,データ行列への変換だと思うんですが。
    しかし,原質問者がそのあたりは一番よく分かっているので,誰よりも適切な関数を書けるでしょうということで。 -- 2006-10-17 (火) 21:39:41
  • このような処理はどこかの本で見たような気がしていたので、適切なキワードがわかれば、すぐに方法があると思ったのですが、一般的ではなかったのですね。確かに汎用化できるかというと、ご指摘のようにケース数によってあいまいな処理を含んでしまいます。例示してくださったものを参考に自分でちょっといじって書いてみます。ありがとうございます。 -- take? 2006-10-18 (水) 01:16:03
  • キーワードということなら、(Rと関係なく恐縮ですが)、Excelのピボットというのが、処理として似ていると思います。例えば、http://www.atmarkit.co.jp/fwin2k/win2ktips/359pivot/pivot.htmlでの実例など、takeさんの挙げられたものにそっくりです。 -- T. I.? 2006-10-18 (水) 12:19:31
  • 各カテゴリに重複するデータがないことが確実な場合は、アドホックなやりかたですが、tapply関数も使えるような気がするのですが、いかがでしょうか?さきの例で使っているデータフレーム df で実験すると、、。
    > # 重複したデータがないことを、lengthで確かめて、、
    > tapply(df$z,list(df$x,df$y),length)     
          A B 
     2001 1 1 
     2002 1 1 
     2003 1 1
    > # sum (または mean)で集計。
    > tapply(df$z,list(df$x,df$y),sum)     
           A B 
     2001  5 1 
     2002 10 3 
     2003 15 5 
    mmk? 2006-10-18 (水) 13:14:04
  • おみごと! -- 2006-10-18 (水) 14:14:28
  • すばらしい!(関数型言語の妙技)。Tips入りしてもよさそう。 -- T. I.? 2006-10-18 (水) 15:03:51
  • お探しの関数は、xtabs関数ではないでしょうか 。-- 2006-10-18 (水) 21:42:50
    > df <- data.frame(year = c(2001, 2002, 2003, 2001, 2002, 2003),
    +   type = c("A", "A", "A", "B", "B", "B"),
    +   value = c(5, 10, 15, 1, 3, 5))
    > xtabs(value~year+type,df)
          type
    year    A  B
      2001  5  1
      2002 10  3
      2003 15  5
  • 重複値があるとどうなるか気になって調べました。加算されてしまうようですね。 -- 2006-10-18 (水) 21:44:38
    > df <- data.frame(year=c(2001,2002,2003,2001,2002,2003,2001),
      type=c("A","A","A","B","B","B","A"),value=c(5,10,15,1,3,5,100))
    > xtabs(value~year+type,df)
          type
    year     A   B
      2001 105   1
      2002  10   3
      2003  15   5
  • これもおみごと!普通,重複はないと思うので(チェックするのが吉),これでよいと思いますね。 -- 2006-10-18 (水) 22:03:04
  • 皆様、ありがとうございます。やはりあるものなのですね。xtabsがそのものぴたりです。tapplyの技は大変関心いたしました、集計に色々応用できそうです。 あと、EXCELも使いこなせば結構使えるもんですねー 大変感謝いたします。 -- take? 2006-10-19 (木) 09:46:58

グラフのカラー・モノクロ線種切り替え

takahashi? (2006-10-17 (火) 16:26:21)

例えばモノクロで印刷される原稿用とプレゼン用で、同じグラフだけどモノクロ出力(この場合、線種で区別)、カラー出力を切り替えたいようなときは、どうするのが一番手っ取り早いかご存知でしょうか?

  • 線種で区別ということは,折れ線グラフの類?
    色の変更と線種の変更はパラメータが違うのですから,そうそう簡単な方法はないのでは?最初から意識して色と同時に線種も考えてグラフを描くようにしておき,印刷時にはグレースケールで印刷(あるいはモノクロコピー)すれば出来上がりの図は線種で識別できるという,はなはだローテクの手法が一番ぴったりなんではないでしょうか。
    どんな図でも,それをモノクロで印刷したときにもちゃんと読み取れるように,最初っから意識して作図しておくのが吉。たとえば,パイチャートやバーグラフなんかは色だけじゃなくパターンも変えておかないといけませんしね。使用する色も水色や黄色なんかを使うとモノクロで印刷したら線が見えないなんてこともありますので,これは経験がものをいうのでしょうね。ちょっと意味合いが違うが,色覚に問題のある人にも識別できる色合いを使うというのも奥ゆかしい心遣い。 -- 2006-10-17 (火) 21:29:43
  • どうもありがとうございます。やっぱりそんな感じですかね。出力デバイスのオプションでモノクロのときはカラー指定をlty、pch、ヒートマップなんかではグレースケールに読み替えてくれるようなスマートなものがあればいいんですけど。 -- takahashi? 2006-10-18 (水) 16:15:20
  • スマートではないですが、ImageMagick?が入っているなら、 -- 谷村 2006-10-23 (月) 15:13:39
    > hist(sqrt(islands), br = 12, col = "lightblue", border = "pink")
    > dev.copy2eps(file="color.eps", family="sans") #カラー
    > dev.copy2eps(file="| convert -colorspace GRAY - bw.eps", family="sans") #グレイ階調
  • convertを通すとベクトル画像ではなくなるみたいです。ベクトル画像のままグレイ階調化するコマンドにbw_convertがありますが、標準入力からデータを読んでくれないので、system()を使う必要があります。 -- 谷村 2006-10-23 (月) 15:24:28
    > hist(sqrt(islands), br = 12, col = "lightblue", border = "pink")
    > dev.copy2eps(file="color.eps", family="sans") #カラー画像に出力
    > system("bw_convert color.eps bw2.eps")
  • ありがとうございます。外部ツールと連動させればいけそうですね。試してみます。 -- takahashi? 2006-10-25 (水) 11:12:12

Rコマンダーでのエクセルファイルの読み込みについて

epep? (2006-10-14 (土) 09:33:30)

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

Rコマンダーのファイルインポートについて教えていただきたいと思います。
最近エクセルファイルの読み込みに対応したということで、エクセルファイルを読み込もうとしているのですが、下記エラーメッセージが出てしまいます。

要求されたパッケージ RODBC をロード中です
以下にエラーsqlTables(channel) : first argument is not an open RODBC channel

データをいろいろ変えたり、エクセルの保存形式を古いものにしたりと、試したのですが、うまくいきませんでした。
どうすれば読み込めるようになるか、教えていただきたいと思います。
よろしくお願い致します。

  • コメントありがとうございます。
    単純に"1"だけを入力したエクセルファイルを、インポートのメニューから指定して読み込んでもエラーメッセージが出るだけで読み込んでもらえませんでした。
    他にも1行目に変数名を入れて、2行目に数字を入れる、更には、一度Rでデータを入力してcsvでエクスポートしたファイルをエクセルで読み込んで、xlsにして保存してもダメでした。-- epep? 2006-10-14 (土) 09:50:34
  • 闇雲に試行錯誤するよりもエラーの意味を吟味した方が効率がよいと思いませんか。sqlTables(channel) : first argument is not an open RODBC channelをキーにGoogleで検索したら53件ヒットしました。これらは読まれましたか? -- 2006-10-14 (土) 22:11:27
  • コメントありがとうございます。すべて英語でしたが読みました。 -- epep? 2006-10-15 (日) 00:08:00
  • アドバイスいただいたことを考え、いろいろ試してみたところ、解決できました。そして私のパソコンだけなのかもしれませんが以下のことが分かりました。
    ・エクセルファイルを保存するフォルダ名に日本語(全角、半角とも)が使ってあるとエラーになる。
     階層の途中に日本語のフォルダがあってもいけない。
    ・エクセルファイルのファイル名に日本語が使ってあるとエラーになる。
    ・エクセルファイル内のシート名に日本語が使ってあるとエラーになる。
    です。csvやtxtだと日本語が使ってあっても問題ないのですが、エクセルファイルの読み込みには対応していないようです。
    ありがとうございました。 -- epep? 2006-10-15 (日) 02:24:44
  • パス名(ディレクトリ名,ファイル名)に日本語が使われていると,問題を起こすプログラムは多いです。あなたのパソコンだけの問題ではないでしょう。 -- 2006-10-15 (日) 09:24:16
  • 教訓:ナガイモノには巻かれてしまいましょう。いくら,日本語対応なんてあっても,なんか問題が生じたときには,日本語を使わないでやってみて,それでうまくいけば,「ああ,日本はまだまだ世界に通用しないんだなあ」と嘆きましょう。 -- 2006-10-15 (日) 21:59:41
  • ざっと, Rcmdrのソースを斜めに読むと, substringは文字で処理, ncharのデフォルトはバイト単位で処理します. パスの処理をする部分を nchar(hoge,type="chars")に置き換えればこの問題は解決する筈です(修正箇所はそんなに多くはありません). 普段Rcmdr(for Windows)の御世話になってる方はトライしてみて下さい. 嘆いていても始まりません. -- なかま 2006-10-16 (月) 10:27:13

オブジェクト名の指定について

dura? (2006-10-12 (木) 16:03:24)

はじめまして。こんにちわ。
オブジェクトの指定について質問させてください。
例えば

> R1
  
    A  B
 A 24  5
 B  5 16

> R2
  
    A  B
 A 22  7
 B  6 15

のようなデータがあるとき

m <- 1
n <- 2
R"m"*R"n"

のような感じで行列演算のオブジェクトを指定できるようなやりかたはありますか?
ちなみにR"m"*R"n"では当然できません。例えばこのような感じで、というつもりで書きました。

  • FAQですね。単語検索で eval(parse( を検索すると,用例もたくさん出てきますので,検索してみましょう(非難しているわけではないです。検索しろといっても,適切な検索語が浮かばなければ,できませんものね)。
    以下のようにすればいかが?
    ところで演算子は"*"でいいんですね。。 -- 2006-10-12 (木) 16:57:01
    > R1 <- matrix(c(24, 5, 5, 16), 2, 2)
    > R2 <- matrix(c(22, 6, 7, 15), 2, 2)
    > m <- 1
    > n <- 2
    > eval(parse(text=paste("R", m, "*", "R", n, sep="")))
         [,1] [,2]
    [1,]  528   35
    [2,]   30  240
  • ばっちりできました。またeval(parse( も検索してみました。便利ですね。非常に勉強になりました。ありがとうございました。 -- dura? 2006-10-12 (木) 23:42:16

GUIの変数ウォッチャー/データフレーム

takahashi? (2006-10-05 (木) 19:04:42)

data.frameのデータをスプレッドシート形式で、目で見て確認できるようなライブラリ(Rcmdrのデータシートみたいなもんです)があるかどうかご存知でしたら、教えて下さい。

また、Matlabをご存知の方に伺いたいのですが、Matlabには変数(構造体というんですかね)の中身を探れるGUIツールがあると思います。
Rでlistが同じように、複雑にネストした構造になりがちだと思いますが、これをGUIで探れるようなライブラリがあるかどうかご存知でしたら教えて下さい。

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

  • 前半については,edit() とか fix() があるようですが? -- 2006-10-05 (木) 19:23:14
  • ありがとうございます。edit.data.frameがちょうど求めていた感じです。盲点でした。 -- takahashi? 2006-10-05 (木) 19:48:19

データを読み込ませた後の分析ができません

4年生? (2006-09-27 (水) 19:23:05)

はじめましてこんにちは.
こないだからRを使おうと頑張っているのですが、分析がうまくできません.
エクセルで乱数を発生させて、54行×4列のデータファイルをつくってCSVファイルで保存しました(data.csv)
その後

dat <- read.csv("data.csv",header=FALSE)
dat

とすると54行×4列のデータが読み込まれます.
次に

summary(dat)

とすると

        V1             V2             V3             V4    
10,006.13: 1   10,038.09: 1   10,058.34: 1   10,100.57: 1  
10,077.88: 1   10,070.11: 1   10,112.10: 1   10,167.06: 1  
10,248.81: 1   10,077.88: 1   10,148.36: 1   10,169.66: 1  
10,264.61: 1   10,173.60: 1   10,171.17: 1   10,281.17: 1  
10,281.25: 1   10,293.74: 1   10,186.77: 1   10,284.54: 1  
10,303.26: 1   10,378.13: 1   10,213.75: 1   10,318.44: 1  
(Other)  :48   (Other)  :48   (Other)  :48   (Other)  :48  

となってしまい,うまく結果が表示されません.
重回帰分析をしようとして

dat <- matrix(c(read.csv("data.csv",header=FALSE)),byrow=TRUE,nc=4)
mreg(dat)

としても

g(dat) :  (subscript) 論理値添え字が長すぎます

となってしまいます.
googleで検索したり色々頑張りましたがうまくいきません・・・
どこに問題があるのでしょうか?
教えてください.よろしくお願いします.

  • すみません,表示したデータは乱数で発生させたものではなくていつかの日経平均の推移です. -- 4年生? 2006-09-27 (水) 19:28:27
  • エクセルで表示を直して,再度挑戦したら分析できました.ご迷惑おかけしました.ありがとうございました. -- 4年生? 2006-09-27 (水) 20:23:14

リストの出力方法を教えてください

Real? (2006-09-25 (月) 14:14:58)

こんにちは。
以下のような数値データが入ったリストを書き込もうとすると

a <- as.list(NULL)
a[[1]] <- c(2.3, 3.5, 4.4, 2.3, 3.5, 4.4, 2.3, 3.5, 4.4)
a[[2]] <- c(2.3, 3.5, 2.3, 3.5, 4.4)
write(a, "aaa.txt", append=T)
以下にエラーcat(list(...), file, sep, fill, labels, append) : 
       argument 1 (type 'list') cannot be handled by 'cat'

と怒られてしまいます。
結果はこのように一番目の要素から順に横に展開したいと思っています。

2.3, 3.5, 4.4, 2.3, 3.5, 4.4, 2.3, 3.5, 4.4
2.3, 3.5, 2.3, 3.5, 4.4

リストを書き込むときはwrite()を使えばよいという情報があったのでこうしたのですが、どうすればよいのでしょうか?
どなたかご指導をお願いします。

  • lapply(a,function(x){write(x, "aaa.txt", append=T,nc=length(x))}) -- takahashi? 2006-09-25 (月) 15:39:19
  • なるほど。リストの要素数が違うときはlapply()を使えばいいんですね。ありがとうございました。 -- Real? 2006-09-25 (月) 15:51:56
  • save(a, file="aaa.txt", ascii=TRUE) 。再度読み込むには load("aaa.txt") 。ascii=TRUE にしないと内部形式で書きこみ。 -- 2006-09-25 (月) 16:43:40
  • R出力の記録を読むと他にも色々な方法が。 -- 2006-09-25 (月) 18:30:59

axis 関数で軸の名前を付けるには?

匿名希望? (2006-09-21 (木) 13:24:21)

x 軸を本来の位置ではない場所に描くとき,plot(...., xaxt="n", xlab="", ...) として,axis(1, pos=foo, ...) としますが,そうすると,x 軸のラベルは付きません。plot 関数の中の xlab="" を取り除くと,本来 x 軸が描かれる位置に相応して x 軸のラベルがつきますが,これは当然意図した場所ではありません。オンラインヘルプを探したのですが,解決法は見つかりませんでした。

  • 補足:par(mgp=c(foo, bar, baz)) の foo を試行錯誤で調整すれば良いようですが,これでは x 軸,y 軸の両方とも同じ相対調整量なので,別々に指定できないとちょっと困ります。 -- 2006-09-21 (木) 13:37:53
  • mtext()はどうですか? -- Akira? 2006-09-21 (木) 13:42:54
  • ありがとうございます。求めていたのは,これです。解決しました。 -- 2006-09-21 (木) 14:12:09

MultiTiff?イメージを取り込む方法

Akira? (2006-09-19 (火) 20:02:59)

rgdalパッケージのreadGDALでtifイメージを取り込むことはできました。
MultiTiff?イメージを同じように取り込むことはできるのでしょうか?
WinXP, R-2.3.1patchです。

  • 自己レスです。rtiffというのがありました。こちらでMultiTiff?をテキストにできるか試してみます。 -- Akira? 2006-09-19 (火) 22:50:40
  • ご報告。rtiffは最大値が1にノーマライズされるようなので、rgdalの方が使いやすそうです。でも、こちらはMultiTiff?はだめっぽい。もう少し調べてみます。
    よくご存知の方、コメントいただけるととてもうれしいです。 -- Akira? 2006-09-20 (水) 21:21:04

コードを用いてグラフの保存

初心者?? (2006-09-19 (火) 10:44:39)

関数の中にいくつかのグラフを表示させるコードが入っています。
コードを使い
R gurahicsに表示されるグラフに名前をつけて保存させることは出来るのでしょうか?

自動化したいと考えてるのですが ご教授ください。

  • どういう類の名前を付けたいのですか。「ファイル名+連番.拡張子」みたいなのでいいのですか?(ファイル名の部分は自分で与えるとして) -- 2006-09-19 (火) 11:10:12
    draw.graph <- function(file.name="")
    {
    	if (file.name != "") pdf(paste(file.name, "%03i.pdf", sep=""), onefile=FALSE)
    	hist(rnorm(1000))
    	plot(rnorm(1000, 1000))
    	barplot(rbinom(10, 1000, 0.25))
    	if (file.name != "") dev.off()
    }
    draw.graph("test")
  • さきほど確認しました。この用に書けばいいのですね、勉強になりました。ありがとうございました。 -- 初心者?? 2006-09-20 (水) 00:02:49

グラフィックをアウトライン化したベクトル画像で出力したい

初心者? (2006-09-17 (日) 19:23:04)

お世話になっています
ポスター制作のため,Rで出力したグラフをdev.copy2epsで出力し,Adobe Illustrator CS2に貼付けているのですが,ギリシャ文字が文字化けしてしまいます(SymbolというフォントがIllustrator側には無いようです).
そこで,アウトライン化したベクトル画像を出力したいのですが,Rではそういったことは可能なのでしょうか?

環境 : R version 2.2.1, 2005-12-20, powerpc-apple-darwin7.9.0

  • ?embedFonts(実際にはgsでフォントを埋め込む)か, epsの中にSymbolを挿入(手で)するか, IllustratorにSymbolフォントを認識させるか. EPSを止めてPDFを渡す(一番良い)が考えられます. -- なかま 2006-09-17 (日) 20:08:23
  • コメントありがとうございます.Illustratorはpdfを読み込めるので,pdfで試してみようと思います. -- 初心者? 2006-09-17 (日) 20:18:51
  • Macintoshでは,いろんな場面でPDFがベストのようです。EPSはたタコです。 -- 2006-09-17 (日) 21:38:38
  • 出来ました。ありがとうございます。 -- DNA? 2006-09-19 (火) 09:19:36

会社(Windows2000)で使用する際に一台で複数ユーザで利用したい

カータン? (2006-09-16 (土) 12:34:25)

こんにちは,現在会社(一般ユーザはuser権限のみ)でRをインストールしたいのですが,admin権限でインストールした後,user権限でそのパソコンにログインしても,Rgui.exeが全く起動しないのですが,そのような物なのでしょうか。ちなみにPowerusersの権限で入っても同じでした。インストール時のユーザ名でないと駄目なのでしょうか。
もし上記のことが可能ならば,パッケージの導入の際にも,admin権限で導入して,その後user権限で使用可能かを現在検討したいのですが,その点もおわかりの方がいらっしゃったらよろしくお願い致します。Powerusers権限ならパッケージ導入も可能なのでしょうか。

  • admin権限の場合は,他の人がインストールした場合でも使用できました。また,poweruserでも自分でインストールした場合は使用できました。 -- カータン? 2006-09-16 (土) 13:19:11
  • Rのパーミッションを変更する -- 2006-09-19 (火) 08:34:39

自前で計算した距離行列をhclust()で使う方法

DNA? (2006-09-15 (金) 17:30:53)

こんにちわ。現在クラスター分析に取り組んでいます。手順としては、

>iris.x <- iris[,-5]
>d.euclid <- dist(iris.x, method="euclid")
>cluster <- hclust(d.euclid, method="average")

以上のように行っています。ただ、こちらの都合でd.euclidを自前で作ったものに差し替えたいのですが、自分で作った行列をhclust()に投げても上手くいきません。
どのようにすればよいのでしょうか?よろしくお願いします。

  • [自前で作ったもの]がどういうものか,どうやって作ったかがわからないので,コメント不能ですね。 -- 2006-09-15 (金) 17:34:32
  • すみません自前でつくったものは、 -- DNA? 2006-09-15 (金) 17:41:30
     A B C
    "A" 1 0.5 0.1
    "B" 0.5 1 0.4
    "C" 0.1 0.4 1
    といった行列をつくって、read.table()で読み込んでいます。
  • 出来ました。ありがとうございます。 -- DNA? 2006-09-15 (金) 18:05:43

データフレームから要素を抽出してリストを作製する

tottie? (2006-09-11 (月) 14:04:05)

こんにちは。以下のようなデータフレームがあり、

NAME	num_1	num_2	num_3	STR	NUM
data_M15  0.322  0.322 -0.040    asd      3
data_G11  0.230  0.507 -0.279    sdf      2
data_B21  0.219  0.160 -0.069    dfg      8
data_C04  0.011  0.292 -0.064    fgh      3
data_J12  0.254  0.238 -0.180    ghj      1
data_P08  0.323  0.320 -0.254    hjk      2
data_M11  0.356  0.330  0.066    jkl      1


一番右列の数字によってSTRを分類したリストを作りたいと思っています。
たとえばリストの一番目の要素にNUM3のSTRであるasdとfgh、
二番目の要素にNUM2のSTRであるsdfとhjkが格納されているリストを作りたいです。

subset()を使用してNUMの値ごとに抽出すればよいとは思うのですが、
どうすればよいのかわかりません。
お手数ですが、ご指導をお願いします。

  • lapply(unique(data$NUM),function(n)data[data$NUM==n,]$STR) factorなのが嫌なら、lapply(unique(data$NUM),function(n)levels(data$STR)[data[data$NUM==n,]$STR]) -- takahashi? 2006-09-11 (月) 14:19:02
  • subset()を使わなくても添え字操作とlapplyの組み合わせで可能なんですね。lapply()の中に関数を定義して引数を渡すとは驚きました。ありがとうございました。 -- tottie? 2006-09-11 (月) 15:05:54

グラフの目盛りを有効数字n桁の小数x10のべき乗で表示したい

fjok? (2006-09-10 (日) 13:28:59)

質問させていただきます.
グラフの軸の表示をデフォルトの1 e-01から有効数字n桁の小数掛ける10のべき乗,という書式に変更する方法をご存知の方はいらっしゃいますでしょうか.よろしくおねがいします.

  • どうしてそう言うことが必要なのか理解できませんが,以下のようにすればよいでしょう。axis 関数を使うわけですね。かくかくの値の所にしかじかのラベルを描いてねとお願いしましょう。どう言うところにどういうラベルを描くかを自動的に行うのは,あなたが考えてください。ちゃんと,汎用解がありますから(実現方法はいろいろだけど) -- 2006-09-10 (日) 19:43:57
    set.seed(123)
    x <- rnorm(1000)*100+600
    hist(x, xaxt="n")
    axis(1, at=3:9*100, labels=sprintf("%.0e", 3:9*100))
    example1.png
  • ありがとうございました.自動的な処理まだはまだできていませんが,できたところまで書いておきます.
    axis関数のパラメータlabelsにどういうラベルを描くかという命令を送ればいいのですね.
    set.seed(123)
    x <- rnorm(1000)*100+600
    hist(x,xaxt=F)
    hoge<-3:9*100
    a<-sub("+","",substr(sprintf("%1.1e",hoge),nchar(sprintf("%1.1e",hoge))-3,
      nchar(sprintf("%1.1e",hoge))),extended=F)#指数部分の取り出し
    b<-sub("+","",substr(sprintf("%1.1e",hoge),1,nchar(sprintf("%1.1e",hoge))-5),
      extended=F)#小数部分の取り出し
    axis(side=1,at=hoge,labels=expression(a^b))
    とでもすればできるのかな,とおもったのですが,これだと最後の行のexpressionの引数のaとbは展開されないのでエラーになってしまうようで,これをなんとか解決しないといけないようです. -- fjok? 2006-09-11 (月) 14:33:05
  • ほんとはもっとちゃんとしたやり方があった気がするけど思い出せない。トリッキーですが、axis(1,hoge,parse(text=paste(a,"^",b)))。あと、hist(x,xaxt="n")が正解。 -- takahashi? 2006-09-11 (月) 15:42:39
  • そんなに難しく考えないでいいのですよ。pretty 関数と sprintf 関数を調べましょう。 -- 2006-09-11 (月) 16:00:25
    set.seed(123)
    x <- rnorm(1000)*10+100
    hist(x,xaxt="n")
    at <- pretty(x)
    labels <- sprintf("%.2e", at) # 有効桁数は自分で(これも引数にするとか)
    axis(side=1,at=at,labels=labels)
    example2.png
  • あとどうでもいいんですが、実数部、指数部の取り出しは、E<-floor(logb(hoge,10));R<-hoge/10^Eでいいと思いますよ。 -- takahashi? 2006-09-11 (月) 16:01:05
  • みなさんありがとうございました.あとは目盛りをmin(), max()等を使って作れば完全に自動で処理できそうですね. -- fjok? 2006-09-11 (月) 21:12:23
  • > あとは目盛りをmin(), max()等を使って作れば
    その必要がないというのが分からないですか? -- 2006-09-11 (月) 21:54:54
  • 以下のようにやってみました.10のべき乗についてはVoltage(V)の括弧の中に入れてしいまいました. -- fjok? 2006-09-12 (火) 01:02:18
    t<-1:100                                   # ある実験データ
    signal <- 1.0e-3*(sin(t/5)+0.1*rnorm(100)) # ある実験データ
    E<-floor(logb(pretty(signal)[length(pretty(signal))],10)) # y軸の目盛りから指数部分を抽出
    plot(signal,type="l",axes=F,xlab="Time(nsec)",
         ylab=parse(text=paste("Voltage(",
         parse(text=paste(10,"^",E)),")")))
    box()
    par(las=1)
    axis(1)
    axis(2,at=pretty(signal),labels=sprintf("%.1f",
       pretty(signal)/10^E*10^floor(logb(pretty(signal)
       [length(pretty(signal))],10))/10^E))
    
  • 画像の質にずいぶん差があるのですね。 -- 2006-09-12 (火) 23:37:13
    alt.png

大量の図のplot、mflowやmfcolは使えない?

ぼう? (2006-09-09 (土) 11:43:42)

おせわになります。
50枚くらいの多数の図をfor文で繰り返したものを、ひとつの画面に表示したいですが、mflow()やmfcol()では「以下にエラーplot.new() : 図の余白が大きすぎます」とでてしまって、表示できません。
pairs()で表示される図のように多数の図を表示するできる関数をどなたかご存知ないでしょうか?

  • 画面に表示しているのですか?グラフィックウインドウが小さすぎるのでしょう。1024×768のモニター画面一杯に広げたら,8行8列でもちゃんと描けましたよ。 -- 2006-09-09 (土) 12:28:07
    > par(mfrow=c(8,8))
    > for (i in 1:60) hist(rnorm(1000)) 
  • ありがとうございます。基礎的なことだったのかも?ですが、表示と出力では違うのですね。さらにmaiを使って余白範囲を狭めたところ、こちらでも1024*768モニタにて70枚くらい余裕で出力できました。 -- ぼう? 2006-09-09 (土) 12:45:47

R CMD BATCHのオプションについて

山嵐? (2006-09-08 (金) 10:23:37)

Rをバッチで実行しようと思っています。
その際、「R CMD BATCH」 + Rコマンドで実行するというのは
理解しているのですが、
オプションとして vanilla --no-save といったものを付けられるように
聞いております。実際にそれらオプションを理解したいとおもい
付けられるオプションの種類と意味をGoogleやR CMD BATCH -h などで調べたのですが、全くわかりませんでした。
R CMD BATCHのオプションについてご存じの方いらっしゃいましたら、オプション詳細についての調べ方やそういった情報が載っているサイトなどご教授いただければと思います。
お手数をおかけします。

  • サイト内検索はしましたか?Windows版RのFAQの2.10に,「詳細については, Rcmd BATCH --help をご利用ください」と書かれてあります。 -- 2006-09-08 (金) 10:43:59
  • R CMD BATCHのオプションという解釈ではなくてR CMDのオプションという解釈が正しいかと。CRANのManualにてIntroductionのInvoking Rあたりにまとまっています。 -- 2006-09-08 (金) 10:45:06
  • http://cran.md.tsukuba.ac.jp/doc/manuals/R-intro.html#Invoking-R -- 山嵐? 2006-09-08 (金) 11:02:16
  • みなさま、ご教授ありがとうございました。疑問がすべて解決いたしました。大変感謝しておりますm(_ _)m -- 山嵐? 2006-09-08 (金) 11:03:05

リストの逆数の和の求め方

inverse? (2006-09-07 (木) 18:25:49)

皆さんにお聞きしたいことがあります。以下のようなリストを作成し、

a <- as.list(NULL)
a <- c(list(1:3), list(4:9))

a1? と a2? ごとに逆数の和を求めたいのですが、どうすればよいでしょうか?
sapplyを使えばいいとは思ったのですが、逆数和を求める関数がないのでどうしようかと思っています。
自分で関数を定義するしかないんでしょうか?
RはVersion 2.3.1、Windows XPを使用しています。
ご回答よろしくお願いします。

  • 以下のようにすれば良いわけですが、もしかして、実際はリスト成分の数がもっと多いのでしょうか。 -- 2006-09-07 (木) 18:30:43
    > c(sum(1/a[[1]]),sum(1/a[[2]]))
    [1] 1.833333 0.995635
  • lapply(lapply(a, function(x) 1/x),sum) -- 2006-09-07 (木) 18:34:13
  • mapply(function(x) sum(1/x), a) こっちの方がいいか。。。 -- 2006-09-07 (木) 18:39:57
  • みなさんご回答ありがとうございます。実際のリスト成分は500ぐらいありますので、mapply()を使えるとはうれしいです。function()にこんな使い方があるとはとても勉強になりました。 -- inverse? 2006-09-07 (木) 19:58:50

MacOS XでのRmapの使用

haruo? (2006-09-07 (木) 15:02:59)

Mac OS X 10.4.6およびR for Mac OS X Version 2.3.1 (2006-06-01)の環境で、Rmapを使用することは可能でしょうか?RmapのWin32のバイナリはhttp://spatial.nhh.no/R/Devel/から入手できますが、Mac用バイナリ・ファイルは入手可能でしょうか?アドバイスを賜りたくよろしくお願い申し上げます。 参考文献:RjpWiki内の「Rmapを使った地図表示」「なんでも掲示板アーカイブ(5) / Rmap」「初級Q&A アーカイブ(3) / Rmapのインストールについて」

> install.packages("Rmap")
download.packages(pkgs, destdir = tmpd, available = available,  中で警告がありました
	 no package 'Rmap' at the repositories
> sessionInfo()
Version 2.3.1 (2006-06-01) 
i386-apple-darwin8.6.1
  • バイナリはありませんが、下の「パッケージrimageのインストール」と同様、ソースパッケージ(.tar.gz)をとってきて展開して、「このコンピュータ上のパッケージディレクトリ」で指定すればインストールできます。ただし、rimageの場合と同様、あらかじめターミナルから外部ライブラリ(gdal, proj, shapelib)をコンパイルしてインストールしておく必要があります。外部ライブラリのコンパイルをやりたくない場合は、MacOSX版でrgdalをインストールにあるようにDarwinPorts?を使って外部ライブラリをバイナリでインストールします。そうした場合には、RmapのインストールはR CMD INSTALL コマンドを用いてターミナルから行う必要があります。 -- 2006-09-07 (木) 16:14:47

異なる区切り幅のヒストグラムについて

crz1? (2006-09-06 (水) 22:48:58)

異なる区切り幅のヒストグラムを描いたところ、
それぞれのバーの横幅が不揃い(指定した区切り間相当の幅)になってしまいます。
このバーの横幅を一定にするにはどうすればいいのでしょうか?

  • 区切り幅というのは,階級幅のことですか?区切り幅を不等間隔にしたのに,それをヒストグラムに描いたときに横幅が不揃いだと文句を言うとは,訳が分かりませんが。横軸の目盛り間隔を変えることになりますよ。どうしてもやりたいなら,各階級の度数を求めて,hist 関数ではなくて barplot 関数を使えば?本質的に barplot 関数は棒グラフなので,棒の横幅は全部同じになりますが。しかしそれはもはやヒストグラムではなく棒グラフ。なんか,当たり前のことをああだこうだ書いたが,質問がそおゆうことなんで。 -- 2006-09-06 (水) 23:13:37
  • 分かりにくくてすみません。ご指摘の通り目盛り間隔を変えたいのです。 -- crz1? 2006-09-06 (水) 23:18:40
  • そうします。ヒストグラムに固執しすぎていました。ありがとうございます。 -- crz1? 2006-09-06 (水) 23:23:00

クラスター図を横向きで描画したい

sally? (2006-08-30 (水) 11:01:13)

こんにちは。いつもWikiにお世話になっています。
さて、今回クラスタ図を描こうとしているのですが、描画を横向きにしたいのです。どうしてもやり方がわかりません。

hc <- hclust(dist(USArrests)^2, "cen")
plot(hc)

これで描画される樹形図を縦ではなく、項目の名前が長いので、横で表示させたいのです。
ちょうど右90度回転させたような図を書くためにはどうしたらよいでしょうか?

ご教授お願い致します。

  • ごめんなさい、今書き方を間違えたので修正しようとしたら、どなたかが修正してくださっていました。ありがとうございます。 -- sally? 2006-08-30 (水) 11:19:21
  • main タイトルだけ text 関数で90度回転させて描き,後で画像ソフトを使って90度回転させればいかがでしょう。 -- 2006-08-30 (水) 13:09:07
  • アドバイスありがとうございます。自動で描画したものをHTMLに出力させたいので、何とか横向きになるようにRの中でやりたいのですが・・。 -- sally? 2006-08-30 (水) 13:24:02
  • 少し,トリッキーですが. 例えばepsで出力したものを, TeXで使ったり, htmlで使う場合, -- なかま 2006-08-30 (水) 14:45:36
    postscript("foo.ps", width = 5, height = 4,
               horizontal = FALSE, onefile = FALSE, paper = "special")
    plot(1:10)
    dev.off()
    等とすると, epsファイルができます, これを例えばスクリプトなんかで,
    # 右回転90度
    @@ -2,10 +2,11 @@
     %%Title: R Graphics Output
     %%Creator: R Software
     %%Pages: (atend)
    -%%BoundingBox: 0 0 360 288
    +%%BoundingBox: 0 0 288 360
     %%EndComments
     %%BeginProlog
    -/bp  { gs gs } def
    +%/bp  { gs gs } def
    +/bp  { gs 0 360.0 translate -90 rotate gs } def
     % begin .ps.prolog
     /gs  { gsave } def
     /gr  { grestore } def
    とか
    # 左回転90度
    @@ -2,10 +2,11 @@
     %%Title: R Graphics Output
     %%Creator: R Software
     %%Pages: (atend)
    -%%BoundingBox: 0 0 360 288
    +%%BoundingBox: 0 0 288 360
     %%EndComments
     %%BeginProlog
    -/bp  { gs gs } def
    +%/bp  { gs gs } def
    +/bp  { gs 288.00 0 translate 90 rotate gs } def
     % begin .ps.prolog
     /gs  { gsave } def
     /gr  { grestore } def
    以上, 「その手でやるんだ(c)沖田艦長」. でした. あとは, convertで? もっとも, ImageMagic?が入っていれば, convertだけでも出来ます.
  • これでいかがでしょう?-- 2006-08-30 (水) 22:27:07
    hc <- hclust(dist(USArrests)^2, "cen")
    par(mar = c(4,1,1,7))
    plot(as.dendrogram(hc), horiz = T)
  • ↑すばらしい!と思ったら左90度回転だった… -- 2006-08-31 (木) 11:08:37
  • わぁ^^答えてくださった皆様、ありがとうございます!!!とっても助かりました。 -- sally? 2006-08-31 (木) 18:24:24

パッケージrimageのインストール

haruo? (2006-08-30 (水) 10:22:28)

finkを使ってfftwとlibjpegをインストールした後、以下の3つの方法により、パッケージrimageのインストールを試みましたが失敗しました。アドバイスを賜りたくよろしくお願い申し上げます。

(方法1)関数 install.packages() でパッケージを指定

> install.packages("rimage")
download.packages(pkgs, destdir = tmpd, available = available,  中で警告がありました
	 no package 'rimage' at the repositories

(方法2)関数 install.packages() で圧縮ファイル(rimage_0.5-7.tar.gz)のパスを指定

> install.packages("/Users/haruo/Desktop/rimage_0.5-7.tar.gz", contriburl = NULL)
以下にエラーgzfile(file, "r") : コネクションを開くことができません
追加情報: Warning message:
圧縮されたファイル 'rimage_0.5-7.tar.gz/DESCRIPTION' を開くことができません

(方法3)メニューの [パッケージとデータ] から [パッケージインストーラ] を選択し、CRAN(ソース)の〔一覧を取得〕 し、rimageを選択して [OK] をクリック

(中略)
checking fftw.h usability... no
checking fftw.h presence... no
checking for fftw.h... no
** Removing '/Library/Frameworks/R.framework/Versions/2.3/Resources/library/rimage'
configure: error: Sorry, can't find fftw header
ERROR: configuration failed for package 'rimage'

<動作環境> Mac OS X 10.4.6; Darwin 8.6.1。R for Mac OS X Version 2.3.1 (2006-06-01)

> sessionInfo()
Version 2.3.1 (2006-06-01)
i386-apple-darwin8.6.1
attached base packages:
[1] "methods"   "stats"     "graphics"  "grDevices" "utils"     "datasets"  "base"    
  • 方法1は、Mac用のバイナリがCRANにないので不可能です。方法2は、圧縮ファイルのパスではなく、展開後のフォルダ自体のパスを指定するべきです。ですが、結果は方法3と同じになるでしょう。 -- 2006-08-30 (水) 12:24:47
  • で、方法3が失敗したのはfinkを使っているためです。現在のMacOS X版Rはfinkの存在を前提としていませんから、パッケージのコンパイルに必要なライブラリはfinkではなく、自分でインストールする必要があります。 -- 2006-08-30 (水) 12:26:01
  • 自分でインストールする方法は、慣れていれば普通にコンパイルしてインストールするだけですが、慣れていないと難しいかもしれません。Darwin Portsがよいかな? -- 2006-08-30 (水) 12:28:01
  • Darwin Portsも結局/swが/optになるだけなので同じですね... -- 2006-08-30 (水) 12:39:19
  • つまり、fftwとかjpeglibといった外部ライブラリは、finkやDarwinPorts?といったパッケージではMacOS X本体への影響を最小限にするために/swとか/optといった、MacOS Xでは使われないフォルダにインストールされてしまうのですが、Rはそれらのフォルダを見ないので、そこにfftwやjpeglibをインストールしてもパッケージのインストールには使われません。Rがみるのは、MacOS Xでも使われる/usr/local/ などです。外部ライブラリをここにインストールしてあげれば問題なくパッケージをインストールできます。 -- 2006-08-30 (水) 13:29:36
  • fftwを/usr/local/にインストールするには、FFTW Installation on MacOSにかいてある通り、FFTW Download Pageからfftw-2.1.5.tar.gzをダウンロードし、解凍してできたフォルダ内でターミナルから "./configure --enable-threads" 次に "make" 最後に"sudo make install" すればよいです。 -- 2006-08-30 (水) 13:29:52
  • libjpegを/usr/local/ にインストールするには、Independent JPEG Group から jpegsrc.v6b.tar.gz をダウンロードし、解凍してできたフォルダ内でターミナルから "./configure" 次に "make" 最後に "sudo make install-lib install-headers" を実行すればよいです。 -- 2006-08-30 (水) 13:30:08
  • 最後にR のパッケージインストーラで、「このコンピュータ上のパッケージディレクトリ」を選択し、「Install」を選んで、でてきたファイル選択ダイアログでは、CRANからダウンロードしたrimage_0.5-7.tar.gzを解凍してできたフォルダの中の適当なフォルダ(srcとかdataとかmanとかでOk)を選択して「開く」をクリックすれば、インストールされるはずです。 -- 2006-08-30 (水) 13:30:18
  • ありがとうございます。正常にインストールされました。 -- haruo? 2006-08-31 (木) 01:15:07
  • もう解決して要るみたいですが、finkでもDarwinPorts?でも、R CMD INSTALLで --configure-args= を使ってやればOKです。参考:MacOSX版でrgdalをインストール 。パッケージシステムを使わずにむやみにtarボールから/usr/localに入れてしまうのはあまり賛成できません。 -- 2006-09-01 (金) 11:06:50
  • 確かにバージョン管理が難しくなってしまうのでオススメできないんですが、R本体もインストール時にtcl/tkなどを/usr/local/にインストールするようにしています。RをfinkやDarwinPorts?に依存しないように、かつR CMDを使わずにできるだけGUIでインストールができるようにするためにはやむを得ないかと思います。本来は、tcl/tkと同様にFFTWやlibjpegについても/usr/local/にインストールするMacOS Xのインストーラが用意されていれば、コマンドラインを全く使わずにすむのですが... -- 2006-09-01 (金) 12:47:32
  • 便乗質問ですが、finkやDarwinPorts?に依存するパッケージを、finkやDarwinPorts?の外部ライブラリを用いてインストールした場合、そのパッケージはDockからRアイコンをクリックする形で起動したRでも使えるのでしょうか?依存するダイナミックリンクライブラリへのパスをRが知らないといけないので、環境変数を設定しないとだめなのでは?と思っているのですが... -- 2006-09-07 (木) 16:18:51

SPSSのvariable labelsを取り出す方法

kam? (2006-08-29 (火) 14:37:57)

はじめまして。SPSSのデータをRで読み込んだ際のラベルの取り出しについて質問させてください。
例えば以下のようなSPSSのsavデータ「test.sav」があった場合

Q1 Q2
1 2
3 3
2 2
4 1
5 2

Q1にSPSSで"国語"、Q2に"算数"とラベル付けしたとして、それをRでdata1として以下のようにして読み込みました。

library(foreign)
data1 <-read.spss("test.sav",use.value.labels=F,to.data.frame = T)

ここで、SPSSで使っていた変数へのラベルをRでも使いたいと思い、取り出す方法を調べていたのですが、分からず困っています。

> attributes(data1)
$names
[1] "Q1" "Q2"
$row.names
[1] "1" "2" "3" "4" "5"
$class
[1] "data.frame"
$variable.labels
    Q1     Q2 
"国語" "英語" 

attributes()を使うと、上のようにvariable.labelsとして読み込んでいるようなのですが、variable.labels(data1)としてもそのようなfunctionは無いというエラーメッセージがでるだけで、"国語"という形で取り出すことができません。
大変初歩的というのは自覚しているのですが、variable.labels、変数名、spssなどのキーワードで、こちらのサイト、リンク先の初心者マニュアル、及びgoogleで調べていても分からず、助言を求めに来た次第です。
使用環境は、Win-XP,RはVersion 2.3.0 です。
実際に扱うSPSSのsavデータは変数が170あり、できればSPSSのラベルをそのままRでも利用したいという状況です。
質問の作法については目を通したのですが、不適切でしたらご指摘ください。よろしくお願いいたします。

  • names(data1) <- data1$variable.labels ではいけませんか。variable.labels という成分は初耳ですが read.spss が付け加える特別なものでしょうか。 -- 2006-08-29 (火) 17:14:37
  • names(data1) <- attr(data1,"variable.labels") としないとうまくいかないかもしれません。 -- 2006-08-29 (火) 17:26:02
  • なるほど、回答ありがとうございます。ちょっと個人的な都合で、実際に確認するのは月曜ということになってしまうのですが、確認次第、また報告いたします。 -- kam? 2006-08-30 (水) 21:12:05
  • names(data1) <- attr(data1,"variable.labels")で、うまくいきました!どうもありがとうございました。names(data1) <- data1$variable.labelsだと、NULLになってしまうようです。 -- kam? 2006-09-04 (月) 10:16:22

日本語のナカグロ(・)の扱い

ニョロニョロ? (2006-08-24 (木) 15:03:05)

おかしなことが生じているので,問題を再現できる最小のデータと共におうかがいします。
test.dat というファイル名で,コードは UTF-8, 改行コードは LF,タブ区切りです。
中身は以下の通り(3行,2列)

変数・あ	変数・ん
1	2
3	6

このデータファイルを読んで,表示したところは何ともなさそうに見えるのですが,

> x <- read.table("test.dat", header=TRUE)
> x
  変数.あ 変数.ん
1       1       2
2       3       6

変数名が変なんです。どうも,ナカグロの変換にともなって,何か変なことが起こっているようです。

> colnames(x)
[1] "変数.あ?0<82>" "変数.ん?0<93>"

コードを EUC-JP にしても,同じようになっています。

  • src/main/character.c の do_makenames のアロケートの問題ですね. 他に影響無いかチェックしてます... orz -- なかま 2006-08-24 (木) 20:52:06
  • ありがとうございました。何も,ナカグロを使わなくても良いので,不便はありません。 -- 2006-08-24 (木) 22:38:39
  • うーん, 我乍ら... 弩級に汚いコードに...(;_;) -- なかま 2006-08-24 (木) 23:12:03
  • R-2.4.0には変なおまけは付かなくなります. -- なかま 2006-08-25 (金) 22:41:25

factor関数の使い方

R初心者? (2006-08-24 (木) 10:48:40)

factor関数でいくつかのカテゴリ変数を1つに纏める方法なんてのはあるのでしょうか?
例えば数量化I類で分析した結果、いくつかのカテゴリ変数が有意ではない、と認められた場合、いくつかのカテゴリ変数をひとつにまとめてしまいたい場合、Microsoft Excelなんかでデータをエディットしてまた再読み込み・・・・なんて手順を踏んでいたんですが、どうせならR上で修正してしまいたいのです。
是非お知恵をお貸し下さい。

  • ? factor は見てみましたか?labels 引数が答えでしょう。 -- 2006-08-24 (木) 11:09:42
    > x <- c(5, 1, 4, 3, 3, 2, 5, 3, 3, 4)
    > factor(x, labels=c(1,1,2,3,4))
     [1] 4 1 3 2 2 1 4 2 2 3
    Levels: 1 1 2 3 4
    > factor(x, labels=c(1,1,1,2,2))
     [1] 2 1 2 1 1 1 2 1 1 2
    Levels: 1 1 1 2 2
  • 趣旨を今一理解できていませんが、いくつかのカテゴリーを統合したいだけのように聞こえますが。例えば次のようなやりかたでは間に合いませんか。-- 2006-08-24 (木) 11:23:21
    > x <- c(1,2,2,3,4,5,5,5,6)
    > factor(x)
    [1] 1 2 2 3 4 5 5 5 6
    Levels: 1 2 3 4 5 6
    > x[x==6] <- 5
    > factor(x)
    [1] 1 2 2 3 4 5 5 5 5
    Levels: 1 2 3 4 5
  • 上の例だと,2を1に統合して,その後を順に詰め合わすという場合に記述量が多くなってしまうし間違いも犯しやすいということがあるかもしれません。 -- 2006-08-24 (木) 11:33:00
  • 早速のレスありがとうございます。・・・後者の方ですね。なるほど。単に要素を置き換えれば良かったんですね。ありがとう御座います。なんか難しく考えていたようです。(思考回路が今だExcelに支配されている!!!Excelでダミー変数を作るのは面倒くさいので・・・・・・汗。) -- R初心者? 2006-08-24 (木) 11:33:41
  • ところで、有意じゃないからといってカテゴリをまとめてはいけません、というのはここでは言わない方がいいのでしょうか? -- 不思議さん? 2006-08-24 (木) 14:18:01
  • その他というカテゴリーに入れるんだと思えば,有意でないカテゴリーをまとめるのは問題ないのでは? -- 2006-08-24 (木) 14:55:48
  • 再解析と言う方針で行こうかと思っていたのですが・・・・。青木先生に質問してみます。 -- R初心者? 2006-08-24 (木) 14:59:49
  • いろんなまとめ方を試して、再解析、再々解析、・・・、再々々・・・解析して、有意になったら発表しよう!じゃ、だめじゃないですか?有意になる、ならないに関わらず、モデルを決めたらそれでいかなきゃいけません。それが許されるなら、カテゴリーをまとめるよりも、有意になるようにデータをけずっちゃった方が手っ取り早いですよ。 -- 摩訶不思議? 2006-08-24 (木) 22:53:17
  • 数量化理論(カテゴリ変数)の場合の作法はちょいと知りませんが、普通の回帰では交互作用を入れて試すところですね。個々は有意でなくても他と組み合わせる(交互作用項)と有意になることもある。 -- 2006-08-24 (木) 23:52:35
  • 数量化の場合は自動化がなかなか難しいが,ステップワイズ変数選択の重回帰分析なんてのは,有効な独立変数を自動的に選択してくれる訳ですから,そんなに非難されるようなことでもないと思います。そのようにして得られた結果をそのまま使うか,ユーザの主体性を反映したモデルに修正して使うかという余地はあるでしょう。

添付ファイル: filemplot1.png 1799件 [詳細] filesupsub.png 1778件 [詳細] fileakima-example2.png 1877件 [詳細] filemulti.png 1866件 [詳細] filecontour-example.png 2072件 [詳細] filemplot3.png 1937件 [詳細] fileuc.png 1961件 [詳細] filedetarame.png 1116件 [詳細] fileclust.png 1975件 [詳細] filemplot2.png 2127件 [詳細] filetest.png 862件 [詳細] fileakima.interp.png 2117件 [詳細] fileexample1.png 1832件 [詳細] filedemo6.png 1973件 [詳細] fileexample2.png 1823件 [詳細] filealt.png 1778件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1569d)