# Single or Multiple Models for Prediction with Categorical Variable

For a project I need to build a model that can predict batch processing time for the steps in a manufacturing process that are performed on a specific tool type. I know that the processing time is a function of the batch size (number of individual pieces of product loaded into the equipment at once) and the step in question.

I assume that the relationship between processing time and batch size is linear and the data indicate that the relationship is quite different based on which step is being run. In a stylized example where there were two steps run on the equipment under investigation, I've assumed an accurate way to describe the relationship between processing time and batch size is $$\mathbf{y=\delta_{0}+\delta_{1}x_{1}+\epsilon}$$ for step 1 and $$\mathbf{y=\alpha_{0}+\alpha_{1}x_{1}+\epsilon}$$ for step 2, where $\mathbf{\delta}\neq\mathbf{\alpha}$. Segregating the data by step and creating a unique predictive model for each step is less efficient than using an indicator variable to represent the step (with the actual data I have more like 10 steps). Using a model that treated the step as an indicator variable, I would need to specify a model $$\mathbf{y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{1}x_{2}+\epsilon}$$ to accurately capture both models stated above within one model that treated the step as an independent variable.

I thought the larger model might produce predictions identical to the two smaller models but I was unsure because the parameters would be estimated with the 'full' data set in the general case and only a portion of the data set in the case of the smaller models. To test this I created a data set in R:

# create variables and define relationships
btchSize <- rnorm(100,12,3)
step <- vector("character",100)
step[1:50] <- "A"
step[51:100] <- "B"
step <- as.factor(step)
y <- vector("numeric",100)
y[1:50] <- 50+25*x1[1:50]+rnorm(50,0,5)
y[51:100] <- 100+15*x1[51:100]+rnorm(50,0,5)

# create linear models
fullDataLm <- lm(y~x1+x2+x1*x2)
aLm <- lm(y[1:50]~x1[1:50])
bLm <- lm(y[51:100]~x1[51:100])

summary(fullDataLm)
summary(aLm)
summary(bLm)


This code gives the following results:

Call:
lm(formula = y ~ x1 + x2 + x1 * x2)

Residuals:
Min      1Q  Median      3Q     Max
-35.332  -9.246   1.227  10.716  26.028

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  49.8611     4.8206  10.343  < 2e-16 ***
x1           24.8566     0.3481  71.410  < 2e-16 ***
x2B          54.9103     7.1373   7.693 1.26e-11 ***
x1:x2B      -10.1553     0.5396 -18.819  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14.07 on 96 degrees of freedom
Multiple R-squared:  0.987, Adjusted R-squared:  0.9866
F-statistic:  2427 on 3 and 96 DF,  p-value: < 2.2e-16

Call:
lm(formula = y[1:50] ~ x1[1:50])

Residuals:
Min      1Q  Median      3Q     Max
-35.332  -8.802  -0.608  11.255  24.448

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  49.8611     4.7467   10.50 4.94e-14 ***
x1[1:50]     24.8566     0.3428   72.52  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.85 on 48 degrees of freedom
Multiple R-squared:  0.991, Adjusted R-squared:  0.9908
F-statistic:  5259 on 1 and 48 DF,  p-value: < 2.2e-16

Call:
lm(formula = y[51:100] ~ x1[51:100])

Residuals:
Min      1Q  Median      3Q     Max
-33.812 -10.049   3.803  10.529  26.028

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 104.7714     5.3427   19.61   <2e-16 ***
x1[51:100]   14.7013     0.4186   35.12   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14.28 on 48 degrees of freedom
Multiple R-squared:  0.9625,    Adjusted R-squared:  0.9618
F-statistic:  1234 on 1 and 48 DF,  p-value: < 2.2e-16


The results indicate that the larger model is identical to the smaller models in terms of the coefficient estimates. I.e. $$\mathbf{\beta_{0}=\delta_{0}} \\ \mathbf{\beta_{1}=\delta_{1}} \\ \mathbf{\beta_{0}+\beta_{2}=\alpha_{0}} \\ \mathbf{\beta_{1}+\beta_{3}=\alpha_{1}} \\$$ My questions are:

(1) Can anyone explain algebraically how the above is true considering that the parameter estimation equation is: $$\mathbf{\hat{\beta}=(X’X)^{-1}X’y}$$ and neither $\mathbf{X}$ nor $\mathbf{y}$ is the same in either of the three regressions. Obviously the indicator variables must cause the parameter estimation equations to end up being equal but I’d like to know how exactly.

(2) Does this result hold for parameter estimation techniques other than OLS?

(3) Which methodology is considered to be correct within the statistics community and why? Since the predictions from the larger model are identical to the predictions from the smaller models, I know that the residuals for the larger model are just the residuals of both smaller models (In my R example the residuals from the 100 observations used in the larger model are just the residuals from the 50 observations used in each of the smaller models). Consequently, if the regression assumptions are satisfied to a larger extent for data within one step then another it seems that the larger model may balance out the data and appear satisfactory in satisfying the regression assumption, while in fact the data lend themselves to regression analysis for one step but for another step are inconsistent with the use of regression.

