化妝品功效評價統(tǒng)計分析中的差異關(guān)系及R語言實現(xiàn)

王福鑫 趙志龍 矯筱蔓

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

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

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

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

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

如果樣本數(shù)在50 以下選擇S-W檢驗，樣本數(shù)在50 以上選擇K-S 檢驗。Q-Q圖可最直觀地顯示樣本數(shù)據(jù)是否符合正態(tài)性分布。例如，有100 名受試者參加功效評價試驗，擬對其進(jìn)行數(shù)據(jù)正態(tài)性檢驗，用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)一步檢驗數(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 來看，說明nordata 不服從標(biāo)準(zhǔn)正態(tài)分布（隨機(jī)數(shù)設(shè)置為平均值為2，sd 為5）。

方差齊性檢驗是檢驗不同樣本的總體方差是否相同的一種方法。常用的檢驗方法有F 檢驗、Bartlet 檢驗和Levene 檢驗。建議首選使用Levene 進(jìn)行樣本齊性檢驗，因其更穩(wěn)健且不依賴總體分布情況，還可用于多組樣本。例如，對3 個實驗組受試者進(jìn)行方差齊性檢驗（t1 對照組、t2 陽性組和t3 試驗組均以隨機(jī)數(shù)為例做樣本），使用car 包中的leveneTest（）函數(shù)進(jìn)行Levene檢驗，程序和結(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)行可視化，觀察三組樣本的離散程度，對比數(shù)據(jù)的方差大小。圖2 中A、B、C分別對應(yīng)t1、t2 和t3 組，由小提琴圖可以發(fā)現(xiàn)，雖然t1 組、t2 組和t3 組數(shù)值范圍不同，但離散成果很相近，這也與使用方差齊性檢驗得到的結(jié)果一致。

2 t檢驗

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

配對樣本t 檢驗，用于分析配對定量數(shù)據(jù)之間的差異對比關(guān)系。與獨立樣本t 檢驗相比，配對樣本t 檢驗要求樣本是配對的。兩組樣本的樣本量要相同，樣本先后的順序要一一對應(yīng)。下面建立兩組數(shù)據(jù)，分別代表用方法A和方法B對皮膚彈性相關(guān)蛋白C進(jìn)行10 次測定，比較兩種檢測蛋白方法是否存在差異，程序與結(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)計對象為兩組以上的獨立樣本則需要運用方差分析。方差分析的目的在于從試驗數(shù)據(jù)中分析出顯著影響因素，以及各個因素間的交互影響，以確定各個因素在該研究中的作用大小。功效評價研究涉及的方差分析大多為單因素方差分析和協(xié)方差分析兩種。

例如，針對某化妝品原料抗氧化能力測試結(jié)果，分析3 個研究組中抗氧化能力指數(shù)（ORAC）是否存在顯著性差異，如果有差異，說明哪組之間的抗氧化能力指數(shù)有差異。首先，對變量的平均值進(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）組平均值之間有差異，但這種差異是否顯著，可使用單因素方差分析來驗證。用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，說明三組ORAC平均值不完全相等。接下來分析哪些組之間的平均值相差較大，使用TukeyHSD（）函數(shù)對方差分析的結(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））

對可視化結(jié)果（見圖4）進(jì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

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

4 卡方檢驗

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

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