数据准备

df <- read.table(file = "D:/Documents/R wd/df.csv", header = T, sep = ",", colClasses = c(year = "character", nitrogen = "character", variety = "character", block = "character")) # 数据导入。
df # 查看数据。
##    year nitrogen variety block   v1   v2  v3   v4
## 1  2020       N1       a     1 1.26 2.14 3.4 4.66
## 2  2020       N1       a     2 1.20 2.90 4.1 5.30
## 3  2020       N1       a     3 1.30 3.00 4.3 5.60
## 4  2020       N1       b     1 1.08 1.72 2.8 3.88
## 5  2020       N1       b     2 1.05 1.65 2.7 3.75
## 6  2020       N1       b     3 1.15 1.35 2.5 3.65
## 7  2020       N2       a     1 1.32 3.78 5.1 6.42
## 8  2020       N2       a     2 1.28 4.32 5.6 6.88
## 9  2020       N2       a     3 1.35 3.95 5.3 6.65
## 10 2020       N2       b     1 1.33 3.47 4.8 6.13
## 11 2020       N2       b     2 1.28 2.72 4.0 5.28
## 12 2020       N2       b     3 1.30 3.90 5.2 6.50
## 13 2021       N1       a     1 1.19 3.61 4.8 5.99
## 14 2021       N1       a     2 1.21 3.29 4.5 5.71
## 15 2021       N1       a     3 1.24 3.26 4.5 5.74
## 16 2021       N1       b     1 1.09 2.71 3.8 4.89
## 17 2021       N1       b     2 1.28 2.32 3.6 4.88
## 18 2021       N1       b     3 1.35 1.95 3.3 4.65
## 19 2021       N2       a     1 1.45 4.35 5.8 7.25
## 20 2021       N2       a     2 1.40 3.80 5.2 6.60
## 21 2021       N2       a     3 1.37 4.23 5.6 6.97
## 22 2021       N2       b     1 1.28 2.72 4.0 5.28
## 23 2021       N2       b     2 1.15 3.35 4.5 5.65
## 24 2021       N2       b     3 1.24 3.46 4.7 5.94

第7章 基本统计分析

7.1 描述性统计分析

7.1.1 方法云集

对于基础安装，你可以使用summary()函数来获取描述性统计量。
summary()函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计。

summary(df[,5:8]) # summary函数计算df数据集5到8列的描述性统计量。
##        v1              v2              v3              v4       
##  Min.   :1.050   Min.   :1.350   Min.   :2.500   Min.   :3.650  
##  1st Qu.:1.198   1st Qu.:2.612   1st Qu.:3.750   1st Qu.:4.888  
##  Median :1.280   Median :3.275   Median :4.500   Median :5.680  
##  Mean   :1.256   Mean   :3.081   Mean   :4.338   Mean   :5.594  
##  3rd Qu.:1.323   3rd Qu.:3.785   3rd Qu.:5.125   3rd Qu.:6.440  
##  Max.   :1.450   Max.   :4.350   Max.   :5.800   Max.   :7.250

sapply(x, FUN, options)
x：数据框（或矩阵）；
FUN为一个任意的函数。如果指定了options，它们将被传递给FUN。

## 通过sapply函数计算描述性统计。 
mystats <- function(x, na.omit=FALSE){
  if(na.omit)
    x <- x[!is.na(x)]
   n <- length(x)
   m <- mean(x)
   s <- sd(x)
   skew <- sum((x-m)^3/s^3)/n
   kurt <- sum((x-m)^4/s^4)/n-3
   return(c(n=n, mean=m, stdev=s, skew=skew, kurtosis=kurt))
} # 定义描述性统计的函数，n为样本量，m为平均值，s为标准差，skew为偏度，kurt为峰度。
sapply(df[5:8], mystats) # sapply函数描述性统计。
##                  v1         v2         v3         v4
## n        24.0000000 24.0000000 24.0000000 24.0000000
## mean      1.2562500  3.0812500  4.3375000  5.5937500
## stdev     0.1018017  0.8771957  0.9430812  1.0149183
## skew     -0.2991286 -0.3609173 -0.3371765 -0.3209242
## kurtosis -0.6567661 -1.0459573 -0.9674555 -0.8823974

扩展学习-偏度和峰度
1 偏度
偏度（Skewness）可以用来度量随机变量概率分布的不对称性。公式如下：

其中是均值，是标准差。
偏度的取值范围为(-∞,+∞)
当偏度<0时，概率分布图左偏。
当偏度=0时，表示数据相对均匀的分布在平均值两侧，不一定是绝对的对称分布。
当偏度>0时，概率分布图右偏。

