邏輯回歸：股票市場數據
library(ISLR) names(Smarket)
dim(Smarket)
summary(Smarket) pairs(Smarket) #pairs（）函數用於返回一個繪圖矩陣，由每個 DataFrame 對應的散點圖組成。
cor(Smarket[,-9]) #cor ()函數可以計算所有預測變最兩兩之間相關系數的矩陣。
解釋：相關的是 Year和Volume ，通過畫圖可以觀察到 Volume 隨時間一直增長，也就是說從 2001 年至 2005 年平均每日股票成交盤在增長。
attach(Smarket) plot(Volume)
邏輯回歸
glm.fits=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomial)#設置family=binomial是請求R執行邏輯斯蒂回歸 summary(glm.fits)
結果分析：這里最小的 值是Lag1的系數。預測變量負的系數表明如果市場昨天的投資回報率是正的，那么今日市場可能不會上漲。
coef(glm.fits) #調用coef ()函數可獲取擬合模型的系數
summary(glm.fits)$coef #查看詳細信息 summary(glm.fits)$coef[,4] glm.probs=predict(glm.fits,type="response") #predict ()函數可用來預測在給定預測變盤值下市場上漲的概率。選項 type =" re~ sponse" 是明確告訴R輸出概率 P(Y=1| X) ，而不必輸出其他信息。 glm.probs[1:10]
contrasts(Direction) #contrasts ()函數表創建了一個啞變量，1 代表 Up (上漲)。
#為預測特定某一天市場上漲還是下跌，必須將預測的概率轉化為類別，Up 或者Down，這里以預測市場上漲概率比0.5 大還是小為依據。
glm.pred=rep("Down",1250)
glm.pred[glm.probs>.5]="Up"
table(glm.pred,Direction)
(507+145)/1250 #計算正確率
mean(glm.pred==Direction) #計算正確率
train=(Year<2005) #先建立一個向致表示2000年至 2004 年之間的觀測，然后由這個向量
定義余下的 2005 年觀測數據集。
解釋：train 是一個布爾向量 (Boolean vector)，train 是個有 1250 個元素的向量，對應數據集中的每條觀測，該向量對數據集2005 年之前的數據設置為TRUE，對2005 年以后的數據標記為 FALSE" 。
Smarket.2005=Smarket[!train,] #2005年的觀測數據 dim(Smarket.2005)
Direction.2005=Direction[!train] #因子型類別變量，Up 或者Down
glm.fits=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomial,subset=train) #建立模型
glm.probs=predict(glm.fits,Smarket.2005,type="response") #預測
glm.pred=rep("Down",252)
glm.pred[glm.probs>.5]="Up"
table(glm.pred,Direction.2005)
mean(glm.pred==Direction.2005) #計算正確率
mean(glm.pred!=Direction.2005) #計算錯誤率
#加如與響應變量無關的預測變量會造成測試錯誤率的增大（因為這樣的預測變量會增大模型方差，但不會相應地降低模型偏差) ，所以去除這樣的預測變量可能會優化模型。所以重新擬合邏輯斯諦回歸模型，只是用 Lag1和 Lag2 兩個預測變量，因為它們在原邏輯斯諦因歸模型中表現出了最佳的預測能力。
glm.fits=glm(Direction~Lag1+Lag2,data=Smarket,family=binomial,subset=train)
glm.probs=predict(glm.fits,Smarket.2005,type="response")
glm.pred=rep("Down",252)
glm.pred[glm.probs>.5]="Up"
table(glm.pred,Direction.2005)
mean(glm.pred==Direction.2005) #計算正確率
106/(106+76) #表明當預測市場上漲時，准確率有 58%
predict(glm.fits,newdata=data.frame(Lag1=c(1.2,1.5),Lag2=c(1.1,-0.8)),type="response") #特定的 Lagl1和Lag2 值下預測投資回報率。
線性判別分析（Linear Discriminant Analysis,LDA）
#在對 Smarket 數據傲線性判別分析 (LDA) ，用 MASS 庫中的 lda ()函數擬合一個 LDA 模型。lda ()函數與 lm() 函數語句是相彤的，也和glm ()類似，除了缺少family 選項。
library(MASS) lda.fit=lda(Direction~Lag1+Lag2,data=Smarket,subset=train) #建立模型 lda.fit
plot(lda.fit) #plot ()函數生成的線性判判別圖像，可以通過對每個訓練觀測計算- 0.642 x Lag1-0.514 x Lag2 獲得。
lda.pred=predict(lda.fit, Smarket.2005) #預測2005年的 names(lda.pred)
lda.class=lda.pred$class table(lda.class,Direction.2005)
mean(lda.class==Direction.2005) #LDA 和邏輯斯諦回歸預測結果幾乎一樣。
二次判別分析(quadratic discriminant analysis,QDA)
#現在對smarket數據擬告 QDA 模型，QDA 在民中可以調用 MASS 庫中的 qda ()函數來實現，其語句與 lda ()一樣。
qda.fit=qda(Direction~Lag1+Lag2,data=Smarket,subset=train) qda.fit
qda.class=predict(qda.fit,Smarket.2005)$class table(qda.class,Direction.2005)
mean(qda.class==Direction.2005)
結果分析：即使沒有用 2005 年始數據來擬合模型，在幾乎 60% 的時間)里， QDA 預測也都是准確的，這表明 QDA假設的二次型比 LDA 和邏輯斯諦回歸的線性假設更接近真實的關系。
邏輯回歸案例：大篷車保險數據的一個應用
對 ISLR庫中的 Caravan (大篷車)數據集運用 KNN 方法。該數據集包括 85 個預測變量，測量了 5822 人的人口特征。響應變最為 Purchase (購買狀態) .表示一個人是否會購買大篷車保險險種。在該數據集中，只有 6% 的人購買了大篷車保險。
dim(Caravan)
attach(Caravan) summary(Purchase)
348/5822
standardized.X=scale(Caravan[,-86]) #標准化數據，第86列為Purchase狀態，不適合標准化 var(Caravan[,1])
var(Caravan[,2])
var(standardized.X[,1]) #命令 standardized.X 將每列均值為0，標准差置為 1 。
var(standardized.X[,2])
test=1:1000
train.X<-standardized.X[-test,] #訓練集
test.X<-standardized.X[test,] #測試集，樣本量為1000
train.Y<-Purchase[-test] #訓練集標簽
test.Y<-Purchase[test] #測試集標簽
glm.fits<-glm(Purchase~.,data=Caravan,family=binomial,subset=-test)
glm.probs<-predict(glm.fits,Caravan[test,],type="response")
glm.pred<-rep("No",1000)
glm.pred[glm.probs>.5]="Yes"
table(glm.pred,test.Y)
glm.pred<-rep("No",1000)
glm.pred[glm.probs>.25]="Yes"
table(glm.pred,test.Y)
11/(22+11)
結果分析：作為比較，對數據擬合一個邏輯斯諦回歸模型。如果用 0.5 作為分類器的預測概率的閾值，則出現一個問題:只有7個測試觀測被預測會購買保險。更糟糕的是，這7個預測都是錯誤的！現在不要求閾值為0.5。如果預測購買概率超過 0.25 時就預測購買的話，結果會好很多，預測有 33 人將會購買保險，而在這些人中正確預測率為33%；這比隨機猜測的五倍還多!
邏輯回歸案例：婚外情數據的應用
當通過一系列連續型和/或類別型預測變量來預測二值型結果變量時，Logistic回歸是一個非常有用的工具。以AER包中的數據框Affairs為例，我們將通過探究婚外情的數據來闡述Logistic回歸的過程。
婚外情數據即著名的“Fair’s Affairs”，取自於1969年《今日心理》（Psychology Today）所做的一個非常有代表性的調查，該數據從601個參與者身上收集了9個變量，包括一年來婚外私通的頻率以及參與者性別、年齡、婚齡、是否有小孩、宗教信仰程度（5分制，1分表示反對，5分表示非常信仰）、學歷、職業（逆向編號的戈登7種分類），還有對婚姻的自我評分（5分制，1表示非常不幸福，5表示非常幸福）。
install.packages("AER")
data(Affairs, package="AER")
summary(Affairs)
table(Affairs$affairs) #統計婚外情次數的人數
結果分析：從這些統計信息可以看到，315/601=52%的調查對象是女性，430/601=72%的人有孩子，樣本年齡的中位數為32歲。對於響應變量，72%的調查對象表示過去一年中沒有婚外情（451/601），而婚外偷情的最多次數為12（占了6%）。
雖然這些婚姻的輕率舉動次數被記錄下來，但此處我們感興趣的是二值型結果（有過一次婚
外情/沒有過婚外情）。按照如下代碼，可將affairs轉化為二值型因子ynaffair：
Affairs$ynaffair[Affairs$affairs > 0] <- 1 #賦值
Affairs$ynaffair[Affairs$affairs == 0] <- 0
Affairs$ynaffair <- factor(Affairs$ynaffair,
levels=c(0,1), #值的結果是0和1
labels=c("No","Yes")) #給二值型因子貼上標簽
table(Affairs$ynaffair)
該二值型因子現可作為Logistic回歸的結果變量：
fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children +
religiousness + education + occupation +rating,
data=Affairs, family=binomial()) #family是設置概率分布
summary(fit.full)
結果分析：從回歸系數的p值（最后一欄）可以看到，性別、是否有孩子、學歷和職業對方程的貢獻都不顯著（后面沒有星星），也就是我們無法拒絕參數為0的假設。去除這些變量重新擬合模型，檢驗新模型是否擬合得好呢？
fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, data=Affairs, family=binomial())
summary(fit.reduced)
結果分析：新模型的每個回歸系數都非常顯著（p<0.05）。由於兩模型嵌套（fit.reduced是fit.full的一個子集），你可以使用anova()函數對它們進行比較，對於廣義線性回歸，可用卡方檢驗。
anova(fit.reduced, fit.full, test="Chisq")
結果分析：結果的卡方值不顯著（p=0.21），表明四個預測變量的新模型與九個完整預測變量的模型擬合程度一樣好。這使得你更加堅信添加性別、孩子、學歷和職業變量不會顯著提高方程的預測精度，因此可以依據更簡單的模型進行解釋。
解釋模型參數：先看看回歸系數
coef(fit.reduced)
對每個觀測值各分解項的值取指數，將結果轉化為原始尺度：
exp(coef(fit.reduced))
結果分析：可以看到婚齡增加一年，婚外情的優勢比將乘以1.106（保持年齡、宗教信仰和婚姻評定不變）；相反，年齡增加一歲，婚外情的的優勢比則乘以0.965。因此，隨着婚齡的增加和年齡、宗教信仰與婚姻評分的降低，婚外情優勢比將上升。
評價預測變量對結果概率的影響
使用predict()函數，可以觀察某個預測變量在各個水平時對結果概率的影響。首先創建一個包含你感興趣預測變量值的虛擬數據集，然后對該數據集使用predict()函數，以預測這些值的結果概率。
現在我們使用該方法評價婚姻評分對婚外情概率的影響。首先，創建一個虛擬數據集，設定年齡、婚齡和宗教信仰為它們的均值，婚姻評分的范圍為1～5。
testdata <- data.frame(rating=c(1, 2, 3, 4, 5), age=mean(Affairs$age),
yearsmarried=mean(Affairs$yearsmarried),
religiousness=mean(Affairs$religiousness))
testdata
testdata$prob <- predict(fit.reduced, newdata=testdata, type="response") #用測試數據集預測相應的概率
#predict：從擬合的穩健廣義線性模型（GLM）對象中獲取預測值，並可選擇估計這些預測值的標准誤差。Newdata：可選地，一個數據框架，在其中尋找用於預測的變量，如果省略，則使用擬合的線性預測器。
testdata
結果分析：從這些結果可以看到，當婚姻評分從1（很不幸福）變為5（非常幸福）時，婚外情概率從0.53降低到了0.15（假定年齡、婚齡和宗教信仰不變）。下面我們再看看年齡的影響：
testdata <- data.frame(rating=mean(Affairs$rating),
age=seq(17, 57, 10), #從17歲到57歲，每隔10歲生成一個
yearsmarried=mean(Affairs$yearsmarried),
religiousness=mean(Affairs$religiousness))
testdata
testdata$prob <- predict(fit.reduced, newdata=testdata, type="response") testdata
結果分析：此處可以看到，當其他變量不變，年齡從17增加到57時，婚外情的概率將從0.34降低到0.11。