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

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

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

plot(base,pch=19,ylim=c(30,180))
abline(lm(y~x,data=base),col="red")

相关视频

回归线（最大程度地减少误差平方和）是红色曲线

Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept)  50.7225    39.3979   1.287     0.23
x             0.4867     0.2747   1.772     0.11

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

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

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

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

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

---

点击标题查阅往期内容

R语言中的多项式回归、局部回归、核平滑和平滑样条回归模型

左右滑动查看更多

01

02

03

04

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

Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.255e+02  5.589e+01  -5.824 0.000394 ***
x            6.569e+00  8.744e-01   7.513 6.84e-05 ***
I(x^2)      -2.203e-02  3.143e-03  -7.011 0.000111 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

现在我们可以证明，对于形式为y =β2x2 +β1x +β0的抛物线，最优值以x⋆=θ= −β1 /2β2来获得。那么θ的自然估计量就是θ= −β 1 /2β2，通过最小化误差平方和。但是如何获得该估计量的方差？通过考虑以下因素，我们自然可以尝试Delta方法

（渐近）方差在这里

是

然后我们通过用未知量的估计值替换未知值来获得此渐近方差的估计量。

theta=-beta[2]/(2*beta[3])
theta
[1]  149.0676
s2=t(dg) %*% sigma %*% dg
s2[,1]
[1,] 94.59911
sqrt(s2)[,1]
[1,] 9.726207

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

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

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

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

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

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

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

optim(par=c(-213,3.5,110,2),logL)
$par
[1] -325.5    6.5  149.0   13.6
$value
[1] 44.3

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

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

我们可以绘制结果

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

在这里达到最大值

opt
$minimum
[1] 149.068

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

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

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

points(vx,vy,pch=19,cex=1.5,col="red")
arrows(min(v1[id]),vy,max(v1[id]),vy,code=3,angle=90,

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

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

loess.smooth(x = newbase$x, y= newbase$y
lines(reg$x,reg$y
for(i in 1:20) simu(TRUE)
lines(loess.smooth(x = base$x, y= base$y, evaluation = 501)

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

hist(V,probability = TRUE
lines(density(V)

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

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

---

---

点击文末“阅读原文”

获取全文完整代码数据资料。

本文选自《R语言多项式线性模型：最大似然估计二次曲线》。

点击标题查阅往期内容

【视频】什么是非线性模型与R语言多项式回归、局部平滑样条、 广义相加GAM分析工资数据|数据分享

MATLAB最小二乘法：线性最小二乘、加权线性最小二乘、稳健最小二乘、非线性最小二乘与剔除异常值效果比较

数据分享|R语言广义线性模型GLM：线性最小二乘、对数变换、泊松、二项式逻辑回归分析冰淇淋销售时间序列数据和模拟

生态学模拟对广义线性混合模型GLMM进行功率（功效、效能、效力）分析power analysis环境监测数据

广义线性模型glm泊松回归的lasso、弹性网络分类预测学生考试成绩数据和交叉验证

有限混合模型聚类FMM、广义线性回归模型GLM混合应用分析威士忌市场和研究专利申请数据

R语言贝叶斯广义线性混合（多层次/水平/嵌套）模型GLMM、逻辑回归分析教育留级影响因素数据

R语言贝叶斯MCMC：GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例

R语言用lme4多层次（混合效应）广义线性模型（GLM），逻辑回归分析教育留级调查数据

R语言广义线性模型GLM、多项式回归和广义可加模型GAM预测泰坦尼克号幸存者

R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）

R语言使用bootstrap和增量法计算广义线性模型（GLM）预测置信区间

R语言广义线性模型(GLMs)算法和零膨胀模型分析

R语言中广义线性模型(GLM)中的分布和连接函数分析

R语言中GLM(广义线性模型)，非线性和异方差可视化分析

R语言中的广义线性模型（GLM）和广义相加模型（GAM）：多元（平滑）回归分析保险资金投资组合信用风险敞口

用广义加性模型GAM进行时间序列分析

R和Python机器学习:广义线性回归glm，样条glm，梯度增强，随机森林和深度学习模型分析

在r语言中使用GAM（广义相加模型）进行电力负荷时间序列分析

用广义加性模型GAM进行时间序列分析

R和Python机器学习:广义线性回归glm，样条glm，梯度增强，随机森林和深度学习模型分析

在r语言中使用GAM（广义相加模型）进行电力负荷时间序列分析