(1)首先安装、载入两个拓展包并导入数据

#拓展包安装

install.packages('table1')

install.packages('lubridate')

library('table1')

library('lubridate')

#载入数据

load('diagnosis.Rdata')

 (2)进行一些分析前处理

#新生成一个分类变量,主要是检查年度的,以便后续展示

data.wide$TestA.year <- as.factor(ifelse(year(data.wide$TestA.date)<2009,'<2009',

ifelse(year(data.wide$TestA.date) == 2009,'=2009','>2009')))

#attach 这个数据框,这样可以告诉R以后的分析都是用这个数据框,就不用重复输入这个数据框的名称了

attach(data.wide)

一、分析一个变量的分布

(一)分类变量

1、给分类变量加上标签

#给分类变量加上label,即告诉R每一个数字对应的类别的含义

data.wide$Status <- factor(data.wide$Status,

levels = c(0,1),

labels = c('No Disease','Disease'))

data.wide$Sex <- factor(data.wide$Sex,

levels = c(0,1),

labels = c('Male','Female'))

 2、描述分类变量的统计量:频数、百分比

#频数

mytable_Sex <- table(Sex)

mytable_Sex

#百分比

prop.table(mytable_Sex)

#百分比的假设检验,H0: 男性的比例为50% , 这里P=0.5表示我们的检验目标为50%

prop.test(mytable_Sex,p=0.5)

> #频数

> mytable_Sex <- table(Sex)

> mytable_Sex

Sex

Male Female

52 48

> #百分比

> prop.table(mytable_Sex)

Sex

Male Female

0.52 0.48

> #百分比的假设检验,H0: 男性的比例为50% , 这里P=0.5表示我们的检验目标为50%

> prop.test(mytable_Sex,p=0.5)

1-sample proportions test with continuity correction

data: mytable_Sex, null probability 0.5

X-squared = 0.09, df = 1, p-value = 0.7642

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

0.4183183 0.6201278

sample estimates:

p

0.52

(二)连续变量

1、描述变量分布的统计量

均值、标准差、中位数、四分位距、极差

#连续变量的统计描述

summary(Age)

> #连续变量的统计描述

> summary(Age)

Min. 1st Qu. Median Mean 3rd Qu. Max.

18.41 33.59 39.63 39.79 45.38 71.71

2、t检验

(当变量服从正态分布时),想要知道某个连续变量的均值,是否等于某个标准值,可以使用t检验

#t 检验,H0:平均年龄为35岁

t.test(Age,mu=35)

> t.test(Age,mu=35)

One Sample t-test

data: Age

t = 4.6678, df = 99, p-value = 9.569e-06

alternative hypothesis: true mean is not equal to 35

95 percent confidence interval:

37.75321 41.82459

sample estimates:

mean of x

39.7889

 t:t统计量df:自由度p-value:P值alternative hypothesis:备择假设95 percent confidence interval:95%置信区间mean of x:样本均值

3、Wilcoxon符号秩检验

(当变量不服从正态分布时) 

#wilcoxon符号秩检验,H0:平均年龄为35

wilcox.test(Age,mu =35)

> #wilcoxon符号秩检验,H0:平均年龄为35

> wilcox.test(Age,mu =35)

Wilcoxon signed rank test with continuity correction

data: Age

V = 3773.5, p-value = 1.778e-05

alternative hypothesis: true location is not equal to 35

二、分析两个变量的关系(组间比较)

1、分类变量&分类变量

例如分析:患病和性别的关系分析方法:列联表和卡方检验(检验)

#分类变量间的卡方检验

#列联表

mytable <- table(Sex,Status)

mytable

#卡方检验

chisq.test(mytable)

> #分类变量间的卡方检验

> #列联表

> mytable <- table(Sex,Status)

> mytable

Status

Sex No Disease Disease

Male 37 15

Female 23 25

> #卡方检验

> chisq.test(mytable)

Pearson's Chi-squared test with Yates' continuity correction

data: mytable

X-squared = 4.6892, df = 1, p-value = 0.03035

2、连续变量&分类变量

(1)(正态分布)&(二分类)

例如:年龄和性别的关系分析方法:t检验

#正态分布&二分类,假设两组组间方差齐

t.test(Age ~ Sex,var.equal = TRUE)

> #正态分布&二分类,假设两组组间方差齐

> t.test(Age ~ Sex,var.equal = TRUE)

Two Sample t-test

data: Age by Sex

t = -2.4941, df = 98, p-value = 0.0143

alternative hypothesis: true difference in means between group Male and group Female is not equal to 0

95 percent confidence interval:

-8.963487 -1.019910

sample estimates:

mean in group Male mean in group Female

37.39288 42.38458

(2)(非正态分布)&(二分类)

例如:某检验指标值与性别的关系分析方法:Mann-Whitney U检验(Wilcoxon秩和检验)注意:检验的是两组的“分布”而不是均值

#非正态分布&二分类