• Seems like you want a Multilevel (eg Hierarchical) Model with the processing time nested in batch size and step. This should (in my limited experience with multilevel models) provide similar results to a multi-stage model, except compacted. As an engineer, I think the multilevel model would be vastly more preferred, as long as you can explain it in layman's terms to the stakeholders. I will leave your question for someone else more knowledgable to answer; just trying to get you on your way. PS- don't confuse Hierarchical Model with Hierarchical Regression... the latter is not of use to you. – TLJ Oct 31 '13 at 22:55

You asked for an algebraic answer to 1, and it is given below. However, a non-algebraic answer would be much better, so that is given first.

The usual technique used to run linear regressions is called Ordinary Least Squares (OLS). The formula you give above is the formula for the OLS estimator. However, the formula you give above is not the definition of the OLS estimator. The definition of the OLS estimator is "that coefficient vector, $\hat{\beta}$, which minimizes the sum of squared residuals." First, verbally . . .

Define the number of observations as $N$, so that each of the $\delta$ and $\alpha$ models have $N$ observations, and the $\beta$ model has $2N$. Consider the $\delta$ and $\alpha$ which solve the first two OLS problems you pose above, and ask, do they (transformed as you have transformed them) solve the third, pooled, $\beta$, OLS problem? Well, think about that third problem's objective function $\sum_{i=1}^{2N} (Y_i-\hat{Y}_i)^2$. We know that $\hat{\alpha}$ minimizes the first $N$ terms (definition of OLS) and has no effect on the last $N$ terms (by inspection). We know that $\hat{\delta}$ minimizes the last $N$ terms (definition of OLS) and has no effect on the first $N$ terms (by inspection). Last, there are "enough" $\beta$ such that the $\alpha$ and $\delta$ may be adjusted independently in the third problem. That is essentially a proof that the OLS estimators have to be exactly the same (again up to your transformation). The proof is not completely formal and explicit, but it's not hard from here to make it so.

To your question 2, yes this result is true much more generally. For any estimator defined by an optimization problem (any "M Estimator," like maximum likelihood for example), there is a result like this. The key requirements are that the parameter spaces of the pooled and separate models are "the same" up to an invertible transformation (like the one you gave to transform $\alpha$ and $\delta$ to $\beta$) and that the objective function of the pooled model can be decomposed (linearly or multiplicatively) into non-interacting parts corresponding (up to an increasing transformation) to the objective functions of the separate problems. This is a lot of models.

The kind of argument I give above is incredibly useful when you are using any kind of M-estimator. For example, it is a famous result for OLS that if you re-scale (change the units of) one of the $X$ variables, that its coefficient will be rescaled by OLS in an exactly offsetting way and that nothing else about the OLS estimator will change. This can be proved by a fairly tedious algebra exercise or by a very brief optimization argument like the one I gave above. The result is not just true for OLS, though. Any M estimator which has the $X$ multiplied against a coefficient will have this magic re-scaling property as well.

Now, to the algebraic demonstration. Changing your notation a little for clarity, start by comparing estimating these two equations by OLS: \begin{align} Y_1 &= X_1\delta + \epsilon_1 \\ Y_2 &= X_1\alpha + \epsilon_2 \end{align} With estimating this equation by OLS: \begin{align} \left(\begin{array}{r} Y_1 \\ Y_2 \end{array} \right) &= \left[\begin{array}{r r} X_1 & 0 \\ X_1 & X_1 \end{array} \right] \left(\begin{array}{r} \gamma_1 \\ \gamma_2 \end{array} \right) + \left(\begin{array}{r} \epsilon_1 \\ \epsilon_2 \end{array} \right) \\ \strut \\ Y &= X\gamma + \epsilon \end{align} My $\gamma_1$ is your $\beta_0$ and $\beta_1$ stacked up and etc. \begin{align} \hat{\gamma} &= (X'X)^{-1}X'Y \\ \strut \\ &= \left[\begin{array}{r r} X_1'X_1 & 0 \\ 2X_1'X_1 & X_1'X_1 \end{array} \right]^{-1} \left(\begin{array}{r} X_1'Y_1 \\ X_1'Y_1 + X_1'Y_2 \end{array} \right) \\ \strut \\ &= \left[\begin{array}{r r} (X_1'X_1)^{-1} & 0 \\ -2(X_1'X_1)^{-1} & (X_1'X_1)^{-1} \end{array} \right] \left(\begin{array}{r} X_1'Y_1 \\ X_1'Y_1 + X_1'Y_2 \end{array} \right) \\ \strut \\ \left( \begin{array}{r} \hat{\gamma}_1 \\ \hat{\gamma}_2 \end{array} \right) &= \left( \begin{array}{r} (X_1'X_1)^{-1}X_1'Y_1 \\ -2(X_1'X_1)^{-1}X_1'Y_1+(X_1'X_1)^{-1}X_1'Y_1 + (X_1'X_1)^{-1}X_1'Y_2 \end{array} \right) \\ \strut \\ \left( \begin{array}{r} \hat{\gamma}_1 \\ \hat{\gamma}_2 \end{array} \right) &= \left( \begin{array}{r} \hat{\delta}_1 \\ \hat{\alpha}_1 - \hat{\delta}_1 \end{array} \right) \end{align}

That's exactly what you got.

• Thank you for the thorough answer. Although I asked for an algebraic answer I really wanted to understand the result intuitively so I especially appreciated the logical argument alongside the algebra. – tjnel Nov 5 '13 at 17:26