2 峰度
峰度（Kurtosis）可以用来度量随机变量概率分布的陡峭程度。公式：

其中是均值，是标准差。
峰度的取值范围为[1,+∞)，完全服从正态分布的数据的峰度值为3，峰度值越大，概率分布图越高尖，峰度值越小，越矮胖。
通常我们将峰度值减去3，也被称为超值峰度（Excess Kurtosis），这样正态分布的峰度值等于0，当峰度值>0，则表示该数据分布与正态分布相比较为高尖，当峰度值<0，则表示该数据分布与正态分布相比较为矮胖。

扩展

Hmisc包中的describe()函数可返回变量和观测的数量、缺失值和唯一值的数目、平均值、分位数，以及五个最大的值和五个最小的值。

library(Hmisc) # 调用Hmisc包，调用前请安装。
describe(df$v1) # 显示结果。
## df$v1 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##       24        0       17    0.994    1.256   0.1169    1.082    1.108 
##      .25      .50      .75      .90      .95 
##    1.197    1.280    1.323    1.364    1.395 
## 
## lowest : 1.05 1.08 1.09 1.15 1.19, highest: 1.33 1.35 1.37 1.40 1.45
##                                                                             
## Value       1.05  1.08  1.09  1.15  1.19  1.20  1.21  1.24  1.26  1.28  1.30
## Frequency      1     1     1     2     1     1     1     2     1     4     2
## Proportion 0.042 0.042 0.042 0.083 0.042 0.042 0.042 0.083 0.042 0.167 0.083
##                                               
## Value       1.32  1.33  1.35  1.37  1.40  1.45
## Frequency      1     1     2     1     1     1
## Proportion 0.042 0.042 0.083 0.042 0.042 0.042

pastecs包中有一个名为stat.desc()的函数，它可以计算种类繁多的描述性统计量。 stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
x是一个数据框或时间序列。若basic=TRUE（默认值），则计算其中所有值、空值、缺失值的数量，以及最小值、最大值、值域，还有总和。若desc=TRUE（同样也是默认值），则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。最后，若norm=TRUE（不是默认的），则返回正态分布统计量，包括偏度和峰度（以及它们的统计显著程度）和Shapiro--Wilk正态检验结果。这里使用了p值来计算平均数的置信区间（默认置信度为0.95）。

library(pastecs) # 调用pastecs包。
stat.desc(df$v1, norm=TRUE, p=0.95) # 返回结果。
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  24.00000000   0.00000000   0.00000000   1.05000000   1.45000000   0.40000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  30.15000000   1.28000000   1.25625000   0.02078019   0.04298709   0.01036359 
##      std.dev     coef.var     skewness     skew.2SE     kurtosis     kurt.2SE 
##   0.10180170   0.08103618  -0.29912856  -0.31669845  -0.65676606  -0.35780260 
##   normtest.W   normtest.p 
##   0.97472260   0.78255589

psych包也拥有一个名为describe()的函数，它可以计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误。

library(psych) # 调用psych包。
describe(df$v1) # 返回结果。
##    vars  n mean  sd median trimmed mad  min  max range skew kurtosis   se
## X1    1 24 1.26 0.1   1.28    1.26 0.1 1.05 1.45   0.4 -0.3    -0.66 0.02

7.1.2 分组计算描述性统计量

aggregate(x, by=list, FUN) 如果有多个分组变量，使用by=list(name1=groupvar1, name2=groupvar2, ... ,groupvarN)这样的语句。

aggregate(df$v1, by=list(df$year, df$nitrogen, df$variety), FUN=mean) # 分组统计。
##   Group.1 Group.2 Group.3        x
## 1    2020      N1       a 1.253333
## 2    2021      N1       a 1.213333
## 3    2020      N2       a 1.316667
## 4    2021      N2       a 1.406667
## 5    2020      N1       b 1.093333
## 6    2021      N1       b 1.240000
## 7    2020      N2       b 1.303333
## 8    2021      N2       b 1.223333

by(data, INDICES, FUN)
data：一个数据框或矩阵；
INDICES：一个因子或因子组成的列表，定义了分组；
FUN是任意函数。

