CODEWITHSUNDEEP
rconstraintslinearmodelsminitab
spencergw
spencergw

Reputation: 299

Constrain linear model in R so that output matches with MINITAB

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?

Upvotes: 1

Views: 69

Answers (1)

Dave2e
Dave2e

Reputation: 24079

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)

Upvotes: 1

Related Questions