The challenge

I would like to create a model that predicts the units sold for any temperature, even outside the range of available data.

I am particular interested how my models will behave in the more extreme cases when it is freezing outside, say the temperature dropped to 0ºC and the prediction for a very hot summer’s day at 35ºC.

Linear regression

I believe the observation \(y_i\) was drawn from a Normal (aka Gaussian) distribution with a mean \(\mu_i\), depending on the temperature \(x_i\) and a constant variance \(\sigma^2\) across all temperatures.

On another day with the same temperature I might have sold a different quantity of ice cream, but over many days with the same temperature the average units of sold ice cream would converge towards \(\mu_i\).

Thus, using a distribution centric notation my model looks like this: \[ \begin{aligned} y_i & \sim \mathcal{N}(\mu_i, \sigma^2), \\ \mathbb{E}[y_i] & =\mu_i = \alpha + \beta x_i \; \text{for all}\; i \end{aligned} \] Note, going forward I will drop the \(\text{for all}\; i\) statement for brevity.

Or, alternatively I might argue that the residuals, the difference between the observations and predicted values, follow a Gaussian distribution with a mean of 0 and variance \(\sigma^2\): \[ \varepsilon_i = y_i - \mu_i \sim \mathcal{N}(0, \sigma^2) \] Furthermore, the equation \[ \mathbb{E}[y_i]=\mu_i = \alpha + \beta x_i \] states my believe that the expected value of \(y_i\) is identical to the parameter \(\mu_i\) of the underlying distribution, while the variance is constant.

The same model written in the classic error terms convention looks like this: \[ \begin{aligned} y_i & = \alpha + \beta x_i + \varepsilon_i\\ \varepsilon_i & \sim \mathcal{N}(0, \sigma^2) \end{aligned} \]

I think the probability distribution centric convention makes it clearer that my observation is just one realisation from a distribution. Furthermore, it emphasises that the parameter of the distribution is modelled linearly.

To model this in R explicitly I use the glm function, specifying the response distribution as Gaussian and the link function from the expected value of the distribution to its parameter as identity. And this is really the trick with a GLMs, it describes the distribution of the observations and how its expected value, often after a smooth transformation, relates to the actual measurements (predictors) in a linear way. The ‘link’ is the inverse function of the original transformation of the data.

That’s all its says on the GLM tin and that all there is.

lin.mod <- glm(units ~ temp, data=icecream, 
              family=gaussian(link="identity"))
display(lin.mod)

## glm(formula = units ~ temp, family = gaussian(link = "identity"), 
##     data = icecream)
##             coef.est coef.se
## (Intercept) -159.47    54.64
## temp          30.09     2.87
## ---
##   n = 12, k = 2
##   residual deviance = 14536.3, null deviance = 174754.9 (difference = 160218.6)
##   overdispersion parameter = 1453.6
##   residual sd is sqrt(overdispersion) = 38.13

\[ \begin{aligned} y_i & \sim \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i & = -159.5 + 30.1 x_i, \; \sigma & = 38.1 \end{aligned} \]

Although the linear model looks fine in the range of temperature observed, it doesn’t make much sense at 0ºC.

The intercept is at -159, which would mean that customers return on average 159 units of ice cream on a freezing day. Well, I don’t think so.

Log-transformed linear regression

Ok, perhaps I can transform the data first. Ideally I would like ensure that the transformed data has only positive values. The first transformation that comes to my mind in those cases is the logarithm.

So, let’s model the ice cream sales on a logarithmic scale. Thus my model changes to: \[ \begin{aligned} \log(y_i) & \sim \mathcal{N}(\mu_i, \sigma^2)\\ \mathbb{E}[\log(y_i)] & =\mu_i = \alpha + \beta x_i \end{aligned} \]

This model implies that I believe the sales follow a log-normal distribution, \(y_i \sim \log\mathcal{N}(\mu_i, \sigma^2)\), which means that I regard higher sales figures as more likely than lower sales figures, as the log-normal distribution is skewed to the right.

Although the model is still linear on a log-scale, I have to remember to transform the predictions back to the original scale (remember \(\mathbb{E}[\log(y_i)] \neq \log(\mathbb{E}[y_i])\)):

\[ \begin{aligned} y_i & \sim \log\mathcal{N}(\mu_i, \sigma^2)\\ \mathbb{E}[y_i] & = \exp(\mu_i + \sigma^2/2) \end{aligned} \]

log.lin.mod <- glm(log(units) ~ temp, data=icecream, 
              family=gaussian(link="identity"))
display(log.lin.mod)

## glm(formula = log(units) ~ temp, family = gaussian(link = "identity"), 
##     data = icecream)
##             coef.est coef.se
## (Intercept) 4.40     0.20   
## temp        0.08     0.01   
## ---
##   n = 12, k = 2
##   residual deviance = 0.2, null deviance = 1.4 (difference = 1.2)
##   overdispersion parameter = 0.0
##   residual sd is sqrt(overdispersion) = 0.14

log.lin.sig <- summary(log.lin.mod)$dispersion
log.lin.pred <- exp(predict(log.lin.mod) + 0.5 * log.lin.sig)
basicPlot()
lines(icecream$temp, log.lin.pred, col="red", lwd=2)
legend(x="topleft", bty="n", lwd=c(2,2), lty=c(NA,1),
       legend=c("observation", "log-transformed LM"),
         col=c("#00526D","red"),  pch=c(1,NA))

