Constrain linear model in R so that output matches with MINITAB

Question

I am trying to recreate in R an example analysis done in a textbook in MINITAB. Below is my code for saving the data and fitting the model.

joint <- c("beveled", "butt", "beveled", "butt", "beveled", "beveled",
           "lap", "beveled", "butt", "lap", "lap", "lap", "butt",
           "lap", "butt", "beveled")
wood <- c("oak", "pine", "walnut", "oak", "oak", "pine", "walnut",
          "walnut", "walnut", "oak", "oak", "pine", "pine", "pine",
          "walnut", "pine")
y <- c(1518, 829, 2571, 1169, 1927, 1348, 1489, 2443, 1263, 1295, 
       1561, 1000, 596, 859, 1029, 1207)

joint <- factor(joint, levels = c("lap", "beveled", "butt"))
wood <- factor(wood, levels = c("walnut", "oak", "pine"))

js_mod <- lm(y ~ joint*wood)

And here is the model summary

summary(js_mod)

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)             1489.0      169.7   8.774 5.03e-05 ***
jointbeveled            1018.0      207.8   4.898  0.00176 ** 
jointbutt               -343.0      207.8  -1.650  0.14289    
woodoak                  -61.0      207.8  -0.293  0.77767    
woodpine                -559.5      207.8  -2.692  0.03100 *  
jointbeveled:woodoak    -723.5      268.3  -2.696  0.03081 *  
jointbutt:woodoak         84.0      293.9   0.286  0.78333    
jointbeveled:woodpine   -670.0      268.3  -2.497  0.04118 *  
jointbutt:woodpine       126.0      268.3   0.470  0.65295    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 169.7 on 7 degrees of freedom
Multiple R-squared:  0.9548,    Adjusted R-squared:  0.9032 
F-statistic:  18.5 on 8 and 7 DF,  p-value: 0.0004713

Here is the MINITAB output given in the textbook

Term               Coef StDev T P
Constant        1375.67 44.22 31.11 0.000
   joint
beveled          460.00 59.63 7.71 0.000
butt            -366.50 63.95 -5.73 0.001
   wood
oak               64.17 63.95 1.00 0.349
pine            -402.50 59.63 -6.75 0.000
   joint* wood
beveled oak     -177.33 85.38 -2.08 0.076
beveled pine    -155.67 82.20 -1.89 0.100
butt oak          95.67 97.07 0.99 0.357
butt pine        105.83 85.38 1.24 0.255

The coefficients are different and I think this must be that a different constraint is used in MINITAB than is used in R. I just don't know what that constraint is or how to implement it in R. I know the analysis must be the same though because the results of the ANOVA tables (not show here) in both cases are the same.

Can anyone help me with that?

Yes. The least squares means are given in the book as well, and when I calculate them using the model in R, the values are the same. In the R output, the intercept is equal to the overall mean where joint is "lap" and wood is "walnut". I don't understand why that wouldn't be the case in MINITAB. — spencergw, May 25 '23 at 22:41
I have MiniTab on my work computer and I cannot reproduce the published results, but I am able to reproduce the R results in MiniTab. I am considering this an error in the book, can you provide a reference to title, author and edition? Maybe someone else would have experience. — Dave2e, May 26 '23 at 02:51
The book is Basic Engineering Data Collection and Analysis by Vardeman and Jobe. The data is on page 64 and the output on page 171. https://www.iastatedigitalpress.com/plugins/books/127/ — spencergw, May 26 '23 at 12:02

Dave2e · Answer 1 · 2023-05-26T14:45:08.610

I am sticking by my comment, that there is an error in the text. Running the following code, will reproduce the Least Squares Means for strength table at the bottom part of page 171.
Running the data in Minitab will also produce the coefficients matching R.

joint <- c("beveled", "butt", "beveled", "butt", "beveled", "beveled",
           "lap", "beveled", "butt", "lap", "lap", "lap", "butt",
           "lap", "butt", "beveled")
wood <- c("oak", "pine", "walnut", "oak", "oak", "pine", "walnut",
          "walnut", "walnut", "oak", "oak", "pine", "pine", "pine",
          "walnut", "pine")
y <- c(1518, 829, 2571, 1169, 1927, 1348, 1489, 2443, 1263, 1295, 
       1561, 1000, 596, 859, 1029, 1207)

joint <- factor(joint, levels = c("lap", "beveled", "butt"))
wood <- factor(wood, levels = c("walnut", "oak", "pine"))

df <- data.frame(joint, wood, y)
js_mod <- lm(y ~ joint*wood, data=df)
summary(js_mod)

library(emmeans)
emmeans(js_mod,  ~ joint, data=df)
emmeans(js_mod,  ~ wood, data=df)
emmeans(js_mod,  ~ joint*wood, data=df)

I had also wondered if it was an error. Thanks for validating me! — spencergw, May 26 '23 at 14:13

Constrain linear model in R so that output matches with MINITAB

1 Answers1