拓端tecdat|R语言多项式线性模型：最大似然估计二次曲线

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

 

“应用线性模型”中，我们打算将一种理论（线性模型理论）应用于具体案例。通常，我会介绍理论的主要观点：假设，主要结果，并进行示范来直观地解释。这里查看一个真实的案例研究，它包含真实数据，2400个观测值，34个变量。

 

这里只有11个观察值，一个简单的线性模型。让我们对这些数据进行线性回归

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

我们可以清楚地看到我们的曲线似乎是凹的，开始时增加，结束时减少，可以进行非参数平滑

1.  
    
2.  
   scatter.smooth(x, y,
3.  
   lpars = list(col = "red")

 

 

我们可以进一步回答“最大数目在哪里吗”，可以建议一个值，找到一个置信区间吗？

我们可以考虑一个二次模型，换句话说，我们的预测将是 抛物线。

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

 

我们可以看到，该模型不仅在视觉上看起来更加符合实际，如果我们看一看回归的结果，该模型也更好

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

现在我们可以证明，对于形式为y =β2x2 +β1x +β0的抛物线，最优值以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

换句话说，如果假设估计量的正态性，我们有以下置信区间

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

 

我们还可以尝试另一种策略 。对于我们的二次模型，我们通常在高斯假设下对数似然函数

1.  
   logL = function(pm){
2.  
    
3.  
   -sum(log(dnorm (base$y-(b0+b1*base$x+b2*base$x^2)) ,0,b3
4.  
   }

在这里，第一个方法是引入θ（其中抛物线的最大值是）作为模型参数之一-例如代替β2

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

如果我们寻找最大似然函数，我们得到

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

这与我们之前的计算是一致的。第二个方法是分析似然函数：我们说在多元参数中，一个比另一个更重要。在给定θ的情况下，其他函数才最大化。从技术上讲，参数θ的轮廓对数似然是

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

我们可以绘制结果

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

 

在这里达到最大值 

1.  
   opt
2.  
   $minimum
3.  
   [1] 149.068

这与我们的计算是一致的。然后，我们可以得到结果的似然比检验。

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

 

然后我们可以在初始图中绘制置信区间

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,

 

就像大多数统计技术一样，这些都是渐近结果，仅凭11个观察结果无法保证其有效性。

另一个解决方案是使用模拟。我们假设观测值是模型，并且是噪声。我们可以将非参数模型（局部平滑）作为模型并假设高斯噪声。为了生成其他样本，我们将观测值保存在x中，另一方面，对于y，我们将使用y +ε，其中ε将根据正态分布绘制

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)

 

给定数据的不对称性，我们再次使用非参数模型。并且我们通过数值计算最大值。我们重复10,000次。

1.  
    
2.  
   hist(V,probability = TRUE
3.  
   lines(density(V)

 

在这里，我们有在10,000个模拟样本上观察到的最大值的经验分布。我们可以通过经验分位数来获得置信区间

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

 

---

最受欢迎的见解

1.R语言多元Logistic逻辑回归 应用案例

2.面板平滑转移回归(PSTR)分析案例实现

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

4.R语言泊松Poisson回归模型分析案例

5.R语言回归中的Hosmer-Lemeshow拟合优度检验

6.r语言中对LASSO回归，Ridge岭回归和Elastic Net模型实现

7.在R语言中实现Logistic逻辑回归

8.python用线性回归预测股票价格

9.R语言如何在生存分析与Cox回归中计算IDI，NRI指标