(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语言基础课程的课程笔记,仅供参考
推荐链接
发表评论