目录
十九：两独立样本非参数检验
###两独立样本非参数检验
#1.正态性检验
mydata <- as.data.frame(mydata)
tapply(mydata[,5],
mydata[,13],
shapiro.test)
#2、两独立样本非参数检验
#wilcox.test(被检验变量~分组变量，paired = T/F)#独立：F
res1 <- wilcox.test(mydata$HDL~mydata$代谢综合征,paired = FALSE)#两独立样本非参数检验的p>0.05,意味着两组间无显著差异
#如果p<0.05,则认为差异显著，需要比较两组的中位数
aggregate(x = list(mydata$HDL),
by = list(mydata$代谢综合征),
FUN = median)
Z <- qnorm(res1$p.value/2)
#多变量的综合演示
#1、正态性检验
k <- NULL
for(i in 2:12){
k <- tapply(mydata[,i],
mydata[,13],
shapiro.test)
if(k$有$p.value<0.05){
l <- paste0('第',i,'组有代谢综合征不服从正态分布')
print(l)}
else{
if(k$无$p.value<0.05){
g <- paste0('第',i,'组无代谢综合征不服从正态分布')
print(g)}
}
}
#以上结果清晰显示，所有变量至少都在一组上非正态分布
#因此选择两独立样本非参数检验
#且数据的表达形式均选择中位数（上四分位数~下四分位数）中位数表达 <- function(x){
paste0(median(x),
'(',
quantile(x,prob = 0.25),
'~',
quantile(x,prob = 0.75),
')')
}res2 <- aggregate(x = list(mydata[,2:12]),
by = list(mydata$代谢综合征),
FUN = 中位数表达)res2 <- t(res2)
res2 <- res2[-1,]
colnames(res2) <- c('无','有')
#3、两独立样本非参数检验
p <- NULL
for(i in 2:12){
pi <- wilcox.test(mydata[,i]~mydata[,13],
paired = FALSE)
p <- rbind(p,pi)
}
p <- p[,3]
View(p)p <- do.call(rbind,p)
colnames(p) <- 'p'
p <- data.frame(p)
p$z <- qnorm(p$p/2)final <- cbind(res2,round(p[,2],3),round(p[,1],3))
colnames(final) <- c('无','有','z','p')从以上的检验结果可以清晰得出：
- 有代谢综合征和无代谢综合征两组在LDL、CHOL、AST三个指标上不存在显著差异，P值全部大于0.05。
- 有代谢综合征和无代谢综合征两组在BMI上存在显著差异，P<0.001， 更进一步比较中位数得出，有代谢综合征组的BMI显著高于无代谢综合征组。
二十：单因素方差分析
#####单因素方差分析
#导入SPSS格式的数据
library(Hmisc)
mydata <- spss.get('NO.24 配套数据.sav',use.value.labels = TRUE)#1、正态性检验
tapply(mydata[,3],
mydata[,2],
shapiro.test)
#2、判断组间的方差齐性
#方差齐性检验使用的是car程序包中的leveneTest函数
library(car)
#leveneTest(被检验变量~分组变量)
leveneTest(mydata$胃半排空时间~mydata$分组)
#莱文方差齐性检验得出的p>0.05,则认为组间方差无显著差异
#3、单因素方差分析
#单因素方差分析使用的函数是oneway.test
res1 <- oneway.test(mydata$胃半排空时间~mydata$分组,
var.equal = TRUE)#方差齐
res1
#以上单因素方差分析得出p<0.05,
#意味着三组间的胃半排空时间的均值存在着显著差异
#则需要继续进行两两比较
#4、两两比较（基于方差齐）
#两两比较需要用到的程序包是PMCMRplus
#PMCMRplus程序包既可以实现方差齐的两两比较
#也可以实现方差不齐的两两比较
library(PMCMRplus)（基于方差齐）
#lsdTest(被检验变量~分组)
mul <- lsdTest(mydata$胃半排空时间~mydata$分组)mul$p.value
#校正(不太懂)
padj <- p.adjust(mul$p.value,method = 'bonferroni')padj
#溃疡样型和动力障碍样型、非特异样型之间的胃半排空存在显著差异
#而动力障碍和非特异性型之间无显著差异，p>0.05.
#溃疡样型和动力障碍样型、非特异样型之间的胃半排空存在显著差异
#而动力障碍和非特异性型之间无显著差异，p>0.05.
#更进一步结合均值
aggregate(x = list(mydata[,3]),
by = list(mydata[,2]),
FUN = mean)
#得出溃疡样型的胃半排空时间显著高于动力障碍样型
#多变量的单因素方差
#1、正态性检验
for(i in 3:8){
print(colnames(mydata)[i])
print(tapply(mydata[,i],
mydata[,2],
shapiro.test))
}
#比较三组间在这"胃半排空时间"、"餐后2小时残留率"、"NO"
#三个变量的差异，可以选择单因素方差分析#2、描述性统计分析
均值标准差 <- function(x){
paste0(round(mean(x),2),'±',round(sd(x),2))
}res2 <- aggregate(x = list(mydata[,3:5]),
by = list(mydata[,2]),
FUN = 均值标准差)
res2 <- t(res2)
res2 <- res2[-1,]
colnames(res2) <- c('溃疡样型','动力障碍样型','非特异性型')
#3、方差齐性检验
for(i in 3:5){
print(colnames(mydata)[i])
print(leveneTest(mydata[,i]~mydata$分组))
}
#这 "胃半排空时间"、"NO"方差齐
#"餐后2小时残留率"方差不齐
#莱文方差齐性检验得出的p>0.05,则认为组间方差无显著差异
#4、单因素方差分析
a <- oneway.test(mydata$胃半排空时间~mydata$分组,
var.equal = TRUE)a$statistic
b <- oneway.test(mydata$NO~mydata$分组,
var.equal = TRUE)
c <- oneway.test(mydata$餐后2小时残留率~mydata$分组,
var.equal = FALSE)r1 <- rbind(a$statistic,b$statistic,c$statistic)
r2 <- rbind(a$p.value,b$p.value,c$p.value)
final <- cbind(res2,r1,r2)
colnames(final)[5] <- 'p'
#以上单因素方差（齐）分析得出p<0.05,
#意味着三组间的胃半排空时间的均值存在着显著差异
#则需要继续进行两两比较#5、两两比较
#lsd是方差齐的两两比较
res5 <- lsdTest(mydata$NO~mydata$分组)
res5
#校正(不太懂)
padj2 <- p.adjust(res5$p.value,method = 'bonferroni')
padj2
#方差不齐的两两比较是Tamhane T2两两比较
tamhaneT2Test(mydata$餐后2小时残留率~mydata$分组)
二十一：K独立样本非参数
######K独立样本非参数检验
#1、导入SPSS格式的数据
library(Hmisc)
mydata <- spss.get('NO.25 配套数据.sav',use.value.labels = TRUE)
#2、正态性检验
tapply(mydata$SS,
mydata$分组,
shapiro.test)
#3、K独立样本非参数检验
#K独立样本非参数检验使用的函数是kruskal.test()
#kruskal.test()
kruskal.test(mydata$SS~mydata$分组)
#以上检验得出p=0.025<0.05,意味着三组间SS存在显著差异#得出了显著差异结论之后，还需要继续进行两两比较
#得出具体每两组之间的差异
#4、两两比较
#K独立样本非参数检验的两两比较会用到wmc函数
source('wmc.txt')
wmc(mydata$SS~mydata$分组,data = mydata,
method = 'holm',sort = FALSE)#以上两两比较得出动力障碍样型和溃疡样型
#之间的SS存在显著差异。
#为了更进一步得出具体什么差异，哪组高：
#哪组低？
#还需要继续比较中位数。
#5、计算中位数
#非参数检验的数据为非正态数据
#因此数据的表达形式是：中位数（下四分位数~上四分位数）
aggregate(x = list(mydata$SS),
by = list(mydata$分组),
FUN = median)
#最终结论是：动力障碍样型显著低于溃疡样型
######多个变量一起进行K独立样本非参数检验
#1、正态性
for(i in 6:8){
print(colnames(mydata)[i])
print(tapply(mydata[,i],
mydata$分组,
shapiro.test))
}
#任意一组非正态，选择K独立
#基于以上结果，选择K独立样本非参数检验比较三组的差异#2、首先进行描述性统计分许
#中位数（上四分位数~下四分位数）中位数描述 <- function(x){
paste0(round(median(x),2),
'(',
round(quantile(x,prob = 0.25),2),
'~',
round(quantile(x,prob = 0.75),2),
')')
}res1 <- aggregate(x = list(mydata[,6:8]),
by = list(mydata$分组),
FUN = 中位数描述)
res1 <- t(res1)res1 <- res1[-1,]
colnames(res1) <- c('溃疡样型','动力障碍样型','非特异性型 ')#3、K独立样本非参数检验
k <- NULL
for(i in 6:8){
print(colnames(mydata)[i])
print(kruskal.test(mydata[,i]~mydata$分组))
l <- kruskal.test(mydata[,i]~mydata$分组)
k <- rbind(k,l)
}
#提取k中的KW卡方
chip <- NULL
for(i in 1:3){
c <- k[i,1]
chip <- rbind(chip,c)
}#提取k中的p
p <- NULL
for(i in 1:3){
pp <- as.numeric(k[i,3])
p <- rbind(p,pp)
}# 将列表转换为向量
values <- unlist(chip)
numbers <- as.numeric(gsub("Kruskal-Wallis chi-squared = ", "", values))
print(numbers)res1 <- cbind(res1,round(numbers,3))
res1 <- cbind(res1,round(p,6))
colnames(res1)[4] <- 'KW卡方'
colnames(res1)[5] <- 'P'
小知识：提取list列表中的数值
在以上列表的chip中，第一列第一行样式为：c(`Kruskal-Wallis chi-squared` = 7.90805836623058)，现在想要提取里面的数值。
方法：
# 将列表转换为向量
values <- unlist(chip)
unlist(k)会合并列表中的所有向量元素到一个单一的向量中。- 当合并后的向量有元素名称时（例如
Kruskal-Wallis chi-squared = 7.90805836623058），这些名称可以很容易地通过names()函数来获取。- 一旦我们有了这些名称，我们可以利用字符串处理函数（例如
gsub()）来提取和操作这些名称中的特定部分，如你所需的数值。如果我们不将列表转换为向量，我们可能需要用更复杂的方式遍历列表的每个元素，并对每个元素单独执行提取操作，这会使代码变得更复杂和冗长。通过将列表转换为向量，我们简化了数据处理的流程。
# 使用正则表达式提取数值
numbers <- as.numeric(gsub("Kruskal-Wallis chi-squared = ", "", values))
gsub("Kruskal-Wallis chi-squared = ", "", values):
gsub()是一个用于字符串替换的函数。在这个函数中:
- 第一个参数
"Kruskal-Wallis chi-squared = "是想要从字符串中替换掉的模式。- 第二个参数
""是你想要将模式替换成的内容。在这里，我们想要将前面提到的模式替换为一个空字符串，也就是删除它。- 第三个参数
values是包含字符串的向量，这是我们想要执行替换操作的数据。因此，这部分的代码的功能是：从
values中的每一个字符串元素中删除"Kruskal-Wallis chi-squared = "这部分内容。as.numeric(...):
这是一个转换函数，用于将数据转换为数字格式。在这里，由于使用
gsub()后我们得到的是一个字符串向量（具体来说，是不再包含“Kruskal-Wallis chi-squared = ”的数字字符串），我们需要将这些数字字符串转换为实际的数值格式。综上，整行代码的目的是从
values向量中删除 “Kruskal-Wallis chi-squared = ” 这部分内容，并将得到的数字字符串转换为数值格式。> print(numbers) [1] 11.732550 7.417471 7.908058
#4、两两比较
for(i in 6:8){
print(colnames(mydata)[i])
print(wmc(mydata[,i]~mydata$分组,data = mydata,
method = 'holm',sort = FALSE))
}
#黄色部分存在显著差异
二十二：秩和检验
#1、比较不同婚姻状况的学历差异
table(mydata$婚姻状况)
#使用K独立样本秩和检验
kruskal.test(mydata$文化~mydata$婚姻状况)
#如果存在显著差异，则继续计算中位数
#分组计算中位数aggregate
#aggregate()
#2、比较男和女在学历等级上的差异
#两组秩和检验，两独立样本的秩和检验
#mydata$文化为factor类型，需要转化成数值型
table(mydata$文化)#从数量高到低进行修改
library(do)
mydata$文化新 <- mydata$文化
mydata$文化新 <- Replace(mydata$文化新,
from = '文盲',
to = 1)
mydata$文化新 <- Replace(mydata$文化新,
from = '小学',
to = 2)
mydata$文化新 <- Replace(mydata$文化新,
from = '初中',
to = 3)
mydata$文化新 <- Replace(mydata$文化新,
from = '高中',
to = 4)
mydata$文化新 <- Replace(mydata$文化新,
from = '本科',
to = 5)
mydata$文化新 <- Replace(mydata$文化新,
from = '研究生及以上',
to = 6)
table(mydata$文化新)
mydata$文化新 <- as.numeric(mydata$文化新)
wilcox.test(mydata$文化新~mydata$性别)
#得出了存在显著差异的结论后，还需要继续比较中位数
aggregate(x=list(mydata$文化新),
by = list(mydata$性别),
FUN = quantile)
#得出男性学历高于女性
二十三：配对样本T检验
###配对样本T检验
#1、导入数据
library(Hmisc)
mydata <- spss.get('NO.27 配套数据.sav',use.value.labels = TRUE)mydata$差值 <- mydata$空腹血糖.治疗前-mydata$空腹血糖.治疗后
apply(mydata[,2:4],2, shapiro.test)
#通过了正态性检验，因此比较前后的
#空腹血糖差异选择的统计学方法是
#配对样本T检验
t.test(mydata$空腹血糖.治疗前,
mydata$空腹血糖.治疗后,
paired = TRUE)
#得出前后差异显著
#3、描述性统计分析
source('均值标准差函数.txt')
res <- apply(mydata[,2:3], 2, 均值标准差函数)
res <- data.frame(res)
二十四：配对样本非参数检验
#1、比较前后疼痛评分的差异
#1、正态性检验
apply(mydata[,2:3],2,shapiro.test)
#评分均为非正态，因此选择独立样本非参数检验
#1.2 配对样本非参数检验
res1 <- wilcox.test(mydata$疼痛评分.治疗前,
mydata$疼痛评分.治疗后,
paired = TRUE)
res1$p.value
Z <- qnorm(res1$p.value/2) #-5.711195
#存在显著差异
#1.3 制作表格
source('中位数函数.txt')
res2 <- apply(mydata[,2:3], 2, 中位数描述)
res2 <- data.frame(res2)
#2、比较严重程度的差异
#严重程度变量属于有序等级变量
#因此，不需要进行正态性检验
#因为是等级变量，直接采用秩和检验
res3 <- wilcox.test(mydata$严重程度.治疗前,
mydata$严重程度.治疗后,
paired = TRUE)
Z2 <- qnorm(res3$p.value/2)
Z2
res4 <- apply(mydata[,4:5], 2, 中位数描述)
res4 <- data.frame(res4)res2 <- t(res2)
res4 <- t(res4)
final <- rbind(res2, res4)Zval <- c(Z, Z2)
Zval
Pval <- c(res1$p.value, res3$p.value)
final <- cbind(final, Zval, Pval)
write.csv(final,'final.csv')从以上的分析结果可以清晰得出：
- 治疗后，疼痛评分显著下降，从原来的42显著降到26，P<0.05；
- 治疗后，严重程度显著下降，从原来的中重度显著降到中度，P<0.05
二十五：K个相关样本非参数检验（肯德尔系数）
####K个相关样本非参检验
#1、friedman.test()(差异显著）
#比较三个状态下的连续型变量的差异
#1.1 正态性检验
apply(mydata1[,2:4],2,shapiro.test)
#其中任意一个非正态，
#则选择K个相关样本非参检验中的friedman检验#1.2friedman检验对数据的要求并整理
#firedman检验要求数据必须是矩阵matrix
mal <- as.matrix(mydata1[,2:4])
#此时矩阵为114个元素，混为一起，所以要告诉软件分成38组，每组三个
dimnames(mal) <- list(1:38,c('术前','术后一周','术后四周'))#1.3 firedman检验
friedman.test(mal)
#得出P<0.05, 意味着三个状态存在显著差异，因此，
#还需要继续进行两两比较，确定具体哪两个状态之间存在显著差异
#1.4 事后多重比较
#PMCRplus程序包是专门用于两两比较的程序包
library(PMCMRplus)
frdAllPairsNemenyiTest(mal,p.adjust = 'bondferroni')
#1.5 表格制作
source('中位数函数.txt')
res1 <- apply(mal,2,中位数描述)
res1 <- data.frame(res1)
res2 <- friedman.test(mal)
res1$chiq <- c(res2$statistic,'NA','NA')
res1$p <- c(res2$p.value,'NA','NA')
write.csv(res1,'res1.csv')
#2、friedman.test()(差异不显著的情况)
#更进一步进行Kendall协同系数检验
ma2 <- as.matrix(mydata2[,2:5])
dimnames(ma2) <- list(1:50,c('评委1','评委2','评委3','评委4'))#2.1 首先进行差异比较
friedman.test(ma2)
#得出P=0.26>0.05,意味着4个评委的结果差异不显著
#即4个评委的结果一致。#接下来，需要选择Kendall协同系数检验计算一致程度
#Kendall协同系数检验需要用到irr程序包中的kendall函数
library(irr)
kendall(ma2)
#kendall's W检验得出P<0.05,意味着四个评委的结论显著一致
#一致程度为0.846，说明四个评委给出的评价结果一致程度较高。
#3、Cochran's Q检验
#Cochran's Q检验用到的是DescTools程序包中的CochranQ函数
library(DescTools)
#CochranQTest()也是要求矩阵matrix形式的数据
ma3 <- as.matrix(mydata3[,2:5])
dimnames(ma3) <- list(1:57,
c('评委1','评委2','评委3','评委4'))
CochranQTest(ma3)
#得出P=0.601>0.05意味着4个评委的结论不存在显著差异。
#假设，存在显著差异，怎么进行事后多重比较？
#两两比较：配对卡方
choose(4,2)
4210
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
