以下内容为阅读医学统计学笔记与R语言所做的笔记。而原文又为孙振球，徐勇勇老师的<<医学统计学>> 第4版的一份笔记，及用R语言的实现。目前还未能找到孙、徐二位老师教材的电子版作为比对，因此需要原文与本人的看法存在出入的地方均以张厚粲<现代心理与教育统计>，r中语法有出入的地方则以<r语言实战>为准。

2 第二章 计量资料的统计描述

2.1 描述统计

• 检测潜在的异常值，检验数据质量

• 助于"理解"数据(了解数据的特点)

• 常用的集中量数(位置度量)：平均值，中位数，四分位数，第三、四分位数，众数，最大值，最小值等

• 常用的离散量数(离散度量)：全距，标准差，方差，四分位间距，变异系数。

例题：

1. 查看全距

#read data，na.strings="NA"means all missing data will change into "NA" as strings
rbc <- read.table("tst.txt",sep="\t", na.strings="NA")
rbc <- as.matrix(rbc)
rbc_q <- c()
#遍历rbc中的列，首尾相接存入rbc中。
for (i in 1:nrow(rbc)){
  rbc_q <- c(rbc_q, rbc[i,])
}
#查看变量类型
typeof(rbc_q)
##[1] "list"

rbc_v <- as.vector(rbc_q)
#omit the NA
rbc_v <- na.omit(rbc_v)
#evaluate range of data
rge<-max(rbc_v)-min(rbc_v)
#range()返回两个极值而非全距
rge
## 2.39

2. 计算频数分布

• 确定组段数
  h=3.49min(s,IQ/1.349) * n^1/3
  s为标准差
  IQ为四分位数

#计算标准差
s<-sd(rbc_v)
## [1] 0.4457298
#计算分位数( 25% 、75% 两个分位数)

iq<-as.numeric(quantile(rbc_v,0.75)-quantile(rbc_v,0.25))
## 0.565

#确定组数
h<-3.49*min(s,iq/1.349)*(length(rbc_v)^(1/3))
## 7.553617

#向上取整
h<-ceiling(h)
## 8

#每组的距离
i<-rge/h

• 频数分布

#设置分组区间左右两侧的值(上下限)
breaks = seq(min(rbc_v), max(rbc_v), length.out = 8)

#按照breaks将rbc切开，计算每个区间内的频次
rbc_v.cut = cut(rbc_v, breaks, right=T,include.lowest=T)

#频次统计
rbc_v.freq = table(rbc_v.cut)


#直方图
## y值为从0到1的频率直方图
hist(rbc_v, right=FALSE, 
      breaks = breaks, labels =TRUE, 
      freq = FALSE, col = "#A8D6FF", 
      border = "white", ylim=c(0,1))
lines(density(rbc_v),col="red",lwd=2)

2.4 描述统计

1. 集中量数

#算术平均值
mean(x)

#几何平均值
exp(mean(log(x)))

几何均值：

几何平均值公式

首先两遍取对数

ln(G) = ln(X1×X2×...×Xn) /n = [ln(X1) + ln(X2) + ... + ln(Xn)] / n

lnxⁿ=nlnx
ln(x^-n)=ln(x)/n

即求完原变量的自然对数的均值后( mean(log(x)) )，再计算以自然常数e为底，ln(G)为幂的值即可，这在R中可以轻松实现。

参考自biye5u.com. R语言中计算几何平均数

#中位数
median(x)
quantile(x)

2. 离散趋势

#极差(全距)
range(x)
max(x)-min(x)

#四分位间距
IQR(x)
quantile(RBC_v, 0.75)-quantile(RBC_v, 0.25)

变异系数(百分之)
sd(x)/mean(x)*100

3. 其他描述统计

summary(x)
#numeric
#min, max, median, mean, 1st-4th quantiles
#factor/charact
#frequency

#众数
mode(x) #为查看数据类型，类似typeof(x)

#可以用table()和sort()来寻找众数
rbc_t <- table(rbc_v) # 计算每个元素的出现的次数
sort(rbc_t, decreasing = TRUE)  # 对计算的次数进行排序，最前边就是众数

#或者结合 which()函数确定众数和其次数
rbc_t <- table(rbc_v)
rbc_t[which(((rbc_t==max(rbc_t))==T))]#选出rbc中最大的数

