机器算法验证 - 为什么输出系数不像线性模型中的真实系数？ - 吾爱随笔录

为什么输出系数不像线性模型中的真实系数？

机器算法验证自习估计

2022-04-08 14:45:32

如何估计香料收获的乘法模型？

为什么输出系数不像线性模型中的真实系数？

这是一个关于通过函数生成值，对其进行噪声处理，然后尝试估计函数参数的故事。我的故事背景发生在沙丘沙漠。

我一直在观察 80 名工人在沙丘沙漠收割香料。我写下了每个工人的三个特征，以及每个工人在整个季节带到粮仓的收获量。这三个特征是家、性和军衔。

突然，一场沙尘暴，我被香料之神俘虏了。神把我带到了沙漠的中央，向我揭示了关于工人带到粮仓的收成负荷的绝对知识。关系具有以下形式：

Harvest = Constant * Home * Sex * Rank * Noise

神向我揭示了这一现象的真正神圣参数以及噪音的本质。噪声是一个随机变量，正态分布，均值为 1，标准差为 5%。

神灵警告我，如果工头进一步压迫工人采收更多的香料，就会发沙虫，这将结束香料的采收。

现在的戏剧是，为了停止对工人的剥削，我必须向公爵解释一切。但公爵不相信香料之神的存在。公爵认为，如果工人受到更大的压迫，他们可以收获更多的产出。说服公爵的唯一方法是他非常重视的计量经济学估计。

我如何估计这个模型的参数？
我怎样才能让公爵相信这种关系是乘法而不是加法的？

这是 80 名工人数据的轨迹：

+------+-----------+--------+---------+
|  Y   |  X1_Home  | X2_Sex | X3_Rank |
+------+-----------+--------+---------+
| 2.82 | Ordos     | M      | Veteran |
| 1.02 | Atreides  | F      | Junior  |
| 4.15 | Ordos     | M      | Elite   |
| 2.78 | Horekonen | M      | Elite   |
| 2.07 | Ordos     | F      | Junior  |
| 3.20 | Ordos     | M      | Veteran |
| 2.16 | Horekonen | M      | Veteran |
| 2.25 | Ordos     | F      | Junior  |
| 2.48 | Ordos     | F      | Veteran |
| 1.10 | Atreides  | F      | Veteran |
| 1.21 | Atreides  | M      | Junior  |
| 1.61 | Horekonen | F      | Junior  |
| 1.07 | Atreides  | F      | Veteran |
| 1.06 | Atreides  | F      | Junior  |
| 1.74 | Atreides  | M      | Elite   |
| 1.14 | Atreides  | F      | Veteran |
| 3.41 | Ordos     | M      | Elite   |
| 1.41 | Atreides  | M      | Veteran |
| 2.59 | Ordos     | F      | Veteran |
| 1.98 | Horekonen | F      | Veteran |
| 2.01 | Horekonen | M      | Junior  |
| 2.98 | Ordos     | M      | Veteran |
| 4.18 | Ordos     | M      | Elite   |
| 1.04 | Atreides  | F      | Veteran |
| 2.77 | Horekonen | M      | Elite   |
| 1.88 | Horekonen | M      | Junior  |
| 2.11 | Horekonen | M      | Junior  |
| 1.47 | Atreides  | F      | Elite   |
| 1.15 | Atreides  | M      | Junior  |
| 1.69 | Atreides  | M      | Elite   |
| 1.47 | Horekonen | F      | Junior  |
| 2.15 | Horekonen | M      | Veteran |
| 1.28 | Atreides  | M      | Veteran |
| 1.91 | Horekonen | F      | Veteran |
| 2.23 | Ordos     | F      | Junior  |
| 2.50 | Horekonen | M      | Elite   |
| 1.75 | Horekonen | F      | Veteran |
| 2.22 | Horekonen | F      | Elite   |
| 2.88 | Ordos     | M      | Junior  |
| 1.62 | Atreides  | M      | Elite   |
| 1.67 | Horekonen | F      | Junior  |
| 2.43 | Ordos     | F      | Veteran |
| 0.92 | Atreides  | F      | Junior  |
| 2.01 | Horekonen | M      | Veteran |
| 1.09 | Atreides  | F      | Veteran |
| 2.12 | Ordos     | F      | Junior  |
| 3.29 | Ordos     | M      | Veteran |
| 2.17 | Horekonen | M      | Veteran |
| 3.17 | Ordos     | F      | Elite   |
| 2.83 | Ordos     | M      | Junior  |
| 1.81 | Atreides  | M      | Elite   |
| 3.20 | Ordos     | F      | Elite   |
| 1.91 | Horekonen | M      | Junior  |
| 0.92 | Atreides  | F      | Junior  |
| 2.32 | Horekonen | F      | Elite   |
| 1.60 | Atreides  | M      | Elite   |
| 1.52 | Atreides  | F      | Elite   |
| 2.40 | Horekonen | F      | Elite   |
| 1.47 | Atreides  | F      | Elite   |
| 1.51 | Horekonen | F      | Junior  |
| 2.58 | Ordos     | M      | Junior  |
| 1.25 | Atreides  | M      | Veteran |
| 2.22 | Horekonen | M      | Veteran |
| 1.22 | Atreides  | M      | Junior  |
| 1.20 | Atreides  | F      | Veteran |
| 1.30 | Atreides  | M      | Veteran |
| 2.50 | Ordos     | F      | Junior  |
| 2.23 | Ordos     | F      | Junior  |
| 3.98 | Ordos     | M      | Elite   |
| 2.26 | Horekonen | F      | Elite   |
| 3.16 | Ordos     | F      | Elite   |
| 1.25 | Atreides  | M      | Junior  |
| 2.20 | Ordos     | F      | Junior  |
| 3.81 | Ordos     | M      | Elite   |
| 1.24 | Atreides  | M      | Veteran |
| 1.66 | Horekonen | F      | Junior  |
| 2.28 | Ordos     | F      | Junior  |
| 2.84 | Ordos     | M      | Veteran |
| 1.01 | Atreides  | F      | Junior  |
| 1.23 | Atreides  | M      | Junior  |
+------+-----------+--------+---------+

