原文鏈接：http://tecdat.cn/?p=27474?

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

考慮我們從實(shí)驗(yàn)、事件等中觀察到一些數(shù)據(jù) y 的情況。我們將觀察結(jié)果 y 解釋為某個(gè)隨機(jī)變量 Y 的實(shí)現(xiàn)：

統(tǒng)計(jì)模型是對(duì)未知參數(shù) θ 的 Y 分布的規(guī)范。通常，觀測(cè)值 y = (y1, . . . , yn) ∈ Rn 是一個(gè)向量，而 Y = (Y1, . . . ., Yn) 是一個(gè)隨機(jī)向量。在這種情況下，統(tǒng)計(jì)模型是 Y1 聯(lián)合分布的規(guī)范， . . , Yn 直到未知參數(shù) θ。

手機(jī)示例

觀察示例： yi = 學(xué)生 i 在講座的前 10 分鐘內(nèi)查看他們的手機(jī)。
模型：

在許多實(shí)驗(yàn)和情況下，觀測(cè)值 Y1, . . . , Yn 不具有相同的分布。Y1 的分布， 。. . , Yn 可能取決于非隨機(jī)量 x1, 。. . , xn 稱為協(xié)變量。

示例：矮小的學(xué)生在數(shù)學(xué)考試中取得更高的分?jǐn)?shù)嗎？

Yi = 平均分，xi = 身高
模型：

?

模型擬合

我們將描述模型擬合的過程如下：
1. 模型規(guī)范——指定觀測(cè)值 Y1, 的分布。. . , Yn 可達(dá)未知參數(shù)。

2. 模型未知參數(shù)的估計(jì)。
3. 推理——這涉及構(gòu)建置信區(qū)間和測(cè)試關(guān)于參數(shù)的假設(shè)。
4. 診斷——檢查模型與數(shù)據(jù)的擬合程度。
“理想”模型應(yīng)該
? 與觀察到的數(shù)據(jù)相當(dāng)吻合。
? 不包含不必要的參數(shù)。
? 易于解釋。

R中的示例

假設(shè)我們有由家庭燃?xì)庀牧亢推骄獠繙囟冉M成的數(shù)據(jù)（見下表）。外部溫度可以用來衡量家庭使用的燃?xì)饬繂幔?br>

我們將 Gas 作為因變量，Yi 和 Temp 作為協(xié)變量 xi。假設(shè)我們使用線性正態(tài)模型來解釋數(shù)據(jù)；其中 Yi 是獨(dú)立的 N(μi, σ2)，其中 μi = β1 + β2xi 對(duì)于 i = 1, . . . , 26. 對(duì)于這個(gè)模型，我們有

?

然后計(jì)算 MLE

使用以下命令。首先，我們合并觀察向量 Y 和設(shè)計(jì)矩陣 X：
?

> X <- cbind(1,dat$Temp)

可以使用以下命令找到：

1. > qr(t(X)%*%X)$rank

2. [1] 2

作為

有滿秩，我們可以計(jì)算它的逆：

> solve(t(X)%*%X)

?

MLE β =

可以計(jì)算如下：
?

> betahat

?

然后繪制模型擬合我們可以使用
?

1. 


2. > lines(x=xs,y=btaht[1]+xs*beaht[2])

此外，殘差平方和 (RSS) 為

1. 


2. > RSS <- t(ehat)%*%ehat

3. > RSS

?

最后，我們可以計(jì)算

:
?

1. > sg2ht <- RSS/(26-2)

2. > sg2at

?

診斷
系數(shù)

衡量線性模型擬合優(yōu)度的一種方法是檢查決定系數(shù)；我們現(xiàn)在解釋。在只有截距項(xiàng)的最簡(jiǎn)單模型中

?

我們有 RSS = ∑in=1(Yi - Y)2。具有更多參數(shù)和大型設(shè)計(jì)矩陣的較大模型將具有較小的 RSS。
?對(duì)于包含截距項(xiàng)的模型，線性模型質(zhì)量的度量是

這稱為決定系數(shù)或 R2 統(tǒng)計(jì)量。請(qǐng)注意，0 ≤ R2 ≤ 1 和 R2 = 1 對(duì)應(yīng)于“完美”模型。

異常值

異常值是不符合其余數(shù)據(jù)的一般模式的觀察結(jié)果。異常值可能由于數(shù)據(jù)記錄錯(cuò)誤、數(shù)據(jù)是兩個(gè)或多個(gè)總體的混合以及模型需要改進(jìn)等原因而發(fā)生。我們將假設(shè)一個(gè)滿秩的設(shè)計(jì)矩陣。

殘差圖

?

杠桿

我們可能對(duì)每個(gè)觀察結(jié)果對(duì)模型擬合的影響程度感興趣。例如，考慮殘差 e 其中
?

杠桿對(duì)應(yīng)于觀察殘差的方差。

庫克距離

衡量觀察影響的另一種方法是考慮它對(duì)估計(jì)量 β 的變化或影響。一種這樣的衡量標(biāo)準(zhǔn)是庫克的距離
?

是不使用第 i 個(gè)觀測(cè)值計(jì)算的估計(jì)量。經(jīng)驗(yàn)法則是查看 Ci 接近 1 的觀測(cè)值。

?庫克距離的解釋

我們現(xiàn)在更詳細(xì)地研究殘差、杠桿和庫克距離。考慮下面給出的 4 個(gè)人工數(shù)據(jù)集，每個(gè)數(shù)據(jù)集都安裝了正常的線性模型。此外，還顯示了殘差與杠桿的關(guān)系圖。

