1 SE, Standard Error #
覈暑レ 蠍語企ゼ 5 豸′.
x <- c(76.2, 76.3, 76.1, 76,3, 76.4)
se <- sd(x)/sqrt(length(x)) #0.6009252
se #1.643168
2,000覈 覈螻 煙 .
set.seed(1000)
x <- rnorm(2000, mean = 70, sd = 10)
summary(x)
蟆郁骸
> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
36.38 63.33 69.89 69.94 76.54 99.19
譴れ姶..
mu <- c()
for(i in 1:100){
mu <- c(mean(sample(x, 5)), mu)
}
se <- sd(mu) / sqrt(5)
mean(mu) #69.82971
se #1.935552
4 One Sample t-test #
2 1覦 2011 襭 蠏 蟆螳 2.1螳 伎. 2012 10覈 覓伎襦 覦 蟆 螳 譟一. 2011螻 2012 るジ螳?
x <- c(3.3, 2.8, 3.0, 2.7, 2.7, 2.0, 1.9, 3.4, 1.4, 1.4)
t.test(x, mu=2.1)
蠏 蟆
> shapiro.test(x)
Shapiro-Wilk normality test
data: x
W = 0.9085, p-value = 0.2711
- 螳
- 譴 0.05手 , p-value螳 0.2711企襦 蠏覓願 蠍郁讌 覈詩. 讀, 蠏覿.
> t.test(x, mu=2.1)
One Sample t-test
data: x
t = 1.5454, df = 9, p-value = 0.1567
alternative hypothesis: true mean is not equal to 2.1
95 percent confidence interval:
1.933026 2.986974
sample estimates:
mean of x
2.46
- 螳
- 譴 0.05手 , p-value螳 0.1567企襦 蠏覓願 蠍郁讌 覈詩.
螳 蟲螳豢 .
> t.test(x)
One Sample t-test
data: x
t = 10.6, df = 9, p-value = 2.269e-06
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
1.93 2.99
sample estimates:
mean of x
2.46
5 Two Sample t-test #
x1 <- c(15,10,13,7,9,8,21,9,14,8)
x2 <- c(15,14,12,8,14,7,16,10,15,12)
蠏 蟆 --> x1, x2螳 譴 0.05 蠏覿.
> shapiro.test(x1)
Shapiro-Wilk normality test
data: x1
W = 0.8666, p-value = 0.09131
> shapiro.test(x2)
Shapiro-Wilk normality test
data: x2
W = 0.9125, p-value = 0.2986
覿一 狩螳?
> var.test(x1, x2)
F test to compare two variances
data: x1 and x2
F = 1.9791, num df = 9, denom df = 9, p-value = 0.3237
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.491579 7.967821
sample estimates:
ratio of variances
1.979094
- 螳
- 蠏覓願: x1螻 x2 覿一 谿願 .
- 襴所: x1螻 x2 覿一 谿願 .
- 譴 0.05手 , p-value螳 0.3237企襦 蠏覓願 蠍郁讌 覈詩. 讀, 覿一 谿願 .
讌 蠏 狩讌 蟆
> t.test(x1, x2, var.equal=T)
Two Sample t-test
data: x1 and x2
t = -0.5331, df = 18, p-value = 0.6005
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-4.446765 2.646765
sample estimates:
mean of x mean of y
11.4 12.3
- 螳
- 蠏覓願: x1螻 x2 蠏 谿願 .
- 襴所: x1螻 x2 蠏 谿願 .
- 譴 0.05手 , p-value螳 0.6005企襦 蠏覓願 蠍郁讌 覈詩. 讀, 讌螳 蠏 谿願 .
豸 蟆 蟆曙 (殊 讀, x1 る 蟆)
> t.test(x1, x2, alternative="less", var.equal=T)
Two Sample t-test
data: x1 and x2
t = -0.5331, df = 18, p-value = 0.3002
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf 2.027436
sample estimates:
mean of x mean of y
11.4 12.3
8 殊覿磯 #
#http://code.google.com/p/sonya/source/browse/trunk/r-project/sample/PlantGrowth.csv
plantGrowth = read.csv("c:\\data\\PlantGrowth.csv")
head(plantGrowth)
boxplot(weight ~ group, data=plantGrowth)
out <- lm(weight ~ group, data=plantGrowth)
summary(out)
anova(out)
> summary(out)
Call:
lm(formula = weight ~ group, data = plantGrowth)
Residuals:
Min 1Q Median 3Q Max
-1.0710 -0.4180 -0.0060 0.2627 1.3690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.0320 0.1971 25.527 <2e-16 ***
grouptrt1 -0.3710 0.2788 -1.331 0.1944
grouptrt2 0.4940 0.2788 1.772 0.0877 .
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 0.6234 on 27 degrees of freedom
Multiple R-squared: 0.2641, Adjusted R-squared: 0.2096
F-statistic: 4.846 on 2 and 27 DF, p-value: 0.01591
> anova(out)
Analysis of Variance Table
Response: weight
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.7663 1.8832 4.8461 0.01591 *
Residuals 27 10.4921 0.3886
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
- 螳
- 蠏覓願: 3蠏碁９ 蠏 螳.
- 襴所: 3蠏碁９ 蠏 るゴ.
- 譴 0.05手 , p-value螳 0.01591企襦 蠏覓願 蠍郁. 讀, 3蠏碁９ 蠏 谿願 .
- 蠏碁, 蠏讌 企 讌 .
par(mfrow=c(2,2))
plot(out)
蠏 -> p-value = 0.4379企襦 蠏覿.
> shapiro.test(resid(out))
Shapiro-Wilk normality test
data: resid(out)
W = 0.9661, p-value = 0.4379
焔一 -> p-value = 0.1714襦 蠏覓願 讌讌. 讀, 焔
#library("lmtest")
> bptest(out)
studentized Breusch-Pagan test
data: out
BP = 3.5273, df = 2, p-value = 0.1714
襴曙
> dwtest(out) #library("lmtest")
Durbin-Watson test
data: out
DW = 2.704, p-value = 0.9502
alternative hypothesis: true autocorrelation is greater than 0
- d糾
- 0 : 蠍一蟯, p-value = 1
- 2 : 襴, p-value = 0
- 4 : 蠍一蟯, p-value = -1
- 旧朱 DW螳 1覲企 蟇磯 3覲企 覃 蠍一蟯 ろ り 朱, 1.5~2.5 伎 蟆曙 襴曙企手
- 襴.
蠏, 焔一, 襴曙, Q-Q Plot 讌煙 誤 覈 蟆 誤. 3螳 蠏碁９ 蠏 谿 蟆 蟆郁骸 谿願 . 企 蠏碁９朱Μ 谿願 覦讌 覲伎.
覦覯1: Dunnett -> 蠏 谿願 譟壱 覲伎譴. (control 觜 觜蟲覯)
install.packages("multcomp")
library("multcomp")
out <- lm(weight ~ group, data=PlantGrowth)
dunnett <- glht(out, linfct=mcp(group="Dunnett")) #蠍一 group plantGrowth$group 企.
summary(dunnett)
plot(dunnett)
95% 襤郁規螳 0
> summary(dunnett)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Dunnett Contrasts
Fit: lm(formula = weight ~ group, data = PlantGrowth)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
trt1 - ctrl == 0 -0.3710 0.2788 -1.331 0.323
trt2 - ctrl == 0 0.4940 0.2788 1.772 0.153
(Adjusted p values reported -- single-step method)
- trt1 - ctrl p-value = 0.323 -> 蠏覓願 讌讌, 蠏 谿
- trt2 - ctrl p-value = 0.153 -> 蠏覓願 讌讌, 蠏 谿
覦覯2: Tukey -> 蠏 谿願 譟壱螻 蠏 谿願 譟壱 覈 覲伎譴.
install.packages("multcomp")
library("multcomp")
out <- lm(weight ~ group, data=PlantGrowth)
tukey <- glht(out, linfct=mcp(group="Tukey")) #蠍一 group PlantGrowth$group 企.
summary(tukey)
plot(tukey)
> summary(tukey)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = weight ~ group, data = plantGrowth)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
trt1 - ctrl == 0 -0.3710 0.2788 -1.331 0.391
trt2 - ctrl == 0 0.4940 0.2788 1.772 0.198
trt2 - trt1 == 0 0.8650 0.2788 3.103 0.012 *
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
(Adjusted p values reported -- single-step method)
- trt1 - ctrl p-value = 0.323 -> 蠏覓願 讌讌, 蠏 谿
- trt2 - ctrl p-value = 0.153 -> 蠏覓願 讌讌, 蠏 谿
- trt2 - trt1 p-value = 0.153 -> 蠏覓願 蠍郁, 蠏 谿
9 伎覿磯 #
#http://code.google.com/p/sonya/source/browse/trunk/r-project/sample/warpbreaks.csv?r=653
warpbreaks = read.csv("c:\\data\\warpbreaks.csv")
伎壱 覲 wool螻 tension
> levels(warpbreaks$wool)
[1] "A" "B"
> levels(warpbreaks$tension)
[1] "L" "M" "H"
- wool -> 螳 襷.
- tension -> 螳 襷讌 . L, M, H 伎伎 .
tension 襯 覦襦 ′.
> warpbreaks$tension = factor(warpbreaks$tension, level = c("L", "M", "H"))
> levels(warpbreaks$tension)
[1] "L" "M" "H"
覿磯
out <- lm(breaks ~ wool*tension, data = warpbreaks)
summary(out)
> summary(out)
Call:
lm(formula = breaks ~ wool * tension, data = warpbreaks)
Residuals:
Min 1Q Median 3Q Max
-19.5556 -6.8889 -0.6667 7.1944 25.4444
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.556 3.647 12.218 2.43e-16 ***
woolB -16.333 5.157 -3.167 0.002677 **
tensionM -20.556 5.157 -3.986 0.000228 ***
tensionH -20.000 5.157 -3.878 0.000320 ***
woolB:tensionM 21.111 7.294 2.895 0.005698 **
woolB:tensionH 10.556 7.294 1.447 0.154327
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 10.94 on 48 degrees of freedom
Multiple R-squared: 0.3778, Adjusted R-squared: 0.3129
F-statistic: 5.828 on 5 and 48 DF, p-value: 0.0002772
- 蠏 谿願 (讌 蠏讌 讌 )
- 蟲語
- woolB:tensionM -> p-value = 0.005698, 譴 0.05 蠏覓願 蠍郁 => 蟲語 誤!!
- woolB:tensionH -> p-value = 0.154327, 譴 0.05 蠏覓願 讌讌
蠏 -> p-value = 0.8162企襦 蠏覿.
> shapiro.test(resid(out))
Shapiro-Wilk normality test
data: resid(out)
W = 0.9869, p-value = 0.8162
焔一 ->
p-value = 0.0006307襦 蠏覓願 蠍郁. 讀, 焔一 . 譬覲 breaks log() sqrt().
#library("lmtest")
> bptest(out)
studentized Breusch-Pagan test
data: out
BP = 21.5744, df = 5, p-value = 0.0006307
襴曙
> dwtest(out)
Durbin-Watson test
data: out
DW = 2.2376, p-value = 0.575
alternative hypothesis: true autocorrelation is greater than 0
- d糾
- 0 : 蠍一蟯, p-value = 1
- 2 : 襴, p-value = 0
- 4 : 蠍一蟯, p-value = -1
- 旧朱 DW螳 1覲企 蟇磯 3覲企 覃 蠍一蟯 ろ り 朱, 1.5~2.5 伎 蟆曙 襴曙企手
- 襴.
譬覲 覲(log) れ 蟆
> out <- lm(log(breaks) ~ wool*tension, data = warpbreaks)
> shapiro.test(resid(out))
Shapiro-Wilk normality test
data: resid(out)
W = 0.9729, p-value = 0.2583
> bptest(out)
studentized Breusch-Pagan test
data: out
BP = 4.8045, df = 5, p-value = 0.4402
> dwtest(out)
Durbin-Watson test
data: out
DW = 2.06, p-value = 0.3167
alternative hypothesis: true autocorrelation is greater than 0
企 譟壱 蠏 るジ螳?
library("multcomp")
out <- lm(log(breaks) ~ wool + tension, data=warpbreaks)
tukey1 <- glht(out, linfct=mcp(wool="Tukey"))
tukey2 <- glht(out, linfct=mcp(tension="Tukey"))
summary(tukey1)
summary(tukey2)
> summary(tukey1)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = log(breaks) ~ wool + tension, data = warpbreaks)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
B - A == 0 -0.1522 0.1063 -1.431 0.159
(Adjusted p values reported -- single-step method)
> summary(tukey2)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = log(breaks) ~ wool + tension, data = warpbreaks)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
M - L == 0 -0.2871 0.1302 -2.205 0.08018 .
H - L == 0 -0.4893 0.1302 -3.758 0.00133 **
H - M == 0 -0.2022 0.1302 -1.553 0.27550
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
(Adjusted p values reported -- single-step method)
![[http]](/moniwiki/imgs/http.png)