6
$\begingroup$

How to estimate multiplicative model of spice harvesting?

Why do output coefficients not resemble true coefficients in a linear model?

This is a story about generating values by function, noising it and then trying to estimate the parameters of the function. The background of my story takes place at Dune desert.

I have been observing 80 workers harvesting spice on the Dune desert. I wrote down three characteristics of each worker and the volume of harvest each worker brought to granary throughout the whole season. These three characteristics were Home, Sex and army Rank.

All of a sudden, there was a sandstorm and I was captured by deity of spice. The deity took me to the middle of the desert and revealed to me the absolute knowledge about the load of harvest the worker brings to the granary. Relationship has the following form:

Harvest = Constant * Home * Sex * Rank * Noise

The deity revealed to me the true divine parameters of the phenomenon as well as the nature of Noise. The Noise is a random variable, normally distributed with mean 1, and standard deviation of 5%.

The deity warned me that if the workers are further oppressed by foremen to harvest more spice, it will send sandworms, which will put an end to spice harvesting.

Now the drama is that to stop worker exploitation I have to explain everything to the duke. But the duke does not believe in the existence the deity of spice. The duke thinks that the workers can harvest more output if oppressed harder. The only way to convince the duke is an econometric estimation which he values a lot.

  1. How can I estimate the parameters of this model?
  2. And how can I convince the duke that the relation is multiplicative and not additive?

Here is the track of the 80 workers data:

+------+-----------+--------+---------+
|  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  |
+------+-----------+--------+---------+

The true divine parameters revealed by deity:

+----------+-----------+-----------+
| 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 |
+----------+-----------+-----------+

The day two. After I received council and comments from the elders and sages of the statistician guild. I perused the path to estimate parameters with logarithms. I prepared the following initial matrix:

+-------+---------+-----------+----------+---------+---------+---------+--------+--------+
|       | 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   |
+-------+---------+-----------+----------+---------+---------+---------+--------+--------+

But the parameters calculated from such a matrix are not only worthless but misleading.

+----------+-----------+-----------+---+--------+---+
| 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 |   |
+----------+-----------+-----------+---+--------+---+

The Asterix sign shows the highest parameter within each variable. The highest values of divine parameters and estimated parameters do not match. For example Ordos has highest parameter within X1 variable. But estimated is Harkonen.

Dear elders and sages of the statistician guild, I turn to you for further council. Is the initial matrix crafted according to the art?

The day four

I submit manuscript in a file which reveals the divine process of harvest values creation. You can play with divine parameters and generate your own harvest values. The manuscript also contains my handcrafted procedure which automatically estimates the divine parameters. The procedure regains the parameters with great precision. There is also a sheet with my unsuccessful trial of estimating parameters through logarithms as advised by the elders of statistics guild.

dune.xlsx

Here are another two manuscripts demonstrating my procedure of regaining parameters with SQL code.

001 create table duna and insert values.sql

002 estymarka.sql

$\endgroup$
9
  • 2
    $\begingroup$ Your question is slightly ambiguous. Are you only asking how a multiplicative model can be estimated, or are you asking how the causal effects can be isolated from these observational data? (Explicitly, it seems the former, but the story implies you want to provide causal information to the Duke to convince him.) Furthermore, in the true data generating process, where is the error? Is it multiplicative as well, or not (see: How to tell the difference between linear and non-linear regression models?)? $\endgroup$ – gung - Reinstate Monica Sep 16 '19 at 16:21
  • 2
    $\begingroup$ This is not a multiplicative model at all: upon taking logarithms, it's a standard additive model with Normal errors. $\endgroup$ – whuber Sep 16 '19 at 17:03
  • 3
    $\begingroup$ You shouldn't have any problems with the deity. As @whuber correctly notes, you take logs and you have a bog-standard additive linear model. Given that that is true, & that the deity knows this well (by virtue of being omniscient), the deity will not give you any grief for providing the correct answer. $\endgroup$ – gung - Reinstate Monica Sep 16 '19 at 18:20
  • 2
    $\begingroup$ The analysis you perform at "worthless but misleading" appears to have nothing to do with your stated question, making the post confusing and occult. Could you restate what you mean by "crafted according to the art"? $\endgroup$ – whuber Sep 18 '19 at 17:50
  • 2
    $\begingroup$ I still have to read through your question but in any case +1 for Dune. $\endgroup$ – Sextus Empiricus Sep 19 '19 at 10:08
