Rで簡単（二元配置分散分析：データには対応なし、繰り返し数は5）
投稿2016.11.17 最終更新日2019.6.5
検定は帰無仮説が大切です！
商品のの主効果：帰無仮説＝カップ麺でもインスタント麺でも評価の母平均は等しい
スープの主効果：帰無仮説＝スープが違っても評価の母平均は等しい
商品とスープの主効果：帰無仮説＝商品とスープの組合せと評価は関係ない
商品 スープ 評価 カップ とんこつ 10 カップ とんこつ 11 カップ とんこつ 11 カップ とんこつ 9 カップ とんこつ 9 カップ 醤油 11 カップ 醤油 8 カップ 醤油 10 カップ 醤油 8 カップ 醤油 8 カップ みそ 7 カップ みそ 3 カップ みそ 5 カップ みそ 2 カップ みそ 3 インスタント とんこつ 11 インスタント とんこつ 11 インスタント とんこつ 10 インスタント とんこつ 9 インスタント とんこつ 10 インスタント 醤油 9 インスタント 醤油 7 インスタント 醤油 9 インスタント 醤油 8 インスタント 醤油 8 インスタント みそ 4 インスタント みそ 4 インスタント みそ 2 インスタント みそ 3 インスタント みそ 3
コピーしてRに読み込みます
x<-read.table("clipboard",header=T) summary(aov(x$評価~x$商品*x$スープ)) Df Sum Sq Mean Sq F value Pr(>F) x$商品 1 1.63 1.63 1.077 0.310 x$スープ 2 231.67 115.83 76.374 3.93e-11 *** x$商品:x$スープ 2 1.67 0.83 0.549 0.584 Residuals 24 36.40 1.52 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
スープによる主効果あり！
次の式でも同じ結果になります
summary(aov(x$評価~x$商品+x$スープ+x$商品:x$スープ))
ベクトルを使用した方法
評価<-c(x$評価) 商品<-factor(c(rep("カップ",15),rep("インスタント",15))) スープ<-factor(rep(c(rep("とんこつ",5),rep("醤油",5),rep("みそ",5)),2)) summary(aov(x$評価~x$商品*x$スープ))
検証するために図を描きましょう！
par(mfrow=c(1,2)) interaction.plot(商品,スープ,評価) interaction.plot(スープ,商品,評価) par(mfrow=c(1,1))
相関行列と共分散行列
相関行列
x1 <- c(93, 89, 115, 90, 75) x2 <- c(121, 136, 121, 161, 125) x3 <- c(101, 115, 118, 122, 106) #データフレームにすること dat <- data.frame(x1, x2, x3) cor(dat) x1 x2 x3 x1 1.0000000 -0.2220449 0.3981155 x2 -0.2220449 1.0000000 0.6767614 x3 0.3981155 0.6767614 1.0000000
散布図
pairs(dat)
共分散行列
これもデータフレームで var(dat) x1 x2 x3 x1 207.80 -54.15 49.80 x2 -54.15 286.20 99.35 x3 49.80 99.35 75.30
参考
install.packages("psych") library(psych) x <- rnorm(100,5,3) y <- rnorm(100,15,8) z <- rnorm(100,2,0.5) dat2 <- data.frame(x,y,z) pairs.panels(dat2)
対応のある一元配置分散分析（繰り返しのない二元配置分散分析）、多重比較
以下のサイトに移転いたしました
生存時間解析の基本
パッケージのインストール
install.packages("survival") library(survival)
survivalのデータセットheartを使用します
head(heart) start, stop, event: Entry and exit time and status for this interval of time age: age-48 years year: year of acceptance (in years after 1 Nov 1967) surgery: prior bypass surgery 1=yes#バイパス手術あり transplant: received transplant 1=yes#移植を受けている id: patient id
バイパス手術あり＆なしの中央値
pfit <- survfit(Surv(stop-start,event) ~ surgery, heart) Call: survfit(formula = Surv(stop - start, event) ~ surgery, data = heart) n events median 0.95LCL 0.95UCL surgery=0 143 66 102 64 280 surgery=1 29 9 897 322 NA
バイパス手術あり＆なしで生存曲線を描く、赤(あり)vs黒(なし)
plot(pfit, col = 1:2)
バイパス手術あり＆なしの有意差検定 log-rank test
rho = 0 （デフォルト）this is the log-rank or Mantel-Haenszel test
rho = 1 it is equivalent to the Peto & Peto modification of the Gehan-Wilcoxon test.
diff_surg <- survdiff(Surv(stop-start,event)~surgery , heart, rho = 0) Call: survdiff(formula = Surv(stop - start, event) ~ surgery, data = heart, rho = 0) N Observed Expected (O-E)^2/E (O-E)^2/V surgery=0 143 66 57.8 1.15 5.17 surgery=1 29 9 17.2 3.87 5.17 Chisq= 5.2 on 1 degrees of freedom, p= 0.02
p値のみを算出
pchisq(diff_surg <- survdiff(Surv(stop-start,event)~surgery , heart, rho = 0)$chisq, 1, lower.tail=FALSE ) [1] 0.02294185
surgeryを1にすれば、全体の曲線と95％信頼区間をプロットする
pfit2 <- survfit(Surv(stop-start,event) ~ 1, heart) plot(pfit2, col = c(1,4,4))
年齢で比較 48歳以上と48歳未満で検定
diff_age1 <- survdiff(Surv(stop-start,event)~ age<0 , heart, rho = 0) Call: survdiff(formula = Surv(stop - start, event) ~ age < 0, data = heart, rho = 0) N Observed Expected (O-E)^2/E (O-E)^2/V age < 0=FALSE 81 39 32.5 1.283 2.29 age < 0=TRUE 91 36 42.5 0.983 2.29 Chisq= 2.3 on 1 degrees of freedom, p= 0.1
年齢で比較 50歳以上と40以下
sub <- subset(heart, age>=2 | age <= -8) diff_age2 <- survdiff(Surv(stop-start,event)~ age>=2 , sub, rho = 0) Call: survdiff(formula = Surv(stop - start, event) ~ age >= 2, data = sub, rho = 0) N Observed Expected (O-E)^2/E (O-E)^2/V age >= 2=FALSE 33 10 18.2 3.67 6.82 age >= 2=TRUE 57 31 22.8 2.92 6.82 Chisq= 6.8 on 1 degrees of freedom, p= 0.009
カテゴリー変数がある場合の層別 カテゴリー変数名をＡとします。Ａには1,2,3,4が含まれており、1と3を比較します
x<-1; y<-3 cat <- subset(heart, Ａ==x | Ａ==y, data=heart ) diff_c <- survdiff(Surv(stop-start,event)~"カテゴリ－変数名" , cat)
COX比例ハザードモデル
cox_age_sur <- coxph(Surv(stop-start,event) ~ age + surgery, data=heart) summary(cox_age_sur) Call: coxph(formula = Surv(stop - start, event) ~ age + surgery, data = heart) n= 172, number of events= 75 coef exp(coef) se(coef) z Pr(>|z|) age 0.02941 1.02985 0.01384 2.125 0.0336 * surgery -0.82104 0.43997 0.35872 -2.289 0.0221 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 exp(coef) exp(-coef) lower .95 upper .95 age 1.03 0.971 1.0023 1.0582 surgery 0.44 2.273 0.2178 0.8887 Concordance= 0.603 (se = 0.039 ) Rsquare= 0.062 (max possible= 0.975 ) Likelihood ratio test= 10.96 on 2 df, p=0.004 Wald test = 9.77 on 2 df, p=0.008 Score (logrank) test = 10.17 on 2 df, p=0.006
生存率の予測
fit_surv <- survfit(cox_age_sur) plot(fit_surv,col = c(2,3,3)) haz50 <- summary(fit_surv) Call: survfit(formula = cox_age_sur) time n.risk n.event survival std.err lower 95% CI upper 95% CI 0.5 172 1 0.995 0.00537 0.9842 1.000 1.0 171 2 0.984 0.00926 0.9659 1.000 2.0 166 3 0.967 0.01314 0.9421 0.994 3.0 160 4 0.945 0.01703 0.9124 0.979 5.0 150 1 0.939 0.01792 0.9050 0.975 6.0 147 2 0.928 0.01962 0.8900 0.967 8.0 144 1 0.922 0.02042 0.8826 0.963 9.0 141 1 0.916 0.02121 0.8750 0.958 10.0 140 1 0.910 0.02197 0.8675 0.954
subset関数で層別化したデータセット
以下のようなデータセットを男女で層別したセットに変更します
id 年齢 性 身長 1 39 女 162.6 2 39 女 161.4 3 41 男 174.1 4 41 女 154.7 5 42 男 184.7 6 43 男 169.3 7 43 男 163.2 8 43 男 164.4 9 44 男 167.8 10 44 男 168.5 11 44 女 154.8 12 45 男 158.5 13 45 男 162.9 14 45 男 161.4 15 45 男 165.9 16 46 男 171.4 17 46 男 172.6 18 47 男 167.6 19 47 男 171.5 20 48 男 166.4 21 48 男 169.2 22 48 男 165.3 23 48 男 167.2 24 48 男 170.2 31 50 女 153.6 32 50 女 163.6 34 51 女 156.6 35 51 女 152.4 36 52 女 146.8 37 52 女 160.7
上記のデータをコピーしてＲにペーストします
data <- read.table("clipboard", header=T)
男性のみのデータセット data_m、女性のみのデータセット data_f
data_m <- subset(data, data$性=="男") data_f <- subset(data, data$性=="女")
男性で身長160cm台のデータセット data_m_160
mh01 <- subset( data, 性=="男" & 身長 >= 160 & 身長 < 170 )
男性で身長165cm以下と170cm以上の群
mh02 <- subset( data, 性=="男" & 身長 <= 165 | 身長 >= 170 )
男性で身長160cm台の人数
それぞれのデータセットのｎ(行数)を確認する
nrow(mh01);nrow(mh02)
身長の差の検定 年齢40未満 vs 年齢48より上
d45 <- subset(data,年齢<45) d48 <- subset(data,年齢>48) t.test(d45$身長,d48$身長) #subset使用せずに、次のように書くことも可能 t.test(data[data[,2]<45,4],data[data[,2]>48,4])