神所揭示的真正神圣参数：

+----------+-----------+-----------+
| variable |   value   | parameter |
+----------+-----------+-----------+
| constant | constant  |       2.0 |
| X1_Home  | Ordos     |       1.4 |
| X1_Home  | Horekonen |       1.0 |
| X1_Home  | Atreides  |       0.6 |
| X2_Sex   | M         |       1.1 |
| X2_Sex   | F         |       0.9 |
| X3_Rank  | Elite     |       1.3 |
| X3_Rank  | Veteran   |       1.0 |
| X3_Rank  | Junior    |       0.9 |
+----------+-----------+-----------+

第二天。在我收到了统计学家公会长老圣贤的议事和意见后。我仔细阅读了用对数估计参数的路径。我准备了以下初始矩阵：

+-------+---------+-----------+----------+---------+---------+---------+--------+--------+
|       | X1_Home |  X1_Home  | X1_Home  | X1_Rank | X1_Rank | X1_Rank | X1_Sex | X1_Sex |
| ln(Y) | Ordos   | Horekonen | Atreides | Elite   | Veteran | Junior  | M      | F      |
+-------+---------+-----------+----------+---------+---------+---------+--------+--------+
| 1.04  | 2.72    | 1.00      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.02  | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 1.42  | 2.72    | 1.00      | 1.00     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 1.02  | 1.00    | 2.72      | 1.00     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.73  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 1.16  | 2.72    | 1.00      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.77  | 1.00    | 2.72      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.81  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.91  | 2.72    | 1.00      | 1.00     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.09  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.19  | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.48  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.06  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.05  | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.55  | 1.00    | 1.00      | 2.72     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.13  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 1.23  | 2.72    | 1.00      | 1.00     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.34  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.95  | 2.72    | 1.00      | 1.00     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.68  | 1.00    | 2.72      | 1.00     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.70  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 1.09  | 2.72    | 1.00      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 1.43  | 2.72    | 1.00      | 1.00     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.04  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 1.02  | 1.00    | 2.72      | 1.00     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.63  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.75  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.38  | 1.00    | 1.00      | 2.72     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 0.14  | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.53  | 1.00    | 1.00      | 2.72     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.39  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.76  | 1.00    | 2.72      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.25  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.65  | 1.00    | 2.72      | 1.00     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.80  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.92  | 1.00    | 2.72      | 1.00     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.56  | 1.00    | 2.72      | 1.00     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.80  | 1.00    | 2.72      | 1.00     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 1.06  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.48  | 1.00    | 1.00      | 2.72     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.51  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.89  | 2.72    | 1.00      | 1.00     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| -0.08 | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.70  | 1.00    | 2.72      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.08  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.75  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 1.19  | 2.72    | 1.00      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.78  | 1.00    | 2.72      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 1.15  | 2.72    | 1.00      | 1.00     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 1.04  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.59  | 1.00    | 1.00      | 2.72     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 1.16  | 2.72    | 1.00      | 1.00     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 0.65  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| -0.08 | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.84  | 1.00    | 2.72      | 1.00     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 0.47  | 1.00    | 1.00      | 2.72     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.42  | 1.00    | 1.00      | 2.72     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 0.87  | 1.00    | 2.72      | 1.00     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 0.38  | 1.00    | 1.00      | 2.72     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 0.41  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.95  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.23  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.80  | 1.00    | 2.72      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.20  | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.18  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 1.00   | 2.72   |
| 0.27  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.92  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.80  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 1.38  | 2.72    | 1.00      | 1.00     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.82  | 1.00    | 2.72      | 1.00     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 1.15  | 2.72    | 1.00      | 1.00     | 2.72    | 1.00    | 1.00    | 1.00   | 2.72   |
| 0.22  | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
| 0.79  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 1.34  | 2.72    | 1.00      | 1.00     | 2.72    | 1.00    | 1.00    | 2.72   | 1.00   |
| 0.22  | 1.00    | 1.00      | 2.72     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.51  | 1.00    | 2.72      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.82  | 2.72    | 1.00      | 1.00     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 1.04  | 2.72    | 1.00      | 1.00     | 1.00    | 2.72    | 1.00    | 2.72   | 1.00   |
| 0.01  | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 1.00   | 2.72   |
| 0.21  | 1.00    | 1.00      | 2.72     | 1.00    | 1.00    | 2.72    | 2.72   | 1.00   |
+-------+---------+-----------+----------+---------+---------+---------+--------+--------+

