#author("2021-09-16T13:32:39+09:00","","")
COLOR(green){SIZE(20){初心者のための R および RjpWiki に関する質問コーナー}}
COLOR(red){SIZE(16){以下の 4 つの「●項目」のどれかをクリック！！}}~
----
COLOR(blue){● [[新規投稿欄:>#sinki]]} --- その前に，「[[投稿における注意事項]]」を読んでから！&br;
COLOR(blue){● [[このページの目次>#mokuji]]} --- 質問への回答・コメントの参照&br;
COLOR(blue){● [[最新のスレッド>#saisin]]} --- 最も最近に書き込まれた「親」記事（最新のコメントが付いた記事ではありません）&br;
COLOR(blue){● [[Ｑ＆Ａ (初級者コース) 過去の目次]]} --- 今までの 17 個の書庫の参照&br;
----
&aname(mokuji);COLOR(green){●● 目次 ●●}
参照は，個々の項目をクリック（質問への回答・コメントの参照も，個々の項目をクリック）~
#contents
----
&aname(sinki);
COLOR(green){●● 新規投稿 ●●}
まずは，「[[投稿における注意事項]]」を読んで，それに従ってください。~
COLOR(red){これに従わない場合には，回答が得られにくい場合があります。}
#article
**パッケージのインストールができない [#g0fbfb35]
>[[もいくんとぅ]] (2021-08-21 (土) 21:38:06)~
~
Rコマンダー、EZR、その他の必要なパッケージがインストールできないので、インストールできる方法をお教えして頂ければ有り難いです。~
~
使用OSは、Mac Book Pro 15inch Mid 2009 Mac 、Mac OS X Lion 10.7.5~
Rのバージョンは、R-3.2.1-snowleopard.pkg~
~
①CRANのホームページからインストーラーをダウンロードしインストール済みです。~
②XQuartzプロジェクトからX11システムをインストール済みです。 ~
③Rというアプリケーションを開き、install.packages("RcmdrPlugin.EZR", dependencies=TRUE)と入力、returnキーを押すが、下記エラーメッセージが表示されてしまいます。~
~
【警告: リポジトリー https://ism//cran.ism.ac.jp/src/contrib に対する牽引にアクセスできません。警告: リポジトリーhttps://ism//cran.ism.jp/bin/ Macosx/contrib/3.2に対する牽引にアクセスできません。警告メッセージ：package ’Rcmdrplugin.EZR'is Not available(for R version3.2.1】~
~
CRANミラーをJapan(Tokyo),Japan(Yonezawa)に変更して試してみましたが、上記同様のメッセージが表示されてしまします。~
~
解決方法をお願い致します。~
//
#comment
**datasetsがありません [#g8a6f17a]
>[[Rビギナー]] (2021-08-20 (金) 01:35:36)~
~
コマンダーにおいて、パッケージ内のデーター→アタッチされたパッケージからデータセットを読み込む を行っても、CarDataとsandwichしかなく、datasetsが表示されません。どうしたらいいでしょうか、お教え願います。~
//
- 初歩的な問題でした。テキスト付属のCDROMからアクセスできました。 -- [[Rビギナー]] &new{2021-09-16 (木) 13:32:39};
#comment
**多変量正規分布の積分について [#ide39545]
>[[よしの]] (2021-07-07 (水) 15:59:48)~
~
多変量正規分布の積分を行いたいのですが，使用できる関数がわかりません．~
~
正規分布（1変量）の場合，~
f_n <- function(x) dnorm(x, myu, sigma)
integrate(f_n, -Inf, Inf)
で計算でき，このような例はいくつかのサイトに示されていました．~
この1変量の場合と同様に，多変量でも積分を行いたいです．~
~
実際に試したいのは，3変量の場合で~
f_mv <- function(x1, x2, x3) dmvnorm(cbind(x1, x2, x3), mean=Mu1, sigma=Si1)~
と，関数を定義して~
これをx1, x2, x3について積分したいのですが，integrateを用いた場合の設定がわかりません．~
~
他の方法でも良いので，なにか情報いただければと思います．~
よろしくお願いします．~
~
RStudio ver 1.3.1093~
//
- cubatureパッケージを使うのが簡単そうですが、数値積分の用途でRを使ったことがないので自信はないですね……。
library(cubature)
library(mvtnorm)
Mu <- c(0, 1, 2)
Sigma <- diag(length(Mu))
f_mv <- function(x) { dmvnorm(x, mean = Mu, sigma = Sigma) }
hcubature(f_mv, c(-3, -2, -1), c(1, 2, 3))
パッケージのヘルプとVignetteを読んで使ってみてください。（わたしは斜め読みしただけなので。）どの程度の精度が必要かにもよるかもしれません。hcubatureはC言語で書かれた関数実体へのラッパーとのことなので、必要があればそちらのコードを確認するべきかと思います。 -- [[しろうと]] &new{2021-07-08 (木) 08:12:13};
- 情報ありがとうございます．いくつかの分布に対して処理を試してみました．多変量関分布のパラメタをデータから推定したものが，きちんと積分して１になるかをチェックしたいというのが目的でしたので，その用途としては十分です．ありがとうございました． -- [[よしの]] &new{2021-07-09 (金) 12:34:07};
#comment
**プロビット分析の仮定 [#sf6abc04]
>[[まるせん]] (2021-06-29 (火) 18:13:37)~
~
被説明変数に0,1の指標、説明変数に４項目程度の連続データでプロビット分析を行っています。データ数は400程度となります。~
Rでglmを使って無事モデルは構築できたのですが、右辺の被説明変数（標準正規分布の引数の部分）が標準正規分布に従っていません。（分布の形も正規分布ではなく、平均、分散もそれぞれ０，１ではありません。）~
そもそも右辺の被説明変数には標準正規分布の仮定を置いていないと思っていますが、私の考えがあっているか教えてください。前提は誤差項の正規性ですよね。~
これ読んどけ、というレポートを教えていただけるだけでもうれしいです。~
~
お願いします。~
//
- 右辺の被説明変数（？）云々のあたりがよく分かりませんが、プロビット分析の理屈自体としては「連続的な変数（群）からカテゴリカルな反応を説明するとき、両者をつなぐS字状の関数関係として正規分布の累積分布を使ってる」というだけではないですか？ -- &new{2021-06-30 (水) 07:40:30};
#comment
**両対数目盛での回帰直線 [#xe0feb6d]
>[[つぐ太郎]] (2021-06-19 (土) 15:47:10)~
~
年齢と死亡率が普通メモリグラフでは曲線ですが、両対数目盛グラフでは直線になることを示したいのですが、その時の回帰直線がうまく書けません。試したプログラムは下記です。~
age <- c(20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85)~
mortality <- c(2, 5, 10, 20, 25, 37, 60, 80, 133, 190, 280, 410, 420, 380)~
plot(mortality ~ age, log = "xy")~
If <- glm(mortality ~ 1 + age)~
cf <- If$coefficients~
curve(cf[1] + cf[2] * x, from = 20, to = 80, col = "blue", cex=1.5, lwd=2, add = T)~
//
- まぁ、出来合いの便利な関数や機能でやる方法もあるでしょうが。。。 -- &new{2021-06-20 (日) 13:37:58};~
「両対数グラフで直線になる」ということは「ageのlogとmortalityのlogが線形関係にある」ということですよね？なので、
rslt <- lm(log10(mortality) ~ log10(age))
summary(rslt)
としてやればよいわけです。
cf <- coef(rslt)
par(mfrow = c(1, 2))
for (l in c("", "xy")) {
plot(mortality ~ age, log = l)
lines(age, 10 ^ (cf[1] + cf[2] * log10(age)))
}
- 明快にありがとうございます -- [[つぐ太郎]] &new{2021-06-22 (火) 23:36:34};
#comment
**Rtoolsのパスが開けません [#x9a472fd]
>[[KF]] (2021-04-10 (土) 20:44:50)~
~
以下を参考にrstanをダウンロードしようとしています。~
http://www.eeso.ges.kyoto-u.ac.jp/emm/materials/bayesian/stan_install~
rstan自体はダウンロード出来たのですが、rstanを使用するためにインストールしたrtoolsのパスを開こうとした際以下のエラーが出ました。~
~
コード~
writeLines('PATH="${RTOOLS40_HOME}\\mingw64\\bin;${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")~
警告メッセージ: ~
file(con, "w") でエラー: コネクションを開くことができません
追加情報: 警告メッセージ:
file(con, "w") で:
ファイル 'C:/Users/*****/OneDrive/??????/.Renviron' を開くことができません: Invalid argument (*****＝ユーザー名)
~
上のリンクでは「RでHomeとして使っているディレクトリ（デフォでgetwd()で表示されるパス）に，.Renviron というファイルが出来ているはず．」と書かれていましたが、そのようなファイルは作成されていないようでした。そのためエラーメッセージが発生したとみられます。~
~
Rtoolsの設定がうまくいっていないせいか、rstanを起動しても以下のエラーメッセージが出ます。~
プロシージャエントリポイント EXTPTR_PTRがダイナミックリンクライブラリC:\Users\*****\R\win-library\4.0\rstan\libs\x64\rstan.dllから見つかりませんでした（*****=ユーザー名）~
~
操作環境は以下の通りです。~
R version 4.0.0 (2020-04-24)~
Platform: x86_64-w64-mingw32/x64 (64-bit)~
Running under: Windows 10 x64 (build 19042)~
~
解決方法をご教示いただけると幸いです。~
よろしくお願いいたします。~
//
- 直接の解決策でなくすみません。~
知り合いのPCをセットアップした際、「起動ディレクトリがOneDriveの配下だと.Rprofileファイルが読み込まれない」という問題にでくわしました。OneDriveまわりのWindowsの仕様は、R関連では邪魔をすることが色々あるようです。（検索してみてください。）管理者権限で起動していれば問題にならない場合もあるとか。~
.Rprofile問題の際は、OneDriveでないユーザホーム直下にファイル実体を置く、みたいな方法で対処した気がしますが、あまり定かではありません。OneDriveが無効な場所で.Renvironを置ける場所は他にないですか？あるいは古典的で泥臭いですが【システムのプロパティ】→【環境変数】でPATHをR外で設定しておくとか。 -- &new{2021-04-11 (日) 14:17:42};
- コメント頂きありがとうございます。OneDriveではない場所に新規のプロジェクトを作成しました。writeLines('PATH="${RTOOLS40_HOME}\\mingw64\\bin;${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")と入力したところ、前回と全く同じエラーが出てしまいました。file(con, "w") でエラー: コネクションを開くことができません 追加情報: 警告メッセージ: file(con, "w") で: ファイル 'C:/Users/*****/OneDrive/??????/.Renviron' を開くことができません: Invalid argument (*****＝ユーザー名) 何が問題なのか全く分からず途方に暮れております。どうぞよろしくお願い致します。 -- [[KF]] &new{2021-04-12 (月) 18:21:14};
- R, Rstudio, Rtoolすべてアンインストールし、再インストールしたらrstanが正常に使えるようになりました。保存場所をOne Driveではないフォルダにしたことが効いたのかもしれません。ありがとうございました。 -- [[KF]] &new{2021-04-12 (月) 20:37:54};
#comment
**read.table()で部分的に文字化けしてしまう [#jc87185b]
>[[Gorou]] (2021-04-07 (水) 20:41:39)~
~
質問です。"UTF-8"ファイルの読み込み時に文字化けしてしまい困っています。~
（Wikiの記法に不慣れな点はお許し願います。)~
~
read.table()での文字化けは非常にありがちな質問ですが、~
~
"read.table()"で読み込んだデータが"「部分的に」文字化け"してしまい、~
~
困っています。~
~
具体的には次のような症状です。~
以下のデータファイルを読み込みたいのですが、~
[data.tsv] (TSV形式、UTF-8、BOMなし）
Age Group Class DBH_mean DBH_sd
31 すべて 全立木 17.9 4.34
31 すべて 生立木 17.9 4.34
31 すべて 残存木 17.9 4.34
36 すべて 全立木 20.0 4.92
36 すべて 生立木 20.0 4.92
36 すべて 残存木 20.0 4.92
41 すべて 全立木 21.9 5.59
41 すべて 生立木 21.9 5.56
41 すべて 残存木 21.9 5.56
~
>df <- read.table("test.tsv",sep="\t",header=TRUE,skip=0, fileEncoding="UTF-8")
>head(df)
~
とすると、~
~
Age Group Class DBH_mean DBH_sd X
1 31 すべて \u5168\u7acb\u6728 17.9 4.34 NA
2 31 すべて \u751f\u7acb\u6728 17.9 4.34 NA
3 31 すべて \u6b8b\u5b58\u6728 17.9 4.34 NA
4 36 すべて \u5168\u7acb\u6728 20.0 4.92 NA
5 36 すべて \u751f\u7acb\u6728 20.0 4.92 NA
6 36 すべて \u6b8b\u5b58\u6728 20.0 4.92 NA
("\"はバックスラッシュ）
~
となり、「Group」のカラムは文字化けしないのに、「Class」のカラムは文字化けしてしまいます。完全に文字化けしているのならば文字コードの間違いだと納得できるのですが、~
このような中途半端な文字化けは、いろいろググって見たのですがさっぱり分かりません。~
~
また、DBH_sd の右に "X"のカラムができてしまうのも不可解です。DBH_sdのカラム、4.34, 4.34, 4.34, ... は数値の直後で改行されています（つまり元のファイルにはヘッダだけのカラムは無い）。~
~
TSVをCSVにしたり、ファイル形式をShift-JISにして fileEncoding='Shift_JIS'としてみても、結果はほぼ同じで途方に暮れております。~
おまじないにと,stringAsFactorsをTRUEにしたりFALSEにしても変わりません(変わるとも思えませんが）。~
~
環境はWindows10 64bit Pro、Rのバージョンは64bitの4.0.4、Rstudioは1.4.1106。症状はRGuiでもRStudioでも同じです。~
~
どなたか解決方法をご教示いただけましたら幸いです。何卒よろしくお願いします。~
//
- とりあえず列が増えるのは、単に行末の余計なタブのせいでしょう。質問者さんが貼られたデータ例にも、行末にタブが入っています。タブのうしろに何もデータがないので、read.table時に自動で列がつくられたものの、値はすべてNAになっています。 ~
~
つぎにClassの列の現象ですが……文字化けというよりは、UTF-8のコードに展開されているというカンジでしょうか。数値をUTF-8として読み取ればもとの「全立木」などの文字列が読み取れますので、データとして「化けている」というのともちょっと違う気がします。~
~
なぜそんなことになってるかについては……難しいですねぇ。可能性として、以下のようないろんなパターンが考えられるかと思います。質問者さんが行末のタブに気づいていらっしゃらなかった（数値後即改行だと思われていた）件もあわせると、最後のような「じつはR外の問題」みたいな可能性も無きにしも非ずではないかな、と。（僭越ですが。）その現象が起こるminimal exampleのファイルをご提示いただければ、どこに問題があるかの切り分けが多少楽になるかと思うのですが。
--その列のデータだけ何か理由があって、read.tableがおかしな挙動をした。
--read.tableはできてるんだけど、それを表示するメソッドprint.data.frameがおかしな表示の仕方をした。（たとえばその列だけ特殊なclass属性がついちゃってて、print時にそれが反映されてしまった。）
--じつはもとのデータもUTF-8のコード直書きみたいな変な形式になってるのだけど、（R外でそれを表示して）確認に使ってるソフトが自動的に文字列として表示してくれてるので気づけていない。Rはそれを愚直に読み込んでるだけ。 -- &new{2021-04-08 (木) 13:45:15};
- exampleファイルです -- [[Gorou]] &new{2021-04-09 (金) 13:09:14};
- 済みません、操作ミスです（削除方法を勉強中です…） -- [[Gorou]] &new{2021-04-09 (金) 13:16:44};
-- 列が増えるのは、余計な半角スペースが入っていたためでした。削除したら解決しました。 -- [[Gorou]] &new{2021-04-09 (金) 13:19:13};
-- オリジナルのデータは、RStudio、秀丸エディタ、notepad++のいずれでもUTF-8として認識、表示されています。 -- [[Gorou]] &new{2021-04-09 (金) 13:22:20};
-- (ubuntuで）nkf -g test.tsv とするとUTF-8が帰ってきますし、"cat test.tsv" でも正しく表示されます。（minimal exampleのファイルの提示方法が分かりません。箇条書きですみませんが、よろしくお願いします。） -- [[Gorou]] &new{2021-04-09 (金) 13:27:07};
-- 度々済みません。根本的な解決策ではないのですが、"read.table"の代わりに"tidyverse"を使ったら文字化けせずに読み込めました。しかし、以前作ったデータとスクリプトの組み合わせが使えないままです。引き続き、情報等ありましたらよろしくお願いします。(Wikiの記法、勉強します…。） -- [[Gorou]] &new{2021-04-09 (金) 14:14:50};
- 確認ですが、このページに添付されたtest_utf8.tsvというファイルは、質問者さんの環境では前記のような問題の再現するファイルなのですか？わたしの環境だと、その「文字化け（と呼んでいらっしゃる現象）」なしでちゃんと読み込めるようです。
> df <- read.table("test_utf8.tsv", sep = "\t", header = TRUE, skip = 0, fileEncoding = "UTF-8")
> head(df)
Age Group Class D_mean H_mean N G V
1 31 すべて 1)間伐前 17.9 12.3 1535 40.95 260.0
2 31 すべて 2)間伐後 17.9 12.3 1535 40.95 260.0
3 36 すべて 1)間伐前 20.0 14.9 1535 51.33 393.1
4 36 すべて 2)間伐後 20.0 14.9 1535 51.33 393.1
5 41 すべて 1)間伐前 21.9 16.6 1535 61.40 517.8
6 41 すべて 2)間伐後 21.9 16.6 1530 61.36 517.6
なおこちらは、R version 4.0.4 on Ubuntu 16.04 LTSをbash上で起動しています。 -- &new{2021-04-09 (金) 17:09:52};
- （最初に提示されたデータと列などが違うので、別のデータなのかな……もしunpublishedで外に出たらまずいデータとかだったら、消すなり差し替えるなりなさってくださいね。） -- &new{2021-04-09 (金) 17:18:14};
- 確認ですが、このページに添付されたtest_utf8.tsvというファイルは、質問者さんの環境では前記のような問題の再現するファイルなのですか？
> はい、その通りです。データを確認していただき、ありがとうございます。ちょっと書き換えましたがほぼ同様のデータです。
理由はさっぱり分からないのですが、Rを4.0.4から4.0.5にアップデートしたら文字化けしなくなりました(RStudioもRGuiも）。
違う点としては4.0.5は管理者として、4.0.4は一般ユーザとしてインストールしたことです。
（管理者だとC:\Program Files\R\R-4.0.5に、一般ユーザだとC:\Users\hogehoge\Documents\R\R-4.0.4にインストールされていました。）
それにしても、元のデータは
> すべて 全立木
であるのに「文字化け」状態だと
> すべて \u5168\u7acb\u6728
と部分的に化けてしまうと言うのはさっぱり訳が分かりません。
いずれにしても、これまでのデータやスクリプトが問題なく使えるので安心しました。
助言をくださったみなさん、どうもありがとうございます。
（ちなみにtidyverseのread_csvも、再確認したところ4.0.5にアップデートする前にはなぜか文字化けしてました。こっちも何が何やら。）
- （返信の編集方法はこれで良かったのかな？みなさん、ありがとうございます。） -- [[Gorou]] &new{2021-04-12 (月) 16:55:33};
#comment
**シミュレーションで確率を見積もりたい [#lba761df]
>[[HO0217]] (2020-12-29 (火) 21:28:39)~
~
「赤札と青札がある。確率1/2でどちらかを引く。~
これを100回繰り返し、赤札を55回以上引く確率を求めよ。」~
という問題を解くため、正規分布に近似して確率を求めると、~
15.87%でした。~
(これは私の計算結果で、答えは知りません。あってるとは思いますが、間違っているかもしれません)~
~
~
近似ではなく、二項分布で計算してみようと思い、~
以下のコードで確率を計算すると、約18.41%(出力は0.1841008)で、~
近似とはやや違う結果となりました。~
tmp_nu<-0~
n_nu<-100~
k_nu<-55~
for(i in k_nu:n_nu){~
tmp_nu<-(
choose(n_nu,i)*(1/2)^(i+(100-i))
)+tmp_nu
}~
tmp_nu~
~
~
次に、シミュレーションで確率を見積もるために以下のコードを実行しました。~
出力は0.13554となりました。~
~
simu_0_nu<-100 #100枚引く。~
simu_1_nu<-1000000 # 何セット試すか。~
result_0_ve<-rep(NA,simu_0_nu)~
result_1_ve<-rep(NA,simu_0_nu)~
result_2_ve<-rep(F,simu_1_nu)~
set.seed(123)~
for(i in 1:simu_1_nu){~
for(j in 1:simu_0_nu){
result_0_ve<-rnorm(simu_0_nu,0,1) #平均0の正規分布は確率1/2で0を越える数値を出す
result_1_ve<-result_0_ve>0 #0を越えたら真
}
result_2_ve[[i]]<-sum(result_1_ve)#真の数を数える
}~
~
sum(result_2_ve>55)~
sum(result_2_ve>55)/simu_1_nu~
~
~
ちなみに、set.seed(123)を~
set.seed(777)にすると0.135447~
set.seed(789)にすると0.136265~
set.seed(111)にすると0.136006~
set.seed(222)にすると0.135414~
となりました。~
どうやら偶然ではなく、私がシミュレーションに使ったコードでは、~
他の計算よりも小さい値になっているようです。~
~
何が原因でシミュレーションと二項分布の計算でこんなにも違ってしまうのでしょうか？~
~
わかる方がいらっしゃいましたら、教えてくださると助かります。~
よろしくお願いいたします。~
//
- 55を含むかどうかの差ではないのですか？ -- &new{2020-12-30 (水) 08:27:58};
> set.seed(1)
> N <- 1000000
> n <- replicate(N, {sum(rnorm(100) > 0)})
> sum(n > 55) / N
[1] 0.135948
> 1 - pbinom(55, 100, 0.5)
[1] 0.1356265
> sum(n >= 55) / N
[1] 0.184649
> 1 - pbinom(54, 100, 0.5)
[1] 0.1841008
- ありがとうございます。全く気が付きませんでした。 -- [[HO0217]] &new{2020-12-30 (水) 09:46:44};
#comment
**所感 [#y895d2fe]
>[[元Rユーザー]] (2020-12-27 (日) 22:33:22)~
~
全盛期とは比較にならないほどの過疎ぶりですね・・・時代は完全にpythonになったのだなと思います。~
//
- ここは初心者に辛辣なコメントが多いのも一因かなぁと。slackのr-wakalangとかは賑わってますよ -- [[Rユーザー]] &new{2021-01-13 (水) 17:29:49};
- 「一時のような辛辣なコメントが減って落ち着いたので、良いんじゃないかなぁ？」と思っていたところです（笑） -- [[現Rユーザ]] &new{2021-01-15 (金) 09:18:03};
#comment
**パネルデータ分析における警告メッセージ（重複エラー？） [#w167b8e1]
>[[R初心者]] (2020-11-08 (日) 22:39:38)~
~
パネルデータでplm分析を行ってますが、下記のメッセージが出て対応に苦慮しております。 ~
~
警告メッセージ: ~
1: pdata.frame(data, index) で: ~
duplicate couples (id-time) in resulting pdata.frame
to find out which, use e.g. table(index(your_pdataframe), useNA = "ifany")
2: is.pbalanced.default(index[[1]], index[[2]]) で: ~
duplicate couples (id-time)
~
恐らく重複エラーなのかもしれませんが、下記を見ると年次などの下線部がエラーの様ですが（？）、削除はしたくなく対応方法を知りたいです。宜しくお願いいたします。~
~
コード PBR 業種コントロール 年次 `1人当たり単体人件費`… 連結子会社数 売上高研究開発費
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 T1301 0.915 0.630 2019 5.87 26 0.0013
2 T1301 0.908 0.618 2018 5.88 24 0.00120
3 T1301 0.872 0.644 2017 6.18 24 0.0011
4 T1301 0.867 0.678 2016 6.15 25 0.00120
//
#comment
**急ぎで知りたいです。 [#k20d893a]
>[[Rビギナー]] (2020-10-28 (水) 01:22:23)~
~
SVM実装の際に、package(kernlab)を使用しています。ksvm関数を使用すると以下のようなエラーが出ます。これはどういう意味でしょうか？~
~
SVM <- ksvm( data_Mic.dummy ~AveMaxTemp7+TP1yp+TN1yp+AveWind7, #モデル式~
data=data_Mic,
kernel ="rbfdot", #カーネル関数の設定
kpar=list(sigma=0.2), #カーネル関数のパラメータ
C=1,
cross=10)
.local(x, ...) でエラー: ~
No Support Vectors found. You may want to change your parameters
追加情報: 警告メッセージ:
.local(x, ...) で: 強制変換により NA が生成されました
//
#comment
**Rで複数の関数のセットをループしたい。 [#u1756c69]
>[[R初心者]] (2020-10-20 (火) 12:55:09)~
~
R初心者で色々わかっていないので、質問の仕方が間違ってたらすみません。~
Rで関数の処理を行いたいのですが、~
~
解析結果のデータは~
~
dds <- 解析結果~
~
に格納しました。~
その大量のデータの中で図示したい５つだけをリストアップし~
dataという下記のようなデータフレームを作りました。~
~
name ID~
1 aaaaa 001~
2 bbbbb 002~
3 ccccc 003~
4 ddddd 004~
5 eeeee 005~
~
処理したい関数が下記のように３つでpng→plotCount→dev.offでワンセットです。~
~
図示したいのが1つであれば~
png("aaaaa.png")~
plotCounts(dds, gene="001", intgroup="condition" )~
dev.off~
~
でデータフレームの一行目のものが図として保存されます。~
ですが連続してdataの二行目、三行目...と自動でループ処理させたくて、~
forやlapplyなどを試したのですが、引数がnameとIDの二つということもあり、うまくプログラムを回せませんでした。~
~
どなたかお詳しい方に上記の関数でpng→plotCount→dev.off→png→plotCount→dev.off→...を自動的に回す方法をご教授頂きたいです。~
~
よろしくお願いします。~
//
- データフレームの中身をループで取得してください
for(i in 1:nrow(dds)){
name_ <- dds$name[i];
ID_ <- dds$ID[i]
}
ループの中でpng, plotCount, dev.offを順次呼び出せば動くと思います。
-- &new{2020-10-23 (金) 11:55:09};
- ありがとうございます。 -- [[R初心者]] &new{2020-11-06 (金) 17:23:08};
#comment
**回帰直線 (RMA) の描き方と統計情報の出し方等について [#bec2344c]
>[[うみねこ]] (2020-07-24 (金) 22:13:58)~
~
R version 3.6.1 (2019-07-05)~
Platform: x86_64-w64-mingw32/x64 (64-bit)~
Running under: Windows 10 x64 (build 18363)~
~
回帰直線の描き方と、その統計情報の出し方について、教えて頂けないでしょうか。~
９列×６行の表があり、その２列目（length）と６列目(weight)にlog変換をして、RMA法を使って回帰直線を描きたいです。gx.rmaというファンクションを見つけましたが、以下のコードではエラーが出て動きません。~
出来たら、95％予測区間や95％信頼区間等についても出し方を知りたいです。~
どうぞご教授よろしくお願いいたします。~
~
data<- read.table('body.txt', header=T)~
row.names(data)=data[,1]~
length <- setNames(log(data[,2]),data[,1])~
weight <- setNames(log(data[,6]),data[,1])~
model <-gx.rma(length,weight)~
model~
//
- 実際の data がどういうものかわからないので，適当にいかのようなデータをつくってやってみたけど，エラーは出ませんでしたが。（質問の仕方がヘタな人がそろっちゃったなあ）
> set.seed = 1234567
> data = as.data.frame(matrix(runif(54), 6, 9))
> data[,1] = letters[1:6]
> row.names(data)=data[,1]
> length <- setNames(log(data[,2]),data[,1])
> weight <- setNames(log(data[,6]),data[,1])
> model <-gx.rma(length,weight)
Reduced Major Axis for length and weight
Means = -1.05 -0.679
SDs = 0.7428 0.8641
Corr = 0.5295
N = 6
Fit = 28%
SE 95% CLs
Slope = 1.163 0.4028 0.1277 <-> 2.199
Intercept = 0.5422 0.518067 -0.7895 <-> 1.874
Accept hypothesis that RMA is (0,1)
> model
$n
[1] 6
$mean
[,1] [,2]
[1,] -1.049833 -0.6789889
$sd
[,1] [,2]
[1,] 0.7428316 0.8640773
$corr
[1] 0.5295494
$a0
[1] 0.5421989
$a1
[1] 1.163221
$sea0
[1] 0.5180675
$sea1
[1] 0.4028333
他の人が追試できるように質問すべし。 -- &new{2020-07-26 (日) 17:41:28};
- すみません。質問が下手で、大変失礼いたしました。データですが、一列目が、人の名前、一行目は身体の長さ、幅、重さ等になっています。再度実行してみましたら、以下のようなエラーが出てしまい、まだ出来ていません。ご教授よろしくお願いいたします。 gx.rma(length, weight) でエラー: 関数 "gx.rma" を見つけることができませんでした -- [[うみねこ]] &new{2020-07-26 (日) 22:19:06};
- > 関数 "gx.rma" を見つけることができませんでした&br;って，library(rgr) やってないんですか？あるいは，install.packages("rgr") さえもしていないとか。 -- &new{2020-07-27 (月) 07:34:13};
- install.packages("rgr") をしていないのが原因でした。確認したところ、gx.rmaのファンクションの記載があるドキュメントの最初のページにrgrと書いてありましたが、スルーしていたようです。無事に上述の結果が出ました！ありがとうございます。次は、この結果を（95％信頼区間と回帰直線）プロットしようとして、plot(model)と打ち込みましたが、xy.coords(x, y, xlabel, ylabel, log) でエラー: 'x' is a list, but does not have components 'x' and 'y'と出て失敗しました。さらに、以下のように打ち込んでみましたが、plot (length, weight,lengthlab=”Length(cm)”,weightlab=”Weight(g)”) エラー: 想定外の入力です in "plot (length, weight,lengthlab=・と出て失敗しました。最小二乗法と同じように作っているつもりですが、plotも出来ません。何度もすみませんが、再度ご教授どうぞよろしくお願いいたします。 -- [[うみねこ]] &new{2020-07-27 (月) 09:44:23};
- > 最小二乗法と同じように作っているつもりですが&br;どの最小二乗法？そもそも，gx.rma の plot method はないと思いますが。あなたが参照している（？）「gx.rmaのファンクションの記載があるドキュメント」に，どのように書いてありますか？書いてなくて，自分で適当に書いてもそれは無理です。lengthlab=”Length(cm)”, weightlab=”Weight(g)”) なんていう引数は plot() にはないのでは？&br;オンラインヘルプを見るとわかりますが，gr.rma の引数に ifplot というのがあるのに気づいていますか？
ifplot
if a x-y plot of the independent variables is required set ifplot = TRUE.
The plot is equi-scaled and the 1:1 line is added. If ifrma = TRUE ifplot will be set to TRUE.
よくわからないまま，適当にプログラムを書いてもだめです。オンラインヘルプくらい読みましょう。それと，エラーメッセージを良く読めば，どこが間違えているか対処法もわかります（わかりやすく書かれているとはいえないでしょうが，それが問題解決の唯一の鍵ですからね） -- &new{2020-07-27 (月) 14:41:47};
- はい。gx.rmaには掲載されていません。if plotは、もちろん気づいていますし、英語の意味も取れますが、それの入力の仕方がわかりません。だから、代表的コマンドの、Plot(x,y)で何とかならないかと、ついでにabline(lm(y,x))も試してみました。おそらくこれらは、ｙ軸方向に残差を取るタイプの一番メジャーな最小二乗法だと思います。オンラインヘルプを読んでも、エラーメッセージを読んでも、読めるけれども、実力がなくて対処が出来ないというのが問題で、初心者打破のため、いろいろ作成して練習していますが、難しいです。 -- [[うみねこ]] &new{2020-07-27 (月) 16:23:18};
- > それの入力の仕方がわかりません&br;Usage の下に，関数の使い方（引数の種類と，その指定法）が書いてあるじゃないですか。
Usage
gx.rma(xx1, xx2, x1lab = NULL, x2lab = NULL, log = FALSE, ifplot = FALSE,
ifrma = FALSE, ifcoeffs = FALSE, ifform = FALSE, iftest = FALSE, ...)
あなたの例の場合だと，
gx.rma(length, weight, x1lab="Length(cm)", x2lab="Weight(g)", ifplot=TRUE, ifrma=TRUE)
のようになる。 -- &new{2020-07-27 (月) 16:27:48};
- ありがとうございます。無事にplotされました。このUsageもですが、もしgx.rma(x,y,xlab=NULL, ylab=NULL...)であったら、もう少し私にはわかりやすいんです。使用例のコードを見ると、余計にUsageが解読し辛くなってしまう位のレベルですが、諦めずに頑張ろうと思います。ありがとうございます。 -- [[うみねこ]] &new{2020-07-27 (月) 17:47:53};
#comment
**Treatment indicator ('Tr') must be a logical variable---i.e., TRUE (1) or FALSE (0) [#bdc9264a]
>[[catdogcat]] (2020-07-24 (金) 18:20:51)~
~
マッチングを試みましたが下記のエラーが出ました~
~
Treatment indicator ('Tr') must be a logical variable---i.e., TRUE (1) or FALSE (0)~
~
比較する群の変数は因子としてあり 0.1 となっているのですが~
解決策がわかりません。~
ご教授お願い致します。~
//
- 実際にあなたが何をやったかが，全くわからない。なのに，解決策を授けられる人がいるだろうか？いや，いない。&br;「群の変数は因子としてあり 0.1」ということは，"must be a logical variable" という条件を満たしていないのではないかな？ -- &new{2020-07-26 (日) 17:46:21};
- 幼稚な質問であり申し訳ございません。確かにlogical variableになってなかったようでした。ありがとうございます。 -- [[catdogcat]] &new{2020-07-26 (日) 22:23:43};
#comment
**Rcppで行列を多次元配列にfor文で代入する [#m2aaf0fb]
> (2020-07-22 (水) 23:44:38)~
~
下記コードをpgm.cppに保存後、Rcpp::sourceCpp("pgm.cpp")でコンパイルし~
test(matrix(c(1,2,3,4),2,2))を実行したのですが1行1列の値が返ってこずエラーでRが落ちます。もし良ければご助言いただけないでしょうか。~
~
~
#include<Rcpp.h>
#define MAX 100
using namespace std;
// [[Rcpp::export]]
int test(Rcpp::IntegerMatrix mtx){
int matrix[MAX][MAX];
for(int i=0; i < mtx.rows(); ++i){
for(int j=0; j < mtx.cols(); ++j){
matrix[i][j] = mtx(i,j);
}}
return matrix[0][0];
}
//
- やってみたけど，確かに1x1の行列は帰ってこないけど，R が堕ちるというようなことはなかったけど？ -- &new{2020-07-26 (日) 17:47:18};
#comment
**くの字に曲がるデータプロットの特異点の算出？ [#b8b95798]
>[[がんば]] (2020-07-22 (水) 08:35:12)~
~
昨年はおせわになりました。実験結果をグラフにすると直線状ではなく「くの字」状に曲がった結果となる場合もあるのですがこの時の特異点をＲで一発算出するような手段がないでしょうか？「データをソートして下限側から徐々に近似曲線を構築するのと上限側から構築するのをそれぞれ行っていって次の点がそれまでの近似曲線のバラツキ範囲外になる点を求める」ような操作をする（人の目視確認の場合はこんな手順だと思います）手段を意図してます。~
動作環境：~
> sessionInfo() ~
R version 4.0.2 (2020-06-22)~
Platform: x86_64-w64-mingw32/x64 (64-bit)~
Running under: Windows 10 x64 (build 18363)~
Matrix products: default~
locale:~
[1] LC_COLLATE=Japanese_Japan.932 LC_CTYPE=Japanese_Japan.932 LC_MONETARY=Japanese_Japan.932~
[4] LC_NUMERIC=C LC_TIME=Japanese_Japan.932 ~
attached base packages:~
[1] stats graphics grDevices utils datasets methods base ~
loaded via a namespace (and not attached):~
[1] compiler_4.0.2 tools_4.0.2~
//
- Googleで"bent line gression R"で検索すると、Rbentパッケージが出てきたけど、rbentfit()が吐くbpがもしかしたら、お望みのものでは -- &new{2020-07-22 (水) 22:13:43};
- http://aoki2.si.gunma-u.ac.jp/R/oresen.html 二本の直線による折線回帰 とかは？ -- &new{2020-07-26 (日) 17:17:55};
- ↑を元に検討してみます。ありがとうございます。 -- [[がんば]] &new{2020-07-29 (水) 08:41:50};
- RbentパッケージについてはhelpのExsamplesをいじってみて検討します。ありがとうごじます。 -- [[がんば]] &new{2020-07-29 (水) 09:01:05};
#comment
**R studio でグラフ化した時に日本語文字化けが起こる [#bacaed2a]
>[[りょう]] (2020-06-30 (火) 21:38:04)~
~
Rによる多変量解析と言う本でRを動かしているのですが、グラフ化すると日本語の部分が□□□のようになります。~
よくわからないのが、"plot", "plotmeans"などでは日本語表示されるのですが、xyplot, histogramだと文字化けすることです。以下がその時のコードです。~
~
>jhk <- read.csv('人事評価結果.csv')~
head(jhk, 3)~
|ID|性別|部署|年代|協調性|技能|知識|ストレス|~
|1|M|A部|中堅|70|65|71|53|~
|2|F|B部|熟練|45|62|70|46|~
・・・・・・~
par(family = "HiraKakuProN-W3") #一応これ試してみました~
library(lattice, gplots)~
histogram(~ストレス, data = jhk, breaks = 20, type='count') #文字化けする~
xyplot(知識~技能|年代+部署, data=jhk) #文字化けする~
plotmeans(協調性~性別, data=jhk, p=0.95, ylim=c(49, 54)) #文字化けしない~
~
・macOS Catalina バージョン10.15.5~
・R studio version 4.0.2 です。この前に使用していたバージョン3.6.3でも同様の状態でした。~
・CSVファイルはUTF-8で保存されています。~
・par(family = "HiraKakuProN-W3") なども試しました。~
~
GoogleやCRAN, マニュアルなどで既に調べたのですがわからず、困っているのでご教授いただけると幸いです。~
//
- ごめんなさい。コードの部分を青くしたり、表形式にしたつもりができていません汗 -- [[りょう]] &new{2020-06-30 (火) 21:39:55};
- お門違いの場所にアップされ，誰も返事をしない，もう見ていないかもしれないが一応書いておこう。&br;histogram, xyplot は lattice のものなので，日本語フォント指定が旨くいかない（なんか方法はあるはずだが，ネットをあさっても...）。&br;回避法としては&br;「画面に描かず PDF ファイルに保存する」とちゃんと描かれるっす。
pdf("test.pdf") # ひつようなら他のパラメータの設定も
histogram(~ストレス, data = jhk, breaks = 20, type='count') #文字化しない
xyplot(知識~技能|年代+部署, data=jhk) #文字化けしない
plotmeans(協調性~性別, data=jhk, p=0.95, ylim=c(49, 54)) #文字化けしない
dev.off()
どんなもんだ？ -- &new{2020-07-14 (火) 14:04:24};
#comment
**ggplot [#bc6d7d37]
>[[ぼう]] (2020-05-14 (木) 12:25:32)~
~
R4.00でggplotを使おうと思っていますが、エラーが出てしまいます。~
~
次のコマンドを実行すると、下のような結果が出てきてしまいます。
> install.packages("ggplot2", dependencies = TRUE)
> library(ggplot2)
> str(iris)
'data.frame': 150 obs. of 5 variables:~
$ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
$ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
$ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
$ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
$ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
> ggplot(data=iris,aes(x=Sepal.Length,y=Sepal.Width))~
ggplot(data = iris, aes(x = Sepal.Length, y = Sepal.Width)) でエラー:
関数 "ggplot" を見つけることができませんでした
//
- ggplot2 がインストールされていない（OS はなんですか。Mac では R4.0 なら，コンパイルが必要なインストールが延々と実行されますが，それがちゃんと出来ているか。Windows ならそうでもないのかな）~
と，思ったが，Windows でもインストールはコンパイルが必要なようですね。で，それでも，ちゃんとインストールできましたが。~
あなたの場合，コンパイルに失敗していると言うことはありませんか。コンパイルが可能な環境が出来ていないのかな。普通は，そんなことないはずなのだけど。 -- &new{2020-05-15 (金) 21:38:52};
- ありがとうございます。Win10を使っています。library(ggplot2)でエラーが出ていないのに、ggplotを実行するとエラーが出てしまうので、原因が分かりません。。。 -- [[ぼう]] &new{2020-05-18 (月) 11:28:18};
- ありがとうございます。Win10を使っています。library(ggplot2)でエラーが出ていないのに、ggplotを実行するとエラーが出てしまうので、原因が分かりません。。。 -- [[ぼう]] &new{2020-05-18 (月) 12:09:38};
#comment
**ccf [#y347335b]
>[[show]] (2020-03-10 (火) 12:58:28)~
~
ccf(rw,rw2)と入力した場合、どちらを基準にしてlagを取っているのでしょうか？~
//
- ? ccf -- "The lag k value returned by ccf(x, y) estimates the correlation between x[t+k] and y[t]." -- &new{2020-03-10 (火) 17:55:37};
- rw2を基準としているということでしょうか。 -- [[show]] &new{2020-03-11 (水) 12:34:59};
- what? rw2? what do you mean?-- &new{2020-03-11 (水) 22:56:56};
- ccf(rw,rw2)の場合、rw２を基準にしてるということでしょうか？ -- [[show]] &new{2020-03-13 (金) 08:33:16};
- x[t+k] and y[t] -- &new{2020-03-13 (金) 10:12:16};
- なぜ自分で確かめもせず，不毛な質問を続けるのかなあ？~
以下のようなテストデータを作ってみる。ベクトルの長さ n = 10
> set.seed(12345)
> n <- 10
> x <- rnorm(n)
> y <- rnorm(n)
ラグ lag = 2，で ccf してみる
> print(ccf(x, y, lag, plot=FALSE), digits=7)
Autocorrelations of series ‘X’, by lag
-2 -1 0 1 2
0.1416173 0.5327757 -0.2716881 0.2595251 -0.1681542
自分で相互相関係数を計算するときに必要なもの。
> mean.x <- mean(x) # x の平均（全部使う）
> mean.y <- mean(y) # y の平均（全部使う）
> ss.x <- var(x)*(n-1) # x の変動（全部使う）
> ss.y <- var(y)*(n-1) # y の変動（全部使う）
相互相関係数を計算する。lag = 2~
使用するデータベクトル~
x[(lag+1):n]-mean.x つまり，lag+1 から n まで，長さ n-lag。平均値 mean.x を引いておく。~
y[1:(n-lag)]-mean.y つまり，1 から n-lag まで，長さ n-lag。平均値 mean.y を引いておく。~
両者掛け合わせて合計し（部分の共変動だ） sum((x[(lag+1):n]-mean.x) * (y[1:(n-lag)]-mean.y) ~
x，y 両方の全体の変動の積の平方根 sqrt(ss.x * ss.y) で割る
> sum((x[(lag+1):n]-mean.x) * (y[1:(n-lag)]-mean.y) / sqrt(ss.x * ss.y)
[1] -0.1681542
上の ccf の lag = 2 の下に書かれている数値 -0.1681542 と一致しているでしょ。~
つまり，ヘルプに書いている 「x[t+k] and y[t]」というのはそいういうこと，です。よ。~
x[(lag+1):n] と y[1:(n-lag)] どちらを基準にしていますか？~
そんなこと，どっちだっていいじゃないか。 -- &new{2020-03-13 (金) 11:11:38};
#comment
**R のアップデート [#rac7c1ae]
>[[私は誰でしょう？]] (2020-03-01 (日) 21:20:52)~
~
今まで，最新の R にするために（新しもの好きなので）~
あ，見ての通り，Mac ユーザです。~
~
curl https://mac.r-project.org/el-capitan/R-3.6-branch/R-3.6-branch-el-capitan-sa-x86_64.tar.gz | sudo tar fvxz - -C /~
~
~
をやっていたのですが，2 ヶ月ほど前から，実行結果が思わしくなくなっています。~
つまるところ，件の URL はアクセス不可ということのようで，アップデートは出来ません。~
~
念の為~
http://mac.r-project.org/el-capitan/R-3.6-branch/R-3.6-branch-el-capitan-signed.pkg~
を使えばアップデートは出来ます。~
//
- サイトの方のWAFとかでブロックされているのか, 単なる設定ミスかはわかりませんが当方でも同じ症状です.
curl -O http://mac.r-project.org/el-capitan/R-3.6-branch/R-3.6-branch-el-capitan-signed.pkg
sudo installer -pkg ./R-3.6-branch-el-capitan-signed.pkg -target /
とかでお茶を濁しては... -- [[何処の誰かは知らないけれど]] &new{2020-03-02 (月) 16:03:50};
- 有難うございます。取りあえずは，それで，しのいでみます。 -- &new{2020-03-02 (月) 22:12:06};
#comment
**連番でファイルを読み込みたいときにファイルに欠番があるときの処理 [#me7fd7f2]
>[[Rビギナー]] (2020-02-26 (水) 10:12:45)~
~
お世話になります．~
下記のように1～3まで連番のファイルを読み込みたいのですが，例えば2番のファイルが無い場合にそのまま1，3のファイルを読み込んでプログラムが止まらないようにするにはどうしたらいいでしょうか．~
try関数で可能かと思いやってみたのですが，だめでした．~
ご教示ください．~
よろしくお願いいたします．~
~
＝＝ここから~
point <- 1:3 # for文を回す数を設定
for (i in point) {
mypath <- file.path(getwd(), paste0(i,"t2-1.csv"))
df <- assign(paste0("data",i), read.csv(mypath, header = T))
df <- data.frame(df)
df
head(df)
df$X <- strptime(df$X,"%Y/%m/%d %H:%M'%S")
head(df)
df$X <- as.POSIXlt(df$X,"%Y/%m/%d %H:%M:%S",tz="GMT0")
head(df)
#https://qiita.com/tktz/items/733c37b1d6102ae52120
library(xts)
df$X <- align.time(df$X-10*60/2,10*60)
head(df)
date = paste0("date",i)
Temp = paste0("Temp",i)
RH = paste0("RH",i)
colnames(df) <- c(date,Temp,RH)
head(df)
mypath <- file.path(getwd(), paste0(i,"t1.csv"))
write.csv(df,mypath)
}
//
- insert 1 sentence. -- &new{2020-02-26 (水) 10:46:10};
mypath <- file.path(getwd(), paste0(i,"t2-1.csv"))
if (! file.exists(mypath)) next
- or -- &new{2020-02-26 (水) 11:09:20};
files = list.files(getwd(), pattern="[0-9]*t2-1.csv")
for (f in files) {
# describe what you want to do
}
- お二人の方，ありがとうございました．おかげさまでうまくいきました！ -- [[Rビギナー]] &new{2020-02-26 (水) 14:57:32};
- お二人の方，ありがとうございました．おかげさまでうまくいきました！ -- [[Rビギナー]] &new{2020-02-26 (水) 15:47:15};
#comment
**時系列 [#h477736a]
>[[マリン]] (2020-02-23 (日) 10:13:04)~
~
２つの時系列のデータ（約８０個）の相関（どのくらい類似性があるか）を~
見たいのですが、方法を教えてください。~
//
- まずはググってみて，それからにしなさい -- &new{2020-02-23 (日) 21:09:21};
- ？ -- [[マリン]] &new{2020-02-24 (月) 07:45:37};
#comment
**固定効果を用いた回帰線 [#s967c99e]
> (2020-02-14 (金) 14:34:30)~
~
Book1というcsvファイルのデータを用いてGLMMによる解析を行ったのですが、固定効果の値を用いた回帰線がうまく引けません。~
csvファイルおよび解析を行ったスクリプトをBook1GLMMというファイル名でアップロードしましたので、ご教授いただけますと幸いです。~
よろしくお願いいたします。~
//
- do! just shown as article '2020-01-30 (木) 12:37:42' the same program isn't it? -- &new{2020-02-14 (金) 17:11:14};
data <- read.csv("Book1.csv")
data
# update.packages
# install.packages("lme4")
library(lme4)
par(family = "HiraKakuProN-W3")
model1 <- glmer(POMN ~ log(siltclayN) + (1 | Place), family = Gamma(link = "log"), data = data)
summary(model1)
place = c("HA", "TN", "OG", "YN", "KG", "WA", "AS") # データフレームではこちらを使っているのでこれでなきゃだめ
place2 = c("半田山", "田野", "小川", "与那", "上賀茂", "西粟倉", "足寄") # 凡例用には別のものを用意する
pch = c(25, 18, 17, 16, 15, 9, 8)
col = c("black", "red", "cadetblue", "magenta", "brown", "coral", "darkorchid")
plot(data$siltclayN, data$POMN, pch = rep(pch, each = 5), col = rep(col, each = 5), cex = 3, cex.axis = 2)
legend("bottomright", legend = place2, pch = pch, col = col, cex = 2)
xx <- seq(min(data$siltclayN), max(data$siltclayN), length = 100)
# yy <- exp(1.1136 + 0.4651 * xx) # そもそも間違い
# lines(xx, yy, lwd = 2)
for (i in seq_along(place)) { # 7 本の曲線を引くんでしょ？
newdata = data.frame(siltclayN=xx, Place=rep(place[i], 100))
yy = predict(model1, newdata=newdata, type="response")
lines(xx, yy, col=col[i])
}
- 7本の曲線は引けました。勉強不足でお恥ずかしいのですが、固定効果の推定値を用いて一本だけ回帰線を引くことはできないのでしょうか？ -- &new{2020-02-14 (金) 18:08:12};
- quit -- &new{2020-02-14 (金) 21:04:44};
#comment
**ガンマ分布 [#f464b927]
> (2020-02-06 (木) 21:14:35)~
~
GLMMなどで、目的変数に用いるデータが連続値で0以上であればガンマ分布を用いると統計の教科書等には書いてあるのですが、連続値で0以上(データの最小値は0)のデータでガンマ分布を用いて解析しようとすると~
「non-positive values not allowed for the 'gamma' family」というエラーが出ます。~
また、同じデータで線形混合効果モデルなら結果は出るのですが、「boundary (singular) fit: see ?isSingular」という警告？が出るせいなのかランダム効果を考慮した回帰線が一本しか出ません。~
解決策などありましたら、お教えいただけますと幸いです。~
//
- show your data and your script(program). -- &new{2020-02-06 (木) 22:07:07};
- model1<-glmer(A~log(B)+(1|Place),family=Gamma(link="log"),data=d) In this script was error 「non-positive values not allowed for the 'gamma' family」. -- &new{2020-02-06 (木) 22:26:52};
- model2<-lmer(A~B+(1|Place),data=d) In this script was 「boundary (singular) fit: see ?isSingular」 -- &new{2020-02-06 (木) 22:28:59};
- The data was uploaded. It is data named data.csv. -- &new{2020-02-06 (木) 22:31:00};
- 'non-positive values not allowed for the 'gamma' family' means 0 is not non-positive, isn't it.
> library(lme4)
> d = read.csv("data.csv")
> model1 <- glmer(A ~ log(B) + (1|Place),f amily=Gamma(link="log"),data=d)
eval(family$initialize, rho) でエラー:
non-positive values not allowed for the 'gamma' family
using d2
> d2 = d[d$A > 0, ]
> model1-2 <- glmer(A~log(B) + (1|Place), family=Gamma(link="log"), data=d2)
> model1-2
no error happend.&br;
can you understand?&br;you must *respect* 'error message'. it warns you are wrong, and he says the reason 'why'.&br;if your data contains 'zero's, it may not be a gamma distribution!!!!! -- &new{2020-02-06 (木) 23:37:25};
- I will try it once. One more thing I want to ask. -- &new{2020-02-07 (金) 00:00:48};
- Now I uploaded another data. When the script of model2 is executed with this data, the message "boundary (singular) fit: see? IsSingular" appears, and only one regression line considering random effects is drawn. What is the cause? -- &new{2020-02-07 (金) 00:03:40};
- you! just , imput '? IsSingular'(without last '?') in console!!!. message says so, ah. -- &new{2020-02-07 (金) 00:11:48};
model2<-lmer(NitMAOMrecover~MAOMCN+(1|Place),data=example)
summary(model2)
place=c("HA","TN","OG","YN","KG","WA","AS")
pch=c(25,18,17,16,15,9,8)
col=c("black","red","cadetblue","magenta","brown","coral","darkorchid")
plot(example$B,example$A,pch=rep(pch,each=5),col=rep(col,each=5),cex=3,cex.axis=2)
legend("bottomright",legend=place,pch=pch,col=col,cex=1.5)
xx<-seq(min(example$B),max(example$B),length=1000)
for(i in seq_along(place)){
yy<-predict(model2,newdata=data.frame(B=xx,Place=rep(place[i],1000)),type="response")
lines(xx,yy,col=col[i],lwd=2)}
- Excuse me. I want to show the script to draw the regression line on the site, but I don't know how to do it. -- &new{2020-02-07 (金) 00:29:15};
- they said. "その前に，「[[投稿における注意事項]]」を読んでから！" and "投稿する前にまず &heart; [[投稿文書の書式:http://www.okadajp.org/RWiki/?%E6%95%B4%E5%BD%A2%E3%83%AB%E3%83%BC%E3%83%AB]] &heart; を読んでください" -- &new{2020-02-07 (金) 08:20:14};
- surprisingly, it's a reckless and messy program.~
Four variables "Sample,Place,A,B" are included in your "example.csv".~
but, your wrote "NitMAOMrecover~MAOMCN+(1|Place),data=example".~
more over, you wrote variables "A", "B" as "example$B" or "example$A"~
a data-frame in lmer function is different from a "new" data-frame in predict function. ~
you have to study a basis of R. -- &new{2020-02-07 (金) 08:23:45};
#comment
**GLMMについて [#tc8ba31b]
>[[ra]] (2020-01-27 (月) 11:58:45)~
~
R初心者です。~
GLMMで解析を行った後、散布図に回帰線を引きたいのですが、predict関数を読み込ませようとすると（オブジェクト"Place"がありません）というエラーが出ます（Placeはランダム効果です）。このエラーを解消し、回帰線を引くにはどうすれば良いでしょうか？~
入力したコマンドを載せておきます。~
よろしくお願いいたします。~
~
model1<-glmer(POMN~log(siltclayN)+(1|Place),family=Gamma(link="log"),data=data)~
plot(data$siltclayN,data$POMN,pch=c(25,18,17,16,15,9,8),cex=3,cex.axis=2)~
legend("bottomright",legend=c("HA","TN","OG","YN","KG","WA","AS"),pch=c(25,18,17,16,15,9,8),cex=2)~
xx<-seq(min(data$siltclayN),max(data$siltclayN),length=100)~
yy<-predict(model1,newdata=data.frame(siltclayN=xx),type="response")~
lines(xx,yy,lwd=2)~
//
- new data set has to include 'siltclayN' *and* 'Place' -- &new{2020-01-27 (月) 13:47:56};
- Sorry,how do I create that data set? -- &new{2020-01-27 (月) 16:26:31};
- just, create -- &new{2020-01-27 (月) 22:11:04};
-
POMN <- c(data$POMN)
siltclayN <- c(data$siltclayN + Place)
dataframe <- data.frame(POMN, siltclayN)
model1 <- glmer(POMN~log(siltclayN) + (1|Place1), family=Gamma(link="log"), data=dataframe)
plot(dataframe$siltclayN, dataframe$POMN, pch=c(25,18,17,16,15,9,8), cex=3, cex.axis=2)
legend("bottomright", legend=c("HA","TN","OG","YN","KG","WA","AS"), pch=c(25,18,17,16,15,9,8),cex=2)
xx <- seq(min(dataframe$siltclayN), max(dataframe$siltclayN), length=100)
yy <- predict(model1, newdata=data.frame(siltclayN=xx), type="response")
-- &new{2020-01-27 (月) 23:14:23};
- What's wrong? -- &new{2020-01-27 (月) 23:16:01};
- newdata=data.frame(siltclayN=xx, Place = ???????) -- &new{2020-01-28 (火) 10:50:09};
- In this command and Place=xx, an error has occurred. Where should I improve? -- &new{2020-01-28 (火) 14:42:23};
- Place is factor isn't it? 'Place = xx' may be invalid!&br;
I'll show you an e xample.
> df = data.frame(x=c(3,2,1,3,5), f=factor(c(1,2,2,3,3)), y=c(2,4,7,5,6))
> fit = lm(y ~ x + f, data=df) # f may not be a factor
> predict(fit) # this prediction uses data-frame df
1 2 3 4 5
2.0 5.4 5.6 5.7 5.3
if you predict with another data-frame, the data-frame must contain 'x' and 'f'
following newdata contains 'x' and 'f'. first 5 elements are same in df(df$x, df$f), last 3 elements is exactly new data
> newdata = data.frame(x=c(3,2,1,3,5, 2,3,5), f=factor(c(1,2,2,3,3, 1,2,2)))
> predict(fit, newdata)
1 2 3 4 5 6 7 8
2.0 5.4 5.6 5.7 5.3 2.2 5.2 4.8 # first five values are identical to previous prediction. last 3 values are new prediction.
do you understand?? -- &new{2020-01-28 (火) 16:28:06};
- Place is a factor. In that case, is Place=c(data$Place) ok? Or should all numerical data and factor data be written in c of x,f and y? -- &new{2020-01-29 (水) 10:40:07};
- I ran it with reference to your example, how can I draw a prediction line from here? -- &new{2020-01-29 (水) 11:37:47};
- if new data-frame is declared as 'newdata=data.frame(siltclayN=xx, Place=pp), pp[i] corresponds to xx[i].&br;examine above example carefullly. predict function must be called as 'predict(fit, newdata=data.frame(x=c(2,3,5), f=c(1,2,2)))' -- &new{2020-01-29 (水) 13:34:36};
- 'prediction line'?? no, no, prediction *prediction plane*. there are *two* !! independent variables(x, f) and one dependent variable(y), so, it's a 3 dimensional data(Actually, if f is a factor and has more than 2 levels, the dimension is more than 3!!). your data is 3 or more dimensinal, too.&br;you know, in case of simple regression analysis y = a*x+b you have one independent variable x, and one dependent variable y, it's a *two* dimensional data, and you can draw two dimensional scatter plot or regression *line*. &br;if you have 3 independent variable model, no one can draw predicted hyperplane(you can't imagine 4 dimensional space, can you??.&br;however, if you draw POMN vs. siltclayN by Place, you can draw 2 dimensional scatter or regression line. -- &new{2020-01-29 (水) 13:41:39};
- your model is POMN = α * log(siltclayN[i]) + β * Place[i]. if you want to predict POMN for new data set, you must prepare new data pair siltclayN[j] and Place[j], j = 1, 2, 3, ... -- &new{2020-01-29 (水) 14:02:09};
- your data-frame is curious, isn't it?&br;
what do you mean
siltclayN <- c(data$siltclayN + Place)
dataframe <- data.frame(POMN, siltclayN)
if Place is a factor variable you cant add to data$siltclayN. 'data$siltclayN + Place' cause an error. -- &new{2020-01-29 (水) 14:47:36};
- upload your data is best way to get solution -- &new{2020-01-29 (水) 14:54:46};
- Did you fall into a deadlock? -- &new{2020-01-29 (水) 22:32:29};
- yes... I don't know how to draw a regression line considering random effects. -- &new{2020-01-30 (木) 10:34:08};
- Can you tell me what to improve specifically? -- &new{2020-01-30 (木) 10:36:21};
- upload your data is best way to get solution -- &new{2020-01-30 (木) 10:46:16};
- Upload now. It is a file called Book1. -- &new{2020-01-30 (木) 10:49:42};
- alright, I was able to upload. -- &new{2020-01-30 (木) 11:31:58};
- Your program is strange overall. for example, pch=c(25,18,17,16,15,9,8) must be pch=rep(c(25,18,17,16,15,9,8), 5). examine each point's coordinates.&br;
you may want a result like this. -- &new{2020-01-30 (木) 12:37:42};
library(lme4)
data = read.csv("Book1.csv")
model1 <- glmer(POMN ~ log(siltclayN) + (1 | Place),
family = Gamma(link = "log"), data = data)
summary(model1)
place = c("HA", "TN", "OG", "YN", "KG", "WA", "AS")
pch = c(25, 18, 17, 16, 15, 9, 8)
col = c("black", "red", "cadetblue", "magenta", "brown", "coral", "darkorchid")
plot(data$siltclayN, data$POMN, pch = rep(pch, each = 5),
col = rep(col, each = 5)) # , cex=3, cex.axis=2)
legend("bottomright", legend = place, pch = pch, col = col) # ,cex=2)
xx <- seq(min(data$siltclayN), max(data$siltclayN), length = 100)
for (i in seq_along(place)) {
yy <- predict(model1,
newdata = data.frame(siltclayN = xx, Place = rep(place[i], 100)),
type = "response")
lines(xx, yy, col = col[i])
}
&ref(glmer.png);
- I wanted to make this diagram. It was very helpful! I didn't have enough knowledge. -- &new{2020-01-30 (木) 13:48:53};
- I'm really thankful to you! -- &new{2020-01-30 (木) 13:50:02};
- I want to draw a regression line with the estimated fixed effects of model1, but I cannot do it well. -- &new{2020-02-13 (木) 21:49:52};
#comment
**文字列からバイト数を指定して切り出す [#dcc3e621]
> (2020-01-26 (日) 22:45:01)~
~
半角・全角が混在する文字列からバイト数を指定して、部分文字列を抽出したいのですが、どのような関数がありますでしょうか。~
例えば~
x <- "あaいbうc"~
という文字列から3-5byte目を抽出したいです。~
（"aい"という部分文字列を抽出したいです。）~
~
substringは文字数を指定する関数のため、byte単位で抽出が不可能なようです。~
~
お手数ですがご教示お願いいたします。~
//
- ? substr / If an input element has declared "bytes" encoding (see Encoding, the subsetting is done in units of bytes not characters. -- &new{2020-01-27 (月) 11:46:10};
- そもそも3-5byte目が「aい」というは正しいですか？ -- [[S]] &new{2020-01-27 (月) 16:29:39};
> x <- "あaいbうc"
> charToRaw(x)
[1] e3 81 82 61 e3 81 84 62 e3 81 86 63
「あ」は3バイト文字なので途中で切られて、3-5byte目は「82 61 e3」です。これは「あ」の3バイト目、1バイト文字の「a」、「い」の1バイト目です。~
全角半角の字数をバイト数で何とかしようとするのは無理があ流のではないでしょうか。
- it's just he wants! -- &new{2020-01-27 (月) 22:13:13};
- euc-jp uses 2 bytes for 1 japanese character
> x <- "あaいbうc" # UTF-8: あ == \xe3\x81\x82
> xx = iconv(x, "utf-8", "euc-jp")
> xx
[1] "\\xa4\\xa2a\\xa4\\xa4b\\xa4\\xa6c" # EUC-JP: あ == \xa4\xa2
> Encoding(xx) = "bytes"
> iconv(substr(xx, 3, 5), "euc-jp", "utf-8")
[1] "aい"
cp932 uses 2 bytes for 1 japanese character, too
> x <- "あaいbうc"
> yy = iconv(x, "utf-8", "cp932")
> yy
[1] "\x82\xa0a\x82\xa2b\x82\xa4c" # CP932: あ == \x82\xa0
> Encoding(yy) = "bytes"
> iconv(substr(yy, 3, 5), "cp932", "utf-8")
[1] "aい"
do you understand? -- &new{2020-01-28 (火) 11:38:51};
- 文字コードを変換してsubstrで切り出すということですね。ご回答ありがとうござました。 -- &new{2020-01-30 (木) 00:44:16};
#comment
**apply系関数でfunに線形回帰lmが適用可能な関数 [#v28d972a]
>[[TY]] (2020-01-07 (火) 15:59:04)~
~
シミュレーションで同じ関数を繰り返し(10^5~10^6回程度)実行する場面を想定しています．~
例えば，以下のsample_dataでREP毎にlm(y~x)を実行したいです．~
for文で実行することはできるのですが，一般にapply系関数の方が実行が速いと聞き，~
探しています．apply系関数に関する知識は，tapply関数にmean関数を適用できる，~
tapply(sample_data$y, sample_data$REP,mean)~
程度のレベルです．~
~
【sample_data】~
y x REP
1 133 3.9 1~
2 81 -2.3 1~
3 104 -1.2 1~
4 89 0.7 1~
5 142 4.0 1~
6 74 -2.9 1~
7 140 3.9 1~
8 170 4.4 1~
9 128 1.6 1~
10 107 1.2 1~
11 49 -4.3 2~
12 79 -2.9 2~
13 77 -3.2 2~
14 125 1.8 2~
15 82 -1.1 2~
16 132 2.6 2~
17 99 -0.0 2~
18 140 2.1 2~
19 140 4.9 2~
20 88 -1.1 2~
//
- 「REP毎にlm(y~x)を実行したい」のか，同じような 10^5 ほどのデータを分析したいのか。教えてあげてもいいけど，そのへんどうなのよ。~
つまり，rep が 1~10^5 くらいの種類があって，それぞれについて x,y の組のデータ（サンプルサイズが 10 くらい）があるんだけど，それぞれについて lm で分析してータセットごとの結果が欲しいのか？&br;結果を待っているより，やってみた方が速いよ（ワラ）。反応があるという保証がない限り，回答しない。&br;それにしても，質問の仕方がへた。 -- &new{2020-01-10 (金) 22:10:35};
- ありがとうございます．出直してきます． -- [[TY]] &new{2020-01-13 (月) 10:59:17};
- 必要なのは，そういう対応ではないと思うが。質問にも答えないし，誠意がない。というか，本当に回答が必要なのか？ -- &new{2020-01-16 (木) 21:03:50};
#comment
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
