#描述性統計分析

myvars <- c("mpg","hp","wt")

head(mtcars[myvars])

# summary函數
# 數值型變量，提供min,max,四分位數（1st QU，Median,3rd Qu）,mean
# 因子向量，邏輯向量，頻數統計

summary(mtcars[myvars])

#sapply函數

#自定義函數,對數值，剔除空值，計算平均數，方差，偏度，峰度，并返回
mystats <- function(x,na.omit=FALSE){
? ? if(na.omit)
? ? ? ? x <- x[!is.na(x)] ?#剔除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))
}
#通過sapply調用自定義函數
myvars <- c("mpg","hp","wt")
sapply(mtcars[myvars], mystats)

#Hmisc包,describe函數
#3 variables 3個變量(3列)，32個觀測值（行）
#獨立變量：mpg,hp,wt
#每個變量下：
#n: 表示總行數
#missiong: 缺失值個數
#unique:非重復值個數
#mean 均值
#百分位數：0.05,0.10，0.2,0.5,0.75,0.9.0.95
#最小的五個值,最大的五個值

library("lattice")
library("survival")
library("Formula")
library("ggplot2")

library("Hmisc") ?#首先完成關聯依賴的加載

myvars <- c("mpg","hp","wt")
describe(mtcars[myvars])

#pastecs包,stat.desc函數

#stat.desc(x,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95)
#當basic=True保持默認值：
#nbr.val ? ? ? ? ? ?所有值的數量
#nbr.null ? ? ? ? ? 空值的數量
#nbr.na ? ? ? ? ? ? 缺失值的數量
#min,max ? ? ? ? ? ?最小值，最大值
#range ? ? ? ? ? ? ?值域
#sum ? ? ? ? ? ? ? ?求和

#當desc=TRUE保持默認值：
#median ? ? ? ? ? ? 中位數
#mean ? ? ? ? ? ? ? 平均數
#SE.mean ? ? ? ? ? ?平均數的標準誤
#CI.mean.0.95 ? ? ? 平均數置信度為95%的置信區間
#var ? ? ? ? ? ? ? ?方差
#std.dev ? ? ? ? ? ?標準差
#coef.var ? ? ? ? ? 變異系數

#當norm=TRUE不保持默認值，默認值為FALSE
#skewness ? ? ? ? ? 偏度
#skew.2SE ? ? ? ? ? 偏度顯著程度
#kurtosis ? ? ? ? ? 峰度
#kurt.2SE ? ? ? ? ? 峰度顯著程度
#normtest.W ? ? ? ? 正態檢驗結果
#normtest.p ? ? ? ? 正態檢驗結果?

#p=0.95 ? ? ? ? ? ? 置信度為0.95

#范例一
install.packages("pastecs")
library(pastecs)

myvars <- c("mpg","hp","wt")
stat.desc(mtcars[myvars],norm = TRUE)

#psych包，describe函數

library("psych")
#加載psych之后，提示Hmisc包中的describe函數被psych包中的describe函數覆蓋；
#最后載入的程序包優先
#采用Hmisc::describe()

myvars <- c("mpg","hp","wt")

#返回參數介紹
# n ? ? ? 非缺失值的數量
# mean ? ?平均數
# sd ? ? ?標準差
# median ?中位數
# trimmed 截尾均值
# mad ? ? 絕對中位差
# min ? ? 最小值
# max ? ? 最大值
# range ? 值域
# skew ? ?偏度
# kurtosis 峰度
# se ? ? ?平均值的標準誤

describe(mtcars[myvars])

#分組計算描述性統計量
myvars <- c("mpg","hp","wt")

#基于mtcars$am分組，對mtcars[mpg,hp,wt]三列，分別求均值
#am 對應分組列，0,1均為組值
#am=mtcars$am,假如不賦值列名，將會被標注為Group.1
aggregate(mtcars[myvars],by=list(am=mtcars$am),mean)