?紅色數(shù)據(jù)點(diǎn)是可疑數(shù)據(jù)點(diǎn)，即高杠桿、高殘差（絕對(duì)值）或兩者兼有。紅線是包括紅色數(shù)據(jù)點(diǎn)的擬合，黑線是不包括紅色數(shù)據(jù)點(diǎn)的擬合。

例子

表中的數(shù)據(jù)顯示了 20 名接受高碳水化合物飲食的男性糖尿病患者從復(fù)合碳水化合物中獲得的總卡路里百分比，以及他們的年齡、體重和卡路里作為蛋白質(zhì)的百分比。

?我們將碳水化合物值作為我們的反應(yīng) Yi，以年齡、體重和蛋白質(zhì)作為協(xié)變量。然后我們用 Yi ～ N(μi, σ2) 擬合正態(tài)線性模型，其中

?然后，要找到 β 的最大似然估計(jì)量，我們需要求解

:

1. > beta.hat <- solve(t(X)%*%X)%*%t(X)%*%y

2. > t(beta.hat)

可以計(jì)算方差的無偏估計(jì)量 RSS/(n - p)，并用于計(jì)算每個(gè)分量的標(biāo)準(zhǔn)差

.
?

> sqrt(diag(sig.sq.hat*solve(t(X)%*%X)))

?

殘差是

> summary(ehat)

?

R2 和它的調(diào)整版本 R2 系數(shù)是

1. 


2. > R2

> 1-(RS/(20-4))/(RS0/(20-1))

我們可以使用 R 中的 lm 函數(shù)檢查這些結(jié)果，如下所示：

> summary(mylm1)

?

?接下來，我們執(zhí)行單向方差分析測(cè)試

> anova(mylm1)

?

?注意輸出的形式——即 Sum Sq 是基于上述模型的差異，即

1. 


2. > sum(mym2$esiuls^2)-sum(mym3$reials^2)

最后，R 還可以為我們生成一些殘差圖，如下所示：

> plot(mylm2)

?

自測(cè)題

Q1) (a) In R fit the normal linear model with:

Based upon the summary of the model, do you think that the model fits the data well? Explain your reasoning using the values reported in the R summary
?(b) Perform a hypothesis test to ascertain whether or not to include the intercept term | use a 5% significance level. Include your code.
(c) Conduct a hypothesis test comparing the models:

E(Y ) = β1 against E(Y ) = β1 + β2x2 + β3x3 + β4x4
as a 5% level. Include your code
(d) By inspecting the leverages and residuals, identify any potential outliers. Name these data points by their index number. Give your reasoning as to why you believe these are potential outliers. You may present up to three plots if necessary.

Q2) We shall now consider a GLM with a Gamma response distribution.
(a) Show that a random variable Y where Y follows a Gamma distribution with probability density function:

?is a member of the exponential family | taking the form a(φ) = φ. State the canonical link function and the variance function in terms of the expected response and the dispersion parameter.

(b) Show that the deviance for a GLM with a Gamma response distribution is

(c) Rewrite (by \hand") the IWLS algorithm ?specifically for the Gamma response and using the link:
?

?This is called the inverse link function.

(d) Write the components of the total score U1; : : : ; Up and the Fisher information matrix for this model.
(e) Given the observations y, what is a sensible initial guess to begin the IWLS algorithm in general?
(f) Manually write an IWLS algorithm to fit a Gamma GLM using your data, mydat, using the inverse link and same linear predictor in Q1a). Use the deviance as the convergence criteria and initial guess of β as (0:5; 0:5; 0:5; 0:5). Present your code and along with your final estimate of β and final deviance.
(g) Based on your IWLS results, compute φbD and φbp and the estimates of var(βb2). In R fit the model again with a Gamma response i.e.
?

> glm(y~x1+x2+x3,family=Gamma,data=mydat)

Note the capital G in Gamma. Verify the results with your IWLS results.
(h) Give a prediction for the response given by the model for x1= 13, x2= 5 x3= 0:255
and give a 91% confidence interval for this prediction. Include your code.
(i) Perform a hypothesis test between this model and another model with the same link and response distribution but with linear predictor η where ηi = β1 + β2xi1 + β3xi2 for i = 1; : : : ; n:
Use a 5% significance level. You may use the deviance function here. Include your code.
(j) Using your IWLS results, manually compute the leverages of the observations for this model | present your code (but not the values) and plot the leverages against the observation index number.
(k) Proceed to investigate diagnostic plots for your Gamma GLM. Identify any potential outliers | give your reasoning. Remove the most suspicious data point
| you must remove 1 and only 1 | and refit the same model. Compare and comment on the change of the model with and without this data point | you may wish to refer to the relative change in the estimated coefficients. You may present up to three plots if necessary.
?

最受歡迎的見解

1.R語言多元Logistic邏輯回歸 應(yīng)用案例

2.面板平滑轉(zhuǎn)移回歸(PSTR)分析案例實(shí)現(xiàn)

3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）

4.R語言泊松Poisson回歸模型分析案例

5.R語言回歸中的Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)

6.r語言中對(duì)LASSO回歸，Ridge嶺回歸和Elastic Net模型實(shí)現(xiàn)

7.在R語言中實(shí)現(xiàn)Logistic邏輯回歸

8.python用線性回歸預(yù)測(cè)股票價(jià)格

9.R語言如何在生存分析與Cox回歸中計(jì)算IDI，NRI指標(biāo)