1-4

在本地建立文件

new directoryempty projectproject1

在官网下载文件：

 点击进入“femaleMiceWeights.csv”，再点击“Raw”,右击“网页另存为”，保存名字“femaleMiceWeights.csv”在本地文件"project1"中

新建R Script文件，命名为“code.R”

在R Script界面中编程：dat <-read.csv("femaleMiceWeights.csv")

计算p-value

dat <-read.csv("femaleMiceWeights.csv")
dat[1:12,2]
mean(dat[13:24,2])-mean(dat[1:12,2])
population <- read.csv("femaleControlsPopulation.csv")

control <- sample(population[,1],12)
mean(control)

n <- 10000
null <- vector("numeric",n)
for(i in 1:n){
control <-sample(population[,1],12)
treatment <- sample(population[,1],12)
null[i] <- mean(treatment)-mean(control)
}
diff <- mean(dat[13:24,2])-mean(dat[1:12,2])
mean(null>diff)

n <- 100
library(rafalib)
mypar(1,1)
plot(0,0,xlim=c(1,30),type="n")
totals <- vector("numeric",11)
for (i in 1:n){
control <- sample(population[,1],12)
treatment <- sample(population[,1],12)
nulldiff <- mean(treatment)-mean(control)
j <- pmax(pmin(round(nulldiff)+6,11),1)
totals[j]<-totals[j]+1
text(j-6,totals[j],pch=15,round(nulldiff,1),cex=0.75)
if(i<15) scan()

}

q-q plot

hist(null)
qqnorm(null)
qqline(null)

pops <- read.csv("mice_pheno.csv")
head(pops)

hf <- pops[pops$Diet=="hf"&pops$Sex=="F",3]
chow <- pops[pops$Diet=="chow"&pops$Sex=="F",3]
mean(hf)-mean(chow)

x <- sample(hf,12)
y <- sample(chow,12)

mean(x)-mean(y)

Ns <- c(3,5,10,25)
B <- 10000
res <- sapply(Ns,function(n){
sapply(1:B,function(j){
mean(sample(hf,n))-mean(sample(chow,n))
})
})

library(rafalib)
for(i in seq(along=Ns)){
title <- paste("Avg=",signif(mean(res[,i]),3))
title <- paste(title,"SD=",signif(sd(res[,i]),3))
qqnorm(res[,i],main=title)
qqline(res[,i])
}

dat<-read.csv("femaleMiceWeights.csv")

dat

control <- dat[1:12,2]

treatment<-dat[12+1:12,2]

diff <- mean(treatment)-mean(control)

print(diff)

t.test(treatment,control)

sd(control)

sd(control)/sqrt(length(control))

se <- sqrt(var(treatment)/length(treatment)+var(control)/length(control))

tstat <- diff/se

1-pnorm(tstat)+pnorm(-tstat)

qqnorm(treatment)

qqline(treatment)

t.test(treatment,control)

dat <- read.csv("mice_pheno.csv")
pop <- dat[dat$Sex=="F" & dat$Diet=="chow",3]
mu <- mean(pop)
mu

N<-30
y <- sample(pop,N)
mean(y)

se <- sd(y)/sqrt(N)
se

Q <- qnorm(1-0.05/2)
interval <- c(mean(y)-Q*se,mean(y)+QQ*se)

plot(mu+c(-7,7),c(1,1),type="n",xlab="weights",
ylab="intervals",ylim=c(1,100))
abline(v=mean(pop))

lines(interval,c(1,1))
for(i in 2:100){
y <- sample(pop,N)
se <- sd(y)/sqrt(N)
interval <- c(mean(y)-Q*se,mean(y)+QQ*se)
color <- ifelse(interval[1]<=mean(pop)&interval[2]>=mean(pop),1,2)
lines(interval,c(i,i),col=color)
}