wilcox.test(TestA ~ Sex)

 

> #非正态分布&二分类

> wilcox.test(TestA ~ Sex)

Wilcoxon rank sum test with continuity correction

data: TestA by Sex

W = 922, p-value = 0.02472

alternative hypothesis: true location shift is not equal to 0

(3)(正态分布)&(多分类变量) 

例如:检查某指标值(log转换后)&和检查年度的关系分析方法:单变量方差分析(ANOVA)

#正态分布&多分类变量

fit <- aov(log(TestA) ~ TestA.year)

model.tables(fit,"means")

summary(fit)

 

> #正态分布&多分类变量

> fit <- aov(log(TestA) ~ TestA.year)

> model.tables(fit,"means")

Tables of means

Grand mean

1.479017

TestA.year

<2009 =2009 >2009

1.523 1.397 1.556

rep 36.000 41.000 23.000

> summary(fit)

Df Sum Sq Mean Sq F value Pr(>F)

TestA.year 2 0.48 0.2415 0.286 0.752

Residuals 97 81.97 0.8451

 (4)(非正态分布)&多分类变量

例如:检查指标值和检查年度的关系分析方法:Kruskal-Wallis检验

#非正态分布&多分类变量

kruskal.test(TestA ~ TestA.year)

 

> #非正态分布&多分类变量

> kruskal.test(TestA ~ TestA.year)

Kruskal-Wallis rank sum test

data: TestA by TestA.year

Kruskal-Wallis chi-squared = 0.225, df = 2, p-value = 0.8936

3、连续变量&连续变量

例如:指标值与年龄

(1)两个连续变量均服从正态分布

Pearson相关系数

#Pearson相关系数

cor(TestA,TestB,method = "pearson")

cor.test(TestA,TestB,method = "pearson")

 

> #Pearson相关系数

> cor(TestA,TestB,method = "pearson")

[1] 0.3480745

> cor.test(TestA,TestB,method = "pearson")

Pearson's product-moment correlation

data: TestA and TestB

t = 3.6756, df = 98, p-value = 0.0003872

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval:

0.1627858 0.5096489

sample estimates:

cor

0.3480745

(2)两个连续变量不均服从正态分布

Spearman相关系数 #spearman相关系数

cor(TestA,TestB,method = "spearman")

> #spearman相关系数

> cor(TestA,TestB,method = "spearman")

[1] 0.3769439

4、小结

要分析的变量数值变量正态分布t检验t检验单因素方差分析非正态分布 wilcoxon 符号秩检验 Mann-Whitney U检验 Kruskal=Wallis检验分类变量二分类z检验卡方检验卡方检验多分类z检验卡方检验卡方检验检验对照表1二分类多分类分类变量

三、描述统计和基本特征表(Table1)

(1)基本特征表 

使用table1包,自动生成基本特征表

#总体特征表

table1(~ Status+Sex+Age+TestA+TestB+TestC+TestA.year,data = data.wide)

#分组特征表

table1(~Sex+Age+TestA+TestB+TestC+TestA.year|Status,data = data.wide)

(2)在特征表中加入单位

# add units to table 1

units(data.wide$Age) <- "years"

table1(~ Sex + Age + TestA + TestB + TestC + TestA.year| Status, data=data.wide)

(3)在特征表中加入更多的统计数据

# add more statistics to table 1

table1(~ Sex + Age + TestA + TestB + TestC + TestA.year| Status, data=data.wide,

render.continuous=c(.="Mean (SD)", .="Median [Min, Max]",.="Median [Q1, Q3]",.="Median [IQR]"))

(4)在特征表中加入P值

# add p value

load('diagnosis.RData')

data.wide$Status <- factor(data.wide$Status,

levels=c(0, 1, 2),

labels=c("No Disease", "Disease", "P-value"))

data.wide$Sex <- factor(data.wide$Sex,

levels=c(0,1),

labels=c("Male", "Female"))

# creat P values

data <- data.wide

outcome <- data.wide$Status

rndr <- function(x, name, ...) {

if (length(x) == 0) {

y <- data[[name]]

s <- rep("", length(render.default(x=y, name=name, ...)))

if (is.numeric(y)) {

p <- t.test(y ~ outcome)$p.value

} else {

p <- chisq.test(table(y, droplevels(outcome)))$p.value

}

s[2] <- sub("<", "<", format.pval(p, digits=3, eps=0.001))

s

} else {

render.default(x=x, name=name, ...)

}

}

rndr.strat <- function(label, n, ...) {

ifelse(n==0, label, render.strat.default(label, n, ...))

}

# make table 1

table1(~ Sex + Age + TestA + TestB + TestC | Status, data=data.wide,

render.continuous=c(.="Mean (SD)"),droplevels=F, render=rndr, render.strat=rndr.strat, overall=F)

#以上内容均来自医咖会R语言基础课程的课程笔记,仅供参考

推荐链接

评论可见,请评论后查看内容,谢谢!!!评论后请刷新页面。