2021-01-25

最新推荐文章于 2023-12-25 08:30:48 发布

qq_37388306

最新推荐文章于 2023-12-25 08:30:48 发布

阅读量299

点赞数

文章标签：统计学

本文链接：https://blog.csdn.net/qq_37388306/article/details/113139166

版权

R语言----三组以上类别做组件比较并计算两两之间的P值

## 自定义函数--生成比较表格--满足正态性--ANOVA；
tab_anova <- function(dat){
  groupNumber <- 3
  dat2 <- data.frame(t1=as.character(1:(groupNumber+1)))
  Star <- 3
  Over <- ncol(dat)
  for(i in c(Star:Over))
  {
    means<-tapply(dat[[i]],dat[[1]],mean)
    means<-sprintf('%.2f',round(means,4))  
    SD<-tapply(dat[[i]],dat[[1]],sd)
    SD<-sprintf('%.2f',round(SD,4)) 
    M.aov<-aov(dat[[i]]~group, data = dat)
    aov.sum<-summary(M.aov)
    pvalue<-aov.sum[[1]]$`Pr(>F)`[1]
    if (pvalue>0.05) {
      aa<-paste(means,'±',SD)
      pvalue<-sprintf('%.2f',round(pvalue,4))
      aa[(groupNumber+1)]<-pvalue
    }else if (pvalue>0.01) {
      out <- LSD.test(M.aov, "group", p.adj="bonferroni")
      marker<-out$groups
      row_name <- row.names(marker)
      newmarker <- data.frame(row_name,marker$groups)
      newmarker <- newmarker[order(newmarker$row_name),]
      pmarker<-newmarker$marker.groups
      aa<-paste(means,'±',SD,pmarker)
      pvalue<-sprintf('%.2f',round(pvalue,4)) 
      aa[(groupNumber+1)]<-paste(pvalue,'*')
    }else {
      out <- LSD.test(M.aov, "group", p.adj="bonferroni")
      marker<-out$groups
      row_name <- row.names(marker)
      newmarker <- data.frame(row_name,marker$groups)
      newmarker <- newmarker[order(newmarker$row_name),]
      pmarker<-newmarker$marker.groups
      aa<-paste(means,'±',SD,pmarker)
      pvalue<-sprintf('%.2f',round(pvalue,4))
      aa[(groupNumber+1)]<-paste(pvalue,'**')
    }
    
    dat2[i-2] <- aa
    names(dat2)[i-2]<-names(dat[i])
  }
  dat3 <- t(dat2)
  return(dat3)
}

## 自定义函数--生成比较表格--不满足正态性；
tab_Krus <- function(dat){
  groupNumber <- 3
  dat2 <- data.frame(t1=as.character(1:(2*groupNumber+1)))
  Star <- 3
  Over <- ncol(dat)
  for(i in c(Star:Over))
  {
    means<-tapply(dat[[i]],dat[[1]],mean)
    means<-sprintf('%.4f',round(means,4))  
    SD<-tapply(dat[[i]],dat[[1]],sd)
    SD<-sprintf('%.2f',round(SD,2))     
    a = kruskal.test(dat[[i]],dat[[1]])
    Pvalue<-a["p.value"]
    b = dunn.test(dat[[i]],dat[[1]], method = "bonferroni",kw = T)
    pvalue = b["P.adjusted"]
    P_value <- c(Pvalue[[1]],pvalue[[1]])
    P_value<-sprintf('%.2f',round(P_value,2)) 
    Results<-paste(means,'±',SD)
    for ( j in c(1:4)) {
      if(P_value[j]>0.05){
        Results[(groupNumber+j)]<-P_value[j]     
      }           
      else if(P_value[j]>0.01){
        Results[(groupNumber+j)]<-paste(P_value[j],'*')
      }
      else {
        Results[(groupNumber+j)]<-paste(P_value[j],'**')
      }
    }

    dat2[i-2] <- Results
    names(dat2)[i-2]<-names(dat[i])
  }
  dat3 <- t(dat2)
  dat3 <- data.frame(dat3)
  names(dat3) <- c("HC","PD","RBD","p.value")
  return(dat3)
}

```{r}
p_star <- function(p_value){
      p_value <- sprintf('%.3f',round(p_value,3))
      if(p_value > 0.05){
        p_value <- p_value    
      } else if(p_value > 0.01){
        p_value <- paste(p_value,'*')
      } else {
        p_value <- paste(p_value,'**')
      }
      return(p_value)
}

两组比较
```r
## 自定义函数--生成比较表格--不满足正态性；
tab_Krus <- function(dat,Star,Over){
  dat2 <- data.frame(t1=as.character(1:3))
  k = 1
  for(i in c(Star:Over))
  {
    means<-tapply(dat[[i]],dat[[1]],mean)
    means<-sprintf('%.3f',round(means,3))  
    SD<-tapply(dat[[i]],dat[[1]],sd)
    SD<-sprintf('%.3f',round(SD,3))     
    a = kruskal.test(dat[[i]],dat[[1]])
    P_value<-a$p.value
    Results<-paste(means,'±',SD)
    Results[3]<-p_star(P_value)  
    dat2[k] <- Results
    names(dat2)[k]<-names(dat[i])
    k = k+1
  }
  dat3 <- t(dat2)
  dat3 <- data.frame(dat3)
  names(dat3) <- c("HC","PD","p.value")
  return(dat3)
}