こういうデータをメタアナリシスで解析する。 testが処理した群、contが対照群 testで病気になった人がtest_posで 全体がtest_total id test_pos test_total cont_pos cont_total 1 test1 5 177 6 144 2 test2 5 132 4 155 3 test3 2 23 33 47 4 test4 6 76 …

さて、よく解析して、p値出して、その通りでしたー、みたいなことを書いてしまいますが、本当に重要なことは、その解析した際に推定された値、変数だったりするわけではないでしょうか。 ということで、より噛み砕こうというのが今回です。 fitted()は、当て…

前回のおさらいから odds ratioの比較に、ただただ区間推定して比較してみたが、vcdパッケージには、woolf検定（層ごとにオッズ比は同じという帰無仮説を検定）をしてくれるパッケージもある。 library(vcd) woolf_test(data2) Woolf-test on Homogeneity of…

コンサルで多次元の分割表を考える機会があったので、メモを残す。 多次元の分割表を考えるうえで陥ってしまいがちなのが、「シンプソンのパラドックス」。要するに、多次元の分割表を勝手に統合したり、分割したりすると、誤った結果を導いてしまうかもしれ…

基本的には、データから１つデータを取り出して、擬似値を計算する。それから、その擬似値の集合（データと同じ分だけできる）から統計量を計算する。 jack.mean { Mean.obs Pseudovalues for (i in 1:length(data)) { Data.minus.one Mean.minus.one Pseudo…

data.frameのカテゴリー数が調べたくなったら、 nlevels(データの列） で調べることができる。

offset = 数字 で棒グラフの基点を設定できる。

http://glmm.wikidot.com/pkg-comparison

Rでモデルを２つ作って，その２つのモデルの当てはまりを検定で（どちらかを帰無仮説として）評価したい場合よく anova(model1, model2, ...) を使います。この時，modelがどのような分布に従うかによって test = ""の中身が変わる。それを知るためには， su…

でも！あまり好まれていないglmmPQLではできた！idを50にしても。 > fglmmPQL fglmmPQL Linear mixed-effects model fit by maximum likelihood Data: dat2 Log-likelihood: NA Fixed: y2 ~ 1 (Intercept) 0.9502553 Random effects: Formula: ~1 | id (Inte…

こちらの続き。 random effectをどのように評価するか？ > leaf disc concleaf) > dataleaf=factor(leaf),disc=factor(disc),conc) > summary(data) leaf disc conc 1:4 1:4 Min. :0.007103 2:4 2:4 1st Qu.:0.084563 3:4 3:4 Median :0.280326 4:4 4:4 Mean…

階層構造と時系列構造のあるデータを解析することになり，忘れかけていた（？）混合モデルをココにメモする。 期待平均平方(expected mean squares)；分散分析的に変量（ランダム）効果を扱おうとする場合，ANOVA表の（調整済みの）平均平方の項にoutputされ…

もしかして，説明変数に相関構造のある場合で，目的変数が順序尺度の場合はGEEなるものを使えばいいのでは…。。

目的変数が３つ以上あるとき，多項ロジット分析を行うが，目的変数が順序尺度のときは累積ロジットモデルなどをよく使う。累積ロジットモデルはある基準のカテゴリー間に確率を算出する。

WinBUGSのインストール。 ベイズ統計データ解析で勉強中。経験ベイズとか。

leafleaf) dataleaf=factor(leaf),disc=factor(disc),conc) ここで，４つの葉っぱからそれぞれ４つの穴を開けて，カルシウム濃度を調べる。で葉っぱごとにどれだけばらつくか（分散）を調べたい場合を考えよう。どのようなときにこの解析を使うかといわれれ…

最近暑さのせいか寝不足だ。そして，明日から盛岡に行きます。 成功を祈る。 ∀はfor all

先日のやり方は間違ってました。 分割表の確率の計算方法が，期待値算出になってました。 正しくは，実際に得られたデータの比率を計算するだけです。 まず，以下の関数を再度定義。 sqr.obsで，先日の例をおさらい。 > data [,1] [,2] [1,] 5 7 [2,] 13 12 …

前回の続き。 > head(data) #またデータを作り直した。 leaf disc conc 1 1 1 0.060311173 2 1 2 0.225865086 3 1 3 0.003394487 4 1 4 0.020026101 5 2 1 0.004372889 6 2 2 0.393164132 > data$leafleaf) > data$leaf [1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4…

非常に参考になるpdfがネットに落ちています。勝手にリンクします。勉強してください。

さて，ここで改めてランダム効果について記述する。 ランダム効果は，一昔前の統計の教科書では「変量効果」などと記述されている。英語では，random effectである。で，どのような時にこの「ランダム効果」を使うのか？以下はこちらやこちらを参考にしまし…

これからますます使われるであろう，ランダム効果について参考までに自分なりの解析法を提示します。 こちらにも似た記述はあります。データはこちらの本を参考にさせていただきました。解析にはR(2.7.1)を用います。 > leaf disc conc dataleaf,disc,conc) …

barnacle z z[1]=0.1 maxt=1000 for(t in 1:(maxt-1)){ W*1 z[t+1] } plot(z,type="l",xlab="generation",ylab="swichpoint",ylim=c(0,1)) text(500,0.1,paste("z=",round(z[1000],4),"equ=",round*2,cex=1.1) } *1:1-z[t])/2*a1+(1-z[t])/2*a0) deltaz *2:a…

イントロでお腹いっぱい。 （正規）線形モデルでは，残差平方和（RSS）は逸脱度。 残差標準誤差は，√（RSS/Df）。新しい関数 termplot() halfplot() MASS中のrlm():ロバスト回帰（外れ値に強くなる） lm()でweight指定で重み付け。

パッケージneetの中の関数multinom、 パッケージVGAMの中の関数vglm を用いる

Rのglm関数などによってデータにモデルを当てはめた際、よく行う（？）summary関数による回帰係数の検定表。 実は今までsummary()で出力されるp値は検定の多重性を考慮していないと思っていた。しかし、summary()で出力される推定値のwald検定のp値は、その…

ordered()

set.seed(整数値) 乱数を発生させる際に、再現性のある乱数を作らせるため、乱数値を（整数で）規定しておける。

interaction.plot(x軸,カテゴリー別,y軸)で書ける。