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

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

?

“應(yīng)用線性模型”中，我們打算將一種理論（線性模型理論）應(yīng)用于具體案例。通常，我會介紹理論的主要觀點：假設(shè)，主要結(jié)果，并進行示范來直觀地解釋。這里查看一個真實的案例研究，它包含真實數(shù)據(jù)，2400個觀測值，34個變量。

?

這里只有11個觀察值，一個簡單的線性模型。讓我們對這些數(shù)據(jù)進行線性回歸

1. 


2. plot(base,pch=19,ylim=c(30,180))

3. abline(lm(y~x,data=base),col="red")

回歸線（最大程度地減少誤差平方和）是紅色曲線

?

1. 


2. Coefficients:

3. Estimate Std. Error t value Pr(>|t|)

4. (Intercept) ?50.7225 ? ?39.3979 ? 1.287 ? ? 0.23

5. x ? ? ? ? ? ? 0.4867 ? ? 0.2747 ? 1.772 ? ? 0.11

我們可以清楚地看到我們的曲線似乎是凹的，開始時增加，結(jié)束時減少，可以進行非參數(shù)平滑

1. 


2. scatter.smooth(x, y,

3. lpars = list(col = "red")

?

?

我們可以進一步回答“最大數(shù)目在哪里嗎”，可以建議一個值，找到一個置信區(qū)間嗎？

我們可以考慮一個二次模型，換句話說，我們的預測將是?拋物線。

1. 


2. lm(y~x+I(x^2),data=base)

?

我們可以看到，該模型不僅在視覺上看起來更加符合實際，如果我們看一看回歸的結(jié)果，該模型也更好

1. 


2. 


3. Coefficients:

4. Estimate Std. Error t value Pr(>|t|)

5. (Intercept) -3.255e+02 ?5.589e+01 ?-5.824 0.000394 ***

6. x ? ? ? ? ? ?6.569e+00 ?8.744e-01 ? 7.513 6.84e-05 ***

7. I(x^2) ? ? ?-2.203e-02 ?3.143e-03 ?-7.011 0.000111 ***

8. ---

9. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

現(xiàn)在我們可以證明，對于形式為y =β2x2 +β1x +β0的拋物線，最優(yōu)值以x?=θ= ?β1 /2β2來獲得。那么θ的自然估計量就是θ= ?β 1 /2β2，通過最小化誤差平方和。但是如何獲得該估計量的方差？通過考慮以下因素，我們自然可以嘗試Delta方法

?

（漸近）方差在這里

?

是

然后我們通過用未知量的估計值替換未知值來獲得此漸近方差的估計量。

1. 


2. theta=-beta[2]/(2*beta[3])

3. theta

4. [1] ?149.0676

5. s2=t(dg) %*% sigma %*% dg

6. s2

7. [,1]

8. [1,] 94.59911

9. sqrt(s2)

10. [,1]

11. [1,] 9.726207

換句話說，如果假設(shè)估計量的正態(tài)性，我們有以下置信區(qū)間

1. 


2. arrows(vx-qt(.975,n-3)*sqrt(s2),vy,

?

我們還可以嘗試另一種策略?。對于我們的二次模型，我們通常在高斯假設(shè)下對數(shù)似然函數(shù)

1. logL = function(pm){

2. 


3. -sum(log(dnorm (base$y-(b0+b1*base$x+b2*base$x^2)) ,0,b3

4. }

在這里，第一個方法是引入θ（其中拋物線的最大值是）作為模型參數(shù)之一-例如代替β2

1. logL = function(pm){

2. 


3. -sum(log(dnorm( base$y-b0+b1 base$x-.5*b1/theta base$x^2) ,0,b3

如果我們尋找最大似然函數(shù)，我們得到

1. optim(par=c(-213,3.5,110,2),logL)

2. $par

3. [1] -325.5 ? ?6.5 ?149.0 ? 13.6

4. 


5. $value

6. [1] 44.3

這與我們之前的計算是一致的。第二個方法是分析似然函數(shù)：我們說在多元參數(shù)中，一個比另一個更重要。在給定θ的情況下，其他函數(shù)才最大化。從技術(shù)上講，參數(shù)θ的輪廓對數(shù)似然是

1. logL=function(theta){

2. 


3. -sum(log(dnorm( base$y- b0+b1*base$x-.5*b1/theta*base$x^2 ,0,b3)

4. 


5. optim(par )$value

我們可以繪制結(jié)果

1. 


2. plot(v1,-v2,type="l",xlim=range(base$x)

?

在這里達到最大值?

1. opt

2. $minimum

3. [1] 149.068

這與我們的計算是一致的。然后，我們可以得到結(jié)果的似然比檢驗。

1. 


2. abline(h=ref,lty=2,col="red")

?

然后我們可以在初始圖中繪制置信區(qū)間

1. 


2. points(vx,vy,pch=19,cex=1.5,col="red")

3. arrows(min(v1[id]),vy,max(v1[id]),vy,code=3,angle=90,

?

就像大多數(shù)統(tǒng)計技術(shù)一樣，這些都是漸近結(jié)果，僅憑11個觀察結(jié)果無法保證其有效性。

另一個解決方案是使用模擬。我們假設(shè)觀測值是模型，并且是噪聲。我們可以將非參數(shù)模型（局部平滑）作為模型并假設(shè)高斯噪聲。為了生成其他樣本，我們將觀測值保存在x中，另一方面，對于y，我們將使用y +ε，其中ε將根據(jù)正態(tài)分布繪制

1. 


2. loess.smooth(x = newbase$x, y= newbase$y

3. lines(reg$x,reg$y

4. 


5. for(i in 1:20) simu(TRUE)

6. lines(loess.smooth(x = base$x, y= base$y, evaluation = 501)

?

給定數(shù)據(jù)的不對稱性，我們再次使用非參數(shù)模型。并且我們通過數(shù)值計算最大值。我們重復10,000次。

1. 


2. hist(V,probability = TRUE

3. lines(density(V)

?

在這里，我們有在10,000個模擬樣本上觀察到的最大值的經(jīng)驗分布。我們可以通過經(jīng)驗分位數(shù)來獲得置信區(qū)間

1. 


2. arrows(quantile(V,.025 ,vy,quantile(V,.975

?

最受歡迎的見解

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

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

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

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

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

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

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

8.python用線性回歸預測股票價格

9.R語言如何在生存分析與Cox回歸中計算IDI，NRI指標

?

?