m==max(m)：m中每个数是不是和最大的数相等，返回为T/F
要首先对rbc进行频次统计，否则就会对取出原数据中最大的数而非频次最大的数。

2.5 正态分布

1. 概率密度函数

概率密度函数probability density。可以求出指定均值/标准差下每个点的概率密度。越高就代表着这个点/区间的概率越密集(大)

#在-10~10区间等分的 100个 数据集x
x <- seq(-10, 10, by = .1)
#创建一个均值是2.5，标准差是0.5正态分布 y
y <- dnorm(x, mean = 2.5, sd = 0.5)
#将 y 中的落在x数据集上的数据画出来
plot(x,y,col="red",pch=20)

2. 概率累计分布函数

累积分布函数Cumulative Distribution ，给出一个正态分布中小于一个给定数字的累计概率

#在-10~10区间等分的 40个 数据集x
x <- seq(-10, 10, by = .5)
#创建一个均值是2.5，标准差是0.5正态分布 y
y <- pnorm(x, mean = 2.5, sd = 0.5)
#将 y 中的落在x数据集上的累计概率画出来
plot(x,y,col="red",pch=20)

为小于该点的概率密度的相加。

3. 正态分位数函数

qnorm()给出某个点在指定均值和标准差下的分位数

和pnorm()相反，pnorm()
为返回该分位数的概率。

4. 生成正态分布函数

rnorm(n，mean，sd)

#设置随机种子，便于重复后续的数据选取
set.seed(50)
#在标准正态分布中随机选取50个数据
y <- rnorm(50)
#对选区的数据绘制频率分布图
hist(y,col="#A8D6FF",labels =TRUE)

为什么设定种子：weixin_30595035。R 语言设定随机数种子。CSDN。

2.6 正态分布检验

方法：

