王福鑫 趙志龍 矯筱蔓

【關(guān)鍵詞】R語(yǔ)言；化妝品功效；統(tǒng)計(jì)分析

2021 年4 月，國(guó)家藥監(jiān)局發(fā)布的《化妝品功效宣稱評(píng)價(jià)規(guī)范》為化妝品及原料提供了功效性評(píng)價(jià)方法的要求和規(guī)范?；瘖y品功效評(píng)價(jià)的研究采用哪種統(tǒng)計(jì)分析工具，對(duì)研究結(jié)論的準(zhǔn)確性、科學(xué)性把控尤為重要。R語(yǔ)言作為一門現(xiàn)代統(tǒng)計(jì)分析工具，不僅在統(tǒng)計(jì)方面有很強(qiáng)大的功能，同時(shí)也具有極強(qiáng)的可視化功能。R語(yǔ)言常用于生物、醫(yī)學(xué)等相關(guān)領(lǐng)域?；瘖y品研究領(lǐng)域中鮮有利用R語(yǔ)言進(jìn)行統(tǒng)計(jì)分析的報(bào)道。在化妝品功效評(píng)價(jià)研究中，R語(yǔ)言可更好地幫助研究人員設(shè)計(jì)試驗(yàn)方案，進(jìn)行數(shù)據(jù)總結(jié)和分析。本文著眼于用R語(yǔ)言實(shí)現(xiàn)化妝品功效評(píng)價(jià)數(shù)據(jù)處理中常見的差異關(guān)系分析，并舉例提供參考代碼。

1 正態(tài)分布和方差齊性

正態(tài)分布是統(tǒng)計(jì)學(xué)中最重要、最常見的分布之一。較常用的差異關(guān)系分析如t 檢驗(yàn)和方差分析均需要對(duì)數(shù)據(jù)分布情況是否為正態(tài)分布和方差齊性進(jìn)行驗(yàn)證。

對(duì)于數(shù)據(jù)的正態(tài)性檢驗(yàn)主要可以通過Q-Q 圖（見圖1）、Pearson 擬合優(yōu)度卡方檢驗(yàn)、K-S擬合優(yōu)度檢驗(yàn)和S-W檢驗(yàn)等。在R 語(yǔ)言中，使用chisq.test（）函數(shù)完成Pearson 擬合優(yōu)度卡方檢驗(yàn)，使用ks.test（）函數(shù)來(lái)進(jìn)行K-S 檢驗(yàn)，使用shapiro.test（）函數(shù)來(lái)完成S-W正態(tài)性檢驗(yàn)。

如果樣本數(shù)在50 以下選擇S-W檢驗(yàn)，樣本數(shù)在50 以上選擇K-S 檢驗(yàn)。Q-Q圖可最直觀地顯示樣本數(shù)據(jù)是否符合正態(tài)性分布。例如，有100 名受試者參加功效評(píng)價(jià)試驗(yàn)，擬對(duì)其進(jìn)行數(shù)據(jù)正態(tài)性檢驗(yàn)，用rnorm 模擬生成隨機(jī)數(shù)作為樣本，利用基礎(chǔ)包中的函數(shù)分析，程序及結(jié)果如下：

> set.seed（12）

> nordata<-rnorm（100，mean=2，sd=5）

> par（pty="s"）

> qqnorm（nordata，pch=1， frame=FALSE）

> qqline（nordata，col="steelblue"，lwd=2）

進(jìn)一步檢驗(yàn)數(shù)據(jù)是否服從標(biāo)準(zhǔn)正態(tài)分布可使用ks.test（）

函數(shù)，程序及結(jié)果如下：

> ks.test（x=nordata，"pnorm"）

## One-sample Kolmogorov-Smirnov test

## data：nordata

## D=0.443，p-value<2.2e-16

## alternative hypothesis：two-sided