但是从这样的矩阵计算的参数不仅没有价值，而且具有误导性。

+----------+-----------+-----------+---+--------+---+
| constant |   value   | estimated |   | divine |   |
+----------+-----------+-----------+---+--------+---+
| constant | constant  |      1.00 |   |   2.00 |   |
| X1_Home  | Ordos     |      1.46 |   |   1.40 | * |
| X1_Home  | Horekonen |      1.62 | * |   1.00 |   |
| X1_Home  | Atreides  |      0.94 |   |   0.60 |   |
| X2_Sex   | M         |      0.82 |   |   1.10 | * |
| X2_Sex   | F         |      1.00 | * |   0.90 |   |
| X3_Rank  | Elite     |      1.00 |   |   1.30 | * |
| X3_Rank  | Veteran   |      1.16 | * |   1.00 |   |
| X3_Rank  | Junior    |      0.62 |   |   0.90 |   |
+----------+-----------+-----------+---+--------+---+

Asterix 符号显示每个变量中的最高参数。神圣参数和估计参数的最高值不匹配。例如鄂尔多斯在 X1 变量中具有最高参数。但估计是哈科宁。

尊敬的统计学家公会的长老和圣人，我请您进一步商议。最初的矩阵是根据艺术制作的吗？

第四天

我在一个文件中提交手稿，该文件揭示了创造收获价值的神圣过程。您可以使用神圣参数并生成自己的收获值。该手稿还包含我手工制作的程序，该程序会自动估计神圣参数。该过程以极高的精度重新获得参数。还有一张纸，上面写着我按照统计行会长老的建议通过对数估计参数的尝试失败了。

沙丘.xlsx

这是另外两篇手稿，展示了我使用 SQL 代码重新获得参数的过程。

001 创建表duna并插入values.sql

第002章

1个回答

第三天

统计公会的长老们发现了神参数的问题。没有单一的解决方案是可能的，因为系统是超定的。我们可以缩放组的不同值，结果将保持不变。

例如，当我们将“常数”系数除以 2 并同时将“家庭”系数乘以 2 时，结果保持不变

Harvest = (Constant/2) * (2*Home) * Sex * Rank * Noise = Constant * Home * Sex * Rank * Noise

