統計解析R たぶんpart3くらい
>>419 これ見て自分も2.0.1にした 確かに問題なくなった 報告ありがとう 間違って機械学習スレに書き込んでしまったので、こちらに再書き込み。 Rmdファイル名に日本語が入っているとRStudioでread系の出力が表示されなくなるという謎のトラップに嵌った。 read.csv関数、readr::read_csv関数、readxl::read_excel関数などを端折ってread系と書いてしまった。 動くんだけどエラーも何も表示されないという現象だった。 適当に日本語のrmdからread_csv試してみたけど動作もメッセージも特に問題ないみたい 自分が条件を勘違いしてるだけかもだけど readrの2.1.0が来てた デフォルトで遅延読込するのが評判悪くて戻したらしい Rstudioで100個のtxtファイル (中身はTsvで3万行2列、1列目は全ファイル共通) を3万行、101列のマトリクスファイルにしようとしています。 left_joinで1列目をキーにして結合するためにデータを100個読み込む際に lf <- list.files(full.names = T) data <- lapply(lf,read.delim) を用いたのですが、 mat <- full_join(data)で結合させようとすると 'full_join' をクラス "list" のオブジェクトに適用できるようなメソッドがありません となってしまいました。 別の関数でこのような結合が可能でしょうか。 >>426 dplyr::bind_cols() があるとわかり(rowsしか見つけられていませんでした)、 こちらでdata2<-dplyr::bind_cols(data[1],data[2]) でエラーが出ないことまで分かったのですが、 今度は 1列目が全残りしてしまいました。。。 >>426 詳細はわかりませんがとりあえず読み込みの部分を以下にすれば結合されたdata.frame(tibble)ができると思います data <- readr::read_tsv(lf) あ、ごめんなさい 勘違いしてたのでいったん428は無視してください bind_colの代わりにdplyr::left_join()で結合してください >>426 確認だけど各ファイルの一列目は言わるインデックスで、どのファイルも全く同じ値で全く同じ並びなの? 例えば、1,2,3,...30,000のような。 >>430 コメントいただきありがとうございます。 data2<-dplyr::left_join(data[1],data[2]) UseMethod("left_join") でエラー: 'left_join' をクラス "list" のオブジェクトに適用できるようなメソッドがありません となってしまいました。 >>431 そのご認識であっています >>432 1列目がすべて同じ値でデータを結合させるキーとして使う必要がなければ、こんな感じ。 最後に1列目をバインドかジョインさせれば、概ね目的が達成できるのではないかと。 require(tidyverse) # サンプルデータの作成 for (i in 1:10) { data.frame( a = seq(from = 1, to = 100), b = rnorm(100) ) %>% readr::write_excel_csv(file = paste0("./sample_data/sample_", i, ".csv")) } # サンプルデータの結合 list.files("./sample_data", full.names = TRUE) %>% purrr::map_dfc(.f = function(x){readr::read_csv(x) %>% dplyr::select(2)}) 列名でデータが識別出来るようにしたければ、もうひと工夫必要だけど。 ファイル名を列の識別に使うとこんな感じ # サンプルデータの読み込みと結合 list.files("./sample_data", full.names = TRUE) %>% purrr::map_df(.f = function(x){readr::read_csv(x) %>% dplyr::mutate(file = basename(x))}) %>% tidyr::pivot_wider(names_from = file, values_from = b) >>433 おお!ありがとうございます。 解釈しながら進めてみます。 >>434 の方は、1列目をキーにleft_joinするのと同じ結果になるので、こちらの方が処理としては汎用性が高いかも。 # サンプルデータの作成 for (i in 0:4) { data.frame( a = seq(from = 1, to = 100), b = rnorm(100) ) %>% readr::write_excel_csv(file = paste0("./sample_data/sample_", i, ".csv")) } # サンプルデータの作成(aの値の範囲を少し変えてある) for (i in 5:9) { data.frame( a = seq(from = 51, to = 150), b = rnorm(100) ) %>% readr::write_excel_csv(file = paste0("./sample_data/sample_", i, ".csv")) } # サンプルデータの読み込みと結合(aが一致しない場合は欠損値NAとなる) list.files("./sample_data", full.names = TRUE) %>% purrr::map_df(.f = function(x){readr::read_csv(x) %>% dplyr::mutate(file = basename(x))}) %>% tidyr::pivot_wider(names_from = file, values_from = b) >>427 一列目がすべて同じとのことなのであれば、427の続きで data2[ , -seq(3,199,2)] とすれば必要な列だけ抽出できると思います あとはマトリクスに変換してください join系でうまく行ってなかったのはリストの扱いがわからなかったからなのですね リストにアクセスする場合、 data[1]でなくdata[[1]]としてデータフレームを取り出します >>437 すみません。本題ではないと思うのですが、 > for (i in 1:10) { + data.frame( + a = seq(from = 1, to = 100), + b = rnorm(100) + ) %>% + readr::write_excel_csv(file = paste0("./sample_data/sample_", i, ".csv")) のところで エラー: Cannot open file for writing:* './sample_data/sample_1.csv' が出てしまいました・・・ >>438 先にディレクトリ(フォルダ)作らないとだめよ。 なお、そのコードはサンプル用のデータファイルを作るだけだから、既にデータファイルがあるなら作らなくても構いません。 >>438 繰り返しになりますが、427にあるようにbind_cols(data)でdata2まで作成できたのであれば、あとは437の式で重複するindex列を消せば抽出が完了します マトリクス形式にするにはas.matrix()を使います エラーが出たという433のコードについては439さんがもう回答してくれていますが念のため補足説明します そのコードは作業ディレクトリにsample_dataという名前のフォルダを作ってから試す必要があります Windowsならエクスプローラからフォルダを作成すれば良いでしょう Rでコードを実行するとそのフォルダの中にテスト用のcsvが作成されます(フォルダを覗いてみてください) そのテスト用csvを使った汎用的な方法が色々と紹介されていますので、どういう挙動をするか試してみてください あとは自分の実データに応用すればOKです cranにR 4.2.0がきてた Windowsもutf8に マルチバイトの苦労が減るかな 利点欠点は知らないけどとりあえず試しに移ってみては Vscodeに移行完了。 他の言語と環境同じなのは楽。簡単に複数のプロセスを走らせることができて便利。 環境変数とかの設定は大変。 RStudioもPositだかに名称変更して多言語対応を前面に出すらしいから色々迷うところだね tidyverseにlubridateが追加されたらしい これは便利かも ALTREPを勉強したノート: https://www.klgrth.io/paste/n3aj8 有効期限の選択肢が最長2日で、2日経つと消えちゃう。 一応HTMLということでアップしたが、HTMLとして表示できない。 ダウンロードしてブラウザーで開くと、多分、表示できると思う。 中身はRというより、殆どC++とJSで、偶に日本語が入るという感じになっている。 facet_wrapで強制的に任意の行×列にする方法がわからない ggh4xのfacet_wrap2はx軸が別の値に変わる致命的なバグがあって駄目だった cowplotとかで強引にやるしかないかな こういう事? [r - What's the difference between facet_wrap() and facet_grid() in ggplot2? - Stack Overflow](https://stackoverflow.com/questions/20457905/whats-the-difference-between-facet-wrap-and-facet-grid-in-ggplot2 ) ```{r} g <- ggplot2::ggplot (ggplot2::mpg, ggplot2::aes (displ, hwy)) g + ggplot2::facet_grid (cyl ~ class) g + ggplot2::facet_wrap (cyl ~ class) ``` 前回のALTREPから離れて配列型のイテレーターに進んでみた。 https://www.klgrth.io/paste/wcdat 今回のファイル形式はRmdにしてみた。 Rcppが使える環境なら動作するんじゃないかと思う 例えばmpgデータならcylでもclassでも常に3*3のサイズにfacetしたい(不足分はスペース) 結局patchworkのplot_layout()が良かった facet_wrap2()よりかなりめんどくさいが仕方ない RのヒープからC++のコンテナーで使うメモリを確保するようにした: https://www.klgrth.io/paste/ampar 今回はサードパーティのライブラリを使っているので、そのままでは コンパイルできないと思う。 * [GitHub - martinus/unordered_dense: A fast & densely stored hashmap and hashset based on robin-hood backward shift deletion](https://github.com/martinus/unordered_dense ) ヘッダーファイル一枚のライブラリなので、Rmdの中にヘッダーファイルを コピペすればコンパイルできるようになると思う。 前回のカウンターの続きで最終回 [zwzg9](https://www.klgrth.io/paste/zwzg9 ) 今回は`ankerl`に加えて`boost`も使っている。`boost`はヘッダー一枚 というわけにはいかないので、面倒かもしれない。使っているのは、 `hash_value`と`hash_combine`という2つの関数だけなので、適当な関数で 差し替えられると思う。`NA`と`NaN`の関係を調べている節は、 intelエンディアンのみの対応で、armエンディアンでは動作しないと思う。 実を言うと、`NaN`がいっぱいあることを知らなかった。1954は何の年だろう? ALTREPを勉強した時にかなりギットハブを徘徊した。そこが勉強のピークだっと 思う。今回の実装では、殆どRの勉強はせずに、知ってる関数だけを使っている。 1. `TYPEOF / ALTREP` 1. `Rf_xlength` 1. `DATAPTR_OR_NULL / DATAPTR / XXX_ELT` 1. `Rf_allocVector` 違うアプリケーション、例えば、評価系だと違う関数セットになって、 また勉強が必要になるかもね。 ファセットの問題は{facetious}のfacet_wrap_strict()でもできた こっちは問題なさそう [facet_wrap_strict](https://github.com/coolbutuseless/facetious )と `facet_wrap`を並べてみる。 ``` {r na_1954, dependson = ""} . = dplyr::mutate (mtcars, cyl = factor (cyl, lev = c (1 : 8, 1000L))) . = ggplot2::ggplot (., ggplot2::aes (mpg, wt)) . = . + ggplot2::geom_point () . + ggplot2::facet_wrap (~ cyl, nrow = 3, ncol = 3) . + facetious::facet_wrap_strict (~ cyl, nrow = 3, ncol = 3) . + ggplot2::facet_wrap (~ cyl, nrow = 3, ncol = 3, drop = F) ``` 神エクセル繋がりということで [A Bayesian probability worksheet](https://terrytao.wordpress.com/2022/10/07/a-bayesian-probability-worksheet/ ) その例で行くなら… cylはいじらずそのままでfacet_wrap_strictするのが求める図かな RのS5クラス使って継承する場合、コンストラクタ中で継承するクラスのコンストラクタを呼ぶことはできないのかな? つまりクラスAとクラスBを定義して、クラスBがクラスAを継承している場合、クラスBのコンストラクタ中でクラスAのコンストラクタを呼ぶ方法 ちょっとググったらS4クラスは出来そうなんだが、S5クラスのは解決策が見当たらんかった 別途自前でクラスAのコンストラクタ本体を定義して、クラスBのコンストラクタからそれを呼ぶしかないか?! dplyrを1.1.0にしてみた .byはまあまあ便利だけどtally()には使えなかった ナンノコッチャと思って調べたら [A Heisenbug lurking in async Python | Hacker News](https://news.ycombinator.com/item?id=34754276 ) 専門用語だったでござる。 [Heisenbug - Wikipedia](https://en.wikipedia.org/wiki/Heisenbug ) こうした言葉を思いつくのは才能なんだろうね。 Rだと簡単にハイゼンバグの例を作れた。 ``` {r a_10850, dependson = ""} uncertainly = with (new.env (), { plus = \(x) \(y) x + y x = 1L plus_1 = plus (x) x = 1000L plus_1 (2) |> print () }); ``` アレレ?ということでデバッグ文を入れてみる。 ``` {r a_29042, dependson = ""} principle = with (new.env (), { plus_debug = \(x) { cat ("Schr?dinger says ", x, "\n", sep = "") \(y) x + y } x = 1L plus_1 = plus_debug (x) x = 1000L plus_1 (2) |> print () }); ``` デバッグ文の有無で挙動が変わる。`uncertainly`の挙動は、バグではなく、 言語仕様だと思った方が良いと思う。 `uncertainly`の挙動を防ぐには関数`force`を使えば良い。 * [How to not fall into R's 'lazy evaluation trap](https://stackoverflow.com/questions/29084193 ) タイトルに"遅延評価"が入っているが、遅延評価自体の問題ではなく、 Rの実装に原因がある。Haskellも遅延評価だが、`uncertainly`のような 挙動が起きれば、Haskellはその存在意義を失う。 ``` ghci import Data.Functor.Identity (Identity (..)) :{ uncertainly :: (Monad f) => f Int uncertainly = do x <- pure 1 plus_1 <- pure (x +) x <- pure 1000 pure $ plus_1 2 :} runIdentity uncertainly flip ($) () uncertainly maybe 0 id uncertainly ``` <pre> 3 3 3 </pre> 言語仕様レベルでのトレードオフかもしれない。 * 局所最適化と大域最悪化 * 整合性とコンパイル時間 * 等など magrittrの例を少し改変 # fns first <- function(x){ message("first") invisible(x) } second <- function(x){ message("second") invisible(x) } # lazy NULL %>% first() %>% second() # eager NULL %!>% first() %!>% second() ビックリパイプは知らなかった。"抜かりなし"だね。ただし、ビックリパイプは ハイゼンバグとは関係ないかな。 ハイゼンバグから離れて、宇宙大戦争について書いてみたい。 `lapply`等のループ系の関数は、最終的にC関数`R_forceAndCall`を呼び出す。 * [apply.c](https://github.com/wch/r-source/blob/trunk/src/main/apply.c ) Rは、"フォースと共にあらんことを"などと呑気なことを言っている場合ではなく、 フォースと共にあらねば死んでしまう。 前回のコード`uncertainly`は変数`x`の使い回しが敗因だが、ループ系の関数 では変数を使い回すしかない。そのために、ループ系の関数ではフォースが必須 になっている。98パーセントぐらいの使用例では、フォース抜きでも動作する (当社調べ)が、残りの2パーセントでコケる。僅か2パーセントでも、 プログラミングをギャンブルにしないためには、穴を塞いでおく必要がある。 上のコード`uncertainly`の挙動を"プロミスの罠"と書くことにする。 ここでの"プロミス"は、JSの"プロミス"ではなくて、Rでの"変数"の 実装方法を指す。 * [6 Functions | Advanced R](https://adv-r.hadley.nz/functions.html ) この記事には次の一節がある。 > You cannot manipulate promises with R code. Promises are like a **quantum state:**: ... 多分、ここでの"量子状態"はハイゼンバグと同じ現象を指しているんだと思う。 観測すると、状態が変化してしまう。 関数`do.call`自体はプロミスの罠と関係しないと思う。 ``` {r a_28697, dependson = ""} do.call (`+`, list (1L, 2L)) ``` Pythonだと次のコードに対応する。 ``` {python a_10786, dependson = ""} (lambda x, y: x + y) (* range (1, 3)) ``` 殆どのプログラミング言語で、関数の引数リストは [一級市民](https://en.wikipedia.org/wiki/First-class_citizen ) でないように思う。一級市民でない代わりに、一級市民の配列からの変換が 用意されている。Pythonでは`*`という関数がビルトインで用意されている。 Rでは引数リストを直接作れるかもしれないが、`do.call`で配列から 変換するのが一般的だと思う。 自分の知る限り、プロミスの罠にハマるのは次のパターンに限られる。 ~~~ {.r} x = "hello" g = f (x) x = "world" g () ~~~ "関数を返す関数"`f`に、"変数"`x`を代入した時にプロミスの罠が可能性が 出てくる。Rの評価戦略はHaskellと同じ [コールバイニード](https://en.wikipedia.org/wiki/Evaluation_strategy ) に分類されている。上の例では、次の場合にプロミスの罠が発生する。 1. 関数`f (x)`の中で引数`x`がニードにならず、 1. 返り値の関数`g`に渡される。 次の例はプロミスの罠が発生する。 ``` {r a_20914, dependson = ""} rude = with (new.env (), { const = \(x) \(...) x x = "hello" g = const (identity (x)) x = "world" g () |> print () }); ``` 関数への代入はニードでないので(by definition)、関数適用の連鎖 `const (identity (x))`の中にはニードがない。そのために、最終的に`g ()`が コールされるまで引数`x`は評価されない。コード`g ()`がコールされた時点で 引数`x`の定義を探しに行くので、プロミスの罠にハマる。 コールバイニードを実現するためには、コンパイラーが内部的に次のように 書き換える必要があると思う。 ``` {r a_1451, dependson = ""} polite = with (new.env (), { const = \(x) \(...) x x_1 = "hello" g = const (identity (x_1)) x_2 = "world" g () |> print () }); ``` [静的単一代入](https://en.wikipedia.org/wiki/Static_single-assignment_form ) と呼ばれる操作と同じだと思う。この書き換えはループ処理では必須になる。 コード`rude`はアカン奴として切り捨てることができるかもしれないが、 ループ処理では変数の多重定義が避けられない。Rでは、何らかの理由で、 この書き換えが省略されているために、プロミスの罠が発生するのかな?と 思っている。 コード`rude`からコード`polite`への書き換えを上とは異なる形で行ってみる。 Rではできないので、Pythonを使う。 ``` {python a_14213, dependson = ""} ##| cache: false # def none (): id = lambda a: a const = lambda x: lambda * y: x def lhs (): x = "hello" g = const (id (x)) x = "world" return g () def rhs (): return (lambda x: (lambda g: (lambda x: g ()) ("world")) (const (id (x)))) ("hello") assert lhs () == rhs () none (); ``` 関数`lhs`から関数`rhs`への書き換え規則は次のようになる。[rule]{#rule} ~~~ {.python} y = f (x) z = g (x, y) rest (x, y, z) == (lambda y: z = g (x, y) rest (x, y, z) ) (f (x)) == (lambda y: (lambda z: rest (x, y, z)) (g (x, y))) (f (x)) ~~~ この書き換え規則を適用した結果、関数`rhs`では、静的単一代入への変換は、 [ラムダ計算](https://en.wikipedia.org/wiki/Lambda_calculus )のアルファ変換 に置き換わる。この意味では、静的単一代入とアルファ変換は同じことになる。 実際のコードでは、分岐、ループ、副作用が入ってくるので、こんな単純は 話では済まないと思うが、極度に単純化すると、静的単一代入とアルファ変換は、 単に、方言の違いということになる。 Rの場合は、書き換え規則[rule](#rule)は等価な変換にならない可能性がある。 特に、プロミスの罠が現れた場合は、挙動が変わると思う。 余談だが、書き換え規則[rule](#rule)をモナドに拡張したものがHaskellの [do-記法](https://en.wikibooks.org/wiki/Haskell/do_notation )になっている。 ~~~ {.python} y <- f (x) z <- g (x, y) rest (x, y, z) == kleisli-extension (lambda y: z <- g (x, y) rest (x, y, z) ) (f (x)) ~~~ モナドが恒等関数の場合に[rule](#rule)に一致するので、do-記法を [rule](#rule)の拡張と言って差し支えないと思う。 おしまい [トランスフォーマー](https://en.wikipedia.org/wiki/Transformer_ (machine_learning_model)) についての記事: オーサーヒアーかな? * [Coinductive guide to inductive transformer heads](https://news.ycombinator.com/item?id=34970877 ) 中を読んでいないが、アテンション機構を [ホップ代数](https://en.wikipedia.org/wiki/Hopf_algebra )で説明、もしくは、 ホップ代数を使ってアテンション機構の類似物を導出するという話のようだ。 この手の救世主的な話はがっかりすることが多いので、過大な期待は禁物かも。 g = do.call(const, list(x)) ではどうだろう RのS5クラス定義(setRefClass)ってsave、loadでシリアライズ、デシリアライズ出来ないのか デシリアライズされたメンバ関数にアクセスできない ## Rにおけるオブジェクト指向への取り組み 趣旨は新規フレームワークR7の紹介だが、既存のフレームワークの俯瞰図 としても使いやすいと思う。 * [What is R7? A New OOP System for R](https://www.jumpingrivers.com/blog/r7-oop-object-oriented-programming-r/ ) ## オブジェクト指向とは? 次の記事を取り上げてみる。 * [Object Oriented Programming Features of Rust](https://doc.rust-lang.org/book/ch17-00-oop.html ) 内容については賛否両論だと思う。しかし、箇条書きできるということが、 オブジェクト指向が単一の概念でないことを示しているように思う。 ## S3はオブジェクト指向か? 比較のために、S3のキーワード`default`を使わないで書く。 ``` {r a_3770, dependson = ""} none = with (new.env (), { to_string = \(x, ...) UseMethod ("to_string", x) as_string = \(x, ...) UseMethod ("as_string", x) to_string.Base = \(x, ...) as_string (x, ...) as_string.Base = \(x, ...) to_string (x, ...) to_string.Derived = \(x, ...) paste0 ("hello ", class (x)) as_string.Derived = as_string.Base structure (identity, class = "Derived") |> as_string () |> print () }); ``` `Base`クラスは、関数`to_string`か関数`as_string`の"どちらか1つを実装せえ"と言っている。 似た内容をPythonで書いてみる。 ``` {python a_8560, dependson = ""} class Base: def to_string (self): return self.as_string () def as_string (self): return self.to_string () class Derived (Base): def to_string (self): return "hello " + type (self).__name__ Derived ().as_string () ``` Haskellで書いてみる。 ``` haskell class Base_ a where to_string :: a -> String to_string = as_string as_string :: a -> String as_string = to_string instance Base_ () where to_string :: () -> String to_string = ("hello " ++) . show as_string () ``` <pre> "hello ()" </pre> ワンセットの関数群をインターフェースと書く。 1. Pythonのクラスでは、 * 構造体の定義とインターフェースの定義と実装が同時に行われる。 1. HaskellのクラスやRのS3では、 * 構造体とインターフェースは別々に定義される。 * インターフェースの実装は既存の"構造体"に対して行われる。 HaskellやRの"構造体"には、通常の構造体に加えて、関数も含まれる。 構造体とインターフェースが別々に定義される点で、RのS3はHaskellのクラスに 近いように思う。 おしまい join_by()が便利で感動した 気象データのマージがすごくシンプルになった 任意の地点・期間における測定値についての風向風速などの気象の影響をみるため近傍アメダスの気象データをマージする際に、便利になったと実感しました 具体的にはマージの際にbetween(時刻,)の条件が追加できるようになったのが大きいです 処理が重いかコードが長いかだった上のマージ処理が、join_by(地点名, between(time, start_time, end_time))で済むようになりました [非等価ジョイン](https://en.wikipedia.org/wiki/Relational_algebra )は 最近になって関数`dplyr::xxx_join`に導入されたらしい。 * Rで非等価結合 (2) * [dplyr 1.1.0: Joins](https://www.tidyverse.org/blog/2023/01/dplyr-1-1-0-joins/ ) [ラッキー](https://www.youtube.com/watch?v=xBJMr1v5Zuw )。 関数`vctrs::vec_locate_matches`には、以前から限定された形の非等価ジョイン のオプションがあったので、実用的な形にするのに時間がかかったのかもしれない。 リスト内包表記が使いたいのでPythonで書く。 ``` {python a_32372, dependson = ""} #| cache: false # def slow_matches (pred, zipper): #{ def go (a, b): #{ a = pandas.DataFrame ({"" : a}).groupby ("").groups b = pandas.DataFrame ({"" : b}).groupby ("").groups key = ((j, k) for j in a.keys () for k in b.keys () if pred (j, k)) val = (zipper (a [j], b [k]) for (j, k) in key) val = zip (* itertools.chain (* val)) return (numpy.fromiter (val, dtype = numpy.int32) for val in val) #} return go #} def none (): #{ a = numpy.random.choice (8, 2 * 8, replace = True) b = a + 4 def go (pred, zipper): #{ (j, k) = slow_matches (pred, zipper) (a, b) out = pandas.DataFrame ({"a" : a [j], "b" : b [k]}) print (out) #} go (lambda a, b: a == b, itertools.product) go (lambda a, b: a == b, lambda a, b: itertools.islice (zip (a, b), 1)) go (lambda a, b: abs (a - b) <= 1, itertools.product) go (lambda a, b: abs (a - b) <= 1, lambda a, b: itertools.islice (zip (a, b), 1)) #} none (); ``` `dplyr::xxx_join`の数々のオプションは理解できていないが、`NA`絡み以外は、 コールバック`pred`と`zipper`でかなりカバーしてるんじゃないかと思う。 ジョインの操作自体は明朗会計だが、実行時性能がボッタクリ価格に転嫁される というプログラミング社会が抱える構造的な問題が反映されている。 おしまい 一応Rスレだしコード一式を載せるならRにしてほしいな 正直、書き方を思い出すためだけに、話題にかこつけて、わざと異なる プログラミング言語で書くこともある。しかし、今回は純粋に、処理の道筋を 単刀直入な形で表せるプログラミング言語を選択した。 [線形回帰](https://en.wikipedia.org/wiki/Linear_regression )を例にとる。 Rで線形回帰はとても書きやすい。しかし、線型回帰をRのイディオムとして 覚えるよりも、線形代数として覚えた方が長い目で見た時に利益になると思う。 同じように、等価ジョインは [プルバック](https://en.wikipedia.org/wiki/Pullback_ (category_theory)#Sets) に対応する。等価ジョインを、Rのイディオムとして覚えるよりも、 単純な`for`ループの形で覚えた方が長い目で見た時に利益になると思う。 リスト内包記法は`for`ループを簡潔に書くための記法になっている。 プルバックの話の続きをアップした。 * [JSFiddle](https://jsfiddle.net/ytkhdpcq/show ) プルバックのカリー化の話になっている。プルバックの話はこれでお終い。 洗練された話は次の記事にある。 * [locally cartesian closed category in nLab](https://ncatlab.org/nlab/show/locally+cartesian+closed+category ) 表計算の立場で見ると、"何が問題なのかわからない"問題かもしれない。 ある規則を満たす表でサマリーするとプルバックのカリー化になる。 4.3.0入れてみた パイプがまたひとつ便利になっていい感じ Windowsでパスが長くても通るようになったのも嬉しい どなたか教えてくだされ 対照実験の医療統計をしたいと思って 傾向スコアマッチングにしようと思ってます 肺がん手術を実施した人たち200症例のうち10人だけ 手術後に合併症を起こしていました この場合 「合併症を起こしたグループ」・・・10人 「合併症を起こさなかったグループ」・・・190人になり 患者の背景で傾向スコアマッチングをキャリパー幅0.2でペア分けすると 各群8人ずつの比較になってしまいます。 これで優位差を出す統計するには数が少なすぎるのであまりにも乱暴な統計になりますか? それとも元々のサンプルが200人から選りすぐりのマッチングさせた8人ずつを選んでいるので8人ずつと少なくても問題ないですか? 共変量のバイアスがより少ないサンプルが得られますが、サイズが8人の二群比較なのは変わらないんじゃないかと思います >>493 つまりサイズが8人ずつの比較では 優位差を言うには弱いと思いますか? >>494 検定はサンプルサイズ込みの判定なので8人であっても二群に有意な差があったと言って差し支えありません 有意水準が甘ければ弱いとは感じますが 一方で、検定は標本のバイアスが除去できていることを保証してくれません スレッドの質をガクッと下げるが、Rの非標準評価で遊んでみた。 * [JSFiddle](https://jsfiddle.net/fo1qjt87/show ) 4.3.0のWindows版はdir()にバグあるね リストがまともに取得できない 4.2.3に戻した tibbleを利用してデータ探索しているんだけど、2つの変数で条件つけて散布図つくるにはどんな手順踏んだらいいのか… 変数1をプロットの色にして変数2をサイズにするとかじゃだめなの? それも一つの方法ですね。ありがとう ただ、目的以外の条件をプロットしたくないのです。 tibbleの構成としては、観測が4時点、変数がたとえば6個(うち1つは個体識別番号で、変数としてはA-Eの5種類あるとしましょう)のような感じで、 散布図のx軸に「観測時点1における全個体の変数A」 y軸に「観測時点4における全個体の変数E」 をプロットしたいのです。 個体識別番号をキューにして、異なる時点の異なる種類の変数を散布図にするとでもいいましょうか… その要件だと以下のような感じでいけると思います library(tidyverse) df <- tibble( time = rep(1:4, each = 5), id = rep(1:5, 4), a = rnorm(20), b = rnorm(20), c = rnorm(20), d = rnorm(20), e = rnorm(20) ) df |> pivot_wider(id_cols = id, names_from = time, values_from = c(a, e)) |> ggplot(aes(a_1, e_4)) + geom_point() ありがとうございます! それをヒントに試してみます。 >>507 要件の理解が違っていなければ散布図作成まで行きますよ time1でのaの値とtime4でのeの値をid毎にプロットしています >>508 失礼しました こちらのビューワの問題で、ブラウザでみたらスクリプトの最後まで表示されました💦 >>508 教えてもらったスクリプトでドンピシャでした、ありがとうございます。 どなたか統計素人の私に教えてください 医療でA群100名とB群100名をいろんな観点から 比較するとします 背景を揃えるために傾向スコア分析で キャリパー0.2で取って30ペア(合計60人)を抽出しました 例えば ①アミラーゼの値とかを統計にかけたいときは、その60人のアミラーゼのデータをまずは正規分布かどうかをShapiro-wilk検定で測って、0.05を上回っていたら「正規分布」と見なして、今度はこれら60人のアミラーゼのデータが等分散かどうかを確認するためにf検定を行なって等分散だったらt検定(Student's t-test)、違ったらWelch's t-testで優位差があるかどうかを調べる。 もしShapiro-wilk検定が0.05未満で非正規分布だったときはMann-whitney's U検定で優位差があるかどうかを調べる。 ↑Q1. この認識であってますか? Q2. この正規分布を計るのはペアを作った後のデータ(n=60)で、正規分布かどうか?等分散かどうか?を見ますか? それともマッチング前のデータ(n=200)で正規分布かどうか?等分散かどうか?を確認するべきですか? Q3. アミラーゼじゃなくて他の連続変数の項目(BMIだったり、血圧だったり、白血球数だったり)の優位差を調べる場合も 全て上の流れで一つ一つの項目ごとで正規分布か?等分散か?など確認していって適宜、該当する算出法を項目ごとに採用して優位差を測る必要がありますか? もし良ければ教えていただけるとありがたいです 素人だと思ううちは傾向スコアには手を出さないほうがいいような気がする >>512 上記の場合はどういうかんじでやるのが 最適解になるのかだけでも教えてもらえないでしょうか? writexlパッケージで保存すると時刻のタイムゾーンは反映されないのね ggraggedパッケージがなかなか便利 facet_grid()を詰めて並べられる R4.3.1にした dir()関連が直ってるようでひと安心 標準パイプの新機能もようやく使える ggplot2のカラーパレットのデフォルトが変わったりしました? 棒グラフの色の割り当てが変わった気がする… >>518 変わってないはず ただggplot2のデフォルトはグループ数に応じて関数で色を割り当ててるのでグループ数が異なると違うパレットのように感じるとは思います 株価データを取得したいがためにRを弄り始めたけど中々難しいね ググってコピペしたコードを一つ一つパッケージと関数の挙動を調べているけど やりたいこと全てができるようになるまでは遠そうだ read.cgi ver 07.5.4 2024/05/19 Walang Kapalit ★ | Donguri System Team 5ちゃんねる