從輸出的p-value<2.2e-16 來(lái)看，說(shuō)明nordata 不服從標(biāo)準(zhǔn)正態(tài)分布（隨機(jī)數(shù)設(shè)置為平均值為2，sd 為5）。

方差齊性檢驗(yàn)是檢驗(yàn)不同樣本的總體方差是否相同的一種方法。常用的檢驗(yàn)方法有F 檢驗(yàn)、Bartlet 檢驗(yàn)和Levene 檢驗(yàn)。建議首選使用Levene 進(jìn)行樣本齊性檢驗(yàn)，因其更穩(wěn)健且不依賴總體分布情況，還可用于多組樣本。例如，對(duì)3 個(gè)實(shí)驗(yàn)組受試者進(jìn)行方差齊性檢驗(yàn)（t1 對(duì)照組、t2 陽(yáng)性組和t3 試驗(yàn)組均以隨機(jī)數(shù)為例做樣本），使用car 包中的leveneTest（）函數(shù)進(jìn)行Levene檢驗(yàn)，程序和結(jié)果如下：

> library（car）

> testdata<-data.frame（x=c（t1，t2，t3））

> group1=c（rep（c“（ A”“， B”“， C”），c（100，100，100））））

> leveneTest（x～group1，data=testdata）

## Levene's Test for Homogeneity of Variance （center=

median）

## Df F value Pr（>F）

## group 2 0.7591 0.469

## 297

使用小提琴圖（利用ggplot2 包）將t1，t2 和t3 組結(jié)果進(jìn)行可視化，觀察三組樣本的離散程度，對(duì)比數(shù)據(jù)的方差大小。圖2 中A、B、C分別對(duì)應(yīng)t1、t2 和t3 組，由小提琴圖可以發(fā)現(xiàn)，雖然t1 組、t2 組和t3 組數(shù)值范圍不同，但離散成果很相近，這也與使用方差齊性檢驗(yàn)得到的結(jié)果一致。

2 t檢驗(yàn)

在研究中，我們常關(guān)注組間差異的比較問題。例如，使用某種具有特定功效的化妝品和之前使用某種化妝品比較，使用效果是否有明顯改善。在樣本整體符合正態(tài)分布且具有方差齊性的條件下，獨(dú)立樣本t 檢驗(yàn)可以很好地實(shí)現(xiàn)分析目的。t 檢驗(yàn)主要分為三種，分別是單樣本t 檢驗(yàn)、配對(duì)樣本t 檢驗(yàn)、獨(dú)立樣本t 檢驗(yàn)。它們本質(zhì)上都是對(duì)比均值，但在不同的分析研究環(huán)境中應(yīng)選擇不同的t 檢驗(yàn)。單樣本t 檢驗(yàn)主要用于比較一組數(shù)據(jù)與一個(gè)特定數(shù)值之間的差異情況，一般可以用于分析整體的態(tài)度傾向。獨(dú)立樣本t 檢驗(yàn)用于分析定類數(shù)據(jù)與定量數(shù)據(jù)之間的差異情況。例如，在兩種背景下（使用和不使用化妝品），比較女性消費(fèi)者的皮膚彈性是否有明顯的差異性。通過兩組數(shù)據(jù)的對(duì)比分析，判斷該化妝品是否會(huì)影響皮膚彈性水平。以數(shù)據(jù)集flex 為例（該數(shù)據(jù)集并非研究結(jié)果，請(qǐng)自行設(shè)置數(shù)據(jù)集），程序和結(jié)果如下：

> t.test（flexindex～use，var.equal=TRUE，data=flex）

## Two Sample t-test

## data：flexindex by use

## t=2.6529，df=187，p-value=0.008667

## alternative hypothesis： true difference in means between

group 0 and group 1 is not equal to 0

## 95 percent confidence interval：

## 72.75612 494.79735

## sample estimates：

## mean in group 0 mean in group 1

## 3055.696 2771.919