唯一重要的神圣值是系数的比率。我们可以在神圣的 R 代码中看到这一点，它只给了我们 $k-1$ 每个大小特征变量的系数 $k$ .

> model <- lm(log(Y)~1+Home+Gender+Rank, data=dune)
> 
> c <- exp(coef(model))
> c
  (Intercept) HomeHorekonen     HomeOrdos       GenderM    RankJunior 
    1.4199111     1.6218721     2.2666767     1.2001998     0.7001083 
  RankVeteran 
    0.7786371 
> 
> #comparing gender
> c['GenderM']       #model
GenderM 
 1.2002 
> 1.1/0.9            #divine
[1] 1.222222
> 
> #comparing homes
> c['HomeHorekonen'] #model
HomeHorekonen 
     1.621872 
> 1/0.6              #divine
[1] 1.666667
> c['HomeOrdos']     #model
HomeOrdos 
 2.266677 
> 1.4/0.6            #divine
[1] 2.333333
> 
> #comparing rank
> c['RankJunior']   #model
RankJunior 
 0.7001083 
> 0.9/1.3           #divine
[1] 0.6923077
> c['RankVeteran']  #model
RankVeteran 
  0.7786371 
> 1/1.3             #divine
[1] 0.7692308
>

比较乘法与线性模型

请注意，有不同的方法可以制作“乘法模型/关系”（确定性部分是乘法的，误差项，或两者兼而有之）。有时人们只是为响应变量的对数计算一个规则的线性模型，但这不是唯一的方法。

当我们与常规线性模型（确定性线性部分 $X\beta$ 和同质误差项 $\epsilon \sim N(0,\sigma^2)$ ：

Y = X β + ϵ

$Y = X\beta + \epsilon$

然后将响应变量的对数建模为线性模型，如下所示：

l o g (Y) = X β + ϵ

$log(Y) = X\beta + \epsilon$

可以转化为：

Y = e^{X β + ϵ}

$Y = e^{X\beta + \epsilon}$

这与

Y = e^{X β} + ϵ

$Y = e^{X\beta} + \epsilon$

因此，误差项在对数变换内部还是外部存在差异。

基本上，您可以根据您是否使用对数/乘法模型以及您是否假设误差项是异构的（包括在对数转换中）或同质的（不包括在对数转换中）来制作以下四种不同的模型。

在 R 中，您可以使用以下方法计算这四种组合：

# copy from your data
dune <- structure(list(  Y =           c(2.82, 1.02, 4.15, 2.78, 2.07, 3.2, 2.16, 2.25, 2.48, 1.1, 1.21, 1.61, 1.07, 1.06, 1.74, 1.14, 3.41, 1.41, 2.59, 1.98, 2.01, 2.98, 4.18, 1.04, 2.77, 1.88, 2.11, 1.47, 1.15, 1.69, 1.47, 2.15, 1.28, 1.91, 2.23, 2.5, 1.75, 2.22, 2.88, 1.62, 1.67, 2.43, 0.92, 2.01, 1.09, 2.12, 3.29, 2.17, 3.17, 2.83, 1.81, 3.2, 1.91, 0.92, 2.32, 1.6, 1.52, 2.4, 1.47, 1.51, 2.58, 1.25, 2.22, 1.22, 1.2, 1.3, 2.5, 2.23, 3.98, 2.26, 3.16, 1.25, 2.2, 3.81, 1.24, 1.66, 2.28, 2.84, 1.01, 1.23), 
                      Home = structure(c(3L, 1L, 3L, 2L, 3L, 3L, 2L, 3L, 3L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 3L, 1L, 3L, 2L, 2L, 3L, 3L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 3L, 2L, 2L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 1L, 3L, 3L, 2L, 3L, 3L, 1L, 3L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 3L, 1L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 3L, 1L, 3L, 3L, 1L, 2L, 3L, 3L, 1L, 1L), .Label = c("Atreides", "Horekonen", "Ordos"), class = "factor"), 
                    Gender = structure(c(2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L), .Label = c("F", "M"), class = "factor"), 
                      Rank = structure(c(3L, 2L, 1L, 1L, 2L, 3L, 3L, 2L, 3L, 3L, 2L, 2L, 3L, 2L, 1L, 3L, 1L, 3L, 3L, 3L, 2L, 3L, 1L, 3L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 3L, 3L, 3L, 2L, 1L, 3L, 1L, 2L, 1L, 2L, 3L, 2L, 3L, 3L, 2L, 3L, 3L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 2L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 3L, 2L, 2L, 3L, 2L, 2L), .Label = c("Elite", "Junior", "Veteran"), class = "factor")), class = "data.frame", row.names = c(NA, -80L))

