anova_note:code01
This is an old revision of the document!
rm(list=ls())
rnorm2 <- function(n,mean,sd) {
mean+sd*scale(rnorm(n))
}
ss <- function(x) {
sum((x-mean(x))^2)
}
set.seed(10)
n <- 30
n.o <- n.p <- n
o <- rnorm(n.o, 100, 10)
p <- rnorm(n.p, 104, 10)
# old way
m.o <- mean(o)
m.p <- mean(p)
df.o <- n.o - 1
df.p <- n.p - 1
diff <- m.o - m.p
pv <- (ss(o)+ss(p))/(df.o+df.p)
se <- sqrt((pv/n.o) + (pv/n.p))
t.cal <- diff/se
t.cal
pt(t.cal, df.o+df.p) * 2
t.out$statistic
t.out$p.value
t.out <- t.test(o,p, var.equal=T)
t.out
#
comb <- list(o = o, p = p)
op <- stack(comb)
head(op)
colnames(op)[1] <- "values"
colnames(op)[2] <- "group"
op$group <- factor(op$group)
head(op)
boxplot(op$values~op$group)
plot(op$values~op$group)
boxplot(op$values~op$group, main="values by group",
yaxt="n", xlab="value", horizontal=TRUE,
col=terrain.colors(2))
abline(v=mean(op$values), col="red", lwd=3)
legend("topleft", inset=.05, title="group",
c("o","p"), fill=terrain.colors(2), horiz=TRUE)
m.tot <- mean(op$values)
m.o <- mean(o)
m.p <- mean(p)
min.x <- min(op$values)
max.x <- max(op$values)
br <- seq(floor(min.x), ceiling(max.x), by = 1)
hist(o, breaks=br,
col=rgb(1,1,1,.5))
abline(v=m.o, col="black", lwd=3)
hist(p, add=T, breaks=br,
col=rgb(.5,1,1,.5))
abline(v=m.p, col="blue", lwd=3)
abline(v=m.tot, col='red', lwd=3)
hist(o, breaks=br,
col=rgb(1,1,1,.5))
abline(v=m.o, col="black", lwd=3)
abline(v=m.tot, col='red', lwd=3)
hist(p, breaks=br,
col=rgb(.5,1,1,.5))
abline(v=m.p, col="blue", lwd=3)
abline(v=m.tot, col='red', lwd=3)
hist(o, breaks=br,
col=rgb(1,1,1,.5))
hist(p, add=T, breaks=br,
col=rgb(.5,1,1,.5))
abline(v=m.tot, col='red', lwd=3)
ss.tot <- ss(op$values)
df.tot <- length(op$values)-1
ss.tot/df.tot
var(op$values)
ss.tot
m.tot
m.o
m.p
ss.o
ss.p
hist(o, breaks=br,
col=rgb(1,1,1,.5))
abline(v=m.o, col="black", lwd=3)
hist(p, add=T, breaks=br,
col=rgb(.5,1,1,.5))
abline(v=m.p, col="blue", lwd=3)
abline(v=m.tot, col='red', lwd=3)
ss.bet <- n.o*(m.tot-m.o)^2 + n.p*(m.tot-m.p)^2
ss.bet
hist(o, breaks=br,
col=rgb(1,1,1,.5))
abline(v=m.o, col="black", lwd=3)
ss.o <- ss(o)
ss.o
hist(p, breaks=br,
col=rgb(.5,1,1,.5))
abline(v=m.p, col="blue", lwd=3)
ss.p <- ss(p)
ss.p
ss.wit <- ss.o+ss.p
ss.wit
ss.bet ss.wit ss.bet+ss.wit ss.tot
df.tot <- length(op$values)-1 df.bet <- nlevels(op$group) - 1 df.wit <- (n.o-1)+(n.p-1) df.tot df.bet df.wit
ms.tot <- ss.tot / df.tot ms.bet <- ss.bet / df.bet ms.wit <- ss.wit / df.wit f.cal <- ms.bet / ms.wit f.cal p.val <- pf(f.cal, df1=df.bet, df2=df.wit, lower.tail = F) p.val summary(aov(op$values~op$group)) t.test(o,p, var.equal = T) diff <- m.o - m.p ssp <- (ss.o + ss.p) / (df.o + df.p) se <- sqrt(ssp/n.o+ssp/n.p) t.cal <- diff/se t.cal p.t.cal <- pt(abs(t.cal), df=df.o+df.p, lower.tail = F)*2 p.t.cal t.cal^2 f.cal df.bet df.wit f.cal
anova_note/code01.1765104920.txt.gz · Last modified: by hkimscil