• 计算综合统计量：动差法(峰度和偏度，三级动差)、虾皮罗-威尔克S-W法、达格斯提诺(D'Agostino)D检验、shapiro-francia W'检验

• 正态分布的拟合优度检验：皮尔逊精确chi-square检验、对数似然比检验、科莫格罗夫-谢米洛夫K-S检验

• 图示：PP图(百分位图)、QQ图(分位数图)、稳定化概率图(SP图)

规则：

SPSS :样本含量3 ≤n ≤5000 时,结果以Shapiro - Wilk (W 检验) 为难,当样本含量n > 5000 结果以Kolmogorov - Smirnov 为准。

SAS 规定:样本含量n ≤2000 时,结果以Shapiro - Wilk (W 检验) 为准,当样本含量n >2000 时,结果以Kolmogorov - Smirnov (D 检验) 为准。

正态检验

函数：

ks.test(x,pnorm(0,1))
#检验x与均值为0，标准差为1的正态分布之间是否有差异。

shaprio.test(x)

#nortest包
lillie.test()#可以实行更精确的Kolmogorov-Smirnov检验。
ad.test()#进行Anderson-Darling正态性检验。
cvm.test()#进行Cramer-von Mises正态性检验。
pearson.test()#进行Pearson卡方正态性检验。
sf.test()#进行Shapiro-Francia正态性检验。

#fBasics包
normalTest()#进行Kolmogorov-Smirnov正态性检验。
ksnormTest()  #进行Kolmogorov-Smirnov正态性检验。
shapiroTest()  #进行Shapiro-Wilk's正态检验。
jarqueberaTest()  #进行jarque-Bera正态性检验。
dagoTest()  #进行D'Agostino正态性检验。
gofnorm()  #采用13种方法进行检验，并输出结果。##并不存在gofnorm这个函数名。

1. 直方图

#确定分组数：
IQ<-quantile(x,0.75)-quantile(x, 0.25)
h=3.49min(s,IQ/1.349)n^1/3^

#直方图
hist(x, right=FALSE, 
      breaks = h, labels =TRUE, 
      freq = FALSE, col = "#A8D6FF", 
      border = "white", ylim=c(0,1))

#正态曲线
lines(density(x),col="red",lwd=2)

2. 概率密度图

直方图的平滑版

#图名
maintxt<-paste("N=",length(x),",","Mean=",round(mean(x),3),",","Sd=",round(s,3))
#作图
plot(density(x),col="red",lwd=2,main = maintxt)

3. QQ-plot

• 只需要确定数据点是否沿着直线(有时也称为Henry’s line)来判断数据是否符合正态。如果点靠近参考线并且在置信区间内，则认为满足了正态性假设。

• 如果qq图所示的非正态分布(系统地偏离参考线)时，通常第一步是对数据进行对数变换， 并重新检查对数变换后的数据是否正态分布。可以应用log()函数进行对数变换。

#载入程集包car
library(car)
#设置种子
set.seed(42)

#随机正态分布的qq图
dat_hist <- data.frame( value = rnorm(12, mean = 165, sd = 5))
qqPlot(dat_hist$value)

library(car)
qqPlot(as.numeric(rbc_v),ylab="RBC", main="RBC QQ-plot")

4. shaprio-wilk检验与k-s检验

shapiro.test(x)

随着观测数的增加，Shapiro-Wilk检验变得非常敏感，甚至对正态性的一个 微小偏差也非常敏感。

样本量较大(n>5000时)使用ks检验效度更高。

第三章 总体均值的估计与假设检验

3.1 推论统计

样本量较小(n<30)时t检验比正态分布要效果更好， 大样本上与正态分布无差异。当自由度df=∞时，t分布曲线为标准正态分布曲线。

1. 假设检验

假设检验：利用反证法，从事情多反面出发(H₀)，通过其发生的概率(p-value)，间接判断要解决的问题(H₁)是否成立。

2. 统计量与标准误

样本：从总体中抽取出来的，能够反映总体特征的少量数据。

标准误：从总体中抽取不同的样本，不同样本之间的均值的标准差、不同样本之间的标准差的标准差这种样本统计量的标准差就叫标准误。反映了抽样统计量的离散程度。(除此之外还有样本均值的均值等，用处暂时不大)

3.2 t分布

当总体形态未知，样本量较大时，可以用均值标准误的无差估计量S_n-1(通常分母为n-1)作为总体标准差。

t分布

t分布的概率密度函数

自由度：自由取值的个数。如 X+Y+Z=1 ，有3个变量，但是能够自由取值的自由两个(这两个可以随意取值，为了使X+Y+Z=1成立，第三个数必须为一个固定的数值)，故其自由度 v=2。自由度的取值一般为：v=n-m。其中n为计算某一统计量是用到的数据个数，m为计算该统计量是用到的其他独立统计量个数。

t分布是一簇曲线，其形态变化与n（即其自由度）大小有关。自由度n越小，t分布曲线越低平；自由度n越大，t分布曲线越接近标准正态分布（u分布） 曲线，当自由度无限大时，t分布就成了正态分布。

1. 概率密度函数dt()

dt(x, df, log=F)

#正态分布曲线，红色
curve(dnorm(x),xlim=c(-5,5),ylim=c(0,0.45),ylab="Student's t Density",col="red",lty=1,lwd=3)
abline(v=0,lwd=1,col="black")#添加中轴线
#df=1时的t分布曲线，绿色
curve(dt(x,1),col="green",lty=2,lwd=3,add=TRUE)
curve(dt(x,2),col="brown",lty=3,add=TRUE)
#df=5时的t分布曲线，蓝色
curve(dt(x,5),col="blue",lty=4,lwd=3,add=TRUE)
curve(dt(x,20),col="dark green",lty=5,add=TRUE)
legend(2,0.38,c("normal","df=1","df=2","df=5","df=20"),
       col=c("red","green","brown","blue","dark green"),lty=1:5)

正态分布与t 分布

当样本量足够大时(df足够大)，比如df=30时，

#正态分布，红色粗线
curve(dnorm(x),xlim=c(-5,5),ylim=c(0,0.45),ylab="Student's t Density",col="red",lty=1,lwd=10)
#df=30,蓝色的细虚线
curve(dt(x,30),col="blue",lty=6,add=TRUE)

样本量较大时的的t分布

可见此时两者已经接近完全重合。

2. 概率累积分布函数pt()

# 单侧 p值，df=5
# P(t => 1.9)t值大于等于1.9的概率
pt(q=1.9, df=5, lower.tail = F)
## [1] 0.05793165

## 双侧 p-value
## 两边对称单侧相加
pt(q=1.9, df=5, lower.tail = F) + pt(q=-1.9, df=5, lower.tail = T)
## 0.1158633

## 对称性：单侧p值*2=双侧p值
pt(q=1.9, df=5, lower.tail = F)*2
## 0.1158633

3. 求置信区间的qt()

# 双尾的95%置信区间的t值(α=0.05，在取值时除以2)
qt(p=0.025, df = 5, lower.tail = T)
## -2.570582

4. 随机生成符合t分布的数据rt()

# 设置种子
set.seed(61)

n <- 10000
# 生成1000个符合df=5的t分布的数据
y_rt <- rt(n, df = 5)
# 做直方图
hist(y_rt, breaks = 100, main = "Randomly t density", freq=FALSE,
col="#A8D6FF",xlim=c(-10,10),ylim=c(0,0.4))
lines(density(y_rt), col="red", lwd=2)

t分布的直方图

5. 参数估计

用已知的样本参数推断总体参数。

• σ未知且n较小(n≤60)时，按照t分布计算。
• σ未知且n较大(n＞60)时，按照标准正态分布(即u分布或Z分布)计算。

6. 假设检验

具体步骤不再赘述。

3.3 t检验的分类与应用

1. t检验与Z检验

t检验基本上有三种类型：

• 单样本t检验(One sample/gourp t-test)
• 配对样本t检验(Paired/matched t-test)
• 两种独立样本t检验(Two sample/gourp t-test)

R的提供了t.test()，可以用于各种t检验。

如果只提供一组数据，则进行单样本t检验，如果提供两组数据， 则进行两样本t检验，主要的相关的参数是：

• paired：TRUE-配对t检验，FALSE-两独立样本t检验。
• var.equal，对于两独立样本t检验，还要注意方差是否齐性，TRUE则进行经典方法，FALSE则采用近似t检验， 将数据取对数处理后进行t检验，t.test()中是采用Welch t检验。
• alternative参数，用来指定检验方式，默认是双侧检验(“two.sided”)，还可以是左侧检验(“less”)和右侧检验(“greater”)。

2. 单样本t检验

t.test(data, mu= , alterntive="two sides"/"lesser"/"greater")

3. 配对t检验

setwd("C:/users/lenovo/desktop")
m <- read.table("1.txt",header=F)
home_df <- as.data.frame(m)
t.test(m$V1,mu=140)
##  t = -2.1367, df = 35, p-value = 0.03969

#左侧检验
t.test(home_df$hb, mu=140, alternative = "less")
#右侧检验
t.test(home_df$hb, mu=140, alternative = "greater")

直接计算出t值后求p值

#直接计算出t值后，使用pt函数计算p值的方式
t.value <- mean(home_df$hb - 140) /
              sd(home_df$hb) * sqrt(nrow(home_df))
#-2.136668
#双侧检验，概率值乘以2
p.value <- pt(t.value, 
              df=nrow(home_df)-1, 
              lower.tail=T)*2
## [1] 0.03969288

4. 两独立样本t检验

两样本含量较小(n1≤60,或(和)n2≤60)，且各自的总体符合正态分布正态分布检验时，再根据两组数据的方差是否一样(方差齐性检验)来采用不同t检验方法。

• 方差齐性检验
  • 使用var.test(x，y)进行方差齐性检验
  • α值一般为0.1

4.1 方差齐性时

#读取数据
bfg_sav<-spss.system.file("ExampleData/SavData4MedSta/Exam03-07.sav")

#将数据形式转换为数据框
bfg_df<-as.data.frame(as.data.set(bfg_sav))
x <- bfg_df$result[which(bfg_df$group==1)]#group为1的数据
y <- bfg_df$result[which(bfg_df$group==2)]

#直接使用t.test()进行计算，注意 var.equal=TRUE
t.test(x,y,var.equal=TRUE)

4.2 方差不齐时

t.teet()中用的是Welch法

而R中还未找到可以实现Cochran&Cox法 和 Satterthwaite法 的函数，但是可以通过自己敲代码实现。如下

cochran's t法

#Cochran&Cox法的检验统计量*t'*的
cochran_test<-function(x, y) {
#计算均值，方差，个案数
mean_x<- mean(x)
mean_y<- mean(y)
dev_x<-var(x)
dev_y<-var(y)
nx<-length(x)
ny<-length(y)

#计算t'值
  t_value <- (mean_x-mean_y)/sqrt(dev_x/nx+dev_y/ny)

#计算P值
  p_value <- pt(cochran_t,df=nx-1,lower.tail=F)*2

#合成列表与命名
  op<-list(t_value, p_value)
  names(op) <- c("cochran' t", "p")

  return(op)
}

satterhwaite法来矫正cochran法的自由度

#Satterthwaite法来矫正Cochran&Cox法的检验统计量*t'*的自由度
sattw_df<-function(x, y){
dev_x<-var(x)
dev_y<-var(y)
nx<-length(x)
ny<-length(y)

  satt_df<-(dev_x/nx+dev_y/ny)^2/
  (((dev_x/nx)^2/(nx-1))+
  ((dev_y/ny)^2/(ny-1)))

  return(satt_df)
}

satt_df<-sattw_df(stdev_x,stdev_y,nx,ny)
#这里注意的是cochran_t>0，我们要计算cochran_t相关的双侧拒绝域面积，需要指定 lower.tail=F，并加倍
pt(cochran_t,satt_df,lower.tail=F)*2
#[1] 0.3427644

3.4 变量转换

常用的变量变换有：

• 对数变换(Logarithm)
• 平方根变换(Square root)
• 平方根反正弦变换(Arcsine, Arcsine-square)
• 倒数变换(Reciprocal)
• 指数变换(Exponential)

数据变换的方法

具体操作

例子：

#install.packages("openintro","e1071","pryr")
#使用mammals数据集测试
library(openintro)
#e1071包中的skewness()计算偏斜度
library(e1071)
#使用pryr包中的绘图对象保存方案
library(pryr)
library("dplyr")
skewness(mammals$body_wt)
#使用log()对body_wt数据进行转换
loged_bd_wt <- log(mammals$body_wt)
#skewness()计算偏斜度
skewness(loged_bd_wt)
##[1] 0.1453599

p1.pryr %<a-% {
hist(mammals$body_wt,breaks = 200,freq=FALSE,
 xlim=c(min(mammals$body_wt),max(mammals$body_wt)),main="Bodt wt. 数据分布")
lines(density(mammals$body_wt), col="red", lwd=1)
}

p2.pryr %<a-% {
  hist(loged_bd_wt,breaks = 200,freq=FALSE,
 xlim=c(min(loged_bd_wt),max(loged_bd_wt)),main="经过Log转换后Bodt wt. 数据分布")
lines(density(loged_bd_wt), col="red", lwd=1)
}

split.screen(c(1, 2))
screen(1)
p1.pryr
screen(2)
p2.pryr

下面是两组样本数据的对数转换示例，在mammalsbrain_wt数据：

#使用log()对brain_wt数据进行转换
loged_br_wt <- log(mammals$brain_wt)
p3.pryr %<a-% {
  plot(mammals$body_wt,mammals$brain_wt,main="body_wt vs brain_wt",col="red",pch=16)
}
p4.pryr %<a-% {
  plot(loged_bd_wt,loged_br_wt,main="body_wt vs brain_wt after both loged",col="red",pch=16)
}
split.screen(c(1, 2))
screen(1)
p3.pryr
screen(2)
p4.pryr

第四章 多个样本均数比较的方差分析

4.1 F分布

4.2 F检验的分类、原理及应用

1. 分类
   按方差分析的基本运算概念，依照所感兴趣的因子数量而可分三大类：

• 单因子方差分析(One-way ANOVA)
• 双因子方差分析(Two-way ANOVA)
• 多因子方差分析(Multi-way ANOVA)

如果依照因子的特性不同而有三种型态，分别是：

• 固定效应方差分析(Fixed-effect ANOVA)
• 随机效应方差分析(Random-effect ANOVA)
• 混合效应方差分析(Mixed-effect ANOVA)

一般的，按照因子数量的分类进行区分使用情况比较常见，这里借用单因素方差分析(One-way ANOVA)介绍方差分析的基本思想。

2. 原理

关于方差分析的原理，可以参考舒华<教育与心理研究中的多因素实验设计>一书。这里不再赘述

顺便说一句，这书还没有看完……柑!

4.3 方差分析的诊断图

4.4 多个样本均数间的多重比较

4.5 多样本方差齐性检验