R语言实战-读书笔记 （第7章 基础统计分析）

***********************************

与导图结合的脚本文件：
创建脚本：文件——新建脚本程序，将以下代码复制粘贴至脚本内，选中右键运行当前或所选代码。

#第7章 基本统计分析#
myvars <- c("mpg", "hp", "wt")#待分析数据集：每加仑汽油行驶英里数（ mpg ）、马力（ hp ）和车重（ wt ）#
head(mtcars[myvars])
#summary()#
summary(mtcars[myvars])
#sapply()#
mystats <- function(x, na.omit=FALSE){
if (na.omit)
x <- x[!is.na(x)]
m <- mean(x)
n <- length(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))
}
myvars <- c("mpg", "hp", "wt")
sapply(mtcars[myvars], mystats)
sapply(mtcars[myvars], mystats,na.omit=TRUE)
# stat.desc()#
library(pastecs)
myvars <- c("mpg", "hp", "wt")
stat.desc(mtcars[myvars])
# describe()#
library(psych)
myvars <- c("mpg", "hp", "wt")
describe(mtcars[myvars])

#分组统计#
#aggregate()#
myvars <- c("mpg", "hp", "wt")
head(mtcars)
aggregate(mtcars[myvars], by=list(am=mtcars$am), mean)
# by=list(name1=groupvar1, name2=groupvar2, ... , nameN=groupvarN) #
aggregate(mtcars[myvars], by=list(am=mtcars$am), sd)
#describeby()#
library(psych)
myvars<-c("mpg","hp","wt")
head(mtcars[myvars])
describeBy(mtcars[myvars],list(am=mtcars$am))

#列联表#
library(vcd)
head(Arthritis)
#一维列联表#
mytable<-with(Arthritis,table(Improved))#with把所有操作都限定在此数据框上#
mytable
prop.table(mytable)*100#转换成百分比#
#二维列联表#
mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
mytable
margin.table(mytable)#求和T#
margin.table(mytable, 1) #Treatment，求和行,R1,R2#
prop.table(mytable, 1) #按照table的第一个变量的边际频数比例O11 /R1 #
margin.table(mytable, 2) #Improved， 求和列，C1,C2#
prop.table(mytable, 2) #按照table的第二个变量的边际频数比例O11 /C1 #
prop.table(mytable) #整个表格频数比例O11 /T#
addmargins(mytable)#频数求和，上表#
addmargins(prop.table(mytable))#所有变量创造边际和#
addmargins(prop.table(mytable, 1), 2)#仅添加各行边际和#
addmargins(prop.table(mytable, 2), 1)#仅添加各列边际和#

#独立性检验#
chisq.test(mytable)#p-value<0.05 存在关系，不独立#
mytable2 <- xtabs(~ Improved + Sex, data = Arthritis)
mytable2
chisq.test(mytable2)#p-value>0.05 独立#
fisher.test(mytable)
mytable3 <- xtabs(~ Treatment + Improved + Sex, data = Arthritis)
mantelhaen.test(mytable3)
#患者接受的治疗与得到的改善在性别的每一水平下并不独立#
#（分性别来看，用药治疗的患者较接受安慰剂的患者有了更多的改善）#

#相关#
states<- state.x77[,1:6]
cov(states)
cor(states)
cor(states, method="spearman")
#偏相关#
library(ggm)
colnames(states)
pcor(c(1,5,2,3,6), cov(states))
#相关检验#
cor.test(states[,3], states[,5])
#alternative hypothesis: true correlation is not equal to 0#
library(psych)
corr.test(states, use="complete")

#方差的齐性检验#
library(MASS)
head(UScrime)
str(UScrime)
#它包含了1960年美国47个州的刑罚制度对犯罪率影响的信息。
结果变量为Prob（监禁的概率）、
U1（14-24岁年龄段城市男性失业率）和U2（35-39岁年龄段城市男性失业率）。
类别型变量So（指示该州是否位于南方的指示变量）将作为分组变量使用。#
UScrime$So<-factor(UScrime$So)#设置因素变量#
bartlett.test(Prob~So,data=UScrime)
install.packages("car")
library(car)
leveneTest(Prob~So,data=UScrime)#方差齐性#
fligner.test(Prob~So,data=UScrime)
#独立样本的 t 检验#
t.test(Prob ~ So, data = UScrime,var.equal=TRUE)#具有显著性差异#
#非独立样本的 t 检验#
sapply(UScrime[c("U1", "U2")], function(x) (c(mean=mean(x), sd=sd(x))))
with(UScrime, t.test(U1, U2, paired = T))
#组间差异非参数，两组比较#
wilcox.test(Prob ~ So, data=UScrime)
with(UScrime, wilcox.test(U1, U2, paired=TRUE))
#组间差异非参数，多组比较#
head(state.x77)
#它包含了美国各州的人口、收入、文盲率、预期寿命、谋杀率和高中毕业率数据。
如果你想比较美国四个地区（东北部、南部、中北部和西部）的文盲率#
states<-data.frame(state.region,state.x77)
head(states)
kruskal.test(Illiteracy~state.region,data=states)