參數(shù)var.equal 用于設(shè)置方差是否齊性，默認(rèn)為FALSE（不齊）。由上面方差齊性檢驗(yàn)的結(jié)論，設(shè)置為TRUE。結(jié)果表明皮膚彈性水平在使用化妝品和不使用化妝品的女性之間有顯著差異（p 值為0.008667）。

配對(duì)樣本t 檢驗(yàn)，用于分析配對(duì)定量數(shù)據(jù)之間的差異對(duì)比關(guān)系。與獨(dú)立樣本t 檢驗(yàn)相比，配對(duì)樣本t 檢驗(yàn)要求樣本是配對(duì)的。兩組樣本的樣本量要相同，樣本先后的順序要一一對(duì)應(yīng)。下面建立兩組數(shù)據(jù)，分別代表用方法A和方法B對(duì)皮膚彈性相關(guān)蛋白C進(jìn)行10 次測(cè)定，比較兩種檢測(cè)蛋白方法是否存在差異，程序與結(jié)果如下：

>x<-c（0.82，0.59，0.61，0.61，0.68，0.88，0.65，0.87，0.99，0.91）

>y<-c（0.68，0.56，0.57，0.39，0.39，0.51，0.65，0.71，0.97，0.61）

>t.test（x，y，paired=TRUE）

## Paired t-test

## data：x and y

## t=3.7119，df=9，p-value=0.004831

## alternative hypothesis： true difference in means is not

equal to 0

## 95 percent confidence interval：

## 0.06131851 0.25268149

## sample estimates：

## mean of the differences

## 0.157

結(jié)果表明，兩種化妝品有顯著差異（p 值=0.004831）。

3 方差分析及多重比較

若統(tǒng)計(jì)對(duì)象為兩組以上的獨(dú)立樣本則需要運(yùn)用方差分析。方差分析的目的在于從試驗(yàn)數(shù)據(jù)中分析出顯著影響因素，以及各個(gè)因素間的交互影響，以確定各個(gè)因素在該研究中的作用大小。功效評(píng)價(jià)研究涉及的方差分析大多為單因素方差分析和協(xié)方差分析兩種。

例如，針對(duì)某化妝品原料抗氧化能力測(cè)試結(jié)果，分析3 個(gè)研究組中抗氧化能力指數(shù)（ORAC）是否存在顯著性差異，如果有差異，說(shuō)明哪組之間的抗氧化能力指數(shù)有差異。首先，對(duì)變量的平均值進(jìn)行可視化，程序如下：

> Testdata<- read. csv（"/Users / Desktop / Testdata. csv"， header=

TRUE）

> library（gplots）

> plotmeans（ORAC～Team， Testdata， xlab= "Team"， ylab="ORAC"，main=""）

從圖3 中可以看出，t1（control）組、t2（positive）組和t3（test）組平均值之間有差異，但這種差異是否顯著，可使用單因素方差分析來(lái)驗(yàn)證。用aov（）函數(shù)進(jìn)行單因素方差分析：

> Testdataaov<-aov（ORAC～Team， Testdata）

> Summary（Testdataaov）

## Df Sum Sq Mean Sq F value Pr（>F）

## Team 2 10.98 5.489 47.36 <2e-16 ***

## Residuals 147 17.04 0.116

## ---

## Signif. codes： 0 ‘*** 0.001 ‘** 0.01 ‘* 0.05 ‘.

0.1‘ 1

由于結(jié)果中p 值遠(yuǎn)小于0.05，說(shuō)明三組ORAC平均值不完全相等。接下來(lái)分析哪些組之間的平均值相差較大，使用TukeyHSD（）函數(shù)對(duì)方差分析的結(jié)果進(jìn)行多重比較，程序如下：library（ggplot2）

tky<-TukeyHSD（Testdataaov）

> tky=as.data.frame（tky$Team）

> tky$pair=rownames（tky）

> ggplot（tky，aes（colour=cut（`p adj`，c（0，0.01，0.05，1），label=c（"p<