The plot looks a little better than the previous linear model and it predicts that I would sell on average 82 ice creams at a 0ºC:

exp(coef(log.lin.mod)[1])

## (Intercept) 
##    81.62131

Although this model makes a little more sense, it appears that is over-predicts sells at lower and higher temperatures.

Furthermore, there is another problem with this model and the previous linear model as well.

The assumed model distributions generate real numbers, but my sales statistics are units and hence always whole numbers. Although the average number of units sold is likely to be a real number any draw from the model distribution should be a whole number.

Poisson regression

The classic approach for count data is the Poisson distribution.

The Poisson distribution has only one parameter, here \(\mu_i\), which is also its expected value. The canonical link function for \(\mu_i\) is the logarithm, which means I have to apply the exponential function to the linear model to get back to the original scale. Here is my model:

\[ \begin{aligned} y_i & \sim \text{Poisson}(\mu_i)\\ \log(\mu_i) & = \alpha + \beta x_i\\ \mathbb{E}[y_i] & =\exp(\alpha + \beta x_i) \end{aligned} \]

Let me say this again, although the expected value of my observation is a real number, the Poisson distribution will generate only whole numbers, in line with the actual sales.

pois.mod <- glm(units ~ temp, data=icecream, 
              family=poisson(link="log"))
display(pois.mod)

## glm(formula = units ~ temp, family = poisson(link = "log"), data = icecream)
##             coef.est coef.se
## (Intercept) 4.54     0.08   
## temp        0.08     0.00   
## ---
##   n = 12, k = 2
##   residual deviance = 60.0, null deviance = 460.1 (difference = 400.1)

pois.pred <- predict(pois.mod, type="response")
basicPlot()
lines(icecream$temp, pois.pred, col="blue", lwd=2)
legend(x="topleft", bty="n", lwd=c(2,2), lty=c(NA,1),
       legend=c("observation", "Poisson (log) GLM"),
         col=c("#00526D","blue"),  pch=c(1,NA))

This looks pretty good. The interpretation of the coefficient should be clear now. I have to use the exp function to make prediction of sales for a given temperature. However, R will do this for me automatically, if I set in the predict statement above type="response".

From the coefficients I can read off that a 0ºC I am expected to sell \(\exp(4.45)=94\) ice creams and that with each one degree increase in temperature the sales are predicted to increase by \(\exp(0.076) - 1 = 7.9\%\).

Note, the exponential function turns the additive scale into a multiplicative one.

So far, so good. My model is in line with my observations. Additionally, it will not predict negative sales and if I would simulate from a Poisson distribution with a mean given by the above model I will always only get whole numbers back.

However, my model will also predict that I should expect to sell over 1000 ice creams if temperature reach 32 Celsius:

predict(pois.mod, newdata=data.frame(temp=32), type="response")

##        1 
## 1056.651

Perhaps the exponential growth in my model looks a little too good to be true. Indeed, I am pretty sure that my market will be saturated at around 800.

Binomial regression

Ok, let’s me think about the problem this way: I have 800 potential sales and I’d like to understand the proportion sold for a given temperature.

This suggests a binomial distribution for the number of successful sales out of 800. The key parameter for the binomial distribution is the probability of success, the probability that someone will buy my ice cream as a function of temperature.

Dividing my sales statistics by 800 would give me a first proxy for the probability of selling all ice cream.

Therefore I need an S-shape curve that maps the sales statistics into probabilities between 0 and 100%.

A canonical choice is the logistic function:

\[\text{logistic}(u) = logit^{-1}(u) = \frac{e^u}{e^u + 1} = \frac{1}{1 + e^{-u}}\]

With that my model can be described as:

\[ \begin{aligned} p_i & \sim \text{B}(n, \mu_i)\\ \mu_i & = \text{logit}^{-1}(\alpha + \beta x_i)\\ \mathbb{E}[p_i] & =\text{logit}(\alpha + \beta x_i) \end{aligned} \]

market.size <- 1000
icecream$opportunity <- market.size - icecream$units
bin.glm <- glm(cbind(units, opportunity) ~ temp, data=icecream, 
    family=binomial(link = "logit"))
display(bin.glm)

## glm(formula = cbind(units, opportunity) ~ temp, family = binomial(link = "logit"), 
##     data = icecream)
##             coef.est coef.se
## (Intercept) -2.85     0.10  
## temp         0.13     0.01  
## ---
##   n = 12, k = 2
##   residual deviance = 75.0, null deviance = 754.1 (difference = 679.1)

bin.pred <- predict(bin.glm, type="response")*market.size
basicPlot()
lines(icecream$temp, bin.pred, col="purple", lwd=2)
legend(x="topleft", bty="n", lwd=c(2,2), lty=c(NA,1),
       legend=c("observation", "Binomial (logit) GLM"),
         col=c("#00526D","purple"),  pch=c(1,NA))

The chart doesn’t look too bad at all!

As the temperature rises higher and higher this model will predict that sales will reach market saturation, while all the other models so far would predict higher and higher sales.

I can predict the sales at 0ºC and 35ºC using the inverse of the logistic function, which is given in R as plogis:

# Sales at 0 Celsius
plogis(coef(bin.glm)[1])*market.size

## (Intercept) 
##    54.78799

# Sales at 35 Celsius
plogis(coef(bin.glm)[1] +  coef(bin.glm)[2]*35)*market.size

## (Intercept) 
##    846.0892