Rの統計解析関数Tips
をテンプレートにして作成
[
トップ
|
Tips紹介
|
初級Q&A
|
R掲示板
|
日本語化掲示板
|
リンク集
]
[
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
COLOR(red){SIZE(25){Rの統計解析関数Tips}}
R に関する Tips で欠けているもの、そしておそらく一番困難(...
#contents
~
*lm 関数を用いた線形モデルの当てはめ結果の要約 (r-help の...
> df1 <- data.frame(x=1:6, y=rep(1:3, 2))
> fit <- lm(y~x, df1)
> Sum <- summary(fit)
# 当てはめオブジェクト fit の要約(普通ユーザーはこの情報...
> Sum
Call:
lm(formula = y ~ x, data = df1)
Residuals:
1 2 3 4 5 6
-0.4286 0.3429 1.1143 -1.1143 -0.3429 0.4286
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2000 0.8177 1.468 0.216
x 0.2286 0.2100 1.089 0.338
Residual standard error: 0.8783 on 4 degrees of freedom
Multiple R-Squared: 0.2286, Adjusted R-squared: 0.03...
F-statistic: 1.185 on 1 and 4 DF, p-value: 0.3375
# 当てはめオブジェクトの要約が持つ情報(リスト形式になっ...
> attributes(Sum)
$names
[1] "call" "terms" "residuals" "coe...
[5] "sigma" "df" "r.squared" "adj...
[9] "fstatistic" "cov.unscaled"
# 回帰係数と関連情報を表形式で得る(summary(fit) が表示す...
> coefficients(Sum) # coefficients 関数は線形回帰当ては...
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2000000 0.8176622 1.467599 0.2161194
x 0.2285714 0.2099563 1.088662 0.3375019
> Sum$coefficients # これでも良い
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2000000 0.8176622 1.467599 0.2161194
x 0.2285714 0.2099563 1.088662 0.3375019
参考までに述べると lm 関数が返す当てはめオブジェクトは実...
> str(fit) # str 関数で当てはめオブジェクトの内容を見る
List of 12
$ coefficients : Named num [1:2] 1.200 0.229
..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
$ residuals : Named num [1:6] -0.429 0.343 1.114 -...
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ effects : Named num [1:6] -4.89898 0.95618 1.2...
..- attr(*, "names")= chr [1:6] "(Intercept)" "x" "" "...
$ rank : int 2
$ fitted.values: Named num [1:6] 1.43 1.66 1.89 2.11 2....
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ assign : int [1:2] 0 1
$ qr :List of 5
..$ qr : num [1:6, 1:2] -2.449 0.408 0.408 0.408 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:6] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:2] "(Intercept)" "x"
.. ..- attr(*, "assign")= int [1:2] 0 1
..$ qraux: num [1:2] 1.41 1.19
..$ pivot: int [1:2] 1 2
..$ tol : num 1e-07
..$ rank : int 2
..- attr(*, "class")= chr "qr"
$ df.residual : int 4
$ xlevels : list()
$ call : language lm(formula = y ~ x, data = df1)
$ terms :Classes 'terms', 'formula' length 3 y ~ x
.. ..- attr(*, "variables")= language list(y, x)
.. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. ..$ : chr "x"
.. ..- attr(*, "term.labels")= chr "x"
.. ..- attr(*, "order")= int 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=length 33 <environment>
.. ..- attr(*, "predvars")= language list(y, x)
$ model :`data.frame': 6 obs. of 2 variables:
..$ y: int [1:6] 1 2 3 1 2 3
..$ x: int [1:6] 1 2 3 4 5 6
..- attr(*, "terms")=Classes 'terms', 'formula' length...
.. .. ..- attr(*, "variables")= language list(y, x)
.. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. .. ..$ : chr "x"
.. .. ..- attr(*, "term.labels")= chr "x"
.. .. ..- attr(*, "order")= int 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=length 33 <environme...
.. .. ..- attr(*, "predvars")= language list(y, x)
- attr(*, "class")= chr "lm"
> str(Sum) # Sum の持つ内部情報の詳細を見る
List of 11
$ call : language lm(formula = y ~ x, data = df1)
$ terms :Classes 'terms', 'formula' length 3 y ~ x
.. ..- attr(*, "variables")= language list(y, x)
.. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. ..$ : chr "x"
.. ..- attr(*, "term.labels")= chr "x"
.. ..- attr(*, "order")= int 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=length 34 <environment>
.. ..- attr(*, "predvars")= language list(y, x)
$ residuals : Named num [1:6] -0.429 0.343 1.114 -...
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ coefficients : num [1:2, 1:4] 1.200 0.229 0.818 0.210...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "(Intercept)" "x"
.. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "...
$ aliased : Named logi [1:2] FALSE FALSE
..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
$ sigma : num 0.878
$ df : int [1:3] 2 4 2
$ r.squared : num 0.229
$ adj.r.squared: num 0.0357
$ fstatistic : Named num [1:3] 1.19 1.00 4.00
..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
$ cov.unscaled : num [1:2, 1:2] 0.8667 -0.2000 -0.2000...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "(Intercept)" "x"
.. ..$ : chr [1:2] "(Intercept)" "x"
- attr(*, "class")= chr "summary.lm"
# 当てはめオブジェクト fit に基づく二次的解析の例(分散分...
> aovfit <- aov(fit)
> aovfit
Call:
aov(formula = fit)
Terms:
x Residuals
Sum of Squares 0.9142857 3.0857143
Deg. of Freedom 1 4
Residual standard error: 0.87831
Estimated effects may be unbalanced
> summary(aovfit)
Df Sum Sq Mean Sq F value Pr(>F)
x 1 0.91429 0.91429 1.1852 0.3375
Residuals 4 3.08571 0.77143
> str(aovfit)
List of 12
$ coefficients : Named num [1:2] 1.200 0.229
..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
$ residuals : Named num [1:6] -0.429 0.343 1.114 -...
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ effects : Named num [1:6] -4.89898 0.95618 1.2...
..- attr(*, "names")= chr [1:6] "(Intercept)" "x" "" ""...
$ rank : int 2
$ fitted.values: Named num [1:6] 1.43 1.66 1.89 2.11 2....
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ assign : int [1:2] 0 1
$ qr :List of 5
..$ qr : num [1:6, 1:2] -2.449 0.408 0.408 0.408 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:6] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:2] "(Intercept)" "x"
.. ..- attr(*, "assign")= int [1:2] 0 1
..$ qraux: num [1:2] 1.41 1.19
..$ pivot: int [1:2] 1 2
..$ tol : num 1e-07
..$ rank : int 2
..- attr(*, "class")= chr "qr"
$ df.residual : int 4
$ xlevels : list()
$ call : language aov(formula = fit)
$ terms :Classes 'terms', 'formula' length 3 y ~ x
.. ..- attr(*, "variables")= language list(y, x)
.. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. ..$ : chr "x"
.. ..- attr(*, "term.labels")= chr "x"
.. ..- attr(*, "order")= int 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=length 35 <environment>
.. ..- attr(*, "predvars")= language list(y, x)
$ model :`data.frame': 6 obs. of 2 variables:
..$ y: int [1:6] 1 2 3 1 2 3
..$ x: int [1:6] 1 2 3 4 5 6
..- attr(*, "terms")=Classes 'terms', 'formula' length...
.. .. ..- attr(*, "variables")= language list(y, x)
.. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. .. ..$ : chr "x"
.. .. ..- attr(*, "term.labels")= chr "x"
.. .. ..- attr(*, "order")= int 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=length 35 <environme...
.. .. ..- attr(*, "predvars")= language list(y, x)
- attr(*, "class")= chr [1:2] "aov" "lm"
* 線形モデル当てはめ関数におけるモデル公式 (2004.01.05) [...
lm (線形回帰)、aov (分散分析)、glm (一般化線形回帰) 等の...
|y ~ ...| 左辺の目的変数を、右辺の説明変数の線形式で線形...
|+ A| 説明変数 A を加える|
|- A |説明変数 A を除く|
|+ 0 または - 1 | 切片(定数)項の除外|
|A:B| 二つの説明変数 A, B の交互作用項 (A と B の積からな...
|A * B | A + B + A:Bと同じ |
|(...)*(...) | 右括弧内の変数と左括弧内の変数の掛け合わせ...
| |ただし A^2, A^2*B 等はそれぞれ A, A*B と解釈|
|(A+B)^2| (A+B)*(A+B) と同じ。A + B + A:B の意味|
| A %in% B | 今一意味がわからない?|
より一般に A:B:C や (A+B+C)^3 といった高次交互作用項も同...
y を目的変数、A, B, C 等を説明変数・因子とすると
| y ~ A| y を A と切片(定数)項で単回帰|
| y ~ A - 1 | 切片項のない、原点を通る単回帰|
| y ~ A + 0 | 同上|
| 0 + y ~ A | 同上|
| y ~ A + B| y を A, B と切片(定数)項で重回帰|
| y ~ A + B - 1| y を A, B で重回帰。切片(定数)項は含まな...
| y ~ A + B + A:B |A,B とそれらの二次交互作用項で重回帰,...
| y ~ (A+B)^2 |同上|
| y ~ (A+B)^2 - 1|A,B とそれらの二次交互作用項で重回帰。...
| y ~ (A+B+C)^2 |A,B,C とそれらのすべての二次交互作用項...
| y ~ A + B + C + A:B + A:C + B:C |同上|
| y ~ (A+B+C)^2 - B:C | y ~ A + B + C + A:B + A:C の意味|
| y ~ (A+B+C)^3 |A,B,C とそれらのすべての二次及び三次交...
| y ~ A + A %in% B | y ~ A + A:B の意味|
モデル式には説明変数の関数を使うことが出来る。
| y ~ A + log(B)| y を A と B の対数値で説明|
| log(y) ~ exp(A) + log(B)| log(y) を exp(A) と log(B) で...
| y ~ A/B| y を A と B の比で説明|
モデル式で特別の意味を持つオペレータ(例えば +,-,*,^)を、...
|y ~ A + I(B+C) | y を A と B+C (B と C の和) で説明|
|y ~ A + I(B*C) | y を A と B*C (B と C の積) で説明、A +...
|y ~ A + I(B^2) | y を A と B^2 (B の自乗値) で説明|
*非線形最小自乗法 nls() 関数 [#iad5cef4]
非線形モデルの非線形最小自乗推定専用の関数で、独自の最適...
- COLOR(magenta){用法}
nls(formula, data = parent.frame(), start, control = nl...
algorithm = "default", trace = FALSE, subset,
weights, na.action)
- COLOR(magenta){引数}
--COLOR(red){formula} 変量とパラメータを含む非線形モデル...
--COLOR(red){data} オプションのデータフレームで、その中で...
--COLOR(red){start} 初期推定値の名前付きリスト、もしくは...
--COLOR(red){control} 制御仕様を決めるオプションのリスト...
--COLOR(red){algorithm} 使用されるアルゴリズムを指定する...
--COLOR(red){trace} 論理値で、繰り返しの過程の記録を表示...
--COLOR(red){subset} オプションのベクトルで、当てはめの過...
--COLOR(red){weights} オプションの数値ベクトルで、(固定さ...
--COLOR(red){na.action} データが NA を含む時に、どうする...
-COLOR(magenta){詳細}
-- nls を人工的な "残差がゼロ" のデータに使ってはならない...
-- nls オブジェクトは典型的には当てはめモデルオブジェクト...
-COLOR(magenta){返り値} 次の成分からなるリスト
--COLOR(red){m} モデルを含む nlsModel オブジェクト
--COLOR(red){data} データ引数として nls に引き渡された公...
**例1 [#yca8db2a]
データフレーム DNase は成分 Run (順序つき因子, 二因子)、c...
data( DNase )
DNase1 <- DNase[ DNase$Run == 1, ] # データの一部を取り...
## 自動スタートモデルを使う。SSlogis はロジスティックモ...
fm1DNase1 <- nls( density ~ SSlogis( log(conc), Asym, xm...
summary( fm1DNase1 ) # 当てはめ結果の要約
Formula: density ~ SSlogis(log(conc), Asym, xmid, scal) ...
Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
Asym 2.34518 0.07815 30.01 2.17e-13 ***
xmid 1.48309 0.08135 18.23 1.22e-10 ***
scal 1.04146 0.03227 32.27 8.51e-14 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1...
Residual standard error: 0.01919 on 13 degrees of freedom
Correlation of Parameter Estimates: 推定パラメータの相関値
Asym xmid
xmid 0.9868
scal 0.9008 0.9063
## 条件付き線形性を使う
fm2DNase1 <- nls( density ~ 1/(1 + exp(( xmid - log(conc...
data = DNase1,
start = list( xmid = 0, scal = 1 ),
alg = "plinear", trace = TRUE )
0.7139315 : 0.000000 1.000000 1.453853 # 当てはめ途中結...
0.1445295 : 1.640243 1.390186 2.461754 # 残差平方和と推...
0.008302151 : 1.620899 1.054228 2.478388
0.004794192 : 1.485226 1.043709 2.347334
0.004789569 : 1.483130 1.041468 2.345218
0.004789569 : 1.483090 1.041455 2.345180
summary( fm2DNase1 ) # 当てはめ結果の要約
Formula: density ~ 1/(1 + exp((xmid - log(conc))/scal)) ...
Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
xmid 1.48309 0.08135 18.23 1.22e-10 ***
scal 1.04145 0.03227 32.27 8.51e-14 ***
.lin 2.34518 0.07815 30.01 2.17e-13 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1...
Residual standard error: 0.01919 on 13 degrees of freedom
Correlation of Parameter Estimates: 推定パラメータの相関値
xmid scal
scal 0.9063
.lin 0.9868 0.9008
## 条件付き線形性を使わない
fm3DNase1 <- nls( density ~ Asym/(1 + exp(( xmid - log(c...
data = DNase1,
start = list( Asym = 3, xmid = 0, scal...
trace = TRUE )
14.32279 : 3 0 1 # 当てはめ途中結果のレポート
0.4542698 : 2.1152456 0.8410193 1.2000640 # 残差平方和...
0.05869602 : 2.446376 1.747516 1.189515
0.005663523 : 2.294087 1.412198 1.020463
0.004791528 : 2.341429 1.479688 1.040758
0.004789569 : 2.345135 1.483047 1.041439
0.004789569 : 2.345179 1.483089 1.041454
summary(fm3DNase1) # 当てはめ結果の要約
Formula: density ~ Asym/(1 + exp((xmid - log(conc))/scal...
Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
Asym 2.34518 0.07815 30.01 2.17e-13 ***
xmid 1.48309 0.08135 18.23 1.22e-10 ***
scal 1.04145 0.03227 32.27 8.51e-14 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1...
Residual standard error: 0.01919 on 13 degrees of freedom
Correlation of Parameter Estimates: 推定パラメータの相関値
Asym xmid
xmid 0.9868
scal 0.9008 0.9063
**例2 重み付き非線形回帰 (非線形モデルを関数で与える) [#...
データ Puromycin は成分 conc (数値), rate (数値), state (...
data(Puromycin)
Treated <- Puromycin[Puromycin$state == "treated", ]
weighted.MM <- function(resp, conc, Vm, K) # 非線形モデ...
{ ## 目的: 白本の p.451 -- RHS にある Michaelis-Menten ...
pred <- (Vm * conc)/(K + conc)
(resp - pred) / sqrt(pred) # sqrt(pred) で割り重みと...
}
# Vm, K がモデルパラメータ、rate および conc はデータフ...
Pur.wt <- nls(~weighted.MM(rate, conc, Vm, K), data = Tr...
start = list(Vm = 200, K = 0.1), trace = T...
112.5978 : 200.0 0.1
17.33824 : 205.67588840 0.04692873
14.6097 : 206.33087396 0.05387279
14.59694 : 206.79883508 0.05457132
14.59690 : 206.83291286 0.05460917
14.59690 : 206.83468191 0.05461109
> summary(Pur.wt)
Formula: 0 ~ weighted.MM(rate, conc, Vm, K) # 関数 weigh...
Parameters:
Estimate Std. Error t value Pr(>|t|)
Vm 2.068e+02 9.225e+00 22.421 7.00e-10 ***
K 5.461e-02 7.979e-03 6.845 4.49e-05 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1...
Residual standard error: 1.208 on 10 degrees of freedom
Correlation of Parameter Estimates:
Vm
K 0.7543
*一次元分布に対する最尤推定。 MASS パッケージ中の fitdist...
-書式: COLOR(magenta){fitdistr(x, densfun, start, ...)}
-引数
-- COLOR(red){x} 数値ベクトル
-- COLOR(red){densfun} 密度関数名を与える文字列か関数(第...
-- COLOR(red){start} 初期値を含む最適化されるべきパラメー...
-- COLOR(red){...} 'densfun' もしくは 'optim' に引き渡さ...
- 詳細
-- 対数尤度の数値微分値を用いた直接的最大化が実行される。
-- 以下の名前指定分布では、もし 'start' が省略されたり、...
-- 再尤推定量が陽に書き下せる時は、それを計算するだけで、...
- 返り値。クラス '"fitdistr"' のオブジェクトで次の二つの...
-- COLOR(red){estimate} 推定パラメータ
-- COLOR(red){sd} 推定標準偏差
**例 [#yfcc8b32]
> set.seed(123) # 乱数種
> x <- rgamma(100, shape = 5, rate = 0.1)
> fitdistr(x, "gamma")
shape rate
6.48592894 0.13649147
(0.89425721) (0.01956555)
> fitdistr(x, dgamma, list(shape = 1, rate = 0.1), lower...
shape rate
2.156761461 0.010000000
(0.277607887) (0.001433778)
> set.seed(123) # 乱数種
> x2 <- rt(250, df = 9) # t 分布に従う乱数...
> fitdistr(x2, "t", df = 9) # 最尤推定による自...
m s # 二つのパラメータ...
-0.01080467 1.04402833
( 0.07222206) ( 0.05433544)
> fitdistr(x2, "t") # 自由度...
m s df # 三つの...
-0.009754149 1.006270714 6.632753559
( 0.071474570) ( 0.077107965) ( 2.716093153)
> mydt <- function(x, m, s, df) dt((x - m)/s, df)/s # 密...
> fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, lower =...
m s # 二つのパラメータの推...
-0.01069635 1.04409435
( 0.07222562) ( 0.05434249)
> set.seed(123) # 乱数種
> x3 <- rweibull(100, shape = 4, scale = 100) # ワイブル...
> fitdistr(x3, "weibull") # ワイブル分布の当てはめ
shape scale # 二つのパラメータの推...
4.0807404 99.9812022
( 0.3127017) ( 2.5816389)
> set.seed(123) # 乱数種
> x4 <- rnegbin(500, mu = 5, theta = 4) # 負の二項分布に...
> fitdistr(x4, "Negative Binomial") # 負の二項分布を当...
size mu # 二つのパラメータ...
4.2178713 4.9439803
(0.5048242) (0.1465538)
Warning messages: # 関連警告
1: NaNs produced in: dgamma(x, shape, scale, log)
2: NaNs produced in: dgamma(x, shape, scale, log)
3: NaNs produced in: dgamma(x, shape, scale, log)
4: bounds can only be used with method L-BFGS-B in:
optim(start, mylogfn, x = x, hessian = TRUE, ...)
5: bounds can only be used with method L-BFGS-B in:
optim(start, myfn, x = x, hessian = TRUE, ...)
6: NaNs produced in: dweibull(x, shape, scale, log)
7: NaNs produced in: dweibull(x, shape, scale, log)
*反復測定分散分析 (Repeated measured ANOVA) [#bda9c928]
**関数 rep.aov [#ee042739]
被験者内要因に関して、Mauchlyの球面性検定、G-G補正、H-F補...
ついでに、Ryan法で多重比較もします。
データ形式はaovに渡すものと同じ。
効果はformulaで、
response~S(subject)+B(a,b)+I(c,d)
とします。
Sのなかには被験者のラベル、Bは被験者間要因、Iは被験者内要...
rep.aov<-function(formula, data=NULL, comment=NULL,...)
{
Terms <- if (missing(data))
terms(formula, c("I", "S", "B"))
else terms(formula, c("I", "S", "B"), data = data)
indI <- attr(Terms, "specials")$I
indB <- attr(Terms, "specials")$B
indS <- attr(Terms, "specials")$S
strI<-if(is.null(indI))NULL
else attr(Terms,"variables")[[1+indI]]
labelsI<-if(is.null(indI))NULL
else sapply(2:length(strI),function(i)deparse(strI[[i]...
strB<-if(is.null(indB))NULL
else attr(Terms,"variables")[[1+indB]]
labelsB<-if(is.null(indB))NULL
else sapply(2:length(strB),function(i)deparse(strB[[i]...
labelsS<-deparse(attr(Terms,"variables")[[1+indS]][[2]])
for(i in rev(labelsI))data<-data[order(data[,i]),]
data<-data[order(data[,labelsS]),]
for(i in rev(labelsB))data<-data[order(data[,i]),]
levelsB<-lapply(labelsB,function(n){data[,n]})
names(levelsB)<-labelsB
levelsI<-lapply(labelsI,function(n){data[,n]})
names(levelsI)<-labelsI
levelsS<-data[,labelsS]
resp<-deparse(formula[[2]])
nc<-prod(sapply(levelsI,function(i){length(levels(i))}))
dmat<-matrix(data[,resp],ncol=nc,byrow=T)
nr<-nrow(dmat)
fm.lm<-NULL;ret.lm<-NULL
if(is.null(indB)){
fm.lm<-dmat~1
ret.lm<-lm(fm.lm)
}else{
lvB<-lapply(labelsB,function(i)levelsB[[i]][1+0:(nr-...
names(lvB)<-labelsB
local({
eval(parse(text=paste(names(lvB),lvB,sep="=")));
fm.lm<<-as.formula(sprintf("dmat ~ %s", paste(labe...
ret.lm<<-lm(fm.lm)
})
}
lvI<-NULL
if(!is.null(indI)){
lvI<-as.data.frame(unique(data[,labelsI]))
colnames(lvI)<-labelsI;rownames(lvI)<-1:dim(lvI)[1]
}
ret.anova<-ret.mauchly<-ret.mcomp<-anvX<-anvM<-list()
terms.anova<-terms(as.formula(sprintf("~ %s",paste(lab...
lv.ta<-attr(terms.anova,"term.labels")
od.ta<-attr(terms.anova,"order")
## execute aov
lv.b<-paste(paste(c(labelsS,labelsB),collapse=":"))
fm.aov<-as.formula(paste(resp,"~",paste(c(labelsI,labe...
"Error(",paste(lv.b,paste(lv.b,lv.t...
collapse=" + "),sep=" + "),")"))
suppressWarnings(ret.aov<-aov(fm.aov,data))
## mauchly's spherical test and G-G & H-F eps, p.val.
for(i in 1:length(lv.ta)){
j<-lv.ta[i]
if(od.ta[i]==1){
anvX[[j]]<-~1
anvM[[j]]<-as.formula(paste("~",j))
}else{
j2<-gsub(":","*",j)
anvX[[j]]<-as.formula(paste("~",j2,"-",j))
anvM[[j]]<-as.formula(paste("~",j2))
}
ret.anova[[j]]<-r<-anova(ret.lm,X=anvX[[j]],M=anvM[[...
if(od.ta[i]==1){
for(k in which(r$Pr<0.05 && r$num>1)){
tmp.aov<-ret.aov[[paste(c(labelsS,labelsB,j),col...
df<-summary(tmp.aov)[[1]]$Df;df<-df[length(df)]
mse<-summary(tmp.aov)[[1]]$"Mean Sq";mse<-mse[le...
dat.tmp<-data[order(data[,j]),]
dat.tmp<-matrix(dat.tmp[,resp],ncol=length(level...
ret.mcomp[[j]]<-pairwise.ryan(dat.tmp,mse,df)
attr(ret.mcomp[[j]],"Effect")<-j
}
}
suppressWarnings( ret.mauchly[[j]]<-mauchly.test(ret...
}
for(i in 1:length(ret.anova)){
hd<-attr(ret.anova[[i]],"heading")
## G-G and H-F epsiron
attr(ret.anova[[i]],"eps")<-c(as.numeric(sub(".*: +"...
## name of within factor
attr(ret.anova[[i]],"within")<-lv.ta[[i]]
}
ret<-list("anova"=ret.anova,
"mauchly"=ret.mauchly,
"mcomp"=ret.mcomp,
"aov"=ret.aov)
attr(ret,"between.factor")<-labelsB
attr(ret,"within.factor")<-labelsI
attr(ret,"subject.factor")<-labelsS
attr(ret,"comment")<-comment
class(ret)<-"rep.anova"
return(ret)
}
**関数 print.rep.anova [#mc2828ed]
rep.anovaの結果をうまいこと出力します。
普通に、
print(rep.aov(...))
とすればいいです。
print.rep.anova<-function(x){
## mauchly's spherical test
ms<-eps<-NULL;for(m in x[[2]])ms<-rbind(ms,c(m$stat,m$...
## construct table
tms1<-tms2<-c();tbl<-NULL
for(r in x[[1]]){
eps<-rbind(eps,attr(r,"eps"))
tms1<-c(tms1,attr(r,"within"))
tr<-r[-dim(r)[1],]
tms2<-c(tms2,sapply(paste(rownames(tr),attr(r,"withi...
tbl<-rbind(tbl,r[-dim(r)[1],-1])
}
mtest<-data.frame(cbind(ms,eps))
names(mtest)<-c("Mauchley's W"," p value"," G-G epsilo...
rownames(mtest)<-tms1
tbl<-data.frame(tbl)
names(tbl)<-c("F stat."," DF^", " DF_", " Pr.(>F)", " ...
rownames(tbl)<-tms2
## rename error term in aov table
tbl.aov<-NULL
tx<-summary(x[[4]])
for(i in 1:length(tx)){
for(j in 1:length(tx[[i]])){
ttx<-tx[[i]][[j]]
if(dim(tx[[i]][[j]])[1]>0){
rownames(ttx)[rownames(ttx)=="Residuals"]<-sub("...
tbl.aov<-rbind(tbl.aov,ttx)
}
}
}
## reformat aov table and add G-G and H-F p.val
gg<-hf<-c();rt<-sapply(strsplit(tms2,":"),function(s)p...
for(i in 1:nrow(tbl.aov)){
rn<-paste(sort(strsplit(sub(" *$","",rownames(tbl.ao...
if(any(rn==rt)){gg<-c(gg,tbl[rn,5]);hf<-c(hf,tbl[rn,...
else{gg<-c(gg,NA);hf<-c(hf,NA)}
}
tbl.aov<-cbind(tbl.aov,"G-G Pr."=gg," H-F Pr."=hf)
## digits...
for(i in 1:4)mtest[,i]<-sprintf("%.4f",mtest[,i])
for(i in c(2,3))tbl.aov[,i]<-ifelse(is.na(tbl.aov[,i])...
for(i in c(4,5,6,7))tbl.aov[,i]<-ifelse(is.na(tbl.aov[...
## output
cat("\n\nRepeated measured ANOVA\n\n")
if(!is.null(cm<-attr(x,"comment")))cat(cm,"\n\n")
cat("Between [",paste(attr(x,"between.factor"),collaps...
cat("Within [",paste(attr(x,"within.factor"),collapse=...
cat("Subject [",attr(x,"subject.factor"),"]\n\n")
cat("Mauchley's test of sphericity\n")
cat("-------------------------------------------------...
print(mtest)
cat("\n")
cat("Analysis of Variance Table\n")
cat("-------------------------------------------------...
print(as.matrix(tbl.aov),na.print="",quote=F,right=T)
cat("\n")
for(r in x[[3]])print(r,comment=paste("Main effect of ...
}
**関数 pairwise.ryan [#a24df7e1]
ライアン法による多重比較
pairwise.ryan<-function(x,MSe=NULL,dfe=NULL,alpha=0.05){
mx<-mean(x)
mc<-apply(x,2,mean)
ns<-dim(x)[1]
nc<-dim(x)[2]
names<-colnames(x)
if(is.null(names))names<-1:nc
if(is.null(dfe))dfe<-length(x)-ns-nc+1
if(is.null(MSe))MSe<-(sum((x-mx)^2)-sum(ns*(mc-mx)^2)-...
ret<-list()
ret$method<-"Ryan's method"
ret$data.name<-deparse(substitute(x))
ts<-matrix(nrow=nc-1,ncol=nc-1,dimnames=list(names[2:n...
ps<-matrix(nrow=nc-1,ncol=nc-1,dimnames=list(names[2:n...
aps<-matrix(nrow=nc-1,ncol=nc-1,dimnames=list(names[2:...
table<-NULL
for(i in 1:(nc-1)){
for(j in (i+1):nc){
d<-abs(mc[i]-mc[j])
t0<-d/(sqrt(MSe)*sqrt(2/ns))
r<-length(which(min(mc[i],mc[j])<=mc&mc<=max(mc[i]...
ap<-2*alpha/(nc*(r-1))
p0<-2*(1-pt(t0,dfe))
ps[i,j-1]<-p0
aps[i,j-1]<-ap
table<-rbind(table,data.frame(paste(names[i],"-",n...
}
}
table<-data.frame(table)
names(table)<-c("pair","r","nominal level","t","p","si...
rownames(table)<-1:dim(table)[1]
ret$p.value<-ps
ret$stat<-ts
ret$adjusted.p.value<-aps
ret$table<-table
ret$means<-mc
attr(ret,"class")<-"pairwise.ryan"
return(ret)
}
**関数 print.pairwise.ryan [#p312c125]
pairwise.ryanの結果をうまく出力します。
print.pairwise.ryan<-function(x,comment=NULL){
cat("\n\tPairwise comparisons using", x$method, "\n\n")
if(!is.null(comment))cat(comment)
cat("data: ", x$data.name, "[",paste(x$means,","),"]\n...
tbl<-x$table
tbl[,3]<-sprintf("%.4f",tbl[,3])
tbl[,4]<-sprintf("%.4f",tbl[,4])
tbl[,5]<-sprintf("%.4f",tbl[,5])
print(tbl)
}
**使用例 [#v3043d61]
言語Rによる分散分析( http://mat.isc.chubu.ac.jp/R/tech.ht...
Anova4 on the web ( http://www.hju.ac.jp/~kiriki/anova4/ ...
但し、G-G補正、H-F補正済みp値については未確認なので、どな...
***one-way within subject [#r05ec03c]
> if(1){
+ RB <- data.frame(
+ a = factor(c(rep(1,8), rep(2,8), re...
+ s = factor(rep(1:8, 4)),
+ result = c(9,7,8,8,12,11,8,13, 6,5,...
+ 10,13,8,13,12,14,14,16, 9,11,13,1...
+ (rep.aov(result ~ I(a) + S(s), RB))
+ }
Repeated measured ANOVA
Between [ ] Within [ a ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
a 0.9751 0.9996 0.9832 1.7774
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. H-F...
Error(s) 7 77.00 11.00 ...
a 3 217.50 72.50 21.4437 0.0000 0.0000 0....
Error(s:a) 21 71.00 3.38 ...
Pairwise comparisons using Ryan's method
Main effect of : a
data: dat.tmp [ 9.5 , 6.5 , 12.5 , 13 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0250 3.2631 0.0037 s.
2 1 - 3 2 0.0250 3.2631 0.0037 s.
3 1 - 4 3 0.0125 3.8070 0.0010 s.
4 2 - 3 3 0.0125 6.5262 0.0000 s.
5 2 - 4 4 0.0083 7.0701 0.0000 s.
6 3 - 4 2 0.0250 0.5439 0.5923 n.s.
***two-way within subject [#race7b7b]
> if(1){
+ RBFpq <- data.frame(
+ a = factor(c(rep(1, 20), rep(2, ...
+ b = factor(rep(c(rep(1, 5), rep(...
+ s = factor(rep(1:5, 8)),
+ result = c(3,3,1,3,5, 4,3,4,5,7,...
+ 3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,...
+ (rep.aov(result ~ I(a,b) + S(s), RBFpq))
+ }
Repeated measured ANOVA
Between [ ] Within [ a , b ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
a 1.0000 NaN 1.0000 1.0000
b 0.6817 0.9618 0.7934 2.0373
a:b 0.2482 0.6036 0.6692 1.3449
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. H...
Error(s) 4 38.15 9.54 ...
a 1 16.90 16.90 8.0958 0.0466 0.0466 ...
Error(s:a) 4 8.35 2.09 ...
b 3 19.70 6.57 6.0383 0.0095 0.0173 ...
Error(s:b) 12 13.05 1.09 ...
a:b 3 30.50 10.17 7.0725 0.0054 0.0169 ...
Error(s:a:b) 12 17.25 1.44 ...
Pairwise comparisons using Ryan's method
Main effect of : b
data: dat.tmp [ 3.5 , 4.4 , 4.9 , 5.4 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0250 1.9298 0.0776 n.s.
2 1 - 3 3 0.0125 3.0019 0.0110 s.
3 1 - 4 4 0.0083 4.0740 0.0015 s.
4 2 - 3 2 0.0250 1.0721 0.3048 n.s.
5 2 - 4 3 0.0125 2.1442 0.0532 n.s.
6 3 - 4 2 0.0250 1.0721 0.3048 n.s.
***two-way mixed: between (1) x within (1) [#x5bb0565]
> if(1){
+ SPFp.q <- data.frame(
+ a = factor(c(rep(1,20), rep(2,2...
+ b = factor(rep(c(rep(1,5), rep(...
+ s = factor(c(rep(1:5, 4), rep(6...
+ result = c(3,3,1,3,5, 4,3,4,5,7...
+ 3,5,2,4,6, 2,6,3,6,4, 3,2,3,6...
+ rep.aov(result ~ I(b) + B(a) + S(s), SPFp.q)
+ }
Repeated measured ANOVA
Between [ a ] Within [ b ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
b 0.8428 0.9502 0.8896 1.3767
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. H...
a 1 16.90 16.90 2.9075 0.1266 ...
Error(s:a) 8 46.50 5.81 ...
b 3 19.70 6.57 5.2013 0.0066 0.0092 ...
b:a 3 30.50 10.17 8.0528 0.0007 0.0012 ...
Error(s:a:b) 24 30.30 1.26 ...
Pairwise comparisons using Ryan's method
Main effect of : b
data: dat.tmp [ 3.5 , 4.4 , 4.9 , 5.4 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0250 1.7911 0.0859 n.s.
2 1 - 3 3 0.0125 2.7861 0.0103 s.
3 1 - 4 4 0.0083 3.7811 0.0009 s.
4 2 - 3 2 0.0250 0.9950 0.3296 n.s.
5 2 - 4 3 0.0125 1.9901 0.0581 n.s.
6 3 - 4 2 0.0250 0.9950 0.3296 n.s.
***three-way mixed : between (2) x within (1) [#ic872892]
> if(1){
+ SPFpq.r <- data.frame(
+ a = factor(c(rep(1, 24), rep(2...
+ b = factor(rep(c(rep(1, 12), r...
+ c = factor(rep(c(rep(1, 4), re...
+ s = factor(c(rep(1:4, 3), rep(...
+ result = c(2,6,5,7, 5,7,9,9, ...
+ 6,6,8,10, 3,6,5,8, 6,7,5,6,
+ 1,2,5,2, 2,1,4,5, 1,5,3,5,
+ 5,3,4,6, 5,6,6,7, 5,5,9,7))
+ (rep.aov(result ~ S(s)+I(c)+B(a,b), SPFpq.r))
+ }
Repeated measured ANOVA
Between [ a , b ] Within [ c ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
c 0.8178 0.3309 0.8459 0.9698
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. ...
a 1 96.33 96.33 15.6216 0.0019 ...
b 1 3.00 3.00 0.4865 0.4988 ...
a:b 1 56.33 56.33 9.1351 0.0106 ...
Error(s:a:b) 12 74.00 6.17 ...
c 2 33.50 16.75 9.5714 0.0009 0.0018 ...
c:a 2 6.17 3.08 1.7619 0.1932 0.1995 ...
c:b 2 24.50 12.25 7.0000 0.0040 0.0067 ...
c:a:b 2 41.17 20.58 11.7619 0.0003 0.0007 ...
Error(s:a:b:c) 24 42.00 1.75 ...
Pairwise comparisons using Ryan's method
Main effect of : c
data: dat.tmp [ 4.875 , 5.5 , 6.875 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0333 1.3363 0.1940 n.s.
2 1 - 3 3 0.0167 4.2762 0.0003 s.
3 2 - 3 2 0.0333 2.9399 0.0072 s.
***three-way mixed : between (1) x within (2) [#nf2f0213]
> if(1){
+ SPFp.qr <- data.frame(
+ a = factor(c(rep(1, 24), rep(2...
+ b = factor(rep(c(rep(1, 12), r...
+ c = factor(rep(c(rep(1, 4), re...
+ s = factor(c(rep(1:4, 6), rep(...
+ result = c(2,6,5,7, 5,7,9,9, ...
+ 6,6,8,10, 3,6,5,8, 6,7,5,6,
+ 1,2,5,2, 2,1,4,5, 1,5,3,5,
+ 5,3,4,6, 5,6,6,7, 5,5,9,7))
+ (r<-rep.aov(result ~ S(s)+I(b,c)+B(a), SPFp.qr))
+ }
Repeated measured ANOVA
Between [ a ] Within [ b , c ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
b 1.0000 NaN 1.0000 1.0000
c 0.6719 0.3700 0.7529 0.9503
b:c 0.8354 0.6378 0.8586 1.1699
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. ...
a 1 96.33 96.33 8.7576 0.0253 ...
Error(s:a) 6 66.00 11.00 ...
b 1 3.00 3.00 2.2500 0.1843 0.1843 ...
b:a 1 56.33 56.33 42.2500 0.0006 0.0006 ...
Error(s:a:b) 6 8.00 1.33 ...
c 2 33.50 16.75 25.1250 0.0001 0.0003 ...
c:a 2 6.17 3.08 4.6250 0.0324 0.0489 ...
Error(s:a:c) 12 8.00 0.67 ...
b:c 2 24.50 12.25 4.3235 0.0385 0.0477 ...
b:c:a 2 41.17 20.58 7.2647 0.0086 0.0128 ...
Error(s:a:b:c) 12 34.00 2.83 ...
Pairwise comparisons using Ryan's method
Main effect of : c
data: dat.tmp [ 4.875 , 5.5 , 6.875 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0333 2.1651 0.0512 n.s.
2 1 - 3 3 0.0167 6.9282 0.0000 s.
3 2 - 3 2 0.0333 4.7631 0.0005 s.
***three-way within subject [#ccafdb19]
> if(1){
+ RBFpqr <- data.frame(
+ a = factor(c(rep(1, 24), rep(2,...
+ b = factor(rep(c(rep(1, 12), re...
+ c = factor(rep(c(rep(1, 4), rep...
+ s = factor(rep(1:4, 12)),
+ result = c(2,6,5,7, 5,7,9,9, 9...
+ 6,6,8,10, 3,6,5,8, 6,7,5,6,
+ 1,2,5,2, 2,1,4,5, 1,5,3,5,
+ 5,3,4,6, 5,6,6,7, 5,5,9,7))
+ (rep.aov(result ~ S(s)+I(a,b,c), RBFpqr))
+ }
Repeated measured ANOVA
Between [ ] Within [ a , b , c ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
a 1.0000 NaN 1.0000 1.0000
b 1.0000 NaN 1.0000 1.0000
c 0.8700 0.8700 0.8850 2.0650
a:b 1.0000 NaN 1.0000 1.0000
a:c 0.5422 0.5422 0.6860 1.0710
b:c 0.9838 0.9838 0.9840 2.8450
a:b:c 0.4831 0.4831 0.6593 0.9736
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr....
Error(s) 3 60.33 20.11 ...
a 1 96.33 96.33 51.0000 0.0057 0.0057...
Error(s:a) 3 5.67 1.89 ...
b 1 3.00 3.00 1.4211 0.3189 0.3189...
Error(s:b) 3 6.33 2.11 ...
c 2 33.50 16.75 60.3000 0.0001 0.0002...
Error(s:c) 6 1.67 0.28 ...
a:b 1 56.33 56.33 101.4000 0.0021 0.0021...
Error(s:a:b) 3 1.67 0.56 ...
a:c 2 6.17 3.08 2.9211 0.1301 0.1625...
Error(s:a:c) 6 6.33 1.06 ...
b:c 2 24.50 12.25 5.3780 0.0459 0.0471...
Error(s:b:c) 6 13.67 2.28 ...
a:b:c 2 41.17 20.58 6.0738 0.0361 0.0670...
Error(s:a:b:c) 6 20.33 3.39 ...
Pairwise comparisons using Ryan's method
Main effect of : c
data: dat.tmp [ 4.875 , 5.5 , 6.875 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0333 3.3541 0.0153 s.
2 1 - 3 3 0.0167 10.7331 0.0000 s.
3 2 - 3 2 0.0333 7.3790 0.0003 s.
**TODO [#s6698f53]
- ライアン法の多重比較で混合要因ANOVAのときに被験者数が異...
- 交互作用が有意のときに下位検定を行いたいんですが、まだ...
- コードが汚いです。
終了行:
COLOR(red){SIZE(25){Rの統計解析関数Tips}}
R に関する Tips で欠けているもの、そしておそらく一番困難(...
#contents
~
*lm 関数を用いた線形モデルの当てはめ結果の要約 (r-help の...
> df1 <- data.frame(x=1:6, y=rep(1:3, 2))
> fit <- lm(y~x, df1)
> Sum <- summary(fit)
# 当てはめオブジェクト fit の要約(普通ユーザーはこの情報...
> Sum
Call:
lm(formula = y ~ x, data = df1)
Residuals:
1 2 3 4 5 6
-0.4286 0.3429 1.1143 -1.1143 -0.3429 0.4286
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2000 0.8177 1.468 0.216
x 0.2286 0.2100 1.089 0.338
Residual standard error: 0.8783 on 4 degrees of freedom
Multiple R-Squared: 0.2286, Adjusted R-squared: 0.03...
F-statistic: 1.185 on 1 and 4 DF, p-value: 0.3375
# 当てはめオブジェクトの要約が持つ情報(リスト形式になっ...
> attributes(Sum)
$names
[1] "call" "terms" "residuals" "coe...
[5] "sigma" "df" "r.squared" "adj...
[9] "fstatistic" "cov.unscaled"
# 回帰係数と関連情報を表形式で得る(summary(fit) が表示す...
> coefficients(Sum) # coefficients 関数は線形回帰当ては...
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2000000 0.8176622 1.467599 0.2161194
x 0.2285714 0.2099563 1.088662 0.3375019
> Sum$coefficients # これでも良い
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2000000 0.8176622 1.467599 0.2161194
x 0.2285714 0.2099563 1.088662 0.3375019
参考までに述べると lm 関数が返す当てはめオブジェクトは実...
> str(fit) # str 関数で当てはめオブジェクトの内容を見る
List of 12
$ coefficients : Named num [1:2] 1.200 0.229
..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
$ residuals : Named num [1:6] -0.429 0.343 1.114 -...
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ effects : Named num [1:6] -4.89898 0.95618 1.2...
..- attr(*, "names")= chr [1:6] "(Intercept)" "x" "" "...
$ rank : int 2
$ fitted.values: Named num [1:6] 1.43 1.66 1.89 2.11 2....
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ assign : int [1:2] 0 1
$ qr :List of 5
..$ qr : num [1:6, 1:2] -2.449 0.408 0.408 0.408 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:6] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:2] "(Intercept)" "x"
.. ..- attr(*, "assign")= int [1:2] 0 1
..$ qraux: num [1:2] 1.41 1.19
..$ pivot: int [1:2] 1 2
..$ tol : num 1e-07
..$ rank : int 2
..- attr(*, "class")= chr "qr"
$ df.residual : int 4
$ xlevels : list()
$ call : language lm(formula = y ~ x, data = df1)
$ terms :Classes 'terms', 'formula' length 3 y ~ x
.. ..- attr(*, "variables")= language list(y, x)
.. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. ..$ : chr "x"
.. ..- attr(*, "term.labels")= chr "x"
.. ..- attr(*, "order")= int 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=length 33 <environment>
.. ..- attr(*, "predvars")= language list(y, x)
$ model :`data.frame': 6 obs. of 2 variables:
..$ y: int [1:6] 1 2 3 1 2 3
..$ x: int [1:6] 1 2 3 4 5 6
..- attr(*, "terms")=Classes 'terms', 'formula' length...
.. .. ..- attr(*, "variables")= language list(y, x)
.. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. .. ..$ : chr "x"
.. .. ..- attr(*, "term.labels")= chr "x"
.. .. ..- attr(*, "order")= int 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=length 33 <environme...
.. .. ..- attr(*, "predvars")= language list(y, x)
- attr(*, "class")= chr "lm"
> str(Sum) # Sum の持つ内部情報の詳細を見る
List of 11
$ call : language lm(formula = y ~ x, data = df1)
$ terms :Classes 'terms', 'formula' length 3 y ~ x
.. ..- attr(*, "variables")= language list(y, x)
.. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. ..$ : chr "x"
.. ..- attr(*, "term.labels")= chr "x"
.. ..- attr(*, "order")= int 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=length 34 <environment>
.. ..- attr(*, "predvars")= language list(y, x)
$ residuals : Named num [1:6] -0.429 0.343 1.114 -...
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ coefficients : num [1:2, 1:4] 1.200 0.229 0.818 0.210...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "(Intercept)" "x"
.. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "...
$ aliased : Named logi [1:2] FALSE FALSE
..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
$ sigma : num 0.878
$ df : int [1:3] 2 4 2
$ r.squared : num 0.229
$ adj.r.squared: num 0.0357
$ fstatistic : Named num [1:3] 1.19 1.00 4.00
..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
$ cov.unscaled : num [1:2, 1:2] 0.8667 -0.2000 -0.2000...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "(Intercept)" "x"
.. ..$ : chr [1:2] "(Intercept)" "x"
- attr(*, "class")= chr "summary.lm"
# 当てはめオブジェクト fit に基づく二次的解析の例(分散分...
> aovfit <- aov(fit)
> aovfit
Call:
aov(formula = fit)
Terms:
x Residuals
Sum of Squares 0.9142857 3.0857143
Deg. of Freedom 1 4
Residual standard error: 0.87831
Estimated effects may be unbalanced
> summary(aovfit)
Df Sum Sq Mean Sq F value Pr(>F)
x 1 0.91429 0.91429 1.1852 0.3375
Residuals 4 3.08571 0.77143
> str(aovfit)
List of 12
$ coefficients : Named num [1:2] 1.200 0.229
..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
$ residuals : Named num [1:6] -0.429 0.343 1.114 -...
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ effects : Named num [1:6] -4.89898 0.95618 1.2...
..- attr(*, "names")= chr [1:6] "(Intercept)" "x" "" ""...
$ rank : int 2
$ fitted.values: Named num [1:6] 1.43 1.66 1.89 2.11 2....
..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
$ assign : int [1:2] 0 1
$ qr :List of 5
..$ qr : num [1:6, 1:2] -2.449 0.408 0.408 0.408 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:6] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:2] "(Intercept)" "x"
.. ..- attr(*, "assign")= int [1:2] 0 1
..$ qraux: num [1:2] 1.41 1.19
..$ pivot: int [1:2] 1 2
..$ tol : num 1e-07
..$ rank : int 2
..- attr(*, "class")= chr "qr"
$ df.residual : int 4
$ xlevels : list()
$ call : language aov(formula = fit)
$ terms :Classes 'terms', 'formula' length 3 y ~ x
.. ..- attr(*, "variables")= language list(y, x)
.. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. ..$ : chr "x"
.. ..- attr(*, "term.labels")= chr "x"
.. ..- attr(*, "order")= int 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=length 35 <environment>
.. ..- attr(*, "predvars")= language list(y, x)
$ model :`data.frame': 6 obs. of 2 variables:
..$ y: int [1:6] 1 2 3 1 2 3
..$ x: int [1:6] 1 2 3 4 5 6
..- attr(*, "terms")=Classes 'terms', 'formula' length...
.. .. ..- attr(*, "variables")= language list(y, x)
.. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:2] "y" "x"
.. .. .. .. ..$ : chr "x"
.. .. ..- attr(*, "term.labels")= chr "x"
.. .. ..- attr(*, "order")= int 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=length 35 <environme...
.. .. ..- attr(*, "predvars")= language list(y, x)
- attr(*, "class")= chr [1:2] "aov" "lm"
* 線形モデル当てはめ関数におけるモデル公式 (2004.01.05) [...
lm (線形回帰)、aov (分散分析)、glm (一般化線形回帰) 等の...
|y ~ ...| 左辺の目的変数を、右辺の説明変数の線形式で線形...
|+ A| 説明変数 A を加える|
|- A |説明変数 A を除く|
|+ 0 または - 1 | 切片(定数)項の除外|
|A:B| 二つの説明変数 A, B の交互作用項 (A と B の積からな...
|A * B | A + B + A:Bと同じ |
|(...)*(...) | 右括弧内の変数と左括弧内の変数の掛け合わせ...
| |ただし A^2, A^2*B 等はそれぞれ A, A*B と解釈|
|(A+B)^2| (A+B)*(A+B) と同じ。A + B + A:B の意味|
| A %in% B | 今一意味がわからない?|
より一般に A:B:C や (A+B+C)^3 といった高次交互作用項も同...
y を目的変数、A, B, C 等を説明変数・因子とすると
| y ~ A| y を A と切片(定数)項で単回帰|
| y ~ A - 1 | 切片項のない、原点を通る単回帰|
| y ~ A + 0 | 同上|
| 0 + y ~ A | 同上|
| y ~ A + B| y を A, B と切片(定数)項で重回帰|
| y ~ A + B - 1| y を A, B で重回帰。切片(定数)項は含まな...
| y ~ A + B + A:B |A,B とそれらの二次交互作用項で重回帰,...
| y ~ (A+B)^2 |同上|
| y ~ (A+B)^2 - 1|A,B とそれらの二次交互作用項で重回帰。...
| y ~ (A+B+C)^2 |A,B,C とそれらのすべての二次交互作用項...
| y ~ A + B + C + A:B + A:C + B:C |同上|
| y ~ (A+B+C)^2 - B:C | y ~ A + B + C + A:B + A:C の意味|
| y ~ (A+B+C)^3 |A,B,C とそれらのすべての二次及び三次交...
| y ~ A + A %in% B | y ~ A + A:B の意味|
モデル式には説明変数の関数を使うことが出来る。
| y ~ A + log(B)| y を A と B の対数値で説明|
| log(y) ~ exp(A) + log(B)| log(y) を exp(A) と log(B) で...
| y ~ A/B| y を A と B の比で説明|
モデル式で特別の意味を持つオペレータ(例えば +,-,*,^)を、...
|y ~ A + I(B+C) | y を A と B+C (B と C の和) で説明|
|y ~ A + I(B*C) | y を A と B*C (B と C の積) で説明、A +...
|y ~ A + I(B^2) | y を A と B^2 (B の自乗値) で説明|
*非線形最小自乗法 nls() 関数 [#iad5cef4]
非線形モデルの非線形最小自乗推定専用の関数で、独自の最適...
- COLOR(magenta){用法}
nls(formula, data = parent.frame(), start, control = nl...
algorithm = "default", trace = FALSE, subset,
weights, na.action)
- COLOR(magenta){引数}
--COLOR(red){formula} 変量とパラメータを含む非線形モデル...
--COLOR(red){data} オプションのデータフレームで、その中で...
--COLOR(red){start} 初期推定値の名前付きリスト、もしくは...
--COLOR(red){control} 制御仕様を決めるオプションのリスト...
--COLOR(red){algorithm} 使用されるアルゴリズムを指定する...
--COLOR(red){trace} 論理値で、繰り返しの過程の記録を表示...
--COLOR(red){subset} オプションのベクトルで、当てはめの過...
--COLOR(red){weights} オプションの数値ベクトルで、(固定さ...
--COLOR(red){na.action} データが NA を含む時に、どうする...
-COLOR(magenta){詳細}
-- nls を人工的な "残差がゼロ" のデータに使ってはならない...
-- nls オブジェクトは典型的には当てはめモデルオブジェクト...
-COLOR(magenta){返り値} 次の成分からなるリスト
--COLOR(red){m} モデルを含む nlsModel オブジェクト
--COLOR(red){data} データ引数として nls に引き渡された公...
**例1 [#yca8db2a]
データフレーム DNase は成分 Run (順序つき因子, 二因子)、c...
data( DNase )
DNase1 <- DNase[ DNase$Run == 1, ] # データの一部を取り...
## 自動スタートモデルを使う。SSlogis はロジスティックモ...
fm1DNase1 <- nls( density ~ SSlogis( log(conc), Asym, xm...
summary( fm1DNase1 ) # 当てはめ結果の要約
Formula: density ~ SSlogis(log(conc), Asym, xmid, scal) ...
Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
Asym 2.34518 0.07815 30.01 2.17e-13 ***
xmid 1.48309 0.08135 18.23 1.22e-10 ***
scal 1.04146 0.03227 32.27 8.51e-14 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1...
Residual standard error: 0.01919 on 13 degrees of freedom
Correlation of Parameter Estimates: 推定パラメータの相関値
Asym xmid
xmid 0.9868
scal 0.9008 0.9063
## 条件付き線形性を使う
fm2DNase1 <- nls( density ~ 1/(1 + exp(( xmid - log(conc...
data = DNase1,
start = list( xmid = 0, scal = 1 ),
alg = "plinear", trace = TRUE )
0.7139315 : 0.000000 1.000000 1.453853 # 当てはめ途中結...
0.1445295 : 1.640243 1.390186 2.461754 # 残差平方和と推...
0.008302151 : 1.620899 1.054228 2.478388
0.004794192 : 1.485226 1.043709 2.347334
0.004789569 : 1.483130 1.041468 2.345218
0.004789569 : 1.483090 1.041455 2.345180
summary( fm2DNase1 ) # 当てはめ結果の要約
Formula: density ~ 1/(1 + exp((xmid - log(conc))/scal)) ...
Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
xmid 1.48309 0.08135 18.23 1.22e-10 ***
scal 1.04145 0.03227 32.27 8.51e-14 ***
.lin 2.34518 0.07815 30.01 2.17e-13 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1...
Residual standard error: 0.01919 on 13 degrees of freedom
Correlation of Parameter Estimates: 推定パラメータの相関値
xmid scal
scal 0.9063
.lin 0.9868 0.9008
## 条件付き線形性を使わない
fm3DNase1 <- nls( density ~ Asym/(1 + exp(( xmid - log(c...
data = DNase1,
start = list( Asym = 3, xmid = 0, scal...
trace = TRUE )
14.32279 : 3 0 1 # 当てはめ途中結果のレポート
0.4542698 : 2.1152456 0.8410193 1.2000640 # 残差平方和...
0.05869602 : 2.446376 1.747516 1.189515
0.005663523 : 2.294087 1.412198 1.020463
0.004791528 : 2.341429 1.479688 1.040758
0.004789569 : 2.345135 1.483047 1.041439
0.004789569 : 2.345179 1.483089 1.041454
summary(fm3DNase1) # 当てはめ結果の要約
Formula: density ~ Asym/(1 + exp((xmid - log(conc))/scal...
Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
Asym 2.34518 0.07815 30.01 2.17e-13 ***
xmid 1.48309 0.08135 18.23 1.22e-10 ***
scal 1.04145 0.03227 32.27 8.51e-14 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1...
Residual standard error: 0.01919 on 13 degrees of freedom
Correlation of Parameter Estimates: 推定パラメータの相関値
Asym xmid
xmid 0.9868
scal 0.9008 0.9063
**例2 重み付き非線形回帰 (非線形モデルを関数で与える) [#...
データ Puromycin は成分 conc (数値), rate (数値), state (...
data(Puromycin)
Treated <- Puromycin[Puromycin$state == "treated", ]
weighted.MM <- function(resp, conc, Vm, K) # 非線形モデ...
{ ## 目的: 白本の p.451 -- RHS にある Michaelis-Menten ...
pred <- (Vm * conc)/(K + conc)
(resp - pred) / sqrt(pred) # sqrt(pred) で割り重みと...
}
# Vm, K がモデルパラメータ、rate および conc はデータフ...
Pur.wt <- nls(~weighted.MM(rate, conc, Vm, K), data = Tr...
start = list(Vm = 200, K = 0.1), trace = T...
112.5978 : 200.0 0.1
17.33824 : 205.67588840 0.04692873
14.6097 : 206.33087396 0.05387279
14.59694 : 206.79883508 0.05457132
14.59690 : 206.83291286 0.05460917
14.59690 : 206.83468191 0.05461109
> summary(Pur.wt)
Formula: 0 ~ weighted.MM(rate, conc, Vm, K) # 関数 weigh...
Parameters:
Estimate Std. Error t value Pr(>|t|)
Vm 2.068e+02 9.225e+00 22.421 7.00e-10 ***
K 5.461e-02 7.979e-03 6.845 4.49e-05 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1...
Residual standard error: 1.208 on 10 degrees of freedom
Correlation of Parameter Estimates:
Vm
K 0.7543
*一次元分布に対する最尤推定。 MASS パッケージ中の fitdist...
-書式: COLOR(magenta){fitdistr(x, densfun, start, ...)}
-引数
-- COLOR(red){x} 数値ベクトル
-- COLOR(red){densfun} 密度関数名を与える文字列か関数(第...
-- COLOR(red){start} 初期値を含む最適化されるべきパラメー...
-- COLOR(red){...} 'densfun' もしくは 'optim' に引き渡さ...
- 詳細
-- 対数尤度の数値微分値を用いた直接的最大化が実行される。
-- 以下の名前指定分布では、もし 'start' が省略されたり、...
-- 再尤推定量が陽に書き下せる時は、それを計算するだけで、...
- 返り値。クラス '"fitdistr"' のオブジェクトで次の二つの...
-- COLOR(red){estimate} 推定パラメータ
-- COLOR(red){sd} 推定標準偏差
**例 [#yfcc8b32]
> set.seed(123) # 乱数種
> x <- rgamma(100, shape = 5, rate = 0.1)
> fitdistr(x, "gamma")
shape rate
6.48592894 0.13649147
(0.89425721) (0.01956555)
> fitdistr(x, dgamma, list(shape = 1, rate = 0.1), lower...
shape rate
2.156761461 0.010000000
(0.277607887) (0.001433778)
> set.seed(123) # 乱数種
> x2 <- rt(250, df = 9) # t 分布に従う乱数...
> fitdistr(x2, "t", df = 9) # 最尤推定による自...
m s # 二つのパラメータ...
-0.01080467 1.04402833
( 0.07222206) ( 0.05433544)
> fitdistr(x2, "t") # 自由度...
m s df # 三つの...
-0.009754149 1.006270714 6.632753559
( 0.071474570) ( 0.077107965) ( 2.716093153)
> mydt <- function(x, m, s, df) dt((x - m)/s, df)/s # 密...
> fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, lower =...
m s # 二つのパラメータの推...
-0.01069635 1.04409435
( 0.07222562) ( 0.05434249)
> set.seed(123) # 乱数種
> x3 <- rweibull(100, shape = 4, scale = 100) # ワイブル...
> fitdistr(x3, "weibull") # ワイブル分布の当てはめ
shape scale # 二つのパラメータの推...
4.0807404 99.9812022
( 0.3127017) ( 2.5816389)
> set.seed(123) # 乱数種
> x4 <- rnegbin(500, mu = 5, theta = 4) # 負の二項分布に...
> fitdistr(x4, "Negative Binomial") # 負の二項分布を当...
size mu # 二つのパラメータ...
4.2178713 4.9439803
(0.5048242) (0.1465538)
Warning messages: # 関連警告
1: NaNs produced in: dgamma(x, shape, scale, log)
2: NaNs produced in: dgamma(x, shape, scale, log)
3: NaNs produced in: dgamma(x, shape, scale, log)
4: bounds can only be used with method L-BFGS-B in:
optim(start, mylogfn, x = x, hessian = TRUE, ...)
5: bounds can only be used with method L-BFGS-B in:
optim(start, myfn, x = x, hessian = TRUE, ...)
6: NaNs produced in: dweibull(x, shape, scale, log)
7: NaNs produced in: dweibull(x, shape, scale, log)
*反復測定分散分析 (Repeated measured ANOVA) [#bda9c928]
**関数 rep.aov [#ee042739]
被験者内要因に関して、Mauchlyの球面性検定、G-G補正、H-F補...
ついでに、Ryan法で多重比較もします。
データ形式はaovに渡すものと同じ。
効果はformulaで、
response~S(subject)+B(a,b)+I(c,d)
とします。
Sのなかには被験者のラベル、Bは被験者間要因、Iは被験者内要...
rep.aov<-function(formula, data=NULL, comment=NULL,...)
{
Terms <- if (missing(data))
terms(formula, c("I", "S", "B"))
else terms(formula, c("I", "S", "B"), data = data)
indI <- attr(Terms, "specials")$I
indB <- attr(Terms, "specials")$B
indS <- attr(Terms, "specials")$S
strI<-if(is.null(indI))NULL
else attr(Terms,"variables")[[1+indI]]
labelsI<-if(is.null(indI))NULL
else sapply(2:length(strI),function(i)deparse(strI[[i]...
strB<-if(is.null(indB))NULL
else attr(Terms,"variables")[[1+indB]]
labelsB<-if(is.null(indB))NULL
else sapply(2:length(strB),function(i)deparse(strB[[i]...
labelsS<-deparse(attr(Terms,"variables")[[1+indS]][[2]])
for(i in rev(labelsI))data<-data[order(data[,i]),]
data<-data[order(data[,labelsS]),]
for(i in rev(labelsB))data<-data[order(data[,i]),]
levelsB<-lapply(labelsB,function(n){data[,n]})
names(levelsB)<-labelsB
levelsI<-lapply(labelsI,function(n){data[,n]})
names(levelsI)<-labelsI
levelsS<-data[,labelsS]
resp<-deparse(formula[[2]])
nc<-prod(sapply(levelsI,function(i){length(levels(i))}))
dmat<-matrix(data[,resp],ncol=nc,byrow=T)
nr<-nrow(dmat)
fm.lm<-NULL;ret.lm<-NULL
if(is.null(indB)){
fm.lm<-dmat~1
ret.lm<-lm(fm.lm)
}else{
lvB<-lapply(labelsB,function(i)levelsB[[i]][1+0:(nr-...
names(lvB)<-labelsB
local({
eval(parse(text=paste(names(lvB),lvB,sep="=")));
fm.lm<<-as.formula(sprintf("dmat ~ %s", paste(labe...
ret.lm<<-lm(fm.lm)
})
}
lvI<-NULL
if(!is.null(indI)){
lvI<-as.data.frame(unique(data[,labelsI]))
colnames(lvI)<-labelsI;rownames(lvI)<-1:dim(lvI)[1]
}
ret.anova<-ret.mauchly<-ret.mcomp<-anvX<-anvM<-list()
terms.anova<-terms(as.formula(sprintf("~ %s",paste(lab...
lv.ta<-attr(terms.anova,"term.labels")
od.ta<-attr(terms.anova,"order")
## execute aov
lv.b<-paste(paste(c(labelsS,labelsB),collapse=":"))
fm.aov<-as.formula(paste(resp,"~",paste(c(labelsI,labe...
"Error(",paste(lv.b,paste(lv.b,lv.t...
collapse=" + "),sep=" + "),")"))
suppressWarnings(ret.aov<-aov(fm.aov,data))
## mauchly's spherical test and G-G & H-F eps, p.val.
for(i in 1:length(lv.ta)){
j<-lv.ta[i]
if(od.ta[i]==1){
anvX[[j]]<-~1
anvM[[j]]<-as.formula(paste("~",j))
}else{
j2<-gsub(":","*",j)
anvX[[j]]<-as.formula(paste("~",j2,"-",j))
anvM[[j]]<-as.formula(paste("~",j2))
}
ret.anova[[j]]<-r<-anova(ret.lm,X=anvX[[j]],M=anvM[[...
if(od.ta[i]==1){
for(k in which(r$Pr<0.05 && r$num>1)){
tmp.aov<-ret.aov[[paste(c(labelsS,labelsB,j),col...
df<-summary(tmp.aov)[[1]]$Df;df<-df[length(df)]
mse<-summary(tmp.aov)[[1]]$"Mean Sq";mse<-mse[le...
dat.tmp<-data[order(data[,j]),]
dat.tmp<-matrix(dat.tmp[,resp],ncol=length(level...
ret.mcomp[[j]]<-pairwise.ryan(dat.tmp,mse,df)
attr(ret.mcomp[[j]],"Effect")<-j
}
}
suppressWarnings( ret.mauchly[[j]]<-mauchly.test(ret...
}
for(i in 1:length(ret.anova)){
hd<-attr(ret.anova[[i]],"heading")
## G-G and H-F epsiron
attr(ret.anova[[i]],"eps")<-c(as.numeric(sub(".*: +"...
## name of within factor
attr(ret.anova[[i]],"within")<-lv.ta[[i]]
}
ret<-list("anova"=ret.anova,
"mauchly"=ret.mauchly,
"mcomp"=ret.mcomp,
"aov"=ret.aov)
attr(ret,"between.factor")<-labelsB
attr(ret,"within.factor")<-labelsI
attr(ret,"subject.factor")<-labelsS
attr(ret,"comment")<-comment
class(ret)<-"rep.anova"
return(ret)
}
**関数 print.rep.anova [#mc2828ed]
rep.anovaの結果をうまいこと出力します。
普通に、
print(rep.aov(...))
とすればいいです。
print.rep.anova<-function(x){
## mauchly's spherical test
ms<-eps<-NULL;for(m in x[[2]])ms<-rbind(ms,c(m$stat,m$...
## construct table
tms1<-tms2<-c();tbl<-NULL
for(r in x[[1]]){
eps<-rbind(eps,attr(r,"eps"))
tms1<-c(tms1,attr(r,"within"))
tr<-r[-dim(r)[1],]
tms2<-c(tms2,sapply(paste(rownames(tr),attr(r,"withi...
tbl<-rbind(tbl,r[-dim(r)[1],-1])
}
mtest<-data.frame(cbind(ms,eps))
names(mtest)<-c("Mauchley's W"," p value"," G-G epsilo...
rownames(mtest)<-tms1
tbl<-data.frame(tbl)
names(tbl)<-c("F stat."," DF^", " DF_", " Pr.(>F)", " ...
rownames(tbl)<-tms2
## rename error term in aov table
tbl.aov<-NULL
tx<-summary(x[[4]])
for(i in 1:length(tx)){
for(j in 1:length(tx[[i]])){
ttx<-tx[[i]][[j]]
if(dim(tx[[i]][[j]])[1]>0){
rownames(ttx)[rownames(ttx)=="Residuals"]<-sub("...
tbl.aov<-rbind(tbl.aov,ttx)
}
}
}
## reformat aov table and add G-G and H-F p.val
gg<-hf<-c();rt<-sapply(strsplit(tms2,":"),function(s)p...
for(i in 1:nrow(tbl.aov)){
rn<-paste(sort(strsplit(sub(" *$","",rownames(tbl.ao...
if(any(rn==rt)){gg<-c(gg,tbl[rn,5]);hf<-c(hf,tbl[rn,...
else{gg<-c(gg,NA);hf<-c(hf,NA)}
}
tbl.aov<-cbind(tbl.aov,"G-G Pr."=gg," H-F Pr."=hf)
## digits...
for(i in 1:4)mtest[,i]<-sprintf("%.4f",mtest[,i])
for(i in c(2,3))tbl.aov[,i]<-ifelse(is.na(tbl.aov[,i])...
for(i in c(4,5,6,7))tbl.aov[,i]<-ifelse(is.na(tbl.aov[...
## output
cat("\n\nRepeated measured ANOVA\n\n")
if(!is.null(cm<-attr(x,"comment")))cat(cm,"\n\n")
cat("Between [",paste(attr(x,"between.factor"),collaps...
cat("Within [",paste(attr(x,"within.factor"),collapse=...
cat("Subject [",attr(x,"subject.factor"),"]\n\n")
cat("Mauchley's test of sphericity\n")
cat("-------------------------------------------------...
print(mtest)
cat("\n")
cat("Analysis of Variance Table\n")
cat("-------------------------------------------------...
print(as.matrix(tbl.aov),na.print="",quote=F,right=T)
cat("\n")
for(r in x[[3]])print(r,comment=paste("Main effect of ...
}
**関数 pairwise.ryan [#a24df7e1]
ライアン法による多重比較
pairwise.ryan<-function(x,MSe=NULL,dfe=NULL,alpha=0.05){
mx<-mean(x)
mc<-apply(x,2,mean)
ns<-dim(x)[1]
nc<-dim(x)[2]
names<-colnames(x)
if(is.null(names))names<-1:nc
if(is.null(dfe))dfe<-length(x)-ns-nc+1
if(is.null(MSe))MSe<-(sum((x-mx)^2)-sum(ns*(mc-mx)^2)-...
ret<-list()
ret$method<-"Ryan's method"
ret$data.name<-deparse(substitute(x))
ts<-matrix(nrow=nc-1,ncol=nc-1,dimnames=list(names[2:n...
ps<-matrix(nrow=nc-1,ncol=nc-1,dimnames=list(names[2:n...
aps<-matrix(nrow=nc-1,ncol=nc-1,dimnames=list(names[2:...
table<-NULL
for(i in 1:(nc-1)){
for(j in (i+1):nc){
d<-abs(mc[i]-mc[j])
t0<-d/(sqrt(MSe)*sqrt(2/ns))
r<-length(which(min(mc[i],mc[j])<=mc&mc<=max(mc[i]...
ap<-2*alpha/(nc*(r-1))
p0<-2*(1-pt(t0,dfe))
ps[i,j-1]<-p0
aps[i,j-1]<-ap
table<-rbind(table,data.frame(paste(names[i],"-",n...
}
}
table<-data.frame(table)
names(table)<-c("pair","r","nominal level","t","p","si...
rownames(table)<-1:dim(table)[1]
ret$p.value<-ps
ret$stat<-ts
ret$adjusted.p.value<-aps
ret$table<-table
ret$means<-mc
attr(ret,"class")<-"pairwise.ryan"
return(ret)
}
**関数 print.pairwise.ryan [#p312c125]
pairwise.ryanの結果をうまく出力します。
print.pairwise.ryan<-function(x,comment=NULL){
cat("\n\tPairwise comparisons using", x$method, "\n\n")
if(!is.null(comment))cat(comment)
cat("data: ", x$data.name, "[",paste(x$means,","),"]\n...
tbl<-x$table
tbl[,3]<-sprintf("%.4f",tbl[,3])
tbl[,4]<-sprintf("%.4f",tbl[,4])
tbl[,5]<-sprintf("%.4f",tbl[,5])
print(tbl)
}
**使用例 [#v3043d61]
言語Rによる分散分析( http://mat.isc.chubu.ac.jp/R/tech.ht...
Anova4 on the web ( http://www.hju.ac.jp/~kiriki/anova4/ ...
但し、G-G補正、H-F補正済みp値については未確認なので、どな...
***one-way within subject [#r05ec03c]
> if(1){
+ RB <- data.frame(
+ a = factor(c(rep(1,8), rep(2,8), re...
+ s = factor(rep(1:8, 4)),
+ result = c(9,7,8,8,12,11,8,13, 6,5,...
+ 10,13,8,13,12,14,14,16, 9,11,13,1...
+ (rep.aov(result ~ I(a) + S(s), RB))
+ }
Repeated measured ANOVA
Between [ ] Within [ a ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
a 0.9751 0.9996 0.9832 1.7774
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. H-F...
Error(s) 7 77.00 11.00 ...
a 3 217.50 72.50 21.4437 0.0000 0.0000 0....
Error(s:a) 21 71.00 3.38 ...
Pairwise comparisons using Ryan's method
Main effect of : a
data: dat.tmp [ 9.5 , 6.5 , 12.5 , 13 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0250 3.2631 0.0037 s.
2 1 - 3 2 0.0250 3.2631 0.0037 s.
3 1 - 4 3 0.0125 3.8070 0.0010 s.
4 2 - 3 3 0.0125 6.5262 0.0000 s.
5 2 - 4 4 0.0083 7.0701 0.0000 s.
6 3 - 4 2 0.0250 0.5439 0.5923 n.s.
***two-way within subject [#race7b7b]
> if(1){
+ RBFpq <- data.frame(
+ a = factor(c(rep(1, 20), rep(2, ...
+ b = factor(rep(c(rep(1, 5), rep(...
+ s = factor(rep(1:5, 8)),
+ result = c(3,3,1,3,5, 4,3,4,5,7,...
+ 3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,...
+ (rep.aov(result ~ I(a,b) + S(s), RBFpq))
+ }
Repeated measured ANOVA
Between [ ] Within [ a , b ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
a 1.0000 NaN 1.0000 1.0000
b 0.6817 0.9618 0.7934 2.0373
a:b 0.2482 0.6036 0.6692 1.3449
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. H...
Error(s) 4 38.15 9.54 ...
a 1 16.90 16.90 8.0958 0.0466 0.0466 ...
Error(s:a) 4 8.35 2.09 ...
b 3 19.70 6.57 6.0383 0.0095 0.0173 ...
Error(s:b) 12 13.05 1.09 ...
a:b 3 30.50 10.17 7.0725 0.0054 0.0169 ...
Error(s:a:b) 12 17.25 1.44 ...
Pairwise comparisons using Ryan's method
Main effect of : b
data: dat.tmp [ 3.5 , 4.4 , 4.9 , 5.4 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0250 1.9298 0.0776 n.s.
2 1 - 3 3 0.0125 3.0019 0.0110 s.
3 1 - 4 4 0.0083 4.0740 0.0015 s.
4 2 - 3 2 0.0250 1.0721 0.3048 n.s.
5 2 - 4 3 0.0125 2.1442 0.0532 n.s.
6 3 - 4 2 0.0250 1.0721 0.3048 n.s.
***two-way mixed: between (1) x within (1) [#x5bb0565]
> if(1){
+ SPFp.q <- data.frame(
+ a = factor(c(rep(1,20), rep(2,2...
+ b = factor(rep(c(rep(1,5), rep(...
+ s = factor(c(rep(1:5, 4), rep(6...
+ result = c(3,3,1,3,5, 4,3,4,5,7...
+ 3,5,2,4,6, 2,6,3,6,4, 3,2,3,6...
+ rep.aov(result ~ I(b) + B(a) + S(s), SPFp.q)
+ }
Repeated measured ANOVA
Between [ a ] Within [ b ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
b 0.8428 0.9502 0.8896 1.3767
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. H...
a 1 16.90 16.90 2.9075 0.1266 ...
Error(s:a) 8 46.50 5.81 ...
b 3 19.70 6.57 5.2013 0.0066 0.0092 ...
b:a 3 30.50 10.17 8.0528 0.0007 0.0012 ...
Error(s:a:b) 24 30.30 1.26 ...
Pairwise comparisons using Ryan's method
Main effect of : b
data: dat.tmp [ 3.5 , 4.4 , 4.9 , 5.4 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0250 1.7911 0.0859 n.s.
2 1 - 3 3 0.0125 2.7861 0.0103 s.
3 1 - 4 4 0.0083 3.7811 0.0009 s.
4 2 - 3 2 0.0250 0.9950 0.3296 n.s.
5 2 - 4 3 0.0125 1.9901 0.0581 n.s.
6 3 - 4 2 0.0250 0.9950 0.3296 n.s.
***three-way mixed : between (2) x within (1) [#ic872892]
> if(1){
+ SPFpq.r <- data.frame(
+ a = factor(c(rep(1, 24), rep(2...
+ b = factor(rep(c(rep(1, 12), r...
+ c = factor(rep(c(rep(1, 4), re...
+ s = factor(c(rep(1:4, 3), rep(...
+ result = c(2,6,5,7, 5,7,9,9, ...
+ 6,6,8,10, 3,6,5,8, 6,7,5,6,
+ 1,2,5,2, 2,1,4,5, 1,5,3,5,
+ 5,3,4,6, 5,6,6,7, 5,5,9,7))
+ (rep.aov(result ~ S(s)+I(c)+B(a,b), SPFpq.r))
+ }
Repeated measured ANOVA
Between [ a , b ] Within [ c ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
c 0.8178 0.3309 0.8459 0.9698
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. ...
a 1 96.33 96.33 15.6216 0.0019 ...
b 1 3.00 3.00 0.4865 0.4988 ...
a:b 1 56.33 56.33 9.1351 0.0106 ...
Error(s:a:b) 12 74.00 6.17 ...
c 2 33.50 16.75 9.5714 0.0009 0.0018 ...
c:a 2 6.17 3.08 1.7619 0.1932 0.1995 ...
c:b 2 24.50 12.25 7.0000 0.0040 0.0067 ...
c:a:b 2 41.17 20.58 11.7619 0.0003 0.0007 ...
Error(s:a:b:c) 24 42.00 1.75 ...
Pairwise comparisons using Ryan's method
Main effect of : c
data: dat.tmp [ 4.875 , 5.5 , 6.875 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0333 1.3363 0.1940 n.s.
2 1 - 3 3 0.0167 4.2762 0.0003 s.
3 2 - 3 2 0.0333 2.9399 0.0072 s.
***three-way mixed : between (1) x within (2) [#nf2f0213]
> if(1){
+ SPFp.qr <- data.frame(
+ a = factor(c(rep(1, 24), rep(2...
+ b = factor(rep(c(rep(1, 12), r...
+ c = factor(rep(c(rep(1, 4), re...
+ s = factor(c(rep(1:4, 6), rep(...
+ result = c(2,6,5,7, 5,7,9,9, ...
+ 6,6,8,10, 3,6,5,8, 6,7,5,6,
+ 1,2,5,2, 2,1,4,5, 1,5,3,5,
+ 5,3,4,6, 5,6,6,7, 5,5,9,7))
+ (r<-rep.aov(result ~ S(s)+I(b,c)+B(a), SPFp.qr))
+ }
Repeated measured ANOVA
Between [ a ] Within [ b , c ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
b 1.0000 NaN 1.0000 1.0000
c 0.6719 0.3700 0.7529 0.9503
b:c 0.8354 0.6378 0.8586 1.1699
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr. ...
a 1 96.33 96.33 8.7576 0.0253 ...
Error(s:a) 6 66.00 11.00 ...
b 1 3.00 3.00 2.2500 0.1843 0.1843 ...
b:a 1 56.33 56.33 42.2500 0.0006 0.0006 ...
Error(s:a:b) 6 8.00 1.33 ...
c 2 33.50 16.75 25.1250 0.0001 0.0003 ...
c:a 2 6.17 3.08 4.6250 0.0324 0.0489 ...
Error(s:a:c) 12 8.00 0.67 ...
b:c 2 24.50 12.25 4.3235 0.0385 0.0477 ...
b:c:a 2 41.17 20.58 7.2647 0.0086 0.0128 ...
Error(s:a:b:c) 12 34.00 2.83 ...
Pairwise comparisons using Ryan's method
Main effect of : c
data: dat.tmp [ 4.875 , 5.5 , 6.875 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0333 2.1651 0.0512 n.s.
2 1 - 3 3 0.0167 6.9282 0.0000 s.
3 2 - 3 2 0.0333 4.7631 0.0005 s.
***three-way within subject [#ccafdb19]
> if(1){
+ RBFpqr <- data.frame(
+ a = factor(c(rep(1, 24), rep(2,...
+ b = factor(rep(c(rep(1, 12), re...
+ c = factor(rep(c(rep(1, 4), rep...
+ s = factor(rep(1:4, 12)),
+ result = c(2,6,5,7, 5,7,9,9, 9...
+ 6,6,8,10, 3,6,5,8, 6,7,5,6,
+ 1,2,5,2, 2,1,4,5, 1,5,3,5,
+ 5,3,4,6, 5,6,6,7, 5,5,9,7))
+ (rep.aov(result ~ S(s)+I(a,b,c), RBFpqr))
+ }
Repeated measured ANOVA
Between [ ] Within [ a , b , c ] Subject [ s ]
Mauchley's test of sphericity
--------------------------------------------------------...
Mauchley's W p value G-G epsilon H-F epsilon
a 1.0000 NaN 1.0000 1.0000
b 1.0000 NaN 1.0000 1.0000
c 0.8700 0.8700 0.8850 2.0650
a:b 1.0000 NaN 1.0000 1.0000
a:c 0.5422 0.5422 0.6860 1.0710
b:c 0.9838 0.9838 0.9840 2.8450
a:b:c 0.4831 0.4831 0.6593 0.9736
Analysis of Variance Table
--------------------------------------------------------...
Df Sum Sq Mean Sq F value Pr(>F) G-G Pr....
Error(s) 3 60.33 20.11 ...
a 1 96.33 96.33 51.0000 0.0057 0.0057...
Error(s:a) 3 5.67 1.89 ...
b 1 3.00 3.00 1.4211 0.3189 0.3189...
Error(s:b) 3 6.33 2.11 ...
c 2 33.50 16.75 60.3000 0.0001 0.0002...
Error(s:c) 6 1.67 0.28 ...
a:b 1 56.33 56.33 101.4000 0.0021 0.0021...
Error(s:a:b) 3 1.67 0.56 ...
a:c 2 6.17 3.08 2.9211 0.1301 0.1625...
Error(s:a:c) 6 6.33 1.06 ...
b:c 2 24.50 12.25 5.3780 0.0459 0.0471...
Error(s:b:c) 6 13.67 2.28 ...
a:b:c 2 41.17 20.58 6.0738 0.0361 0.0670...
Error(s:a:b:c) 6 20.33 3.39 ...
Pairwise comparisons using Ryan's method
Main effect of : c
data: dat.tmp [ 4.875 , 5.5 , 6.875 , ]
pair r nominal level t p sig.
1 1 - 2 2 0.0333 3.3541 0.0153 s.
2 1 - 3 3 0.0167 10.7331 0.0000 s.
3 2 - 3 2 0.0333 7.3790 0.0003 s.
**TODO [#s6698f53]
- ライアン法の多重比較で混合要因ANOVAのときに被験者数が異...
- 交互作用が有意のときに下位検定を行いたいんですが、まだ...
- コードが汚いです。
ページ名:
WWW を検索
OKADAJP.ORG を検索