# Regression imputation of missing data based on OLS effects

Let's say we have a two-way with interaction experiment with missing data. Being the dataset:

> dt
X1 X2 X3
1   1  1 58
2   1  1 45
3   1  1 47
4   1  2 61
5   1  2 65
6   2  1 31
7   2  1 35
8   2  1 29
9   2  2 43
10  2  2 51
11  2  2 49
12  3  1 45
13  3  1 55
14  3  2 31
15  3  2 37
16  3  2 37
17  4  1 78
18  4  1 83
19  4  2 36
20  4  2 34
21  4  2 34
> dput(dt)
structure(list(X1 = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L), .Label = c("1",
"2", "3", "4"), class = "factor"), X2 = structure(c(1L, 1L, 1L,
2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L,
2L, 2L), .Label = c("1", "2"), class = "factor"), X3 = structure(c(58,
45, 47, 61, 65, 31, 35, 29, 43, 51, 49, 45, 55, 31, 37, 37, 78,
83, 36, 34, 34), .Dim = c(21L, 1L))), .Names = c("X1", "X2",
"X3"), row.names = c(NA, -21L), class = "data.frame")


I then perform an analysis of variance and estimate the means by the least squares method (lsmeans). So far so good...

But I have a quantitative X1, and want to do a cubic regression analysis for this factor. The problems begin:

• If I treat this case as a complete case scenario, the regression coefficients do not match the lsmeans means (Since I have four levels of the factor, the $R^2$ of the fitted model should be 1, with the curve containing the four means along its extension). In other words, the curve does not pass through the four lsmeans (one for each level of X1).
• If I input the missing data in the dataset based on the effects estimated by OLS (using the overparametrized design matrix, which gives me the mean effects of each factor and interaction) I can make the regression curve to match the lsmeans means. But now I have different sums of squares for the model and, most important, different degrees of freedom for the error estimate (due to the imputation of the missing values).

I am running this regression as a post-hoc analysis for the analysis of variance. What is the right approach? Won't the second one lead me to a type I error?

• What I am saying is that the fitted cubic curve passes through the four means (four levels of X1). Jul 5, 2015 at 1:21
• By definition, the $R^2$ is the $SSR/SST$, T standing for Treatment, since we are fitting the model to the means, trying to explain treatment effects with the model. In the case we have four treatment levels and fitting a cubic regression curve to those levels we have the sum of squares of the regression model equals to the treatment sum of squares, resulting in a $R^2=1$. Jul 5, 2015 at 1:38
• I understand what you are saying. This is a post-hoc analysis, made after the ANOVA. I am fitting to the means, testing the model SS with the treatments'. Besides that, all this has nothing to do with my question at all. You are now just trying to make an useless point to the subject of the question. Jul 5, 2015 at 1:53

Here's what I would do with these data.

First, I'd rename the variables to something sensible. This is statistics, not math, and even if the names I chose are wrong, they are a far cry better than X1 - X3. Also, create a numeric version of trt (nee X1) for later use

> names(dt) = c("trt", "grp", "resp")
> dt = transform(dt, trtnum = as.numeric(trt))


Fit an ANOVA model (including interaction), and obtain the predictions for each treatment, by group:

> dt.lm = lm(resp ~ trt * grp, data = dt)

> library(lsmeans)
> dt.lsm = lsmeans(dt.lm, ~ trt | grp)

> dt.lsm
grp = 1:
trt   lsmean       SE df lower.CL upper.CL
1   50.00000 2.524241 13 44.54671 55.45329
2   31.66667 2.524241 13 26.21337 37.11996
3   50.00000 3.091552 13 43.32111 56.67889
4   80.50000 3.091552 13 73.82111 87.17889

grp = 2:
trt   lsmean       SE df lower.CL upper.CL
1   63.00000 3.091552 13 56.32111 69.67889
2   47.66667 2.524241 13 42.21337 53.11996
3   35.00000 2.524241 13 29.54671 40.45329
4   34.66667 2.524241 13 29.21337 40.11996

Confidence level used: 0.95


Display these results graphically

> lsmip(dt.lsm, grp ~ trt)


We apparently have different trends for the two groups. Examine the polynomial contrasts

> contrast(dt.lsm, "poly")
grp = 1:
contrast    estimate        SE df t.ratio p.value
linear    109.833333 12.621207 13   8.702  <.0001
quadratic  48.833333  5.644375 13   8.652  <.0001
cubic     -24.500000 12.621207 13  -1.941  0.0742

grp = 2:
contrast    estimate        SE df t.ratio p.value
linear    -97.666667 12.494358 13  -7.817  <.0001
quadratic  15.000000  5.354725 13   2.801  0.0150
cubic       9.666667 11.428989 13   0.846  0.4130


These results suggest that we need a quadratic fit but may be able to do without the cubic trends. So let's try a quadratic model:

> dt.poly = lm(resp ~ (trtnum + I(trtnum^2)) * grp, data = dt)
> anova(dt.lm, dt.poly)
Analysis of Variance Table

Model 1: resp ~ trt * grp
Model 2: resp ~ (trtnum + I(trtnum^2)) * grp
Res.Df   RSS Df Sum of Sq      F Pr(>F)
1     13 248.5
2     15 334.2 -2   -85.705 2.2418 0.1457


The comparison of the two models is non-significant, so I'll choose the quadratic model as being more parsimonious. Here's the plot

> lsmip(dt.poly, grp ~ trtnum, at = list(trtnum = seq(1, 4, by = .5)))


Obtain the equations of the two parabolas

> summary(dt.poly)

Call:
lm(formula = resp ~ (trtnum + I(trtnum^2)) * grp, data = dt)

Residuals:
Min      1Q  Median      3Q     Max
-5.6067 -3.2520 -0.1951  1.2927  9.4100

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)       88.83000    7.99653  11.109 1.23e-08 ***
trtnum           -52.50833    7.49433  -7.006 4.23e-06 ***
I(trtnum^2)       12.69833    1.49887   8.472 4.22e-07 ***
grp2              -0.04951   11.82550  -0.004 0.996715
trtnum:grp2       23.62622   10.61939   2.225 0.041862 *
I(trtnum^2):grp2  -8.88939    2.08080  -4.272 0.000668 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.72 on 15 degrees of freedom
Multiple R-squared:  0.928, Adjusted R-squared:  0.9041
F-statistic: 38.69 on 5 and 15 DF,  p-value: 4.723e-08


Thus, the fitted parabolas are:

• Group 1: $\hat y = 88.83 - 52.51 t + 12.70 t^2$
• Group 2: $\hat y = (88.83 - .05) + (-52.51 + 23.63) t + (12.70 - 8.89) t^2 = 88.78 - 31.89 t + 3.81 t^2$
• Thank you for your post. Unfortunately it still does not answer my question. You made your decision for the quadratic curve based on the command contrast(dt.lsm, "poly"). But then fitted the curve considering it a complete case scenario using lm(resp ~ (trtnum + I(trtnum^2)) * grp, data = dt). The estimates of the two commands are not the same, and that is exactly my question. Jul 6, 2015 at 20:37
• The cubic trend will pass through the 4 means, but not through the 4 lsmeans... Jul 6, 2015 at 20:39
• I see that the data are imbalanced, but there are no incomplete cases. All 8 factor combinations have data. If that weren't true, those two quadratic curves would not be estimable. The model with interaction is mandatory, as there is a very strong interaction between the two factors. So any reasonable analysis would include 8 lsmeans, not 4; nor does a cubic trend interpolating the 4 marginal least-squares means. I don't think I can in good conscience offer the answers you want; sorry. Jul 7, 2015 at 1:07