#基于mtcars$am分組，對mtcars[mpg,hp,wt]三列，分別求標準差
#am 對應分組列，0,1均為組值
#aggregate() 僅允許在每次調用中使用平均數，標準差這樣的單返回值函數，無法一次返回若干個統計量
aggregate(mtcars[myvars],by=list(am=mtcars$am),sd)

#使用by分組計算描述性統計量
dstats <- function(x)sapply(x,mystats) #嵌套函數定義

myvars <- c("mpg","hp","wt")

by(mtcars[myvars],mtcars$am,dstats) #規避aggregate僅能返回單個返回值問題

#doby包中的summaryBy()分組計算概率統計量
library(doBy)

#summaryBy(formula,data=dataframe,FUN=function)
#var1+var2+var3~groupvar1+groupvar2+...
#~左側的變量是需要分析的數值型變量
#~右側的變量是類別性的分組變量
summaryBy(mpg+hp+wt~am,data=mtcars,FUN = mystats)
? ??
#psych包中describeBy
library(psych)
myvars <- c("mpg","hp","wt")

#基于分組，分別統計各列的各種統計指標
#統計指標：
#n,mean,sd,median,trimmed,mad,min,max,range,skew,kurtosis,se
#不允許指定任意函數，普適性較低
describeBy(mtcars[myvars],list(am=mtcars$am))

#針對類別性變量的圖表

#數據源
library("vcd")
head(Arthritis)

#一維列聯表

#table函數默認忽略缺失值NA,如頻數統計時需要視為一個有效的類別，請useNA="ifany"
mytable <- with(Arthritis,
? ? ? ? ? ? ? ? table(Improved)
? ? ? ? ? ? ? ? ) #table 函數,頻數統計

prop.table(mytable) #基于頻數統計，轉成成比例值,0-1間的小數形式表示
prop.table(mytable)*100 #基于頻數統計，轉成成比例值,1-100間的整數形式表示

#二維列聯表

#mytable <- xtabs(~A+B,data=mydata)
#mydata 是一個矩陣或者數據框
#交叉分類的變量應出現在公式的右側，即~符號的右方，以+作為分隔符
#A的值 行變量 第一個變量
#B的值 列變量 第二個變量
mytable <- xtabs(~Treatment+Improved,data=Arthritis) #生成二維頻數表格

margin.table(mytable,1) #基于二維頻數統計，二次加工，按照行，生成邊際頻數（按行求和）
margin.table(mytable,2) #基于二維頻數統計，二次加工，按照列，生成邊際頻數（按列求和）

prop.table(mytable,1)#基于二維頻數統計，二次加工，按照行，生成當前頻數/邊際頻數（按行求和）
prop.table(mytable,2)#基于二維頻數統計，二次加工，按照列，生成當前頻數/邊際頻數（按列求和）

prop.table(mytable)#基于二維頻數表格，不指定行維度，列維度時，所有元素得到的比例總和為1

addmargins(mytable) #基于二維頻數表格，追加一列，保存基于行求和的值；追加一行，保存基于列求和的值

addmargins(prop.table(mytable))#基于百分比表格，追加一列，保存基于行求和的%值；追加一行，保存基于列求和的%值
addmargins(prop.table(mytable,1),2)#基于行百分比的表格，追加一列，保存行的和值
addmargins(prop.table(mytable,2),1)#基于列百分比的表格，追加一行，保存列的和值

#使用CrossTable生成二維列聯表
library(gmodels)

#結果解釋(不同隔斷中，相同行位置，或者相同列位置值相加)

#第一行
#整數部分剝離
#相當于addmargins(mytable)),返回二維頻數邊際統計表
#行方向：29+7+7=43 13+7+21=41
#列方向：29+13=42,7+7=14,7+21=28
#42+14+28=43+41=84

#第二行
#Chi-square contribution

#第三行
#相當于addmargins(prob.table(mytable,1),2)
#行方向: 0.674+0.163+0.163=1,0.317+0.171+0.512=1

