R語言實戰讀書筆記(七)基本統計分析

summary()

sapply(x,fun,options):對數據框或矩陣中的每個向量進行統計

mean

sd：標準差

var：方差

min:

max:

median:

length:

range:

quantile:

vars <- c("mpg", "hp", "wt")
head(mtcars[vars])

summary(mtcars[vars])

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))
}

sapply(mtcars[vars], mystats)

#如結果，分佈右偏(0.61)，較正態分佈稍平(-0.37)

#Hmisc包中的describe可返回變量和觀測的數量，缺失值和惟一值的數目、平均值、分位數及五個最大的值和最小的值

library(Hmisc)
describe(mtcars[vars])

#pastects包中有一個stat.desc函數，能夠計算各類描述性統計量

library(pastecs)
stat.desc(mtcars[vars])

#psych包的describe能夠計算描述性統計量

library(psych)
describe(mtcars[vars])

#aggregate()：分組獲取描述性統計量

aggregate(mtcars[vars], by = list(am = mtcars$am), mean)
aggregate(mtcars[vars], by = list(am = mtcars$am), sd)

#by(data,indices,fun)

dstats <- function(x)(c(mean=mean(x), sd=sd(x)))
by(mtcars[vars], mtcars$am, dstats)

#doBy的summaryBy能夠分組計算描述性統計量

library(doBy)
summaryBy(mpg + hp + wt ~ am, data = mtcars, FUN = mystats)

#psych中的desribe.by能夠分組計算描述性統計量，但不容許指定任何函數。

library(psych)
describe.by(mtcars[vars], mtcars$am)

#使用reshape包

library(reshape)
dstats <- function(x) (c(n = length(x), mean = mean(x), sd = sd(x)))

#id.vars，按什麼分組，measure看什麼變量？
dfm <- melt(mtcars, measure.vars = c("mpg", "hp", "wt"), id.vars = c("am", "cyl"))
cast(dfm, am + cyl + variable ~ ., dstats)

table：使用N個類別型變量建立一個N維列聯表

xtabs：根據一個公式和一個矩陣或數據框建立一個N維列聯表

prop.table：將表中條目表示爲分數形式

margin.table：計算表中條目的和

addmargins：將結果放入表中

ftable：建立一個緊湊的平鋪式列聯表

#一維列聯表

mytable <- with(Arthritis, table(Improved))
mytable
prop.table(mytable)

#二維列聯表

mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
mytable

#計算邊緣頻數，1表示第一個變量
margin.table(mytable, 1)

#計算邊緣比例，1表示第一個變量
prop.table(mytable, 1)
margin.table(mytable, 2)
prop.table(mytable, 2)
prop.table(mytable)

#每行和每列最後加sum
addmargins(mytable)
admargins(prop.table(mytable))

#僅添加各行的和
addmargins(prop.table(mytable, 1), 2)

#僅添加各列的和
addmargins(prop.table(mytable, 2, 1)

#能夠用gmodels的CrossTable生成二維列聯表

library(gmodels)
CrossTable(Arthritis$Treatment, Arthritis$Improved)

#獨立性檢驗

library(vcd)
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
chisq.test(mytable)

#以上代碼結果p值爲0.001，存在某種關係
mytable <- xtabs(~Improved+Sex, data=Arthritis)
chisq.test(mytable)

#以上代碼結果p值爲0.08，不要緊

#如下代碼結果p值爲0.001，有關係

mytable <- xtabs(~Treatment+Improved, data=Arthritis)
fisher.test(mytable)

#如下代碼p值爲0.0006，假設不存在三階交互做用

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

#如下爲相關性度量，值越大相關性越強

library(vcd)
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
assocstats(mytable)

#cor函數：計算Pearson相關係數，Spearman相關係數，Kendall相關係數

#cov函數：計算協方差

states <- state.x77[, 1:6]

#如下代碼計算方差和協方差
cov(states)

#如下代碼計算Pearson積差相關係數
cor(states)

#如下代碼計算Spearman等級相關係數
cor(states, method="spearman")

#能夠看出收入和高中畢業率之間存在強正相關，文盲率和預期壽命之間存在很強的負相關

#設定行列，能夠得出更直觀的結果

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

#相關性檢驗

cor.test(states[, 3], states[, 5])

#一次檢驗多個

library(psych)
corr.test(states, use = "complete")

#獨立樣本t檢驗

library(MASS)
t.test(Prob ~ So, data=UScrime).

#非獨立樣本t檢驗

library(MASS)
sapply(UScrime[c("U1", "U2")], function(x) (c(mean = mean(x),
sd = sd(x))))
with(UScrime, t.test(U1, U2, paired = TRUE))

#組間差別的非參數檢驗，wilcox.test適用於兩組成對數據和沒法保證正態性假設的情境

with(UScrime, by(Prob, So, median))

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

#paired設爲TRUE，配對檢驗

sapply(UScrime[c("U1", "U2")], median)
with(UScrime, wilcox.test(U1, U2, paired = TRUE))

#多於兩組的檢驗

states <- as.data.frame(cbind(state.region, state.x77))kruskal.test(Illiteracy ~ state.region, data=states)