4
+50
$\begingroup$

The day three

The elders of the statistics guild have discovered a problem in the divine parameters. There is no single solution possible because the system is over-determined. We can scale the different values of the groups and the results will remain true.

For instance when we divide the 'constant' coefficient by two and at the same time multiply the 'home' coefficients by two then the result remains unchanged

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

The only divine values that matter are the ratio's of coefficients. We can see this in the divine R-code which only gives us $k-1$ coefficients for each characteristic variable of size $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
> 

Comparing multiplicative versus linear model

Note that there are different ways to make a 'multiplicative model/relation' (is the deterministic part multiplicative, the error term, or both). Sometimes people just compute a regular linear model for the logarithm of the response variable, however this is not the only way.

When we compare with a regular linear model (a deterministic linear part $X\beta$ and a homeogeneous error term $\epsilon \sim N(0,\sigma^2)$:

$$Y = X\beta + \epsilon$$

then modelling the logarithm of the response variable as such a linear model would look like this:

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

which can be transformed into:

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

and this is different from

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

So there is a difference whether the error term is inside or outside the logarithmic transformation.

Basically you can make the following four different models based on whether you use a logarithmic/multiplicative model and whether you assume the error term to be heterogeneous (include in the log transformation) or homogeneous (not included in the log transformation).

In R you can compute these four combinations by using:

# 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))

The difference in the use of homogeneous or heterogeneous error terms is not so much important for the predicted values. Except the error estimates for the predicted values will be different (and given the increase of residuals for larger values of $Y$ the use of heterogeneous error terms would not be so bad).

More important is that there is a difference whether you use a multiplicative model or an additive model:

example

you can see that the linear model (on the right) has a bias and low and high values will be underestimated and middle values will be overestimated.

You can also note that this curve on the right image seems to be like a nice function. And this might make you wonder whether there could not be some function to adapt it and make the fit better. For instance use $Y = f(X\beta) + \epsilon$. And indeed what you call 'multiplicative model' is just like an additive model with an exponential function.

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

with the relationships between the coefficients $d$ and $c$ as $d = exp(c)$

See the correspondence of the coefficients for the three different implementations of the model of the first category:

> 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 
$\endgroup$
13
  • $\begingroup$ Bravo! I have not mentioned it but the deity told me that it uses median of Y in place of constant. So alternatively we might assume that we model is Y/(medianY) = Home * Sex * Rank (with out constant). What you have estimated looks promising. $\endgroup$ – Przemyslaw Remin Sep 19 '19 at 11:02
  • 1
    $\begingroup$ The elders also forgot to mention something. It is not just about the constant coefficient, but also the parameter coefficients can be scaled relatively to each other. $$\text{Y/(medianY) = (Home/a/b) * (Sex*a) * (Rank*b)}$$ will hold for any values of $a$ and $b$. So the divine is not so easy to capture and we may never know what values the deity has been using. We will have to remain humble and just accept that ratio's are the only thing that we can capture and for the rest we need to trust whatever the deity tells us. $\endgroup$ – Sextus Empiricus Sep 19 '19 at 11:06
  • $\begingroup$ In addition to this problem with ratio's you have been expression your parameters without taking the exponent. When you perform a linear model with the logarithm $$\text{log(Y) = log(Home) + log(Sex) + log(Rank)}$$ then this means that you are computing the logarithm of the values $Home$, $Sex$ and $Rank$. And since you seem to be adding into it some minus sign (although I can't see exactly where that discrepancy arises) you end up with low values in your model when the divine values are high. $\endgroup$ – Sextus Empiricus Sep 19 '19 at 11:13
  • $\begingroup$ I got the lesson: "The only divine values that matter are the ratio's of coefficients". There is still shade of doubt. You have thrown out "one of a kind" of each variable. For example you have thrown out Atreides from Home variable. Does it mean the coefficient for Atreides is equal 1? $\endgroup$ – Przemyslaw Remin Sep 19 '19 at 11:17
  • $\begingroup$ It is arbitrary and depends on the setup of the computation. It is actually amazing how you got to compute your values. I wonder how you did that. It must have been under influence of quite some spice, since only the deity might help you with a computation that results into output like that. $\endgroup$ – Sextus Empiricus Sep 19 '19 at 11:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.