myfun <- function(x) {
  c(mean = mean(x), sd = sd(x))
} # 定义函数myfun。
by(df$v1, df$year, myfun) # 对df数据的v1以year进行平均值和标准差的统计。
## df$year: 2020
##       mean         sd 
## 1.24166667 0.09943781 
## ------------------------------------------------------------ 
## df$year: 2021
##      mean        sd 
## 1.2708333 0.1063834
by(df[,5:8], df[,1:3], summary) # 对df数据的5到8列进行描述性统计，以前3列分分组依据。
## year: 2020
## nitrogen: N1
## variety: a
##        v1              v2             v3              v4       
##  Min.   :1.200   Min.   :2.14   Min.   :3.400   Min.   :4.660  
##  1st Qu.:1.230   1st Qu.:2.52   1st Qu.:3.750   1st Qu.:4.980  
##  Median :1.260   Median :2.90   Median :4.100   Median :5.300  
##  Mean   :1.253   Mean   :2.68   Mean   :3.933   Mean   :5.187  
##  3rd Qu.:1.280   3rd Qu.:2.95   3rd Qu.:4.200   3rd Qu.:5.450  
##  Max.   :1.300   Max.   :3.00   Max.   :4.300   Max.   :5.600  
## ------------------------------------------------------------ 
## year: 2021
## nitrogen: N1
## variety: a
##        v1              v2              v3             v4       
##  Min.   :1.190   Min.   :3.260   Min.   :4.50   Min.   :5.710  
##  1st Qu.:1.200   1st Qu.:3.275   1st Qu.:4.50   1st Qu.:5.725  
##  Median :1.210   Median :3.290   Median :4.50   Median :5.740  
##  Mean   :1.213   Mean   :3.387   Mean   :4.60   Mean   :5.813  
##  3rd Qu.:1.225   3rd Qu.:3.450   3rd Qu.:4.65   3rd Qu.:5.865  
##  Max.   :1.240   Max.   :3.610   Max.   :4.80   Max.   :5.990  
## ------------------------------------------------------------ 
## year: 2020
## nitrogen: N2
## variety: a
##        v1              v2              v3              v4       
##  Min.   :1.280   Min.   :3.780   Min.   :5.100   Min.   :6.420  
##  1st Qu.:1.300   1st Qu.:3.865   1st Qu.:5.200   1st Qu.:6.535  
##  Median :1.320   Median :3.950   Median :5.300   Median :6.650  
##  Mean   :1.317   Mean   :4.017   Mean   :5.333   Mean   :6.650  
##  3rd Qu.:1.335   3rd Qu.:4.135   3rd Qu.:5.450   3rd Qu.:6.765  
##  Max.   :1.350   Max.   :4.320   Max.   :5.600   Max.   :6.880  
## ------------------------------------------------------------ 
## year: 2021
## nitrogen: N2
## variety: a
##        v1              v2              v3              v4       
##  Min.   :1.370   Min.   :3.800   Min.   :5.200   Min.   :6.600  
##  1st Qu.:1.385   1st Qu.:4.015   1st Qu.:5.400   1st Qu.:6.785  
##  Median :1.400   Median :4.230   Median :5.600   Median :6.970  
##  Mean   :1.407   Mean   :4.127   Mean   :5.533   Mean   :6.940  
##  3rd Qu.:1.425   3rd Qu.:4.290   3rd Qu.:5.700   3rd Qu.:7.110  
##  Max.   :1.450   Max.   :4.350   Max.   :5.800   Max.   :7.250  
## ------------------------------------------------------------ 
## year: 2020
## nitrogen: N1
## variety: b
##        v1              v2              v3              v4       
##  Min.   :1.050   Min.   :1.350   Min.   :2.500   Min.   :3.650  
##  1st Qu.:1.065   1st Qu.:1.500   1st Qu.:2.600   1st Qu.:3.700  
##  Median :1.080   Median :1.650   Median :2.700   Median :3.750  
##  Mean   :1.093   Mean   :1.573   Mean   :2.667   Mean   :3.760  
##  3rd Qu.:1.115   3rd Qu.:1.685   3rd Qu.:2.750   3rd Qu.:3.815  
##  Max.   :1.150   Max.   :1.720   Max.   :2.800   Max.   :3.880  
## ------------------------------------------------------------ 
## year: 2021
## nitrogen: N1
## variety: b
##        v1              v2              v3              v4       
##  Min.   :1.090   Min.   :1.950   Min.   :3.300   Min.   :4.650  
##  1st Qu.:1.185   1st Qu.:2.135   1st Qu.:3.450   1st Qu.:4.765  
##  Median :1.280   Median :2.320   Median :3.600   Median :4.880  
##  Mean   :1.240   Mean   :2.327   Mean   :3.567   Mean   :4.807  
##  3rd Qu.:1.315   3rd Qu.:2.515   3rd Qu.:3.700   3rd Qu.:4.885  
##  Max.   :1.350   Max.   :2.710   Max.   :3.800   Max.   :4.890  
## ------------------------------------------------------------ 
## year: 2020
## nitrogen: N2
## variety: b
##        v1              v2              v3              v4       
##  Min.   :1.280   Min.   :2.720   Min.   :4.000   Min.   :5.280  
##  1st Qu.:1.290   1st Qu.:3.095   1st Qu.:4.400   1st Qu.:5.705  
##  Median :1.300   Median :3.470   Median :4.800   Median :6.130  
##  Mean   :1.303   Mean   :3.363   Mean   :4.667   Mean   :5.970  
##  3rd Qu.:1.315   3rd Qu.:3.685   3rd Qu.:5.000   3rd Qu.:6.315  
##  Max.   :1.330   Max.   :3.900   Max.   :5.200   Max.   :6.500  
## ------------------------------------------------------------ 
## year: 2021
## nitrogen: N2
## variety: b
##        v1              v2              v3             v4       
##  Min.   :1.150   Min.   :2.720   Min.   :4.00   Min.   :5.280  
##  1st Qu.:1.195   1st Qu.:3.035   1st Qu.:4.25   1st Qu.:5.465  
##  Median :1.240   Median :3.350   Median :4.50   Median :5.650  
##  Mean   :1.223   Mean   :3.177   Mean   :4.40   Mean   :5.623  
##  3rd Qu.:1.260   3rd Qu.:3.405   3rd Qu.:4.60   3rd Qu.:5.795  
##  Max.   :1.280   Max.   :3.460   Max.   :4.70   Max.   :5.940