# stuff to compute the models in various ways
#
# logarithmic deterministic model 
# homogeneous error terms
modelglm <- glm(Y ~ 1+Home+Gender+Rank, family=gaussian(link="log"), data=dune)

modelnls1 <- nls(Y ~ a * c(1,b1,b2)[Home] * c(1,c1)[Gender] * c(1,d1,d2)[Rank], 
             start = c(a=1,b1=1,b2=1,c1=1,d1=1,d2=1), data=dune)

modelnls2 <- nls(Y ~ exp(a + c(0,b1,b2)[Home] + c(0,c1)[Gender] + c(0,d1,d2)[Rank]), 
                start = c(a=1,b1=1,b2=1,c1=1,d1=1,d2=1), data=dune)

# logarithmic deterministic model 
# heterogeneous error terms
modellm <- lm(log(Y)~1+Home+Gender+Rank, data=dune)

# linear deterministic model 
# heterogeneous error terms
modelquasi <- glm(Y~1+Home+Gender+Rank, family=quasi(link="identity", variance="mu"), data=dune)

# linear deterministic model 
# homogeneous error terms
modelind <- lm(Y~1+Home+Gender+Rank, data=dune)

#### stuff to create the plots below
    plot(exp(predict(modellm)),dune$Y, 
     ylab = "observed values", xlab="estimated mean",
     cex=0.7,pch=1,col=rgb(0,0,0,0.5),bg=rgb(0,0,0,0.5))
lines(c(0,10),c(0,10))
title(expression(Y == exp(X * beta + epsilon)))
plot((predict(modelind)),dune$Y, 
     ylab = "observed values", xlab="estimated mean",
     cex=0.7,pch=1,col=rgb(0,0,0,0.5),bg=rgb(0,0,0,0.5),log="")
lines(seq(0.1,10,0.1),seq(0.1,10,0.1))
title(expression(Y == X * beta + epsilon))

使用同质或异质误差项的差异对于预测值并不那么重要。除了预测值的误差估计会不同（并且考虑到较大值的残差增加 $Y$ 使用异构误差项不会那么糟糕）。

更重要的是，使用乘法模型还是加法模型是有区别的：

你可以看到线性模型（右边）有一个偏差，低值和高值会被低估，中间值会被高估。

您还可以注意到，右图中的这条曲线似乎是一个不错的函数。这可能会让你想知道是否不能有一些功能来适应它并使其更适合。例如使用 $Y = f(X\beta) + \epsilon$ . 实际上，您所说的“乘法模型”就像具有指数函数的加法模型。

Y = e^{c_{0} + c_{H o m e} + c_{S e x} + c_{r a n k}} + ϵ = d_{0} + d_{H o m e} + d_{S e x} + d_{r a n k} + ϵ

$Y = e^{c_0 + c_{Home} + c_{Sex} + c_{rank}} + \epsilon = d_0 + d_{Home} + d_{Sex} + d_{rank} + \epsilon$

与系数之间的关系 $d$ 和 $c$ 作为 $d = exp(c)$

查看第一类模型的三种不同实现的系数对应关系：

> coefficients(modelnls1)
        a        b1        b2        c1        d1        d2 
1.4118260 1.6168696 2.2726165 1.2082349 0.7037927 0.7817258 
> exp(coefficients(modelnls2))
        a        b1        b2        c1        d1        d2 
1.4118251 1.6168700 2.2726173 1.2082353 0.7037931 0.7817260 
> exp(coefficients(modelglm))
  (Intercept) HomeHorekonen     HomeOrdos       GenderM    RankJunior 
    1.4118257     1.6168696     2.2726169     1.2082349     0.7037929 
  RankVeteran 
    0.7817259

其它你可能感兴趣的问题

上一篇为什么只有当它们具有相同的域时才能添加两个随机变量？下一篇损失函数是否一定要添加到观察中？