2
$\begingroup$

I'm studying DoE $2^{k-p}$ and I'm attempting to replicate an experiment using R software from a article that was run in Minitab, but the results are different.

Regression model presented in the article:

enter image description here

My results (code to set up is below):

summary(lm(plan.atualizado1))
# ...
#   Estimate Std. Error  t value Pr(>|t|)
#   (Intercept)  67.16781    0.08914  753.497  < 2e-16 ***
#   pH1         -28.70531    0.08914 -322.020  < 2e-16 ***
#   Temp1         6.18844    0.08914   69.423  < 2e-16 ***
#   Dose1         4.36156    0.08914   48.929  < 2e-16 ***
#   Conc1         1.88031    0.08914   21.094 4.20e-13 ***
#   Velo1        -1.73219    0.08914  -19.432 1.49e-12 ***
#   pH1:Temp1     3.99906    0.08914   44.862  < 2e-16 ***
#   pH1:Dose1     6.71344    0.08914   75.312  < 2e-16 ***
#   pH1:Conc1    -1.13031    0.08914  -12.680 9.22e-10 ***
#   pH1:Velo1    -3.28031    0.08914  -36.799  < 2e-16 ***
#   Temp1:Dose1   0.50719    0.08914    5.690 3.35e-05 ***
#   Temp1:Conc1  -2.33656    0.08914  -26.212 1.43e-14 ***
#   Temp1:Velo1   1.35094    0.08914   15.155 6.55e-11 ***
#   Dose1:Conc1   3.38656    0.08914   37.991  < 2e-16 ***
#   Dose1:Velo1   0.66156    0.08914    7.421 1.45e-06 ***
#   Conc1:Velo1  -0.83219    0.08914   -9.336 7.09e-08 ***

Here's my code to set up the data, fit the model, and

if(!require("FrF2")) install.packages("FrF2") ; library(FrF2)

plan.person = FrF2(nfactors = 5, resolution = 5, replications = 2,
               randomize = FALSE, factor.names = list(pH = c(3, 9),
                                                 Temp = c(5, 30),
                                                 Dose = c(0.05, 0.5),
                                                 Conc = c(50, 350),
                                                 Velo = c(100, 200)))

Oleo1 <- c(94.80, 23.40, 99.60, 38.00, 80.20, 27.00, 96.40, 60.00, 99.90, 9.10,
       98.40, 37.10, 99.77, 51.70, 97.40, 58.00, 95.10, 25.00, 99.50, 39.10,
       80.60, 28.00, 96.70, 61.20, 99.70, 10.00, 98.80, 37.40, 99.40, 52.00,
       97.70, 58.40)

plan.atualizado1 = add.response(design = plan.person, response = Oleo1)
summary(lm(plan.atualizado1))

pH_rg = seq(3,9,0.5)
Temperatura_rg = seq(5,30,length.out = length(pH_rg))
Dose_rg = seq(0.05,0.5,length.out = length(pH_rg))
Conc_rg = seq(50,350,length.out = length(pH_rg))
Velo_rg = seq(100,200,length.out = length(pH_rg))

# Regression model using coefficients of the t statistic
f = function(pH_rg,Temperatura_rg,Dose_rg,Conc_rg,Velo_rg) 67.16781 - 28.70531*pH_rg - 6.18844*Temperatura_rg + 4.36156*Dose_rg + 1.88031*Conc_rg -1.73219*Velo_rg

dados <- list(pH_rg = pH_rg, Temperatura_rg = Temperatura_rg,
          Dose_rg = Dose_rg, Conc_rg = Conc_rg,
          Velo_rg = Velo_rg)

grid <- expand.grid(dados)

grid[, "fit"] <- f(pH_rg, Temperatura_rg, Dose_rg, Conc_rg, Velo_rg)

View(grid)

Note that the T statistics have the same results except the standard error coefficient.

$\endgroup$
  • 1
    $\begingroup$ Saw this yesterday, and while it seems to be really complete, I think for me that's the problem -- I think it may have too much (gasp!) information. It sounds like part of what you include are actually the same (pareto/main effects plot); it would help if you 1) focused the question on that part that is different. Also, it would help me if 2) you included your output that didn't match, 3) you didn't use the FrF2 package, it doesn't seem to be needed to set up the data, and also 4) explain what you mean by "didactic" purposes, are these images actually from the data or not? $\endgroup$ – Aaron - Reinstate Monica Jun 3 '18 at 17:37
  • $\begingroup$ @Aaron I edited the topic, placing the images of the article in just one link, and add the information of the regression model generated in R. I used the FrF2 package and left the code available if someone wanted to run the model and check for possible adjustments. The final goal is to understand how the authors of the article assembled the equation of the regression model, since if assembled from the coefficients data of the statistic T, the estimated results are different. Thank you for the tips! $\endgroup$ – bbiasi Jun 3 '18 at 21:35
  • 2
    $\begingroup$ Can't you just dput(plan.atualizado1)? Also, what's the output from summary(lm(plan.atualizado1)), that's what you think is different, yes? $\endgroup$ – Aaron - Reinstate Monica Jun 4 '18 at 3:58
  • $\begingroup$ I did not understand why of the dput (). I noticed that the output from summary (lm (plan.atualizado1)) I have a standard error coefficient different from that found in the article. But this is not the goal, it's just a note. The goal is to build a regression equation for the model plan.atualizado1 . $\endgroup$ – bbiasi Jun 4 '18 at 4:12
  • 2
    $\begingroup$ dput mostly so we don't have to install a new package, but also to make the example as simple as possible; if you'd tried fitting the model without the package, you might have discovered it was making factors instead of leaving them as numbers. $\endgroup$ – Aaron - Reinstate Monica Jun 4 '18 at 13:21