扩展

doBy包中summaryBy()函数分组计算描述性统计量。
summaryBy(formula, data=dataframe, FUN=function)
formula可接受的格式var1 + var2 + var3 +....+ varN ~ groupvar1 + groupvar2 + groupvar3 +....+ groupvarN.
在~左侧的变量是需要分析的数值型变量，而右侧的变量是类别型的分组变量。function可为任何内建或用户自编的R函数。

myfun <- function(x) {
  c(mean = mean(x), sd = sd(x))
} # 定义函数myfun。
library(doBy) # 调用doBy包。
summaryBy(v1 + v2 ~ year + nitrogen + variety, data = df, FUN = myfun) # 返回结果。
##   year nitrogen variety  v1.mean      v1.sd  v2.mean     v2.sd
## 1 2020       N1       a 1.253333 0.05033223 2.680000 0.4703190
## 2 2020       N1       b 1.093333 0.05131601 1.573333 0.1965536
## 3 2020       N2       a 1.316667 0.03511885 4.016667 0.2761038
## 4 2020       N2       b 1.303333 0.02516611 3.363333 0.5971879
## 5 2021       N1       a 1.213333 0.02516611 3.386667 0.1939931
## 6 2021       N1       b 1.240000 0.13453624 2.326667 0.3800439
## 7 2021       N2       a 1.406667 0.04041452 4.126667 0.2891943
## 8 2021       N2       b 1.223333 0.06658328 3.176667 0.3992910

psych包中的describe.by()函数分组计算概述统计量。 一个以上的分组变量时可以list(groupvar1, groupvar2, ... , groupvarN)

