如何計算描述統(tǒng)計量

1. summary()函數(shù)能颁。

最基本的就是summary函數(shù)杂瘸。可以輸出 最小值劲装，最大值胧沫，四分位數(shù)，數(shù)值型變量的均值占业，以及因子向量和邏輯型向量的頻數(shù)統(tǒng)計绒怨。
以下用mtcars數(shù)據(jù)集來舉例

head(mtcars)
str(mtcars)
myvars <- c("mpg", "hp", "wt")
summary(mtcars[myvars])

2. sapply()函數(shù)。

使用sapply函數(shù)調用任意函數(shù)來實現(xiàn)谦疾。
使用方法：sapply(x, FUN, options)
定義一個mystats函數(shù)南蹂，后通過sapply調用。

mystats <- function(x, na.omit=FALSE){
  if (na.omit) 
    x <- x[!is.na(x)] #去掉缺失值
  m <- mean(x) #均值
  n <- length(x) #個數(shù)
  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來調用自己定義好的函數(shù)mystats
sapply(mtcars[myvars], mystats)

3. pastecs包中的stat.desc()函數(shù)念恍。推薦使用這個六剥。

stat.desc()這個函數(shù)是比較好用的計算描述統(tǒng)計量的函數(shù)，可以輸出的結果比較多峰伙。
用法：stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
basic=TRUE疗疟，計算的基礎指標包括：數(shù)據(jù)個數(shù)(nbr.val)、空值的個數(shù)(nbr.null)瞳氓、缺失值的個數(shù)(nbr.na)策彤、最小值(Min)、最大值(Max)匣摘、值域(range店诗，即max-min)、所有非缺失值的總和(Sum)
desc=TRUE音榜，計算：中位數(shù)(median)庞瘸、平均數(shù)平均數(shù)(Mean)、平均數(shù)的標準誤(SE.Mean)赠叼、平均數(shù)置信度為95%的置信區(qū)間(CI.Mean)擦囊、方差(Var)、標準差(std.dev)和變異系數(shù)(coef.var,其定義為標準差除以平均數(shù))嘴办。
平均值置信度（95%）指的是在95%的置信度下計算出的平均值的允許誤差,可以用平均值+或-這個數(shù)來計算置信區(qū)間的上限和下限.
norm=TRUE霜第，計算偏度和峰度（以及他們的統(tǒng)計顯著程度），s-w正態(tài)檢驗的結果户辞。偏度系數(shù)skewness,其顯著性判別依據(jù)(如果skew.2SE>1泌类，則偏度顯著不等于零)、峰度系數(shù)kurtosis、其顯著性依據(jù)(與skew.2SE相同)刃榨、正態(tài)性的Shapiro-Wilk檢驗(Normaltest.W)及其概率(normtest.p)

install.packages("pastecs")
library(pastecs)
stat.desc(mtcars[myvars],basic=TRUE, desc=TRUE, norm=TRUE, p=0.95)

如何分組計算描述統(tǒng)計量

1. aggregate()函數(shù)弹砚。
   使用這個函數(shù)，可以計算平均數(shù)枢希、標準差等返回值桌吃。但是每次只能計算一種。

data <-aggregate(df$decision.RESP, by=list(Subject=df$Subject,Intention=df$intention,Outcome=df$outcome), FUN=mean, na.rm=TRUE)
aggregate(mtcars[myvars], by=list(am=mtcars$am), mean)
aggregate(mtcars[myvars], by=list(am=mtcars$am), sd)

doBy包中的summaryBy()函數(shù) 推薦使用這個苞轿。和stat.desc搭配使用茅诱。

可以同時返回若干個統(tǒng)計量。
用法：summaryBy(formula, data=dataframe, FUN=function)
formula的格式：var1 + var2 + var3 + ... + varN ~ groupvar1 + groupvar2 + ... + groupvarN
～左側的變量是需要分析的數(shù)值型變量搬卒，右側的變量是類別型的分組變量瑟俭。

function可以是任何內建或者用戶自編的R函數(shù)。

install.packages("doBy")
library(doBy)
# 根據(jù)am的不同水平來 分別計算 am不同水平下 mpg hp wt的描述統(tǒng)計值契邀。調用的函數(shù)是stat.desc摆寄。
# 這樣可以輸出stat.desc()函數(shù)的很多描述統(tǒng)計值。
summaryBy(mpg+hp+wt~am, data=mtcars, FUN=stat.desc)

直方圖

par(mfrow=c(2,2)) # 畫4幅圖

par(mfrow=c(1,2))實現(xiàn)一頁多圖的功能
其中坯门，通過設定函數(shù)par()的各個參數(shù)來調整我們的圖形
mfrow=c(2,2) 就是畫4幅圖微饥，
mfrow=c(3,5),是畫15幅圖，
例如 par(mfrow=c(2,3)) 一個圖版顯示2行古戴，3列
作者：巴拉巴拉_9515
鏈接：http://www.reibang.com/p/e5f86fd8d033
來源：簡書
著作權歸作者所有欠橘。商業(yè)轉載請聯(lián)系作者獲得授權，非商業(yè)轉載請注明出處现恼。

hist(mtcars$mpg)

hist(mtcars$mpg, 
     breaks=12,# 指定組的數(shù)量简软，也就是直方圖中柱形的數(shù)量。
     # col="red", # 指定顏色
     xlab="Miles Per Gallon", #指定x軸的標題
     main="Colored histogram with 12 bins") #指定圖形的標題

hist(mtcars$mpg,
     freq=FALSE, # 根據(jù)概率密度而不是頻數(shù)繪制圖形
     breaks=12,
     # col="red",
     xlab="Miles Per Gallon",
     main="Histogram, rug plot, density curve")
rug(jitter(mtcars$mpg,amount=0.01)) 
#軸須圖(rug plot)：實際數(shù)據(jù)值的一種一維呈現(xiàn)方式述暂。
# jitter 將向每個數(shù)據(jù)點添加一個小的隨機值（一個±amount之間的均勻分布隨機數(shù)），以避免重疊的點產生影響建炫。
lines(density(mtcars$mpg), col="blue", lwd=2) # 添加一條藍色畦韭、寬度2的 密度曲線。

x <- mtcars$mpg 
h<-hist(x, 
        breaks=12,
        # col="red",
        xlab="Miles Per Gallon",
        main="Histogram with normal curve and box") 
xfit<-seq(min(x), max(x), length=40) 
# dnorm是密度函數(shù)肛跌，也就是正態(tài)分布曲線函數(shù) 
# dnorm(x, mean = 0, sd = 1, log = FALSE)
yfit<-dnorm(xfit, mean=mean(x), sd=sd(x))
# h$mids[1:2] 含義是 h這一個list里面的mids變量的第一個和第二個數(shù)據(jù)艺配，這里表示的是，橫坐標每個柱形的中間點的坐標衍慎。
# diff(h$mids[1:2]) 第一個和第二個坐標的差值
yfit <- yfit*diff(h$mids[1:2])*length(x) 
lines(xfit, yfit, col="blue", lwd=2)
box() #將圖形圍繞起來的盒型

h[["mids"]]

箱線圖

用法：boxplot(formula, data=dataframe)

par(mfrow=c(1,1))
boxplot(mtcars$mpg, main="Box plot", ylab="Miles per Gallon")
a<-boxplot.stats(mtcars$mpg)

a[["stats"]] 畫圖所需的數(shù)據(jù)：最小值转唉、下四分位、中位數(shù)稳捆、上四分位赠法、最大值。
a[["n"]] 數(shù)據(jù)的個數(shù)
a[["out"]] 離群點乔夯，是 范圍±1.5IQR以外的值砖织。IQR表示四分位距款侵，即上四分位數(shù)與下四分位數(shù)的差值。
默認情況下侧纯，箱線圖兩條須的延伸極限不會超過 盒型各端±1.5IQR的范圍新锈。

boxplot(mpg~cyl, #y~a，將為類別型變量a的每個水平并列的生成數(shù)值型變量y的箱線圖眶熬。
        # 如果有兩個變量妹笆，則寫為 y~a*b 表示為變量a和b所有水平的兩兩組合生成數(shù)值型變量y的箱線圖。
        data=mtcars,
        main="Car Mileage Data",
        xlab="Number of Cylinders",
        ylab="Miles Per Gallon")

mtcars[c("cyl","am")]<-lapply(mtcars[c("cyl","am")],as.factor)
boxplot(mpg~cyl*am,
        data=mtcars,
        main="Car Mileage Data",
        xlab="Number of Cylinders",
        ylab="Miles Per Gallon")

boxplot(mpg~am*cyl,
        # 畫圖的順序是娜氏，cyl的一個水平上拳缠，am的兩個水平都畫完。在畫cyl的下一個水平
        data=mtcars,
        col=c("gold","darkgreen"),
        main="MPG Distribution by Auto Type",
        xlab="Auto Type",
        ylab="Miles Per Gallon")

小提琴圖

是箱線圖和核密度圖的結合體牍白。
用法：vioplot(x1, x2, ... , names=, col=)

install.packages("vioplot")
library(vioplot)
x1 <- mtcars$mpg[mtcars$cyl==4]
x2 <- mtcars$mpg[mtcars$cyl==6]
x3 <- mtcars$mpg[mtcars$cyl==8]
vioplot(x1, x2, x3,
        names=c("4 cyl", "6 cyl", "8 cyl"),
        col="gold")
# 圖形中白點是中位數(shù)脊凰，黑色盒型的范圍是上下四分位點，細黑線表示須茂腥。外部形狀即為核密度估計狸涌。
title("Violin Plots of Miles Per Gallon", 
      ylab="Miles Per Gallon",
      xlab="Number of Cylinders")

點圖

用法：dotchart(x, labels=)

dotchart(mtcars$mpg,
         labels=row.names(mtcars),# labels是由每個點的標簽組成的向量
         cex=.7, # 指定標簽的大小，包括labels最岗，main和xlab
         main="Gas Mileage for Car Models",
         xlab="Miles Per Gallon")

x <- mtcars[order(mtcars$mpg),]
x$cyl <- factor(x$cyl)
x$color[x$cyl==4] <- "red"
x$color[x$cyl==6] <- "blue"
x$color[x$cyl==8] <- "darkgreen"
dotchart(x$mpg, 
         labels = row.names(x), 
         cex=.7, 
         groups = x$cyl,# 通過groups來選定一個因子帕胆，用以指定x中元素的分組方式。
         gcolor = "black",# 控制不同組 標簽的顏色般渡。這里是指 4 6 8三個數(shù)字的顏色懒豹。
         color = x$color,# 點和標簽的顏色來自向量color
         pch=19, #在R語言中，點的樣式由pch的取值決定驯用。 
         main = "Gas Mileage for Car Models\ngrouped by cylinder", 
         xlab = "Miles Per Gallon")

頻數(shù)表和列聯(lián)表

install.packages("vcd")
library(vcd)
head(Arthritis)

table()

table()函數(shù)默認忽略缺失值
使用N個類別型變量（因子）創(chuàng)建一個N維列聯(lián)表
用法：table(var1, var2, ..., varN)

一維列聯(lián)表

mytable<-table(Arthritis$Improved) 
mytable
mytable <- with(Arthritis,
                table(Improved)
                )
mytable

prop.table()將頻數(shù)轉化為比例

prop.table(mytable) #prop.table()將頻數(shù)轉化為比例
prop.table(mytable)*100 #prop.table()將頻數(shù)轉化為比例脸秽，之后乘以100，變成百分比

二維列聯(lián)表

mytable <- table(A, B)
A是行變量蝴乔，B是列變量

mytable<-table(Arthritis$Improved,Arthritis$Treatment) 
x<-prop.table(mytable) 
# 不加數(shù)字的話记餐，分母是 Improved和Treatment加起來的總人數(shù)。
# 得到的所有比例數(shù)字加起來等于1薇正。也就是sum(x)等于1
x<-prop.table(mytable,1) 
# 1表示根據(jù)table()語句中的第1個變量(這里指Arthritis$Improved)來計算比例：
# 分別計算 Improved的三個水平的比例片酝，分母是 Improved的三個水平的總數(shù)。
# 每一行的比例加起來等于1挖腰。apply(x,1,sum)等于1,1,1
x<-prop.table(mytable,2) 
# 2表示根據(jù)table()語句中的第2個變量(這里指Arthritis$Treatment)來計算比例：
# 分別計算 Treatment的兩個水平的比例雕沿，分母是 Treatment 的兩個水平的總數(shù)。
# 每一列的比例加起來等于1猴仑。apply(x,2,sum)等于1,1

mytable <- xtabs(~ A + B, data=mydata)

A是行變量审轮，B是列變量

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

margin.table 生成頻數(shù)

margin.table(mytable, 1)
# 1表示根據(jù)table()語句中的第1個變量(這里指Arthritis$Treatment)來計算頻數(shù)
# 也就是計算 Treatment 的兩個水平各自的頻數(shù)。
margin.table(mytable, 2)
# 1表示根據(jù) table()語句中的第2個變量(這里指Arthritis$Improved)來計算頻數(shù)
# 也就是計算 Improved 的三個水平各自的頻數(shù)。
margin.table(mytable)
# 不加數(shù)字的話断国，計算得到的是總人數(shù)贤姆。

addmargins()

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

CrossTable() 生成二維列聯(lián)表 推薦使用這個∥瘸模可以生成多個參數(shù)霞捡。

install.packages("gmodels")
library(gmodels)
CrossTable(Arthritis$Treatment, Arthritis$Improved)
# 得到一個表格，包含數(shù)據(jù)總數(shù)薄疚、每種條件下的數(shù)據(jù)個數(shù)碧信、以及比例。
# 這里比例包括 三種街夭，分母分別是每一行的總數(shù)砰碴、每一列的總數(shù)、數(shù)據(jù)總數(shù)板丽。

多維列聯(lián)表

ftable()函數(shù)可以 以緊湊的方式輸出多維列聯(lián)表呈枉。推薦使用這個。

mytable <- xtabs(~ Treatment+Sex+Improved, data=Arthritis)
mytable
prop.table(mytable, c(1, 2))
# c(1,2)代表 Treatment 和 Sex埃碱，計算這兩個變量2*2比例的列聯(lián)表猖辫。
# 但是要根據(jù) Improved 的不同水平來分別計算。
ftable(prop.table(mytable, c(1, 2)))
ftable(addmargins(prop.table(mytable, c(1, 2)),1))
ftable(addmargins(prop.table(mytable, c(1, 2)),2))

相關系數(shù)和顯著性的計算

cor()函數(shù)用來計算相關系數(shù);

用法：cor(x, use= , method= )
method有三種砚殿，pearson spearman kendall啃憎，默認是 pearson

cov()計算協(xié)方差。

states<- state.x77[,1:6]
cov(states) #協(xié)方差
cor(states) #相關系數(shù)
cor(states, method="spearman") # spearman等級相關
x <- states[,c("Population", "Income", "Illiteracy", "HS Grad")]
y <- states[,c("Life Exp", "Murder")]
cor(x,y) #計算得到4*2的相關矩陣

cor.test() 用來對相關系數(shù)的顯著性進行檢驗

但是一次只能檢驗一種相關關系似炎。
用法：cor.test(x, y, alternative = , method = )
alternative 取值有："two.side" "less"或"greater"辛萍。默認是 two.side，表示雙尾檢驗（總體相關系數(shù)不等于!=0）羡藐。
less 表示贩毕，假設為 總體相關系數(shù)<0；greater 表示仆嗦，假設為 總體相關系數(shù)>0
method有三種辉阶，pearson spearman kendall，默認是 pearson
cor.test(states[,3], states[,5])

如果要檢驗多種相關關系欧啤，需要使用 psych包中的corr.test()

用法：corr.test(x, y = NULL, use = "pairwise",method="pearson",adjust="holm",alpha=.05,ci=TRUE,minlength=5)

install.packages("psych")
library(psych)
corr.test(states, use="complete")
# 輸出的結果里，對角線上方的數(shù)據(jù) 是經過多重比較矯正的結果启上。下方的是沒有矯正的p值邢隧。
# use的取值 可以是 pairwise 和 complete，分別表示對缺失值執(zhí)行成對刪除 或 行刪除冈在。

pcor() 計算偏相關

用法：pcor(u, S)
u是數(shù)值向量倒慧，前兩個數(shù)值為要計算相關系數(shù)的變量下標，其余的變量為要排除影響的變量下標。

install.packages("ggm")
library(ggm)
colnames(states)
pcor(c(1,5,2,3,6), cov(states)) #偏相關 等價于 pcor(c("Population","Murder","Income","Illiteracy","HS Grad"), cov(states))
# 在本例中纫谅，1,5對應的變量("Population","Murder")是計算相關的變量炫贤，2,3,6 ("Income","Illiteracy","HS Grad")是需要排除影響的變量。

pcor.test 偏相關系數(shù)的顯著性檢驗

用法：pcor.test(r, q, n)
r是偏相關系數(shù)付秕，q是要控制的變量數(shù)(以數(shù)值表示位置)兰珍，n為樣本大小。
pcor.test(pcor(c(1,5,2,3,6), cov(states)),3,50)# 輸出結果為 t值询吴、df值和p值

t檢驗

獨立樣本t檢驗

使用方法：t.test(y ~ x, data) 或t.test(y1, y2)
默認參數(shù)是掠河，假設方差不齊的t檢驗，雙側檢驗猛计。
如果想要輸出方差齊的t檢驗唠摹，需要加入 var.equal=TRUE

library(MASS) #使用MASS包中的數(shù)據(jù)集來做練習
w<-t.test(Prob ~ So, data=UScrime)
w<-t.test(Prob ~ So, data=UScrime,var.equal=TRUE)

獨立樣本的非參數(shù)檢驗

wilcoxon秩和檢驗，也被稱作為 mann-whitney u檢驗

使用方法：wilcox.test(y ~ x, data) 或 wilcox.test(y1, y2)

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

配對樣本t檢驗

不需要方差齊性的假設
使用方法: t.test(y1, y2, paired=TRUE)

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

非獨立(配對)樣本的非參數(shù)檢驗

wilcoxon符號秩和檢驗

使用方法：wilcox.test(y1, y2, paired=TRUE)

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

回歸

對于OLS回歸奉瘤，通過使得 預測誤差（殘差）平方和最小 和 對影響變量的解釋力度（R平方）最大勾拉，可獲得模型參數(shù)。

簡單線性回歸

lm()函數(shù)

用法：myfit <- lm(formula, data)

formula的形式為：Y~ X1 + X2 + ... + Xk
~ 左邊為響應變量(因變量)盗温，~右邊是各個預測變量(自變量)藕赞，預測變量之間用+符號隔開。
:符號表示變量的交互項肌访。例如 y~x+z+x:z
*符號表示所有可能交互項的簡潔方式找默。例如代碼 y~xzw 等價于 y ~ x + z + w + x:z + x:w + z:w + x:z:w
^符號表示交互項達到某個次數(shù)。代碼 y~(x+z+w)^2 就代表交互項就到2階交互吼驶。等價于 y ~ x + z + w + x:z + x:w + z:w
.符號表示 除了因變量外的所有變量惩激，這里不包含交互作用。如果一個數(shù)據(jù)框包含變量 x y z w蟹演。那么 代碼 y~. 等價于 y ~ x + z + w
-符號表示從等式中移除某個變量风钻。例如 y~(x+z+w)^2-x:w 等價于 y ~ x + z + w + x:z + z:w

-1 表示 刪除截距項。y~x-1 擬合y在x上的回歸酒请，并強制直線通過原點骡技。截距為0。

summary() 展示擬合模型的詳細結果
coefficients() 列出擬合模型的模型參數(shù)（截距項和斜率β）
confint() 計算模型參數(shù)的置信區(qū)間（默認是95%）
fitted() 列出擬合模型的預測值
residuals() 列出擬合模型的殘差：預測值-實際值
anova() 1.生成一個擬合模型的方差分析表 或2. 比較兩個或多個擬合模型的方差分析表羞反。

vcov() 列出擬合模型的協(xié)方差矩陣
AIC() 輸出赤池信息統(tǒng)計量
plot() 生成評價擬合模型的診斷圖
predict()用擬合模型對新的數(shù)據(jù)集預測響應變量值

fit <- lm(weight ~ height, data=women)
summary(fit)

結果里面的兩個參數(shù) Multiple R-squared 和 Adjusted R-squared 的值：
R-squared 表示模型可以解釋預測變量的方差布朦，也就是模型預測變量（預測值）解釋響應變量（真實值）的程度。也是實際值和預測值之間相關系數(shù)的平方昼窗。Adjusted R-squared 與之類似是趴，但是考慮了模型的參數(shù)數(shù)目。

women$weight
fitted(fit)
residuals(fit)
plot(women$height,women$weight,
     xlab="Height (in inches)",
     ylab="Weight (in pounds)")
abline(fit)

par(mfrow=c(2,2)) #兩行兩列的圖展示在一起
plot(fit) #輸出四張圖 可以檢驗正態(tài)性 方差齊性 和異常值

線性回歸的前提條件

1. 正態(tài)性

1.1 使用performance包 check_normality()

library(performance)
x<-check_normality(fit)
plot(x,type = "qq") #畫qq圖
plot(x, type = "pp") #畫pp圖
plot(x) #畫殘差的正態(tài)分布圖

1.2 不使用包 直接畫標準化殘差 qq圖

qqnorm(residuals(fit))
qqline(residuals(fit))

1.3 使用car包里的qqPlot() 畫t分布下的學生化/標準化殘差 qq圖

library(car)
qqPlot(fit)

1.4 直接畫直方圖

hist(residuals(fit))

2. 檢驗自變量和因變量是否是線性關系

library(car)
crPlots(fit)

3. 檢驗方差齊性

3.1 直接畫預測值和殘差圖 殘差沒有標準化

plot(fitted(fit),residuals(fit)) 
abline(h=0,lty=2)

3.2 使用performance包 check_heteroscedasticity()

x<-check_heteroscedasticity(fit)
plot(x)

3.3 使用car包中的 ncvTest()來做檢驗澄惊，同時用spreadLevelPlot()畫圖

library(car)
ncvTest(fit)
spreadLevelPlot(fit)

4. 誤差的獨立性

檢驗誤差的序列相關性（看到一個資料說唆途，僅僅時序類數(shù)據(jù)需要進行該診斷）富雅。

使用car包中的durbinWatsonTest()函數(shù)來進行 Durbin-Watson 檢驗。

durbinWatsonTest(fit)

線性模型假設的綜合檢驗

gvlma包中的gvlma()函數(shù)

gvlma()函數(shù)只適用于 lm()模型

install.packages("gvlma")
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

回歸中的多重共線性問題

可用統(tǒng)計量 VIF(Variance Inflation Factor,方差膨脹因子) 來檢測
vif>5 就認為存在嚴重的共線性肛搬。<5可以忽略

1. 可用car包中的vif()函數(shù)來計算得到

library(car)
vif(fit)
# sqrt(vif(fit)) > 2 # problem?

2. 可用performance包中的check_collinearity()函數(shù)來計算得到

check_collinearity(fit)

離群點檢測

car包中的outlierTest()

outlierTest()函數(shù)可以求得最大標準化殘差絕對值 Bonferroni 調整后的p值
該函數(shù)只根據(jù)單個最大（正或負）殘差值的顯著性來判斷是否有離群點没佑。若不顯著則說明數(shù)據(jù)集中沒有離群點；若顯著温赔，則你必須刪除該離群點蛤奢，然后在檢測是否還有其他的離群點存在。

library(car)
outlierTest(fit)

強影響點檢測

cutoff <- 4/(nrow(states)-length(fit$coefficients)-2) 
cutoff<-1
plot(fit, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")

1. performance包中的check_outliers()

check_outliers(fit)

二項式回歸

fit2 <- lm(weight ~ height + I(height^2), data=women)
summary(fit2)
fitted(fit2)
plot(women$height,women$weight,
     xlab="Height (in inches)",
     ylab="Weight (in lbs)")
lines(women$height,fitted(fit2))
par(mfrow=c(2,2)) 
plot(fit2)

多元線性回歸

多元就是指 多個自變量让腹。

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
cor(states) #做多元回歸之前远剩，最好先計算一下變量之間的相關性。
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
summary(fit)
confint(fit)

檢驗誤差的序列相關性（看到一個資料說骇窍，僅僅時序類數(shù)據(jù)需要進行該診斷）瓜晤。

使用car包中的durbinWatsonTest()函數(shù)來進行 Durbin-Watson 檢驗。

library(car)
durbinWatsonTest(fit)
crPlots(fit)

有顯著交互項的多元線性回歸

fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
summary(fit)
install.packages("effects")
library(effects) 
# 用effects包的 effect()函數(shù)來畫圖展示 交互項的結果
# 用法：plot(effect(term,mod,,xlevels),multiline=TRUE)
# term是要畫的項腹纳，mod是lm()擬合的模型痢掠，xlevels是一個list，用來指定變量要設定的常量值嘲恍。
# multiline=TRUE 表示添加多條相應直線足画。multiline=FALSE 表示多個xlevels分開畫多圖。
plot(effect("hp:wt", fit,, list(wt=c(2.2,3.2,4.2))), multiline=TRUE)
plot(effect("hp:wt", fit,, list(wt=c(2.2,3.2,4.2))), multiline=FALSE)

模型比較

1. 用anova()函數(shù)進行模型比較

注意用anova()函數(shù)比較模型時佃牛，模型需要是嵌套模型淹辞。就比如下面的例子中，一個模型fit1包含了fit2俘侠。

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
fit2 <- lm(Murder ~ Population + Illiteracy, data=states)
anova(fit2, fit1)
# 結果顯示p>0.05象缀，并不顯著。說明兩個模型fit1和fit2的擬合效果一樣好爷速。所以可以從模型中去掉"Income", "Frost"兩個變量央星。

2. 用AIC(Akaike Information Criterion，赤池信息準則)來比較模型

AIC考慮了模型統(tǒng)計擬合度以及用來擬合的參數(shù)數(shù)目惫东。AIC越小莉给，代表模型擬合的越好，說明模型用較少的參數(shù)獲得了足夠的擬合度廉沮。

AIC方法不需要模型之間是嵌套關系颓遏。

fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
fit2 <- lm(Murder ~ Population + Illiteracy, data=states)
AIC(fit1,fit2)

怎么選擇“最佳”模型

如何確定在模型中保留哪些變量 模型是最優(yōu)的呢。anova()一般用于有2滞时，3個模型比較的時候叁幢。如果模型特別多，該怎么辦呢漂洋？

1. 逐步回歸：向前逐步回歸遥皂、向后逐步回歸、向前向后逐步回歸

使用MASS包中的stepAIC()函數(shù)

依據(jù)的是 精確AIC準則
下面是一個向后回歸的例子

library(MASS)
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
stepAIC(fit, direction="backward")

2. 全子集回歸

全子集回歸中刽漂，所有可能的模型都會被檢驗演训。

使用leaps包中的regsubsets()函數(shù)

可以通過R平方，調整R平方或Mallows Cp 統(tǒng)計量等準則來選擇“最佳”模型

install.packages("leaps")
library(leaps) 
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")
# 結果可以用car包中的subsets()函數(shù)繪制
library(car)
subsets(leaps, statistic="cp",
        main="Cp Plot for All Subsets Regression")
abline(1,1,lty=2,col="red")

交叉驗證

通過交叉驗證法贝咙，可以評價回歸方程的泛化能力样悟。

交叉驗證：將一定比例的數(shù)據(jù)挑選出來作為訓練樣本（訓練集），另外的樣本作為保留樣本（測試集）庭猩。先在訓練樣本上獲取回歸方程窟她，然后在保留樣本上做預測。
在k重交叉驗證中蔼水，樣本被分為k個子樣本震糖，輪流k-1個子樣本組合作為訓練集，另外1個樣本作為測試集趴腋。
這樣就可以得到k個預測方程吊说，將k個預測集的預測表現(xiàn)結果求平均值。

1. R平方的k重交叉驗證：bootstrap包中crossval()函數(shù)

方差分析

以下的概念中：大寫字母表示組別因子优炬，小寫表示定量變量颁井。
單因素方差分析ANOVA：y~A
含單個協(xié)變量的單因素協(xié)方差分析ANCOVA：y~x+A
雙因素ANOVA：y~AB
含兩個協(xié)變量的兩因素協(xié)方差分析ANCOVA：y~x1+x2+AB
隨機化區(qū)組：y~B+A（B是區(qū)組因子）
單因素組內ANOVA（單因素重復測量anova）：y~A+Error(Subject/A)
含單個組內因子（W）和單個組間因子（B）的重復測量anova：y~B*W+Error(Subject/W)
均衡設計（balanced design）是指 每個實驗條件下的被試量相等。否則就是非均衡設計（unbalanced design）蠢护。
有多個（大于一個）因變量的設計被稱為多元方差分析（MANOVA）雅宾，若有協(xié)變量，則被稱為 多元協(xié)方差分析（MANCOVA）葵硕。

表達式的順序 在以下兩種情況下眉抬，會對結果產生影響：

1. 因子不止一個，且是非均衡設計贬芥；
2. 存在協(xié)變量吐辙。

• 類型I方法（序貫型） 計算ANOVA效應：
  因為R中默認 使用類型I方法: y~A+B+A:B
  A對y的影響；A不做調整
  控制A時蘸劈，B對y的影響昏苏；B根據(jù)A調整
  控制A和B的主效應時，A與B的交互效應威沫。A:B交互項根據(jù)A和B調整贤惯。
• 類型II（分層型）：
  效應根據(jù)同水平或低水平的效應做調整。A根據(jù)B調整棒掠，B根據(jù)A調整孵构，A:B交互項同時根據(jù)A和B調整。
• 類型III（邊界型）：SPSS 和 SAS 默認使用的是類型III烟很。所以如果要得到相同的結果颈墅，需要使用類型III蜡镶。aov()函數(shù)默認使用的是類型I的方法，若要使用類型III恤筛，需要用car包中的Anova()函數(shù)官还，具體可參考 help(Anova,package="car")。
  每個效應根據(jù)模型其他各效應做相應調整毒坛。A根據(jù)B和A:B做調整望伦，A:B交互項根據(jù)A和B調整。

樣本大小越不平衡煎殷，效應項的順序對結果的影響越大屯伞。一般來說，越基礎性的效應越需要放在表達式前面豪直。具體來講：
首先是協(xié)變量劣摇，然后是主效應，接著是雙因素的交互項弓乙，再接著是三因素的交互項饵撑，以此類推。
對于主效應唆貌，越基礎性的效應越應放在表達式前面滑潘。

aov(formula, data=dataframe)

單因素方差分析 one-way anova

單因素方差分析 對應的 非參數(shù)檢驗 數(shù)據(jù)之間獨立

Kruskal-Wallis檢驗

使用方法：kruskal.test(y ~ A, data)

states <- data.frame(state.region, state.x77)
kruskal.test(Illiteracy ~ state.region, data=states)
# 通過r語言實戰(zhàn)第二版書中 作者寫的wmc函數(shù) 做事后檢驗的多重比較
source("http://www.statmethods.net/RiA/wmc.txt")
wmc(Illiteracy ~ state.region, data=states, method="holm")
help(p.adjust) #查看其它多重比較矯正的方法
# p.adjust.methods 
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
#   "fdr", "none")

單因素重復測量方差分析 對應的 非參數(shù)檢驗 數(shù)據(jù)之間不獨立

Friedman檢驗

使用方法：friedman.test(y ~ A | B, data)

RoundingTimes <-
  matrix(c(5.40, 5.50, 5.55,
           5.85, 5.70, 5.75,
           5.20, 5.60, 5.50,
           5.55, 5.50, 5.40,
           5.90, 5.85, 5.70,
           5.45, 5.55, 5.60,
           5.40, 5.40, 5.35,
           5.45, 5.50, 5.35,
           5.25, 5.15, 5.00,
           5.85, 5.80, 5.70,
           5.25, 5.20, 5.10,
           5.65, 5.55, 5.45,
           5.60, 5.35, 5.45,
           5.05, 5.00, 4.95,
           5.50, 5.50, 5.40,
           5.45, 5.55, 5.50,
           5.55, 5.55, 5.35,
           5.45, 5.50, 5.55,
           5.50, 5.45, 5.25,
           5.65, 5.60, 5.40,
           5.70, 5.65, 5.55,
           6.30, 6.30, 6.25),
         nrow = 22,
         byrow = TRUE,
         dimnames = list(1 : 22,
                         c("Round Out", "Narrow Angle", "Wide Angle")))
friedman.test(RoundingTimes)
## => strong evidence against the null that the methods are equivalent
##    with respect to speed
wb <- aggregate(warpbreaks$breaks,
                by = list(w = warpbreaks$wool,
                          t = warpbreaks$tension),
                FUN = mean)
wb
friedman.test(wb$x, wb$w, wb$t)
friedman.test(x ~ w | t, data = wb)