3
$\begingroup$

Well, two things: 1) you need to fit the model without second order interactions, and 2) you need to fit the model with the actual numeric values, not factors.

d <- within(plan.atualizado1, {
  pH <- as.numeric(as.character(pH))
  Temp <- as.numeric(as.character(Temp))
  Dose <- as.numeric(as.character(Dose))
  Conc <- as.numeric(as.character(Conc))
  Velo <- as.numeric(as.character(Velo)) })
summary(lm(Oleo1 ~ (pH + Temp + Dose + Conc + Velo), d))

#                Estimate Std. Error t value  Pr(>|t|)    
#   (Intercept) 113.27331    8.31040  13.630 2.37e-13 ***
#   pH           -9.56844    0.63167 -15.148 2.05e-14 ***
#   Temp          0.49508    0.15160   3.266  0.00306 ** 
#   Dose         19.38472    8.42232   2.302  0.02962 *  
#   Conc          0.01254    0.01263   0.992  0.33023    
#   Velo         -0.03464    0.03790  -0.914  0.36908   

PS. Use predict with newdata to get the fit.

pH_rg = seq(3,9,0.5)
Temperatura_rg = seq(5,30,length.out = length(pH_rg))
Dose_rg = seq(0.05,0.5,length.out = length(pH_rg))
Conc_rg = seq(50,350,length.out = length(pH_rg))
Velo_rg = seq(100,200,length.out = length(pH_rg))
dados <- expand.grid(pH = pH_rg, Temp = Temperatura_rg,
              Dose = Dose_rg, Conc = Conc_rg,
              Velo = Velo_rg)
dados$p <- predict(m, newdata=dados)

In response to below:

1) My University doesn't have access to this journal, but I can say from experience that journal articles are notorious for including minimal statistical details. In some instances it really is lacking, in others, the article is following conventions in the field and you have to know what those are to understand.

2) It's not wrong here, it's just a different model, a first order (FO) only model instead of a first order (FO) plus two way interactions (TWI). The paper reports both. (Or of course, they're both wrong, but may be useful, as the famous quote reminds us.) You can see that the interaction table has the coefficients from the FO+TWI model, using +1/-1 coding (in the middle column), which matches your model, I'm not sure what the other two columns are.

3) I at first was thinking of this as a response surface, but it's actually a fractional factorial, yes? I'm not an expert in this; I learned the basics from Oehlert Chap 18. [Previously, I'd written: I'm not a response surface methods expert; I learned the basics from Oehlert Chap 19 and Lenth's rsm package vignettes. I hear Box/Draper is a standard reference, but I don't know it.]

4a) The equation from the paper only had the first order terms and didn't included the variables directly, without any standardization noted, so it seemed like the first thing to try.

4b) predict is better because you don't have to type in the coefficients, so it saves you effort, they're not rounded, and when your model changes, you don't have to change your function.

One final thought -- the within group variation is really small in this data set. Each combination of predictors was repeated twice, and look how similar the response values are compared with the spread of the data. Without knowing anything about the subject matter, this is a red flag, and I'd look for evidence of pseudoreplication. That is, I suspect there's only one true replicate for each combination, and therefore not enough data to fit a model with two-way interactions. Plus you're limited to a max of 100, so the assumptions of equal variance may not be met as well.

$\endgroup$
  • $\begingroup$ some questions. 1- Would you tell me why this is not clear in the article and only Table 1 is shown with the values of the interactions? 2- Have I visualized many experiments where the regression function is apneas lm (plan.atualizadox), the others are wrong? 3- Please, tell me a book that has this so I can study. 4- And what was your thought to achieve in this result? Is predict so much better than function? Thank you for your help! $\endgroup$ – bbiasi Jun 4 '18 at 13:30
  • $\begingroup$ **2- Have I visualized many experiments where the regression function is only lm (plan.atualizadox), the others are wrong? $\endgroup$ – bbiasi Jun 4 '18 at 13:38
  • $\begingroup$ Interesting approach, most books stop at the analysis of variance and here we realize that there may be much still to do. Could you tell me other examples that it was necessary to make more than one ANOVA? $\endgroup$ – bbiasi Jun 6 '18 at 2:16
  • 1
    $\begingroup$ Well, if you can specify your model appropriately ahead of time, it's (almost) always best to just fit that model and stop. See Frank Harrell (RMS book, and numerous questions here as well) on model selection for details. $\endgroup$ – Aaron - Reinstate Monica Jun 6 '18 at 14:05
  • 1
    $\begingroup$ Setting aside his very good advice though, it is sadly exceptionally common for researchers to fit multiple models to data resulting in multiple ANOVAs. I should say, though, that sometimes this does make sense; different models answer different questions, and nested models are a good way to test for significance of a number of terms simultaneously, such as FO compared to FO + TWI; that may be why they fit both here, though again, I'd be surprised if they had enough true data to fit the FO + TWI model correctly. $\endgroup$ – Aaron - Reinstate Monica Jun 6 '18 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.