library(psych) # 调用psych包。
describe.by(df, df$year) # 返回结果，计算df数据集5到8列的描述性统计量，分组变量是year。
## 
##  Descriptive statistics by group 
## group: 2020
##           vars  n mean   sd median trimmed  mad  min  max range  skew kurtosis
## year*        1 12 1.00 0.00   1.00    1.00 0.00 1.00 1.00  0.00   NaN      NaN
## nitrogen*    2 12 1.50 0.52   1.50    1.50 0.74 1.00 2.00  1.00  0.00    -2.16
## variety*     3 12 1.50 0.52   1.50    1.50 0.74 1.00 2.00  1.00  0.00    -2.16
## block*       4 12 2.00 0.85   2.00    2.00 1.48 1.00 3.00  2.00  0.00    -1.74
## v1           5 12 1.24 0.10   1.28    1.25 0.07 1.05 1.35  0.30 -0.78    -0.98
## v2           6 12 2.91 1.01   2.95    2.92 1.32 1.35 4.32  2.97 -0.18    -1.59
## v3           7 12 4.15 1.09   4.20    4.17 1.41 2.50 5.60  3.10 -0.23    -1.60
## v4           8 12 5.39 1.18   5.45    5.42 1.50 3.65 6.88  3.23 -0.28    -1.60
##             se
## year*     0.00
## nitrogen* 0.15
## variety*  0.15
## block*    0.25
## v1        0.03
## v2        0.29
## v3        0.32
## v4        0.34
## ------------------------------------------------------------ 
## group: 2021
##           vars  n mean   sd median trimmed  mad  min  max range  skew kurtosis
## year*        1 12 1.00 0.00   1.00    1.00 0.00 1.00 1.00  0.00   NaN      NaN
## nitrogen*    2 12 1.50 0.52   1.50    1.50 0.74 1.00 2.00  1.00  0.00    -2.16
## variety*     3 12 1.50 0.52   1.50    1.50 0.74 1.00 2.00  1.00  0.00    -2.16
## block*       4 12 2.00 0.85   2.00    2.00 1.48 1.00 3.00  2.00  0.00    -1.74
## v1           5 12 1.27 0.11   1.26    1.27 0.12 1.09 1.45  0.36  0.06    -1.19
## v2           6 12 3.25 0.73   3.32    3.28 0.80 1.95 4.35  2.40 -0.19    -1.11
## v3           7 12 4.53 0.77   4.50    4.52 0.89 3.30 5.80  2.50  0.09    -1.20
## v4           8 12 5.80 0.82   5.72    5.76 0.95 4.65 7.25  2.60  0.31    -1.20
##             se
## year*     0.00
## nitrogen* 0.15
## variety*  0.15
## block*    0.25
## v1        0.03
## v2        0.21
## v3        0.22
## v4        0.24
describe.by(df$v1, list(df$year, df$nitrogen, df$variety)) # 返回结果，计算df数据v1的描述性统计量，分组变量是year,nitrogen,variety。
## 
##  Descriptive statistics by group 
## : 2020
## : N1
## : a
##    vars n mean   sd median trimmed  mad min max range  skew kurtosis   se
## X1    1 3 1.25 0.05   1.26    1.25 0.06 1.2 1.3   0.1 -0.13    -2.33 0.03
## ------------------------------------------------------------ 
## : 2021
## : N1
## : a
##    vars n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 3 1.21 0.03   1.21    1.21 0.03 1.19 1.24  0.05 0.13    -2.33 0.01
## ------------------------------------------------------------ 
## : 2020
## : N2
## : a
##    vars n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 3 1.32 0.04   1.32    1.32 0.04 1.28 1.35  0.07 -0.09    -2.33 0.02
## ------------------------------------------------------------ 
## : 2021
## : N2
## : a
##    vars n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 3 1.41 0.04    1.4    1.41 0.04 1.37 1.45  0.08 0.16    -2.33 0.02
## ------------------------------------------------------------ 
## : 2020
## : N1
## : b
##    vars n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 3 1.09 0.05   1.08    1.09 0.04 1.05 1.15   0.1 0.24    -2.33 0.03
## ------------------------------------------------------------ 
## : 2021
## : N1
## : b
##    vars n mean   sd median trimmed mad  min  max range  skew kurtosis   se
## X1    1 3 1.24 0.13   1.28    1.24 0.1 1.09 1.35  0.26 -0.27    -2.33 0.08
## ------------------------------------------------------------ 
## : 2020
## : N2
## : b
##    vars n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 3  1.3 0.03    1.3     1.3 0.03 1.28 1.33  0.05 0.13    -2.33 0.01
## ------------------------------------------------------------ 
## : 2021
## : N2
## : b
##    vars n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 3 1.22 0.07   1.24    1.22 0.06 1.15 1.28  0.13 -0.23    -2.33 0.04

7.1.3 结果的可视化

df1 <- aggregate(df$v1, by = list(df$year, df$nitrogen, df$variety),mean) # 分组汇总df数据框v1列。
df1 # 显示结果
##   Group.1 Group.2 Group.3        x
## 1    2020      N1       a 1.253333
## 2    2021      N1       a 1.213333
## 3    2020      N2       a 1.316667
## 4    2021      N2       a 1.406667
## 5    2020      N1       b 1.093333
## 6    2021      N1       b 1.240000
## 7    2020      N2       b 1.303333
## 8    2021      N2       b 1.223333

plot(density(df1$x)) # 简单显示峰度和偏度。

image.png

boxplot(df$v1) # 单变量箱线图。

image.png

boxplot(df$v1 ~ df$year) # 按分组显示变量箱线图。

image.png

参考资料：

1. 《R语言实战》（中文版），人民邮电出版社，2013.
2. 偏度和峰度，https://zhuanlan.zhihu.com/p/84614017