模型的解釋力和預(yù)測力的區(qū)別
解釋力的標(biāo)準(zhǔn)：R方，調(diào)整后R方
常用的預(yù)測力標(biāo)準(zhǔn)：Cp，AIC，BIC，AUC等等
案例：1986年美國各大核心棒球球隊(duì)隊(duì)員的表現(xiàn)數(shù)據(jù)和次年的收入
數(shù)據(jù)包：ISLR，Hitters數(shù)據(jù)
fix (Hitters)會出現(xiàn)bug，可以使用Hitters = ISLR::Hitters語句
首先進(jìn)行數(shù)據(jù)清洗，去掉空白的數(shù)據(jù)，判斷是否空白（NA）的函數(shù)：
is.na()
找出那些行是NA：which（is.na(Hitters[,19])）#Summary顯示只有19列salary存在空白數(shù)據(jù)
命名為drop = which (is.na (Hitters[, 19]))
新的沒有NA的數(shù)據(jù)：newhitters = Hitters[-drop,]
新數(shù)據(jù)有多少行：dim (newhitters)
先對球員工資進(jìn)行一個(gè)探索性分析，做出頻數(shù)直方圖
ggplot (newhitters, aes (x = Salary))+ #在ggplot中，+號類似于ps里的疊加圖層
? geom_histogram (bins = 40, fill = 'red')+
? ggtitle ("Histogram of Salary")+
? ylab ("Frequency")#y軸名稱
如果要改成頻率直方圖：
ggplot (newhitters) +
geom+histogram (aes (x = Salary, y = ..density..), bins = 40, fill = "blue") +
ggtitle ("HIstogram of Salary")
直方圖明顯右斜（頭部效應(yīng)），需要轉(zhuǎn)化成正態(tài)分布的數(shù)據(jù)，處理方法：取對數(shù)
newhitters = newhitters %>% mutate (lnsalary = log (Salary))
制作相應(yīng)的兩兩相關(guān)散點(diǎn)圖，可以看出由于數(shù)據(jù)數(shù)量較多速度很慢
當(dāng)變量數(shù)量多的時(shí)候，可以選擇相關(guān)性熱力圖
所需程序包：dplyr，ggplot2，tidyverse和reshape2，代碼如下：
plotData = melt (cor(newhitters[sapply (newhitters, is.numeric)]))#將非數(shù)變量剔除，melt函數(shù)將數(shù)據(jù)轉(zhuǎn)換成可以用于制圖的數(shù)據(jù)
ggplot (plotData,
? ? ? ? aes (x = Var1, y = Var2, fill = value))+
? geom_tile()+
? ylab("")+
? xlab("")+
? scale_fill_gradient(low = "#56B1F7", high = '#132B43')+
? guides (fill = guide_legend (title = 'correlation'))
說明大多數(shù)球員表現(xiàn)得數(shù)據(jù)是和來年的薪水相關(guān)的
數(shù)據(jù)預(yù)測能力的判斷：預(yù)測誤差1/n(∑n, i = 1((yi-y'i)^2)))越小越好
兩種測量預(yù)測誤差的方法：直接法與間接法
直接法：
再隨機(jī)收集一批新數(shù)據(jù)，康康這個(gè)模型在新數(shù)據(jù)里的預(yù)測情況
不可以用生成模型的數(shù)據(jù)
如果收集新數(shù)據(jù)的成本可能太高，將已有的數(shù)據(jù)分成兩個(gè)部分
或者cross validation交叉驗(yàn)證
間接法：
找一個(gè)類似于R方的統(tǒng)計(jì)量，用來估算預(yù)測誤差
常用的間接統(tǒng)計(jì)量：Mallow's Cp or BIC
Cp = 1/n(RSS + 2d*^σ^2)
n是數(shù)據(jù)的個(gè)數(shù)，RSS是殘差平方和，^σ模型均方誤（mean square error）
Cp的性質(zhì)：
隨著數(shù)據(jù)量增大，Cp會無線趨近于真實(shí)的預(yù)測誤差
Cp越小，模型誤差越小，也就是預(yù)測精度越高
BIC： Bayesian InformationCriterion 貝葉斯信息標(biāo)準(zhǔn)
貝葉斯方法：不存在頻率學(xué)派生成的客觀概率，所謂概率是主觀臆斷和后驗(yàn)數(shù)據(jù)的信息之和
BIC = 1/n (RSS + k * ln(n)*^σ^2)
k是要計(jì)算的參數(shù)個(gè)數(shù)，對于自變量個(gè)數(shù)加一，因?yàn)槌嗣總€(gè)自變量前面的系數(shù)以外，還需要計(jì)算一個(gè)截距項(xiàng)。BIC越小，預(yù)測能力越好
對模型進(jìn)行選擇，首先對原始數(shù)據(jù)進(jìn)行處理，刪掉salary列改為lnsalary
共有20個(gè)變量，其中潛在有19個(gè)自變量
需要使用程序包：leaps
fit1 = regsubsets (lnsalary~., newhitters, nvmax = 20, method = 'exhaustive')
如果需要包含二次項(xiàng)，.處改為.*.
Nvmax參數(shù)：最多使用的變量個(gè)數(shù)
Exhaustive：還有forward和backward參數(shù)，推薦exaustive
觀察回歸結(jié)果：summary (fit1)$which 標(biāo)準(zhǔn)：殘差平方和小
左側(cè)表頭意味分別選取1/2/3/4……個(gè)自變量的時(shí)候，那一個(gè)或幾個(gè)自變量的相關(guān)性好
此時(shí)已獲取10個(gè)局部最優(yōu)的模型，接下來需要尋找全局最優(yōu)的模型。
本次涉及的全部代碼：
library (GGally)
library (dplyr)
library (ggplot2)
library (tidyverse)
library (reshape2)
Hitters = ISLR::Hitters
str (Hitters)
drop = which (is.na (Hitters[, 19]))
newhitters = Hitters[-drop,]
dim (newhitters)
ggplot (newhitters, aes (x = Salary))+
? geom_histogram (bins = 40, fill = 'red')+
? ggtitle ("Histogram of Salary")+
? ylab ("Frequency")
newhitters = newhitters %>% mutate (lnsalary = log (Salary))
newhitters %>% select (lnsalary, AtBat, Hits, HmRun, Runs, RBI) %>% ggpairs()
plotData = melt (cor(newhitters[sapply (newhitters, is.numeric)]))
ggplot (plotData,
? ? ? ? aes (x = Var1, y = Var2, fill = value))+
? geom_tile()+
? ylab("")+
? xlab("")+
? scale_fill_gradient(low = "#56B1F7", high = '#132B43')+
? guides (fill = guide_legend (title = 'correlation'))
newhitters[, 19] = NULL
fit1 = regsubsets (lnsalary~., newhitters, nvmax = 20, method = 'exhaustive')
summary (fit1)
summary (fit1)$which