#第四行
#相當于addmargins(prob.table(mytable,2),1)
#列方向: 0.690+0.310=1,0.500+0.500=1,0.250+0.750=1

#第五行
#相當于addmargins(prob.table(mytable))
#行方向：0.345+0.083+0.083=0.512,0.155+0.083+0.250=0.488
#列方向: 0.345+0.155=0.500,0.083+0.083=0.167,0.083+0.250=0.333
#0.512+0.488=0.500+0.167+0.333=1

CrossTable(Arthritis$Treatment,Arthritis$Improved)?

#多維列聯表

#原生成二維列聯表的函數，直接推廣到多維

#生成三維數據源，頻數表
#Treatment 可理解成原始行維度，滿足塊條件下的統計
#sex ? ? ? 可理解成原始列維度，滿足塊條件下的統計
#Improved ?只能是第三維度，塊標識
mytable <- xtabs(~Treatment+Sex+Improved,data=Arthritis)

ftable(mytable) #區別于xtabs的多維表格輸出格式，這個更緊湊些

#margin.table 單維度
margin.table(mytable,1) #根據變量1(行變量)， 多個塊計算邊際值，橫向求和
margin.table(mytable,2) #根據變量2(列變量)， 多個塊計算邊際值，縱向求和
margin.table(mytable,3) #根據變量3(塊變量)， 單塊計算邊際值，塊內求和

#margin.table 超過單維度
margin.table(mytable,c(1,3)) #相當于忽略第二個維度，僅根據第1，第3的維度的條件，統計頻數

ftable(prop.table(mytable,c(1,2))) #第三個維度仍為塊維度，比較難，沒太懂
ftable(addmargins(prop.table(mytable,c(1,2)),3)) #結合上面一行理解

ftable(addmargins(prop.table(mytable,c(1,2)),3))*100 #結合上一行理解

#獨立性檢驗
#顯著性檢驗評估了是否存在充分的證據以拒絕變量間相互獨立的原假設

#卡方獨立性檢驗
library(vcd)
mytable <- xtabs(~Treatment+Improved,data=Arthritis)?

#卡方假設檢驗
#用于檢驗原假設

#范例一
#p-value p值，p值落入小概率范圍(p<0.01)則拒絕原假設?
#df 自由度
chisq.test(mytable)

#范例二
#p-value p值，p值未落入小概率范圍(p>0.01)則接受原假設，兩者相互獨立?
mytable <- xtabs(~Treatment+Sex,data=Arthritis)
chisq.test(mytable)

#fisher 精確檢驗
mytable <- xtabs(~Treatment+Improved,data=Arthritis)
fisher.test(mytable) #可以用于任意行數大于等于2的二維列聯表上使用，但是2行，2列不行

#Cochran-Mantel-Haenszel檢驗
#其原假設：兩個名義變量在第三個變量的每一層中都是條件獨立的

#范例一
mytable <- xtabs(~Treatment+Improved+Sex,data=Arthritis)
mantelhaen.test(mytable)

#相關性的度量
#如果拒絕了相互獨立的原假設，開始轉向衡量相關性強弱和相關性度量

#二維列聯表的相關性度量
#phi-coefficient: ?fi相關系數
#Contigency Coeff: 列聯相關系數
#Cramer's V ? ? ?: V相關系數
#較大值意味著較強的相關性
library(vcd)
mytable <- xtabs (~Treatment+Improved,data=Arthritis)
assocstats(mytable)

#結果的可視化

#相關
#相關系數可以用來描述定量變量之間的關系。相關系數的符號表明關系的方向（正相關或負相關），
#其值的大小表示關系的強弱程度(完全不相關時0，完全相關時為1)

#Pearson，積差相關系數，衡量了兩個定量變量之間的線性相關程度
#spearman，等級相關系數則衡量了分級定序變量之間的相關程度
#Kendall's Tau相關系數也是一種非參數的等級相關度量

