相關係數(Correlation) 描述兩個變數X、Y之間的線性相關 Example: data1中的身高及體重 如何量化這樣的線性關係呢? Correlation! Linear correlation!
相關係數(Correlation) By definition, the correlation between X and Y is Its estimate, Pearson’s correlation coefficient
相關係數(Correlation) r>o: positively correlated r<0: negatively correlated r=0: no linear correlation r=0不代表、Y之間沒有關係 ,有可能只是他們之間的關 係不是線性的 →畫圖還是必要的
相關係數(Correlation) R程式:cor(x,y,method = c("pearson", "kendall", "spearman")) ) x: 數值向量或是矩陣 y: 數值向量,當x是矩陣的時候,可以不需輸入
相關係數(Correlation) 若想進一步檢定 vs. 檢定統計量 95% confidence interval:
相關係數(Correlation) R程式:cor.test(x, y, alternative = c("two.sided", "less", "greater"), method = c("pearson", "kendall", "spearman"), exact = NULL, conf.level = 0.95, continuity = FALSE, ...) x: 數值向量 y: 數值向量 exact: T或F,表示是否計算exact p-value continuity: 是否需要進行連續校正 所以身高與體重有統計顯著的正相關
Practice 請畫出在Surgical data中,liver與clot的散佈圖。 請問由圖中,可以看出liver與clot的關係嗎? Q: 除了看相關性的強度,能不能看彼此如何影 響?Regression!
Linear Regression 血壓是否和體重有線性相關; 該線性關係如何描述; 如何描述血壓和體重、性別、等等的關係。 Y: response variable, dependent variable (say, bp) X: covariate, explanatory variable, independent variable (say, weight)
Linear Regression Q: how does X affect Y? Can we fit a line in the scatter plot? In fact, we should say , where is called error, is normal with zero mean and variance 2.
Linear regression using R R程式:lm(formula, data, ...) formula: y~x,其中y是response,x是covariate 3.943=70.8432/17.9663
Linear regression Confidence interval of and ? Use t-distribution with df=n-2 Testing if the coefficient =0? If =0? Use t with df=n-2 An increase of 1kg in Weight leads to an increase of 0.7167 in Bp. If someone weighs 70kg, then his/her bp is estimated by 70.84 + 0.7270 = 121.24 --- interpolation
Practice 想知道在Surgical data中,clot如何影響liver,請 建立liver與clot之迴歸模式。 如何解釋此模型呢? 請問clot對liver的影響是顯著的嗎?
Homework 想知道在Surgical data中,enzyme如何影響 SVtime,請建立enzyme與SVtime之迴歸模式。 如何解釋此模型呢? 請問enzyme對SVtime的影響是顯著的嗎?
How good is the regression? 14 How good does the line explain all the variation in y? How good does the fitted correlation of (X,Y) explain Y? 因為 定義判斷係數(coefficient of determination): Pearson’s correlation coefficient In simple linear regression, SSTO SSE SSR deviation of fitted values around the grand mean total deviation in responses around the grand mean deviation of observations around fitted line percentage of variation explained by regression line
Example 15 R2=0.4149
AVOVA table of regression 16 SSR SSE
Practice 17 在Surgical data中,模式為liver~clot 請問在此模型中,判斷係數為多少
Diagnostics 基本假設:殘差平均為0,相差變異數相同,殘 差之間不相關 看殘差和index的關係(應該要沒關係) 18 基本假設:殘差平均為0,相差變異數相同,殘 差之間不相關 看看殘差的分佈情況 看殘差和index的關係(應該要沒關係) 殘差應該要與解釋變數無關 殘差應該要和fitted value無關
Diagnostics 19 If… From minus to positive! Model may not be proper. Time effect? (If x=time) Randomly scattered around zero! Constant var有問題;若X值大則var大;試試加別的X或是weighted LS? Linearity 有問題試試polynomial 或transform X?
Example 20
Diagnostics in R 21
Diagnostics plots to examine 6→fitted model→2→3→1→4→5 22 plots to examine The linear effect of each predictor: or Constant variance: Independence of samples: or Normality assumption: Q-Q plot Other important predictors? Say : Are there outliers: , scatter plot, … If Yes, examine if it is true outlier, or gross error. If Yes, more data near this point. If No, delete the data point before regression analysis. 6→fitted model→2→3→1→4→5
Multiple linear regression 23 Extension of SLR, including more than one predictors in the model Linear? Linear? Difference?
Multiple linear regression 24 Model: : regression coefficients : observed data are independent In matrix form
Multiple linear regression 25 哪些term可以放到X中呢? Predictors: 如例子中的weight, age, sex Transformations of predictors Polynomials: and Dummy variables and factors Interactions and other combinations of predictors:
Example 26
Inference of regression coefficients 27 和SLR時一樣,用最小平方法 satisfy Gauss-Markov Thm
Inference of regression coefficients 28 和在SLR中相同,我們想要估計 的confidence interval, 或是進行檢定,需要先估計出 Recall, in SLR H is called hat matrix SST=SSE+SSR There are p-1 covariates in the regression model. There are n observations and p parameters.
Inference of regression coefficients 29 想要知道整個模式fit如何: Under , E(MSR)= ; otherwise E(MSR)> Define , with df=(p-1,n-p) 在H0之下, ,所以如果F偏離1太遠,我們就傾 向拒絕H0 H1 是什麼呢?
Inference of regression coefficients 30 若是針對某個 ,想知道 是否和 有線性關係 在H0之下, 所以拒絕H0 ,如果 你可以由此推出 的confident interval嗎?
Example 31
Practice 32 在Surgical data中想知道影響存活時間(SVtime )的因素,將存活時間取自然對數。有興趣的 因素為clot、prog、enzyme與age 請寫下此迴歸模式 請問prog的係數為0嗎? 請問此模式顯著嗎?
Homework 在bodyfat資料中,共包含4個變項(Y、X1、X2 、X3) 33 在bodyfat資料中,共包含4個變項(Y、X1、X2 、X3) 請分別畫出Y與X1、X2、X3的散佈圖,請問Y和X1 、X2、X3有線性關係嗎? 請分別檢定X1、X2、X3的迴歸係數是否為0 請問此模式是顯著的嗎?
Inference of regression coefficient 34 Estimation Least square estimator Normal assumption for interval estimate Testing For overall model , F-test For single , t-test
Example 35
Example 36
Practice 37 在Surgical data中想知道影響存活時間(SVtime )的因素,將存活時間取自然對數。有興趣的 因素為clot、prog、enzyme與age 請寫下此迴歸模式 請問此模式之adjusted R2 為多少? 請問prog的係數為0嗎? 請問此模式顯著嗎?
Simultaneous tests for partial coefficients 38 To test several parameters simultaneously, it is equivalent to “compare” two regression models, one contains all covariates and the other contains less covariates Use “extra sums of squares” to “distinguish” the two models
Extra sums of squares For two regression models 39 For two regression models model A: model B: Their SSEs will be different The difference is defined as extra sum of squares Similarly, the extra sum of squares 在已有X1的情況下,模式中增加X2的影響
在已有X1、X2的情況下,模式中增加X3的影響 What if using SSR? 40 From thus Similarly, the extra sum of squares And 在已有X1、X2的情況下,模式中增加X3的影響
Example 41 SSE decreases by 320.77; SSR increases by 320.77; Extra sums of squares, SSR(X4 | X1,X2,X3).
Use Extra Sums of Squares to test a partial sets of coefficients 42 Test statistic is Ex: p-value=0.0083 It can be written using extra sum of squares When only 1 coefficient is considered in the test (as in this case), 9.44=(3.073)2; F*=t2!!!
Practice 43 在Surgical data中想知道影響存活時間(SVtime )的因素,將存活時間取自然對數。模式A內包 含clot、prog、enzyme與gender。若再加入性別 與prog、enzyme之交互作用,請問交互作用是 否應列入模式中考慮?
例子 2008年台灣各縣市與澎湖縣的流浪狗資料 各縣市自1999~2008年的流浪狗處理數字,以及 各縣市在2008年的其他指標數字 city 縣市名稱 farmArea 各縣市耕地佔總面積的比例 captured 流浪狗累積補抓數目 divorced 各縣市離婚者所佔的比例 adoptedR 被捕抓之流浪狗被認 養的比例 unemployed 各縣市失業率 killedR 被捕抓之流浪狗被安 樂死之比例 crimeR 各縣市每10萬人刑事案件數目 unknownR 被捕抓之流浪狗被狀 況不明之比例 oldR 各縣市老人福利金額佔年度支出的比 例 population 各縣市於2008年的人 口數 computerR 各縣市平均每100個家庭的電腦數目 graduate 各縣市研究所畢業者 的人數
例子 以adoptedR(認養比例)為應變數 先刪除兩個變數:city(縣市名稱)以及與應變數 adoptedR(認養比例)有高度相關的unknownR( 狀況不明比例)變數 將所有的解釋變數放進模型中 model0=lm(adoptedR~captured+killedR+population+ graduate+farmArea+divorced+unemployed+crimeR+ oldR+computerR, data=dogs2, x=T)
Stepwise regression逐步迴歸 方法一:使用step函數搭配AIC指標進行逐步迴 歸變數篩選 step語法: step(lm物件,direction=“both”,k=2) direction可以選的值為forward, backward, both,其 中both是指任何解釋變數被加入模型後,仍有可能 在稍候被刪除;或是被刪除後,仍有可能在稍候被 加入。 k=2是使用模型的AIC作為篩選標準,若k=log(n),n 為樣本數,則是使用BIC為準則。
Stepwise regression逐步迴歸 方法二:使用step函數搭配BIC指標,進行逐步 迴歸變數篩選 summary(step(model0, k=log(nrow(dogs)), method="both"))
Stepwise regression逐步迴歸 方法三:使用leaps套件的regsubsets函數 regsubsets語法: regsubsets(X,y, nbest=k1, nvmax=k2, method) 其中X為包含所有解釋變數的矩陣;y為應變數向量 ;nbest=k1指定所有解釋變數數目相同的候選模型 中,都要挑出k1個最佳模型;nvmax=k2指定候選 模型中最多包含k2個解釋變數。mothod選項的值可 以是”forward”, “backword”, “exhausitive”(即all possible所有可能)與”seqrep”
Stepwise regression逐步迴歸 (A)Forward selection library(leaps) out.forward=regsubsets(as.matrix(dogs2[- 2]),y=dogs2$adoptedR, nbest=1, method="forward") s.forwd=summary(out.forward) 將候選模型的R2(rsq), SSE(rss), R2(adj), Cp, BIC值列出 round(cbind(s.forwd$which, rsq=s.forwd$rsq, adjr2=s.forwd$adjr2, rss=s.forwd$rss, cp=s.forwd$cp, bic=s.forwd$bic),2)
Stepwise regression逐步迴歸 以上計算結果中,每一個橫列代表一個候選模 型,各直行底下的1代表該直行的解釋變數有出 現在某個橫列的候選模型中,0則表示沒有出現 。 依據最小BIC值-21.8,forward selection的最佳模 型為模型編號3:解釋變數包含:killedR, unemployed,與crimeR,其R2值約為77.53%
Stepwise regression逐步迴歸 (B)Backward selection out.backward=regsubsets(as.matrix(dogs2[- 2]),y=dogs2$adoptedR, nbest=1, method="backward") s.back=summary(out.backward) #將候選模型的R2(rsq), SSE(rss), R2(adj), Cp, BIC值列 出 round(cbind(s.back$which, rsq=s.back$rsq, adjr2=s.back$adjr2, rss=s.back$rss, cp=s.back$cp, bic=s.back$bic),2)
Stepwise regression逐步迴歸 依據最小BIC值-25.74,backward selection的最佳 模型為模型編號6:解釋變數包含:captured, graduate, farmArea, crimeR, oldR, computerR,其 R2值約為87.43%
All possible subset selection (c)所有可能模型選取法 All possible subset selection out.all=regsubsets(as.matrix(dogs2[- 2]),y=dogs2$adoptedR, nbest=1, method="exhaustive") s.all=summary(out.all) #將候選模型的R2(rsq), SSE(rss), R2(adj), Cp, BIC值列 出 round(cbind(s.all$which, rsq=s.all$rsq, adjr2=s.all$adjr2, rss=s.all$rss, cp=s.all$cp, bic=s.all$bic),2)
All possible subset selection 依據最小BIC值-25.74,all possible subset selection的最佳模型為模型編號6:解釋變數包 含:captured, graduate, farmArea, crimeR, oldR, computerR,其R2值約為87.43% 可以搭配identify函數畫出all possible 篩選法的 Cp圖,並即時點選最佳的模型。Cp圖的X座標 是各候選模型的迴歸係數數目,Y座標是相對的 Cp值。
All possible subset selection q=as.vector(rowSums(s.all$which)) #迴歸係數數目 q plot(q, s.all$cp, xlim=c(0,8),ylim=c(0,18)) abline(0, b=1) identify(q, s.all$cp) 可以看出候選模型編號4的Cp直最靠近45度斜線, 具有最佳的Cp值。
Stepwise regression逐步迴歸 (D)Sequential replacement逐次替換法 out.all=regsubsets(as.matrix(dogs2[- 2]),y=dogs2$adoptedR, nbest=1, method="seqrep") s.step=summary(out.step) #將候選模型的R2(rsq), SSE(rss), R2(adj), Cp, BIC值列 出 round(cbind(s.step$which, rsq=s.step$rsq, adjr2=s.step$adjr2, rss=s.step$rss, cp=s.step$cp, bic=s.step$bic),2)
Stepwise regression逐步迴歸 依據最小BIC值-24.06,逐次替換法選出的最佳 模型為模型編號7:解釋變數包含:captured, graduate, farmArea, divorced, crimeR, oldR, computerR,其R2值約為84%
Outliers-Leverage hii值可協助偵測相對於解釋變數x’s的離群值:若 hii的值大於2p/n,則第i個觀察值可能是離群值 ,其中p-1為解釋變數數目。 hatvalues函數可以算出模型的槓桿值 程式: 2*7/23 (hii=hatvalues(model1)) which(as.vector(hii)>2*7/23) 結果顯示,第22個觀察值可能是相對於解釋變 數的離群值。
Outliers-Cook’s Distance Cook’s D指標值di若大於1(Cook and Weisberg, 1982),則第i個觀察值可能是離群值。 cooks.distance函數可算出di值 程式: cooks.distance(model1) which(as.vector(cooks.distance(model1))>1) 從Cook’s D指標來看,沒有任何di值大於1,沒有 離群值。
Outliers-T化殘差值 程式 (root.MSE=summary(model1)$sigma) #MSE^0.5 hii=hatvalues(model1) student.residual=resid/(root.MSE*sqrt(1-hii)) which(as.vector(abs(student.residual))>2.5) 從殘差值可知,第1、3、9、16、20個觀察值 的t化殘差大於2.5,有可能是離群值。
Outliers-T化殘差值 也可畫出t化殘差vs.Fitted values,並使用identify 函數來即時點選出可能的離群值 plot(fitted(model1),student.residual) abline(h=0) identify(fitted(model1),student.residual)
Outliers-T化去點殘差值 若沒有離群值存在,則此一統計量會服從t(n-p-1)分佈 (p=迴歸係數數目,包含beta0) 流浪狗資料,n=23, p=7, 因此n-p-1=23-7-1=15 若t化去點殘差絕對值大於t(0.95,15)則該觀察值可能是 離群值。 此外,由於離群值的偵測往往是一次針對所有觀察值來 檢查,因此一般建議用Bonferroni校正的α*值: α*/2= α/2n n為樣本數。 程式 jackknife.residual=rstudent(model1)
Outliers-T化去點殘差值 (1)先用α=0.1,不做Bonferronni校正: qt(0.95,15) which(as.vector(abs(jackknife.residual))>qt(0.95,15)) 指出第9、15、17、20個觀察值可能是離群值。 (2)採用Bonferronni校正 0.1/(2*23) qt(0.1/(2*23),15,lower.tail=F) which(as.vector(abs(jackknife.residual))>3.354188) 採用Bonferronni校正後,沒有任何殘差值超過界限,但要 注意的是, α值經過校正後,由於與樣本數n成反比,得 出的查表值會更大,導致偵測結果趨於保守,離群值需 有很大的殘差值才會被偵測出來。
Outliers-T化去點殘差值 畫出t化去點殘差vs.Fitted values,並使用identify 函數即時點選出可能的離群值 plot(fitted(model1),jackknife.residual) abline(h=0) identify(fitted(model1),jackknife.residual)
Influential observations 影響點 影響點的偵測可由槓桿值(leverages)、Cook’s D、DFBETAS、DFFITS等指標偵測出來 Leverages(hii) 2*7/23 (hii=hatvalues(model1)) which(as.vector(hii)>2*7/23) 結果顯示,第22個觀察值可能是影響點。
Influential observations 影響點 Cook’s D Cook’s D衡量每一個觀察值被移除後,對於迴歸係數估 計值的影響是否顯著。 通常以F(0.5, p, n-p)查表值當作比較的門檻值 cooks.distance(model1) qf(0.5,7,16) which(as.vector(cooks.distance(model1)))>qf(0.5,7,16) 所有觀察值的Cook’s D都沒有超過0.9457994,所以沒有 相對於迴歸係數估計值變化的影響點。
Influential observations 影響點 DFBETAS與Cook’s D一樣是以迴歸係數估計變化 量的大小當作影響點的偵測指標,若DFBETAS值 大於2,則相對的觀察值可能是影響點 dfbetas(model1) which(dfbetas(model1)>2) DFBETAS指標找不到任何影響點。
Influential observations 影響點 DFFITS衡量觀察值被移除後,對於應變數估計值 的影響,若DFFITS值大於2,則相對的觀察值可 能是影響點 dffits(model1) which(dffits(model1)>2) 第20個觀察值為影響點
Influential observations 影響點 可用滑鼠在圖上點出紅色圓,圈選所屬的觀察 直編號 influencePlot(model1) 點選出15、20兩個影響點
共線性 解釋變數之間如果存在嚴重的共線性問題,則某些解釋 變數的VIF值應該會很大。 一般判斷標準是若VIF值大於10,則可能有共線性問題 。可用car套件的VIF函數來計算模型各解釋變數的VIF library(car) vif(model1) mean(vif(model1)) 從VIF值發現,雖然解釋變數graduate, oldR, computerR的 VIF值超過平均VIF,但是沒有任何解釋變數的VIF值超過 10,故這個模型的共線性不是很嚴重。