swtj ryy ks?

[color=#008000][b]      17级[/b][u][/u]    
张yu答案[/color]
一、    mnsy（40分）
从正态分布N(0,σ)中随机抽取样本，根据样本估计总体均数的95%可信区间，重复1000 次，有多少次能估计到总体0，均数。不同的样本含量N对结果有什么影响，不同的标准差σ对结果有什么影响，不同的可信度又会对结果有什么影响？（1）进行模拟实验，填入下面表格（30分）；
（2）总结规律（10分）。
1000次模拟中估计到总体均数的次数
模拟参数    标准差σ=1    标准差σ=10
    N=10    N=20    N=50    N=10    N=20    N=50
可信度=90%    367    296    182    405    289    191
可信度=95%    436    331    201    446    315    210
可信度=99%    563    392    270    561    422    292

总结结论如下：
当样本量和标准差不变时，随着可信度的增加，落入可信区间内的次数增多；当样本量和可信度不变时，随着标准差的增加落入可信区间内次数基本不变；当可信度和标准差不变时，随着样本量的增加落入可信区间内的次数减少。

答：set.seed(17200115)
ttest<-function(n1,mean1,sd1,alpha)
{
  count<-0
  for(i in 1:1000)
  {
    d1<-rnorm(n1,mean=mean1,sd=sd1)
    xb<-mean(d1)
   
    tmp<-sd1/sqrt(n1)*qnorm(1-alpha/2);df<-n1
   
    a<-xb-tmp
    b<-xb+tmp
   
   
    if(a>mean1)
    {
      count<-count+1
    }
    if(b<mean1)
    {
      count<-count+1
    }    
  }
  return(count)
}

ttest(10,0,1,0.9)
ttest(20,0,1,0.9)
ttest(50,0,1,0.9)
ttest(10,0,1,0.95)
ttest(20,0,1,0.95)
ttest(50,0,1,0.95)
ttest(10,0,1,0.99)
ttest(20,0,1,0.99)
ttest(50,0,1,0.99)

ttest(10,0,10,0.9)
ttest(20,0,10,0.9)
ttest(50,0,10,0.9)
ttest(10,0,10,0.95)
ttest(20,0,10,0.95)
ttest(50,0,10,0.95)
ttest(10,0,10,0.99)
ttest(20,0,10,0.99)
ttest(50,0,10,0.99)

二、数据分析（60分）
请载入R语言自带的数据集iris，该数据集描述了三种不同的鸢尾属植物对应的萼片和花瓣的长度和宽度。其中，Species变量代表植物种类，Sepal代表萼片，Petal代表花瓣，Length为长度，Width为宽度。请按照下述要求进行分析。
（1）    根据不同的Species种类，分别使用均数±标准差描述不同种类植物Sepal.Length和Sepal.Width两个属性，使用中位数（四分位数间距）描述不同种类植物Petal.Length和Petal.Width两个属性（10分）。

（2）计算植物属性的4个变量间的Pearson相关系数矩阵和相关系数假设检验的P值矩阵，无需分组（10分）。
根据数据所得Pearson相关系数矩阵如下

（3）请自编变量转换函数，要求实现如下功能：将定量数据根据固定界值划分为定性数据， [0~2]定义为1，(2~4]定义为2，(4~6]定义为3, (6~8]定义为4。将编写的函数填写在试卷上并对Petal.Length变量做变换，对变换后的变量做频数统计（10分）。

（4）绘制Sepal.Length在不同Species种类下的分布的箱式图和散点图，要求：1）在一张图中绘制，散点图直接绘制在箱式图上；
2）箱式图不同种类填充色不同，并添加图例说明不同颜色与Species种类对应关系；
3）散点为蓝色，且排列方式为随机扰动，而不是一串排列（15分）。
结果如下

（5）使用三种不同的分类算法，基于植物属性的4个变量数据将其分为3类，并展示分类结果与真实标签比较，使用准确率指标进行评价。分类算法不限（15分）。

答：##第二大题##

##1##
iris <- datasets::iris
iris1<-iris[which(iris$Species=='setosa'),]
iris2<-iris[which(iris$Species=='versicolor'),]
iris3<-iris[which(iris$Species=='virginica'),]

paste(round(mean(iris1$Sepal.Length),2),"±",round(sd(iris1$Sepal.Length),2))

paste(round(mean(iris2$Sepal.Length),2),"±",round(sd(iris2$Sepal.Length),2))

paste(round(mean(iris3$Sepal.Length),2),"±",round(sd(iris3$Sepal.Length),2))

paste(round(mean(iris1$Sepal.Width),2),"±",round(sd(iris1$Sepal.Width),2))

paste(round(mean(iris2$Sepal.Width),2),"±",round(sd(iris2$Sepal.Width),2))

paste(round(mean(iris3$Sepal.Width),2),"±",round(sd(iris3$Sepal.Width),2))

a<-quantile(iris1$Sepal.Length,0.25)
b<-quantile(iris1$Sepal.Length,0.75)
paste(median(iris1$Sepal.Length),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris2$Sepal.Length,0.25)
b<-quantile(iris2$Sepal.Length,0.75)
paste(median(iris2$Sepal.Length),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris3$Sepal.Length,0.25)
b<-quantile(iris3$Sepal.Length,0.75)
paste(median(iris3$Sepal.Length),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris1$Sepal.Width,0.25)
b<-quantile(iris1$Sepal.Width,0.75)
paste(median(iris1$Sepal.Width),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris2$Sepal.Width,0.25)
b<-quantile(iris2$Sepal.Width,0.75)
paste(median(iris2$Sepal.Width),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris3$Sepal.Width,0.25)
b<-quantile(iris3$Sepal.Width,0.75)
paste(median(iris3$Sepal.Width),"(",a[[1]],"-",b[[1]],")")

##2##
install.packages("psych")
library(psych)

data1<-iris[,c(1,2,3,4)]
data1
a<-cor(data1)
corr.test(data1)
 

##3##

zhang<-function(a){
    for(i in 1:length(a)){
        if(a[i]>=0 & a[i]<=2){
            a[i]=1
        }else if(a[i]>2 & a[i]<=4){
            a[i]=2
        }else if(a[i]>4 & a[i]<=6){
            a[i]=3
        }else{
            a[i]=4
        }
    }
    return(a)
}

zhang(iris$Petal.Length)

freq<-table(zhang(iris$Petal.Length))
freq

##4##
boxplot(iris$Sepal.Length~iris$Species,col=c(2,3,5),xlab="x",ylab="y")
par(new=TRUE)
plot(iris$Sepal.Length,col="blue",xlab="x",ylab="y")

legend("topright",                                    #图例位置为右上角
legend=c("setosa","versicolor","virginica"),          #图例内容
col=c(2,3,5),                                 #图例颜色
lty=1,lwd=2)                                          #图例大小

##5##

#加载数据
iris <- datasets::iris
iris2 <- iris[,-5]
species_labels <- iris[,5]
library(colorspace) # 颜色包
species_col <- rev(rainbow_hcl(3))[as.numeric(species_labels)]
#绘制 SPLOM:
pairs(iris2, col = species_col,
      lower.panel = NULL,
       cex.labels=2, pch=19, cex = 1.2)
 
# 添加图例
par(xpd = TRUE)
legend(x = 0.05, y = 0.4, cex = 2,
   legend = as.character(levels(species_labels)),
    fill = unique(species_col))
par(xpd = NA)

par(las = 1, mar = c(4.5, 3, 3, 2) + 0.1, cex = .8)
MASS::parcoord(iris2, col = species_col, var.label = TRUE, lwd = 2)
 
# 添加标题
title("Parallel coordinates plot of the Iris data")
# 添加图例
par(xpd = TRUE)
legend(x = 1.75, y = -.25, cex = 1,
   legend = as.character(levels(species_labels)),
    fill = unique(species_col), horiz = TRUE)

data1<-iris[,c(1,2,3,4)]
data1
kmeans(data1)

[color=#99CC00]孙chuanrui答案[/color]
[color=#008000]
#模拟试验[/color]

set.seed(17200130)
ttest <- function(num,n,alpha,mean,sd){
  count <- 0
  for (i in 1:num) {
    x1 <- rnorm(n,mean,sd)
    x2 <- rnorm(n,mean,sd)
    test <- t.test(x1,x2,var.equal = T)
    p <- test$p.value
    if( p< alpha ){count <- count +1}
  }
  
  return(count)
  
}
#第一行
for (alpha in c(0.05,0.10)) {
  for (n in c(10,15,20)) {
    print(paste("alpha = ",alpha,"N = ",n," ",ttest(1000,n,alpha,50,5)))
  }#i
}#j
#第二行
for (alpha in c(0.05,0.10)) {
  for (n in c(10,15,20)) {
    print(paste("alpha = ",alpha,"N = ",n," ",ttest(1000,n,alpha,50,8)))
  }#i
}#j
#第三行
for (alpha in c(0.05,0.10)) {
  for (n in c(10,15,20)) {
    print(paste("alpha = ",alpha,"N = ",n," ",ttest(1000,n,alpha,70,10)))
  }#i
}#j

[color=#008000]
二、[/color]
#实际数据分析
data2 <- read.csv("C:\\Users\\阿瑞\\Desktop\\test_data.csv")
#(1)绘制年龄分布的柱形图和核密度曲线

par(mfcol = c(2,1))
plot(as.factor(data2$age),main = "年龄分布柱形图",xlab = "年龄（岁）",ylab = "人数",col = "#99CCCC")

plot(density(data2$age),main = "年龄分布核密度曲线图",xlab = "年龄（岁）",col = "#003366")
#(2)自定义函数计算定量资料的四分位数
fun <- function(n,data2){
  p50 <- round(quantile(data2[,c(n)],na.rm = T)[3],2)
  p25 <- round(quantile(data2[,c(n)],na.rm = T)[2],2)
  p75 <- round(quantile(data2[,c(n)],na.rm = T)[4],2)
  p <- cbind(p50,p25,p75)
  
  per <- paste0(p[1],"(",p[2],"-",p[3],")")
  name <- colnames(data2)[n]
  x <- data.frame(name,per)
  return(x)
}

q2 <- rbind(fun(3,data2) ,fun(6,data2) ,fun(9,data2) )
write.csv(q2,"C:\\Users\\阿瑞\\Desktop\\q2.csv")
#q3评价不同性别之间差异是否显著
library(tableone)
data2$sex <- as.factor(data2$sex)

data2$education<-as.factor(data2$education)
data2$smoke<-as.factor(data2$smoke)
data2$drink<-as.factor(data2$drink)
data2$exercise<-as.factor(data2$exercise)
vars <- c("age","education","bmi","smoke","drink","exercise")
table1 <- CreateTableOne(vars = vars,strata = c("sex") ,smd = T ,addOverall = T,data = data2)
table2 <- print(table1,nonnormal = c(),showAllLevels = T)
write.csv(table2,file = "C:\\Users\\阿瑞\\Desktop\\q3.csv")

#q4使用适当的广义线性模型评估肺癌风险与基线因素的关联强度，总结哪些因素关联性较强。

#使用单因素logistic回归，评价肺癌风险与基因因素之间的关联强度
gen_logistic <- function(var){
  in_formula <- as.formula(paste("lung_ca~",var)) 
  p <- glm(in_formula,family=binomial(link=logit),data=data2)
  beta <- summary(p)$coefficients[2,1]
  LCI <- summary(p)$coefficients[2,1] -summary(p)$coefficients[2,2]*1.96 
  UCI <- summary(p)$coefficients[2,1] +summary(p)$coefficients[2,2]*1.96 
  OR <- exp(beta) 
  OR_LCI <- exp(LCI)
  OR_UCI <- exp(UCI)
  p_value <- summary(p)$coefficients[2,4]
  name <- var
  data_var <- data.frame(name,OR,OR_LCI,OR_UCI,p_value)
  return(data_var)
}
#对年龄进行分层
data2$ageu[data2$age >= 40 &data2$age < 50] <- 1
data2$ageu[data2$age >= 50 &data2$age < 60] <- 2
data2$ageu[data2$age >= 60 &data2$age < 70] <- 3
data2$ageu[data2$age >= 70 &data2$age < 80] <- 4
#对bmi进行分层
data2$bmiu[data2$bmi <18] <- 1
data2$bmiu[data2$bmi >=18 &data2$bmi <24] <- 2
data2$bmiu[data2$bmi >=24] <- 3

q4 <- rbind(gen_logistic("ageu"),gen_logistic("sex"),gen_logistic("education"),
      gen_logistic("bmiu"),gen_logistic("family_ca"),gen_logistic("smoke"),
      gen_logistic("packyr"),gen_logistic("respdis"),gen_logistic("secsmoke"),
      gen_logistic("exposure"),gen_logistic("drink"),gen_logistic("exercise"))  
write.csv(q4,file = "C:\\Users\\阿瑞\\Desktop\\q4.csv")
#多因素分析
mul <- glm(lung_ca ~ageu+family_ca+sex+bmiu+smoke+education+respdis,data = data2)
m <- summary(mul)$coefficients[-1,]
beta <- m[,1]
LCI <- beta -m[,2]*1.96 
UCI <- beta +m[,2]*1.96 
OR <- exp(beta) 
OR_LCI <- exp(LCI)
OR_UCI <- exp(UCI)
p <- m[,4]
model1 <- rbind(beta,LCI,UCI,OR,OR_LCI,OR_UCI,p)
model1 <- t(model1)
write.csv(model1,file = "C:\\Users\\阿瑞\\Desktop\\q4_1.csv")

[color=#008000]

王yuanyuan[/color]

simu1<-function(alpha,sigma,n){
  count<-0
  for (i in 1:1000) {
    sample<-rnorm(10000,0,sigma)
    data<-sample(sample,n)
    t<-t.test(data,conf.level =alpha)
    low<-t$conf.int[1]
    up<-t$conf.int[2]
    if(low<0 & up>0) {
     count<-count+1
    }
  }
  return(count)
}

simu1(0.9,1,10)
simu1(0.9,1,20)
simu1(0.9,1,50)

simu1(0.9,10,10)
simu1(0.9,10,20)
simu1(0.9,10,50)

simu1(0.95,1,10)
simu1(0.95,1,20)
simu1(0.95,1,50)

simu1(0.95,10,10)
simu1(0.95,10,20)
simu1(0.95,10,50)

simu1(0.99,1,10)
simu1(0.99,1,20)
simu1(0.99,1,50)

simu1(0.99,10,10)
simu1(0.99,10,20)
simu1(0.99,10,50)

第二题见别tie

###(1)

attach(iris)

describe<-function(x,parametric=TRUE,print=TRUE){

   a<-which(colnames(iris) == x)

   if(parametric){

    df_mean<-aggregate(iris[,a],by=list(Species),FUN=mean)

    df_sd<-aggregate(iris[,a],by=list(Species),FUN=sd)

  }else{

    median<-aggregate(iris[,a],by=list(Species),FUN=median)

    df_q1<-aggregate(iris[,a],by=list(Species),FUN=quantile,probs=0.25)

    df_q3<-aggregate(iris[,a],by=list(Species),FUN=quantile,probs=0.75)

  }

  if(parametric&print){

    cat(paste0(mean,"±",sd))

  }else if(!parametric&print){

    cat(paste0(median,"(",q1,",",q3,")"))

  }

   return()

}

describe(Sepal.Length)

describe(Sepal.Width)

if(1){

  mean

  a1<-cat(paste0(mean[1,2],"±",sd[1,2]))

  a2<-cat(paste0(mean[2,2],"±",sd[2,2]))

  a3<-cat(paste0(mean[3,2],"±",sd[3,2]))

  res<-rbind(a1,a2,a3)

}else if(!parametric&print){

  cat(paste0(median,"(",q1,",",q3,")"))

}

###(2)

install.packages("psych")

library(psych)

data<-iris[,-5]

cor(data)

corr.test(data)

####(3)

trans<-function(var){

  iris$cat[0<=var&var<=2]<-1

  iris$cat[2<var&var<=4]<-2

  iris$cat[4<var&var<=6]<-3

  iris$cat[6<var&var<=8]<-4

  iris1<<-iris

  return(iris1)

}

trans(Petal.Length)

table(iris1$cat)

####(4)

install.packages("ggplot2")

library(ggplot2)

ggplot(data=iris,aes(x=Species,y=Sepal.Length,colour=Species))+geom_boxplot()+geom_jitter(colour="blue")

####(5)

##第一种：层次聚类

d<-dist(data)

fit.average<-hclust(d,method = "average")

plot(fit.average)

clusters<-cutree(fit.average,k=3)

data2<-cbind(iris,clusters)

data2

table(clusters)

##2.逻辑回归

train<-sample(1:nrow(iris),0.7*nrow(iris))

df.train<-iris[train,]

df.validate<-iris[-train,]

fit.logit<-glm(Species~.,data=df.train,family = binomial())

logit.pred<-predict(fit.logit,df.validate)

logit.perf<-table(df.validate$Species,logit.pred,dnn=c("Actual","Predicted"))

##3.随机森林

install.packages("randomForest")

library(randomForest)

train<-sample(1:nrow(iris),0.7*nrow(iris))

df.train<-iris[train,]

df.validate<-iris[-train,]

fit.forest<-randomForest(Species~.,data=df.train,importance=TRUE)

forest.pred<-predict(fit.forest,df.validate)

forest.perf<-table(df.validate$Species,forest.pred,dnn=c("Actual","Predicted"))

forest.perf

###决策树

install.packages("party")

library(party) #导入party包

myctree<-ctree(df.train$Species~.,df.train)

pre.forest<-predict(myctree,df.validate)

jueceshu<-table(pre.forest,df.validate$Species)

###准确率

##决策树

acc<-(jueceshu[1,1]+jueceshu[2,2])/(jueceshu[1,1]+jueceshu[2,2]+jueceshu[1,2]+jueceshu[2,1])

• scr 期中
• 模拟实验

从相同的正态分布总体N(50,52)中随机抽两个样本，样本含量均为10。对两个样本进行t检验，检验水准为0.05，以比较其均数。重复上述步骤1000次，有多少次拒绝H0? 改变检验水准，例如0.10，结果如何？改变样本含量，结果如何？改变均数和标准差结果如何？完成表格并下结论总结。

结论：

（1）当均数，方差和样本量不变时，检验水准增加，1000次模拟试验中拒绝的次数增加；

（2）当均数，方差和检验水准不变时，样本量增加，1000次模拟试验中拒绝的次数基本不变；

（2）当检验水准不变时， 1000次模拟试验中拒绝的次数与正态分布的均数、标准差和样本量无关；

set.seed(17200122)
Ttest<-function(n,mean,sd,alpha){
  count<-0
  for(i in 1:1000){   #i只是一个指针 后面不一定要出现
    x1<-rnorm(n,mean,sd)
    x2<-rnorm(n,mean,sd)
    ttest<-t.test(x1,x2,var.equal=TRUE)
    P<-ttest$p.value
    if(P<alpha) {count<-count+1}
  }
  return(count)
}
#mean=50 sd=5 alpha=0.05 n=c(10,15,20)

for(alpha in c(0.05,0.1)){
  for(n in c(10,15,20)){
    print(Ttest(n,50,5,alpha))
  }
}
for(alpha in c(0.05,0.1)){
  for(n in c(10,15,20)){
    print(Ttest(n,50,8,alpha))
  }
}
for(alpha in c(0.05,0.1)){
  for(n in c(10,15,20)){
    print(Ttest(n,70,10,alpha))
  }
}

• wyy期末
• 模拟实验（40分）

从正态分布N(0,σ)中随机抽取样本，根据样本估计总体均数的95%可信区间，重复1000 次，有多少次能估计到总体均数。不同的样本含量N对结果有什么影响，不同的标准差σ对结果有什么影响，不同的可信度又会对结果有什么影响？（1）进行模拟实验，填入下面表格（30分）；

（2）总结规律（10分）。

结论：有结果表格可知，从正态分布N(0,σ)中随机抽取样本，根据样本估计总体均数的95%可信区间，重复1000 次

1. 可信区间估计到总体参数的次数不会随着样本含量以及标准差的增加而增加，在不同的标准差以及样本量下，可信区间估计到总体参数的次数大致不变，都在（可信度*1000）次左右。
2. 可信度的改变会影响可信区间估计到总体参数的次数。可信度越大，可信区间估计到总体参数的次数也就随之增大，可信区间估计到总体参数的次数都在（可信度*1000）次左右。

set.seed(17200122)
Test<-function(n,sd,Z){
  count<-0
  for(i in 1:1000){
    x<-rnorm(n,0,sd) #mean +- Z*se  sd/根号n sqrt CI1 CI2 1.64 1.96 2.58
    Mean<-mean(x)
    Sd<-sd(x)
    Se<-Sd/sqrt(n)
    CI1<-Mean-Z*Se
    CI2<-Mean+Z*Se
    if(CI1<0&CI2>0){count<-count+1}
  }
return(count)
}
Test(10,1,1.96)

for(sd in c(1,10)){
  for(n in c(10,20,50)){
    print(Test(n,sd,1.64))
  }
}
for(sd in c(1,10)){
  for(n in c(10,20,50)){
    print(Test(n,sd,1.96))
  }
}
for(sd in c(1,10)){
  for(n in c(10,20,50)){
    print(Test(n,sd,2.58))
  }
}

期中 scr

#17200130
#数据录入
data2 <- read.csv("C:\\Users\\闃跨憺\\Desktop\\test_data.csv")
data2 <- read.csv("C:/Users/lenovo/Desktop/test_data.csv")
#(1)请绘制年龄age变量分布的柱形图和核密度曲线图（要求在同一个图中绘制）

par(mfcol = c(2,1)) #col 列 row 行 一页多图 一页两幅图 mfrow=c(2,2) 4张图 两行两列
plot(as.factor(data2$age),main = "年龄分布柱形图",xlab = "年龄",ylab = "人数",col = "#99CCCC")
#col颜色类别 xlab ylab 横纵坐标轴标题 main 图表标题
plot(density(data2$age),main = "年龄分布核密度曲线图",xlab = "年龄",col = "#003366")

#(2)请编制一个自定义函数（function），
#要求：能够计算定量资料的四分位数，并以“P50（P25-P75）”的形式输出
#输出体重指数bmi和吸烟强度packyr以及年龄age的结果 age 第三列 bmi 第六列 吸烟强度 第九列

quantile(data2,na.rm = T) 
#na.rm = T 允许有遗漏值和NaN 输出为 0% 25% 50% 75% 100% 对应quantile()[2] [3] [4] 一列数中的第二位 第三位 第四位

fun <- function(n,data2){
  p50 <- round(quantile(data2[,c(n)],na.rm = T)[3],2) 
  #round(timedata,2) 保留两位小数 quantile 分位数函数 data2[,c(n)] c()赋值 c(n) data2第n列
  p25 <- round(quantile(data2[,c(n)],na.rm = T)[2],2) 
  p75 <- round(quantile(data2[,c(n)],na.rm = T)[4],2)
  #三个数据由纵向变成横向 再按照格式组合起来  再写一个列名  再转换为数据框  返回数据框
  p <- cbind(p50,p25,p75)#根据列进行合并 横向  rbind 根据行进行合并 纵向
  per <- paste0(p[1],"(",p[2],"-",p[3],")") #任意数量的参数组合起来 paste0 以空字符串连接字符
  name <- colnames(data2)[n]
  #colnames() and rownames()函数可以修改矩阵或数据框的行名和列名，但不可修改矩阵的变量名。
  #colnames(data)[n]="xxx"修改data第n列列名为xxx
  x <- data.frame(name,per)#转换为数据框
  return(x)
}
#rbind 根据行进行合并 纵向
q2 <- rbind(fun(3,data2) ,fun(6,data2) ,fun(9,data2) )
write.csv(q2,"C:/Users/lenovo/Desktop/q2.csv") #导出数据
#三线图 居中 无边框 加入上下框线宽度为1.5 最后一条线是标题行下面的线 
#（鼠标选中标题行的下面一行，进入边框和底纹的面板，设置线宽为0.5，右边预览中选择上框线）

#（3）年龄age，教育程度education，体重指数bmi，吸烟smoke，饮酒dink，锻炼exercise，
#上述哪几个因素在不同性别!间差异显著？请使用适当统计学方法评价。
install.packages("tableone")
library(tableone) #TableOne是生成统计表的工具，常用于生成论文中的表格
data2$sex <- as.factor(data2$sex)#data2中的sex变量提取

data2$education<-as.factor(data2$education)
data2$smoke<-as.factor(data2$smoke)
data2$drink<-as.factor(data2$drink)
data2$exercise<-as.factor(data2$exercise)
vars <- c("age","education","bmi","smoke","drink","exercise") #赋值给vars 向量
table1 <- CreateTableOne(vars = vars,strata = c("sex") ,smd = T ,addOverall = T,data = data2)
#以不同组(sex变量)分层创建Table 1
#创建table one的函数非常简单，CreatTableOne()函数只需要指出需描述的变量（即前面的vars变量列表），strata参数说明按照trt变量分层即可。
#注意的是如果前面没有指定分类变量类型，tableone会以数值型变量处理你的变量，这也就是需要提前指定分类变量的原因。
table2 <- print(table1,nonnormal = c(),showAllLevels = T)
#指定变量分析方法
#默认情况下，tableone使用正态分布方法分析资料，因此会出现“(mean (sd))”的描述，
#两组计量资料用oneway.test函数，t检验分析，分类资料用 chisq.test函数，卡方检验分析，默认有矫正卡方，精确性检验用fisher.test函数Fisher检验等等。
#我们可以指定变量是进行正态或非正态检验方法。
#summary(table2)
write.csv(table2,file = "C:/Users/lenovo/Desktop/q3.csv")
#特征    Overall    男    女    P值 一一对应 level    Overall    1    2    p    n是合计 三线表
#总结：对年龄和BMI使用t检验比较在不同性别之间的差异性，结果显示，年龄在不同性别之间的差异具有统计学意义；
#对教育程度，吸烟，饮酒和锻炼使用列联表的卡方检验，比较在不同性别之间的差异，
#结果显示，教育程度，吸烟，饮酒的不同性别之间的差异具有统计学意义。

data2
#（4）使用适当的广义线性模型评估肺癌风险与基线因素的关联强度，总结哪些因素关联性较强
#先进行单因素分析 后进行多因素分析 肺癌lung_ca 与12个变量
#挨个进行单因素分析 进行循环
gen_logistic <- function(var){
  in_formula <- as.formula(paste("lung_ca~",var)) #formula拟合公式 family 指定分布族
  p <- glm(in_formula,family=binomial(link=logit),data=data2)#正态分布gaussian(link=identity) 二项分布 binomial(link=logit)
  beta <- summary(p)$coefficients[2,1]
  LCI <- summary(p)$coefficients[2,1] -summary(p)$coefficients[2,2]*1.96 
  UCI <- summary(p)$coefficients[2,1] +summary(p)$coefficients[2,2]*1.96 
  OR <- exp(beta) 
  OR_LCI <- exp(LCI)
  OR_UCI <- exp(UCI)
  p_value <- summary(p)$coefficients[2,4]
  #P判断是否具有统计学意义 OR 首先在可信区间内 其次与1比较危险因素还是保护因素
  name <- var
  data_var <- data.frame(name,OR,OR_LCI,OR_UCI,p_value)
  return(data_var)
}
gen_logistic("age")
#？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？
data2$ageu[data2$age >= 40 &data2$age < 50] <- 1
data2$ageu[data2$age >= 50 &data2$age < 60] <- 2
data2$ageu[data2$age >= 60 &data2$age < 70] <- 3
data2$ageu[data2$age >= 70 &data2$age < 80] <- 4

data2$bmiu[data2$bmi <18] <- 1
data2$bmiu[data2$bmi >=18 &data2$bmi <24] <- 2
data2$bmiu[data2$bmi >=24] <- 3

q4 <- rbind(gen_logistic("ageu"),gen_logistic("sex"),gen_logistic("education"),
      gen_logistic("bmiu"),gen_logistic("family_ca"),gen_logistic("smoke"),
      gen_logistic("packyr"),gen_logistic("respdis"),gen_logistic("secsmoke"),
      gen_logistic("exposure"),gen_logistic("drink"),gen_logistic("exercise"))  
write.csv(q4,file = "C:\\Users\\闃跨憺\\Desktop\\q4.csv")
#多因素分析  单因素分析有意义的变量#有交互作用的怎么办？？？？？？？？？？？？？？？？？？？？？？
mul <- glm(lung_ca ~ageu+family_ca+sex+bmiu+smoke+education+respdis,data = data2)
m <- summary(mul)$coefficients[-1,]#？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？
beta <- m[,1] #截距
LCI <- beta -m[,2]*1.96 
UCI <- beta +m[,2]*1.96 
OR <- exp(beta) 
OR_LCI <- exp(LCI)
OR_UCI <- exp(UCI)
p <- m[,4]
model1 <- rbind(beta,LCI,UCI,OR,OR_LCI,OR_UCI,p)
model1 <- t(model1)#转置
write.csv(model1,file = "C:\\Users\\闃跨憺\\Desktop\\q4_1.csv")
 

期末 zy

##第一大题##

set.seed(17200115)
ttest<-function(n1,mean1,sd1,alpha)
{
  count<-0
  for(i in 1:1000)
  {
    d1<-rnorm(n1,mean=mean1,sd=sd1)
    xb<-mean(d1)
   
    tmp<-sd1/sqrt(n1)*qnorm(1-alpha/2);df<-n1
   
    a<-xb-tmp
    b<-xb+tmp
   
   
    if(a>mean1)
    {
      count<-count+1
    }
    if(b<mean1)
    {
      count<-count+1
    }    
  }
  return(count)
}

ttest(10,0,1,0.9)
ttest(20,0,1,0.9)
ttest(50,0,1,0.9)
ttest(10,0,1,0.95)
ttest(20,0,1,0.95)
ttest(50,0,1,0.95)
ttest(10,0,1,0.99)
ttest(20,0,1,0.99)
ttest(50,0,1,0.99)

ttest(10,0,10,0.9)
ttest(20,0,10,0.9)
ttest(50,0,10,0.9)
ttest(10,0,10,0.95)
ttest(20,0,10,0.95)
ttest(50,0,10,0.95)
ttest(10,0,10,0.99)
ttest(20,0,10,0.99)
ttest(50,0,10,0.99)

##第二大题##

##1##
#根据不同的Species种类，分别使用均数±标准差描述不同种类植物Sepal.Length和Sepal.Width两个属性，
#使用中位数（四分位数间距）描述不同种类植物Petal.Length和Petal.Width两个属性
iris <- datasets::iris
iris1<-iris[which(iris$Species=='setosa'),]
iris2<-iris[which(iris$Species=='versicolor'),]
iris3<-iris[which(iris$Species=='virginica'),]

paste(round(mean(iris1$Sepal.Length),2),"±",round(sd(iris1$Sepal.Length),2))

paste(round(mean(iris2$Sepal.Length),2),"±",round(sd(iris2$Sepal.Length),2))

paste(round(mean(iris3$Sepal.Length),2),"±",round(sd(iris3$Sepal.Length),2))

paste(round(mean(iris1$Sepal.Width),2),"±",round(sd(iris1$Sepal.Width),2))

paste(round(mean(iris2$Sepal.Width),2),"±",round(sd(iris2$Sepal.Width),2))

paste(round(mean(iris3$Sepal.Width),2),"±",round(sd(iris3$Sepal.Width),2))

a<-quantile(iris1$Sepal.Length,0.25)
b<-quantile(iris1$Sepal.Length,0.75)
paste(median(iris1$Sepal.Length),"(",a[[1]],"-",b[[1]],")")#？？？？？？？？？？？？？？？？？

a<-quantile(iris2$Sepal.Length,0.25)
b<-quantile(iris2$Sepal.Length,0.75)
paste(median(iris2$Sepal.Length),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris3$Sepal.Length,0.25)
b<-quantile(iris3$Sepal.Length,0.75)
paste(median(iris3$Sepal.Length),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris1$Sepal.Width,0.25)
b<-quantile(iris1$Sepal.Width,0.75)
paste(median(iris1$Sepal.Width),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris2$Sepal.Width,0.25)
b<-quantile(iris2$Sepal.Width,0.75)
paste(median(iris2$Sepal.Width),"(",a[[1]],"-",b[[1]],")")

a<-quantile(iris3$Sepal.Width,0.25)
b<-quantile(iris3$Sepal.Width,0.75)
paste(median(iris3$Sepal.Width),"(",a[[1]],"-",b[[1]],")")

##2##计算植物属性的4个变量间的Pearson相关系数矩阵和相关系数假设检验的P值矩阵，无需分组
install.packages("psych")
library(psych)

data1<-iris[,c(1,2,3,4)]
data1
a<-cor(data1)
corr.test(data1)
 

##3##请自编变量转换函数，要求实现如下功能：
#将定量数据根据固定界值划分为定性数据， [0~2]定义为1，(2~4]定义为2，(4~6]定义为3, (6~8]定义为4。
#将编写的函数填写在试卷上并对Petal.Length变量做变换，对变换后的变量做频数统计

zhang<-function(a){
    for(i in 1:length(a)){ #i作为指针进行循环  下面也用到了
        if(a[i]>=0 & a[i]<=2){
            a[i]=1
        }else if(a[i]>2 & a[i]<=4){
            a[i]=2
        }else if(a[i]>4 & a[i]<=6){
            a[i]=3
        }else{
            a[i]=4
        }
    }
    return(a)
}

zhang(iris$Petal.Length)

freq<-table(zhang(iris$Petal.Length))#table 统计频数
freq

##4##绘制Sepal.Length在不同Species种类下的分布的箱式图和散点图，
#要求：1）在一张图中绘制，散点图直接绘制在箱式图上；
#2）箱式图不同种类填充色不同，并添加图例说明不同颜色与Species种类对应关系；
#3）散点为蓝色，且排列方式为随机扰动，而不是一串排列 ？？？？？？？？？？？？？？？？？如果不是随机扰动呢？？？？？
boxplot(iris$Sepal.Length~iris$Species,col=c(2,3,5),xlab="x",ylab="y")
par(new=TRUE)
plot(iris$Sepal.Length,col="blue",xlab="x",ylab="y")

legend("topright",                                    #图例位置为右上角
legend=c("setosa","versicolor","virginica"),          #图例内容
col=c(2,3,5),                                 #图例颜色
lty=1,lwd=2)                                          #图例大小

##5##使用三种不同的分类算法，？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？
#基于植物属性的4个变量数据将其分为3类，
#并展示分类结果与真实标签比较，使用准确率指标进行评价。分类算法不限

#加载数据
iris <- datasets::iris
iris2 <- iris[,-5]#？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？
species_labels <- iris[,5]
library(colorspace) # 颜色包
species_col <- rev(rainbow_hcl(3))[as.numeric(species_labels)]
#绘制 SPLOM:
pairs(iris2, col = species_col,
      lower.panel = NULL,
       cex.labels=2, pch=19, cex = 1.2)
 
# 添加图例
par(xpd = TRUE)
legend(x = 0.05, y = 0.4, cex = 2,
   legend = as.character(levels(species_labels)),
    fill = unique(species_col))
par(xpd = NA)

par(las = 1, mar = c(4.5, 3, 3, 2) + 0.1, cex = .8)
MASS::parcoord(iris2, col = species_col, var.label = TRUE, lwd = 2)
 
# 添加标题
title("Parallel coordinates plot of the Iris data")
# 添加图例
par(xpd = TRUE)
legend(x = 1.75, y = -.25, cex = 1,
   legend = as.character(levels(species_labels)),
    fill = unique(species_col), horiz = TRUE)

data1<-iris[,c(1,2,3,4)]
data1
kmeans(data1)
 

1、单因素分析

变量

OR

OR_LCI

OR_UCI

p_value

ageu

2.31

1.95

2.73

<0.01

sex

0.60

0.47

0.77

<0.01

education

0.73

0.53

1.00

0.0510

bmiu

0.96

0.76

1.22

0.7377

family_ca

1.43

1.00

2.05

0.0483

smoke

1.90

1.44

2.51

<0.01

packyr

1.02

1.01

1.04

<0.01

respdis

1.77

1.29

2.42

<0.01

secsmoke

1.33

1.03

1.72

0.0302

exposure

0.97

0.69

1.37

0.8616

drink

1.22

0.92

1.62

0.1705

exercise

1.08

0.84

1.38

0.5657

单因素分析结果显示，年龄，肺癌家族史，吸烟，呼吸系统疾病史，被动吸烟是肺癌发生的危险因素，与肺癌之间的关联有统计学意义；性别是肺癌的保护因素，女性患肺癌的风险低于男性，教育程度是肺癌的保护因素，教育程度越高，肺癌发生的风险越低。

2、多因素分析

考虑到吸烟和饮酒之间的交互作用，以及是否吸烟与吸烟多少之间的交互作用，考虑纳入以下变量进行多因素分析得：

变量

OR

OR_LCI

OR_UCI

p

ageu

1.21

1.17

1.25

<0.01

family_ca

1.00

0.91

1.11

0.943597

sex

0.94

0.88

1.01

0.090639

bmiu

0.97

0.92

1.03

0.310501

smoke1

1.11

1.02

1.20

0.013562

education2

0.97

0.90

1.04

0.363829

education3

1.04

0.93

1.15

0.490749

respdis

1.11

1.02

1.21

0.015094

多因素分析结果显示，在校正了性别等因素后，年龄，吸烟，呼吸系统疾病史是肺癌的危险因素。

遇到年龄 bmi  一般分析都会把这两个变量做个等级

coefficient（-1，）提取模型系数  删除第一行截距项

先做单因素分析 有意义的变量纳入多因素模型 再去检验

单有多有√  单有多无：间接关联 解释变量之间存在共线性  单无多有：单因素分析存在混杂因素的影响 多因素分析消除了混杂因素的影响 所以有意义呈现了出来

绘制Sepal.Length在不同Species种类下的分布的箱式图和散点图，

要求：1）在一张图中绘制，散点图直接绘制在箱式图上；

2）箱式图不同种类填充色不同，并添加图例说明不同颜色与Species种类对应关系；

3）散点为蓝色，且排列方式为随机扰动，而不是一串排列

install.packages("ggplot2")
library(ggplot2)
ggplot(data=iris,aes(x=Species,y=Sepal.Length,colour=Species))+geom_boxplot()+geom_jitter(colour="blue")#随机扰动
ggplot(data=iris,aes(x=Species,y=Sepal.Length,colour=Species))+geom_boxplot()+geom_point(colour="blue")#不是随机扰动

分组变量赋值给颜色属性

ggplot(data = df, mapping = aes(x = x, y = y, colour = z)) + geom_point(size = 3)

分组变量赋值给形状属性

shape = z

geom_jitter(width = 0, height = 0.3, alpha = 0.5)