0.01"，"p<0.05"，"Non-Sig"））））

+theme_bw（base_size=16）

+geom_hline（yintercept=0， lty="11"，colour="grey30"，size=1）

+ geom_errorbar（aes（pair， ymin=lwr， ymax=upr）， width=0.2， size=

1）

+geom_point（aes（pair，diff），size=2）

+labs（colour=""）+theme（axis.text.x=element_text（size=14））

對(duì)可視化結(jié)果（見圖4）進(jìn)行程序驗(yàn)證：

> tukeyHSD（Testdataaov）

## Tukey multiple comparisons of means

## 95% family-wise confidence level

## Fit：aov（formula=ORAC～Team，data=Testdata）

## $Team

## diff lwr upr p adj

## Positive-Control 0.648 0.48679895 0.809201

0.0000000

## Test-Control 0.204 0.04279895 0.365201

0.0089521

## Test-Positive - 0.444 - 0.60520105 - 0.282799

0.0000000

上述計(jì)算結(jié)果顯示，在顯著性水平為0.05 的情況下，3 組ORAC值兩兩之間均有顯著差異，與可視化結(jié)果一致。

4 卡方檢驗(yàn)

t 檢驗(yàn)和方差分析主要研究定類和定量數(shù)據(jù)的差異性，卡方檢驗(yàn)則主要用于分析定類數(shù)據(jù)與定量數(shù)據(jù)之間的差異。若卡方值越大，二者偏差程度越大；反之，二者偏差越小。簡(jiǎn)言之，卡方分析是用來(lái)判斷兩個(gè)樣本間的差異程度，從而推斷兩個(gè)變量之間有沒有關(guān)系。常見的卡方分析是2×2 列聯(lián)表形式。例如，評(píng)價(jià)使用某種化妝品對(duì)皮膚彈性是否有影響。在R語(yǔ)言中，使用chisq.test（）函數(shù)即可實(shí)現(xiàn)卡方檢驗(yàn)，可以返回卡方值和對(duì)應(yīng)的p 值，同時(shí)還可以計(jì)算自由度。但是其對(duì)數(shù)據(jù)集的格式有一定的要求，如下為實(shí)際的程序代碼：

> Useful<-c（33，42）

> Useless<-c（17，8）

> data<- data. frame（Useful， Useless， row. names=c（"Control"，

"Test"））

> chisq.test（data）

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

correction

## data：data

## X-squared=3.4133，df=1，p-value=0.06467> data

## Useful Useless

## Control 33 17

## Test 42 8

t 檢驗(yàn)和方差分析都屬于參數(shù)統(tǒng)計(jì)。非參數(shù)統(tǒng)計(jì)也能計(jì)算得出結(jié)果，只因?yàn)榉菂?shù)統(tǒng)計(jì)沒有利用數(shù)據(jù)的全部信息，從而檢驗(yàn)功效往往較低。因此，在條件允許的情況下盡可能采用參數(shù)統(tǒng)計(jì)，使寶貴的原始數(shù)據(jù)得到充分利用。數(shù)據(jù)的正態(tài)分布和方差齊性在實(shí)際應(yīng)用過程中，即使不完全滿足這些前提條件，在很多時(shí)候統(tǒng)計(jì)效果也是可以被接受的，只有在嚴(yán)重違背這些前提條件導(dǎo)致統(tǒng)計(jì)結(jié)果不可信時(shí)，才采用非參數(shù)統(tǒng)計(jì)。用R語(yǔ)言能清楚、便捷、可視化地實(shí)現(xiàn)化妝品功效評(píng)價(jià)數(shù)據(jù)差異關(guān)系的統(tǒng)計(jì)分析功能，方便化妝品研究人員更好地設(shè)計(jì)研究、分析數(shù)據(jù)，迅速找出差異點(diǎn)并進(jìn)行下一步研究。