 示例 评价某药物耐受性及安全性的I期临床实验中，将符合纳入标准的30名健康自愿者随机分为3组，每组10名，各组注射剂量分别为0.5U、1U、2U，观察49小时部分凝血活酶时间（s）。不同剂量组的部分凝血活酶时间是否不同？ 示例来源：李康,贺佳等.医学统计学（第6版）.北京:人民卫生出版社,2013. 很明显，剂量可以看做是研究因素，有3个水平，这种完全随机设计的资料首先考虑的就是方差分析。如果方差分析有统计学意义，我们可能还会进一步通过两两比较考察到底是那两组有统计学差异。此外，本例中分组因素为有序变量，我们可能还想知道部分凝血活酶时间是否会随着给药剂量的增加而呈现某种变化趋势，此时我们可以进行趋势检验，检验是否满足线性、二次、三次等多项式变化。对于分组变量和因变量都是分类资料的趋势检验，我们曾在《线性趋势检验》中做过介绍，当前示例是分组变量有序、结局变量为连续的资料，我们可以采用趋势方差检验，即在方差分析中采用多项式对比检验。当然，方差分析也有自己的适用条件：独立性、正态性（各水平因变量服从正态分布，准确地说应该是模型残差服从正态分布）、方差齐性（各水平的总体具有相同的方差），需要进行评估。 【1】数据导入。不同类型的数据载入涉及到不同的程序包，相同数据类型的导入也有不同的包，excel数据的载入可参见《为什么是R？》，SPSS数据的载入可参见《正态分布的检验》、《方差齐性检验》，STATA数据的载入可参见《描述性统计分析》。本例我们使用的是RStudio，导入可以使用菜单操作 Files >> Import Dataset >> From Excel… library(readxl) fdata <- read_excel("D:/Temp/fdata.xlsx") library(stats) Fa<-aov(time~dose,data=fdata) summary(Fa) 结果显示可以认为3个给药剂量的部分凝血活酶时间不同（F=6.524，P=0.00488<0.05）。 我们也可以采用拟合线性模型的方式来进行方差分析。不论是t检验、方差分析，还是非参数检验，其实都是线性回归的特殊形式而已，感兴趣的可以去阅读Jonas Kristoffer Lindeløv的文章：Common statistical tests are linear models (or: how to teach stats)，我们贴出文章的一张总结表： 线性回归用到函数lm，lm {stats}：used to fit linear models. It can be used to carry out regression, single stratum analysis of variance and analysis of covariance (although aov may provide a more convenient interface for these). Usage:lm(formula, data, subset, weights, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,singular.ok = TRUE, contrasts = NULL, offset, ...) lrfit<-lm(time~dose,data=fdata) summary(lrfit) 【3】事后比较。方差分析后如果发现各组之间存在显著的统计学意义，接下来我们可能想知道哪几个组间会存在差异，这就涉及到两两比较（或者叫多重比较）。多重比较的方法有很多，比如SPSS中提供的两两多重比较方法如下： 这些在R中都不是事儿，因为R有包，“包”治百病！比如： pairwise.t.test {stats}：Calculate pairwise comparisons between group levels with corrections for multiple testing. pairwise.t.test Usage pairwise.t.test(x, g, p.adjust.method = p.adjust.methods,pool.sd = !paired, paired = FALSE,alternative = c("two.sided", "less", "greater"),...) TukeyHSD Usage TukeyHSD(x, which, ordered = FALSE, conf.level = 0.95, ...) PostHocTest {DescTools}：A convenience wrapper for computing post-hoc test after having calculated an ANOVA.仅用于R3.6.3及以上版本。 TukeyHSD Usage PostHocTest(x, which = NULL,method = c("hsd", "bonferroni", "lsd", "scheffe", "newmankeuls", "duncan"),conf.level = 0.95, ordered = FALSE, ...) glht {multcomp}：General linear hypotheses and multiple comparisons for parametric models, including generalized linear models, linear mixed effects models, and survival models. glht Usage glht(model, linfct, ...) pairwise.t.test(fdata$time, fdata$dose, p.adj = "bonf") ##p.adj方法有"holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr", "none"，当前命令以bonferroni进行校正 plot(TukeyHSD(Fa, "dose")) library(DescTools) PostHocTest(Fa,,method = c("lsd")) glht {multcomp}函数可用于多种模型，如aov、lm、glm等，但需要将要比较的变量指定为因子变量。命令清单与结果如下，具体结果就不再赘述了。 Dose<-factor(fdata$dose) #采用glht函数需要首先将要比较的变量指定为因子变量 Fc<-aov(fdata$time~Dose) library(multcomp) summary(glht(Fc, linfct=mcp(Dose="Tukey")))
上述命令清单也可以用比对系数的形式进行，命令如下： contrast<-rbind("0.5U vs 1U" = c(1,-1,0),"0.5U vs 2U" = c(1,0,-1),"1U vs 2U" = c(0,1,-1)) multic<-glht(Fc, linfct=mcp(Dose=contrast)) summary(multic) Dose<-factor(fdata$dose) lrfitc<-lm(fdata$time~Dose) summary(glht(lrfitc, linfct=mcp(Dose="Tukey")))
【4】趋势检验。给药剂量0.5U、1U、2U为有序变量，我们可以考察随着给药剂量的增加，部分凝血活酶时间呈二次项变化趋势还是呈线性趋势。在R中要进行趋势方差分析，需要先将分析的变量设置为有序因子，然后采用正交多项式的对照，就可以进行趋势分析（线性、二次、三次等），但比较遗憾的是只能用于等距水平的有序因子。本例实际上不等距，因此仅做演示。在SPSS中，一般线性模型中的单变量分析（Univariate）的Polynomial也提供了同样的等距有序因子的趋势分析，而One-Way ANOVA的Polynomial比对可进行等距和不等距的因子分析。 Dose_ord<-factor(fdata$dose,order=TRUE,levels=c("0.5U", "1U", "2U")) #设置有序因子 lrfit_ord<-lm(fdata$time~Dose_ord,contrasts = "contr.poly") #contr.poly用于趋势分析（线性、二次、三次等）和等距水平的有序因子 summary(lrfit_ord) plot(Dose_ord, fdata$time) #用图形展示变化趋势，或者使用命令：boxplot(time~dose,data=fdata) 结果显示二次项(quadratic)结果的系数（Dose_ord.Q）具有统计学意义，表明二次项的系数不为0，即随着给药剂量的增加，部分凝血活酶时间呈二次项变化趋势，系数估计值为-2.8332<0，表明二次项开口向下，大体呈先升后降的趋势，后面的趋势图也证实的确如此。但需要说明的是本例0.5U、1U、2U并不等距，
统计分析方法都有自己的适用条件，方差分析要求独立、正态和方差齐同，最后我们评估一下下正态性和方差齐性。 library(stats) by(fdata$time,fdata$dose,shapiro.test) #对fdata中的time变量按dose分组进行shapiro检验
library(car) library(carData) leveneTest(time~dose,data=fdata,center=mean) #levene方差齐性检验，默认以中位值进行检测
 2020.05.11 |
|
|
来自： Memo_Cleon > 《待分类》
