一、 模拟实验
从相同的正态分布总体N(50,52)中随机抽两个样本，样本含量均为10。对两个样本进行成组t检验，检验水准为0.05，以比较其均数。重复上述步骤1000次，有多少次拒绝H0? 改变检验水准，例如0.10，结果如何？改变样本含量，结果如何？完成表格并下结论总结。
#第一题
ttest<- function(n,alpha,mean1,mean2,sd1,sd2)
{
set.seed(18200135)
count <- 0
for(i in 1:1000)
{
x <- rnorm(n,mean1,sd1)
y <- rnorm(n,mean2,sd2)
if(t.test(x,y,var.equal = T)[[3]]<alpha)
{
count <- count + 1
}
}
return(count)
}
ttest(10,0.05,50,50,5,5)
ttest(20,0.05,50,50,5,5)
ttest(50,0.05,50,50,5,5)
ttest(10,0.10,50,50,5,5)
ttest(20,0.10,50,50,5,5)
ttest(50,0.10,50,50,5,5)
ttest(10,0.20,50,50,5,5)
ttest(20,0.20,50,50,5,5)
ttest(50,0.20,50,50,5,5)
#参考
for(n in c(10,20,50)){
for(alpha in c(0.05,0.1,0.2)){
num=0
for(i in 1:1000){
x = rnorm(n,mean=50,sd=5)
y = rnorm(n,mean=50,sd=5)
test=t.test(x,y,var.equal = T)
if(test$p.value<=alpha)num=num+1
}
print(num)#print(paste(n,alpha,num))
}
}
结论：随着检验水准的增加，拒绝H0的次数增加；随着样本量的增加，拒绝H0的次数会逐渐趋向于1000×检验水准次.
实际数据分析
肺癌是恶性肿瘤中发病率和死亡率均为第一的疾病，是对人群健康和生命威胁最大的恶性肿瘤之一。大量研究表明，肺癌的发生与多种外界环境因素均有着紧密的关联，如主动吸烟、被动吸烟（二手烟）、肺癌家族史等。医学研究中，探索肺癌发生的相关风险因素、预测肺癌风险具有着重要意义。在某地一次回顾性病例对照研究中，纳入肺癌患者500例，健康对照人群500例，收集年龄、性别、吸烟等相关因素，拟调查肺癌的相关风险因素，数据集见test_data.csv。变量编码字典如下：
（1）请绘制年龄age变量分布的柱形图和核密度曲线图（要求在同一个图中绘制）
（2）请编制一个自定义函数（function），要求：能够计算定量资料的四分位数，并以“P50（P25-P75）”的形式输出，如年龄age: 60 (53-65)岁。然后输出体重指数bmi和吸烟强度packyr的结果。
体重指数bmi：23.87511478(22.21646536525.64291739)40)
吸烟强度packyr：25(10
（3）年龄age，教育程度education，体重指数bmi，吸烟smoke，上述哪几个因素在不同性别间差异显著？请使用适当统计学方法评价。
（4）使用适当的广义线性模型评估肺癌风险与基线因素的关联强度，总结哪些因素关联性较强。
setwd("C:/Users/Administrator/Desktop")
datas <- read.csv("test_data.csv",header = T)
colnames(test)
#第一问
a <- density(datas$age)
plot(a)
par(new=TRUE)
barplot(table(datas$age))
#参考
test_data <- read.csv("test_data.csv")
head(test_data)
dim(test_data)
hist(test_data$age)
lines(density(test_data$age),lwd=2,col="red")
hist(test_data$age,freq=F)
lines(density(test_data$age),lwd=2,col="red")
#第二问 请编制一个自定义函数（function），
#要求：能够计算定量资料的四分位数，并以“P50（P25-P75）”的形式输出
#，如年龄age: 60 (53-65)岁。然后输出体重指数bmi和吸烟强度packyr的结果。
quartile<-function(var){
q1<-quantile(var,0.25,na.rm=T)
q2<-quantile(var,0.5,na.rm=T)
q3<-quantile(var,0.75,na.rm=T)
Q<-cat(paste0(q2,"(",q1,"~",q3,")"))
}
quartile(datas$bmi)
quartile(datas$packyr)
#参考
describe=function(x){
q50=median(x,na.rm=T)
q25=quantile(x,0.5,na.rm=T)
q75=quantile(x,0.75,na.rm=T)
stats=paste0(round(q50,2),"(",round(q25,2),"-",round(q75,2),")")
return(stats)
}
describe(test_data$bmi)
#第三问 年龄age，教育程度education，体重指数bmi，
#吸烟smoke，上述哪几个因素在不同性别间差异显著？请使用适当统计学方法评价。
t.test(age~sex,data=datas)
wilcox.test(education~sex,data=datas)
t.test(bmi~sex,data=datas)
wilcox.test(smoke~sex,data=datas)
#参考
t.test(test_data$age~test_data$sex)#wilcox.test(test_data$age~test_data$sex)
chisq.test(table(test_data$education,test_data$sex))#fisher.test(table(test_data$education,test_data$sex))
t.test(test_data$bmi~test_data$sex)#wilcox.test(test_data$bmi~test_data$sex)
chisq.test(table(test_data$smoke,test_data$sex))#fisher.test(table(test_data$smoke,test_data$sex))
#第四问
#univariate logistic
logistic<-function(x){
fml<-as.formula(paste0("lung_ca~",x))
glog<-glm(fml,family = binomial(link = "logit"),data=datas)
gsum<-summary(glog)
OR<-round(exp(gsum$coefficients[2,1]),3)
Z<-round(gsum$coefficients[2,3],3)
P<-round(gsum$coefficients[2,4],3)
L<-round(exp(confint(glog)[2,1]),3)
U<-round(exp(confint(glog)[2,2]),3)
CI<-paste0(L,sep='-',U)
logresult<-data.frame('variable'=x,
'OR'=OR,
'CI'=CI,
'Z'=Z,
'P'=P)
return(logresult)
}
var<-colnames(datas)[3:14]
univar<-lapply(var,logistic)
#参考
test_data$packyr[is.na(test_data$packyr)]=0
test_data$bmi[is.na(test_data$bmi)]=mean(test_data$bmi,na.rm=T)
test_data1=test_data[,-1]
model=glm(lung_ca ~ .,data=test_data1,family=binomial)
summary(model)
step(model,direction = "backward")
apply(test_data1,2,function(x) sum(is.na(x)))
model=glm(lung_ca ~ age+sex+packyr+respdis+exposure+drink,data=test_data1,family=binomial)