#協方差
states <- state.x77[,1:6]
cov(states)

cor(states) #協方差，默認pearson系數
cor(states,method="spearman") #spearman相關系數

#相關系數矩陣默認為一個方陣，所有變量之間兩兩計算相關，也可以計算非方形的相關矩陣

x <- states[,c("Population","Income","Illiteracy","HS Grad")]
y <- states[,c("Life Exp","Murder")]
cor(x,y)

#偏相關
#控制一個或多個定量變量時，另外兩個定量變量之間的相互關系

library(ggm)
colnames(states) #獲取列名

#pcor(u,S),偏相關系數常用于社會科學的研究中
#u是一個數值向量，前兩個數值表示要計算相關系數的變量下標，其余的數值為
#條件變量(即要排除影響的變量)，S為變量的協方差陣
pcor(c(1,5,2,3,6),cov(states)) #控制其他變量，求得Poputation與Murder之間的相關關系
#返回0.346

#相關性的顯著性檢驗

#cor.test(x,y,alternative=,method=),每次只能檢驗一種相關關系
#x,y 要檢驗相關性的變量
#alternative 則用來指定進行雙側檢驗或者單側檢驗，
#greater,當研究的總體的相關系數大于0時
#less，當研究的總體的相關系數小于0時
#two.side 總體相關系數不等于0

#原假設：檢驗了預期壽命和磨砂率的Pearson相關系數為0的原假設
# p-value = 1.258e-08,拒絕原假設，即預期壽命和磨砂率之間的總體相關度不為0
cor.test(states[,3],states[,5])

#psych包 corr.test函數
library(psych)

#Correlation matrix ?相關矩陣，對稱陣
#Probability values ?顯著性檢驗
corr.test(states,use="complete")

#獨立樣本的t檢驗
#一個針對兩組的獨立樣本t檢驗可以用于檢驗兩個總體的均值相等的假設
t.test(y~x,data)

#范例一

library(MASS)

#原假設：南方各州和非南方各州擁有相同監禁概率的假設
#p-value = 0.0006506，p<0.01,落入拒絕區域
#正態化變換，Y/1-Y,log(Y/1-Y),arcsin(Y),arcsin(sqrt(Y))
t.test(Prob~So,data=UScrime)

#非獨立樣本的t檢驗
t.test(y1,y2,paired = TRUE)

library(MASS)
sapply(UScrime[c("U1","U2")],function(x)(c(mean=mean(x),sd=sd(x))))

#p-value < 2.2e-16，落入拒絕域
#可以拒絕年長和年輕男性的平均失業率相同的假設
with(UScrime,
? ? ?t.test(U1,U2,paired=TRUE)
? ? ?)

#如果想在多于兩個的組之間進行比較，假設數據是從正態總體中獨立抽樣而得的，那你
#可以使用方差分析(ANOVA),如果數據無法滿足t檢驗或者ANOVA的參數檢驗，可以轉而使用
#非參數方法

#wilcoxon秩和檢驗，是非獨立樣本t檢驗的一種非參數替代方法
#試用于兩組成對數據和無法保證正態性假設的情景

#范例一
with(UScrime,
? ? ?by(Prob,So,median)
? ? )

wilcox.test(Prob~So,data=UScrime)

#多于兩組的比較

#kruskal-wallis檢驗
#kruskal.test(y~A,data) y是一個數值型結果變量，A是一個擁有兩個或更多水平的分組變量

#范例一

states <- data.frame(state.region,state.x77)

#p-value =4.726e-05,落入拒絕域，拒絕文盲率相同的假設
kruskal.test(Illiteracy~state.region,data=states)

#范例二 InternetOpenUrl failed: '無法與服務器建立連接'
source("http://www.statmethods.net/RiA/wmc.txt")
#執行失敗， InternetOpenUrl failed: '無法與服務器建立連接'
?

總結

以上是生活随笔為你收集整理的读书笔记_第七章的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。