Fitting response surface using rsm package in R - Lack of fit test is missing

Question

I am trying to fit an rsm model to a data set with three factors, to try to find optimum parameters of a simple empirical environmental model for disease risk prediction, or influence of a change in each predictor. I am not really sure how to make a reproducible example so I will try to describe what I am doing.
My data looks like:

    rh_thresh temp_thres hours       auc
           93         11    11 0.6198718
           87          5    15 0.6256410

And is stored:

    Data are stored in coded form using these coding formulas ...
t ~ (temp_thres - 10)
rh ~ (rh_thresh - 90)
h ~ (hours - 12)

My independent variable is area under the curve of empirical ROC. Values are ranging from 0.6 to 0.82. I have fitted second order model using SO function from rsm package.
Problem is that I dont get lack of fit test. Summary of the model fit:

Call:
rsm(formula = auc ~ SO(t, rh, h), data = cd_data)

               Estimate  Std. Error  t value  Pr(>|t|)    
(Intercept)  0.73647659  0.00219900 334.9145 < 2.2e-16 ***
t            0.00465510  0.00059000   7.8900 1.127e-14 ***
rh          -0.01554212  0.00059000 -26.3424 < 2.2e-16 ***
h           -0.01402361  0.00059000 -23.7687 < 2.2e-16 ***
t:rh        -0.00144076  0.00015847  -9.0917 < 2.2e-16 ***
t:h         -0.00152461  0.00015847  -9.6208 < 2.2e-16 ***
rh:h        -0.00150469  0.00015847  -9.4950 < 2.2e-16 ***
t^2          0.00046880  0.00018059   2.5958  0.009628 ** 
rh^2        -0.00167358  0.00018059  -9.2671 < 2.2e-16 ***
h^2         -0.00190715  0.00018059 -10.5604 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.6838,    Adjusted R-squared:  0.6799 
F-statistic: 172.8 on 9 and 719 DF,  p-value: < 2.2e-16

Analysis of Variance Table

Response: auc
               Df  Sum Sq  Mean Sq F value    Pr(>F)
FO(t, rh, h)    3 0.88339 0.294462 361.896 < 2.2e-16
TWI(t, rh, h)   3 0.21592 0.071975  88.458 < 2.2e-16
PQ(t, rh, h)    3 0.16610 0.055367  68.046 < 2.2e-16
Residuals     719 0.58502 0.000814                  
Lack of fit   719 0.58502 0.000814                  
Pure error      0 0.00000                           

Stationary point of response surface:
         t         rh          h 
-7.3837482 -1.3846032 -0.1790263 

Stationary point in original units:
temp_thres  rh_thresh      hours 
  2.616252  88.615397  11.820974 

Eigenanalysis:
eigen() decomposition
$`values`
[1]  0.0007978626 -0.0010301864 -0.0028796105

$vectors
         [,1]        [,2]       [,3]
t   0.9540713 -0.02300125 -0.2986952
rh -0.2143994 -0.74880338 -0.6271574
h  -0.2092386  0.66239296 -0.7193433

I guess the problem is over-saturation, but how do I report this? I do not have replicates.

Answer 1

I don't see a terrible problem here. You don't have replications, and that is why there is no estimate of pure error or test of lack of fit. But a response-surface model is just a regression model, and people fit regression equations to unreplicated data all the time. They just use the residual error as the denominator of tests (as is done in the output that is shown).

A separate issue is whether it is a good model. You should do a lot of residual plots (see any regression textbook) to explore this. You might also try a third-order model to see if it makes a significant improvement to the fit. Something like:

mod3 = lm(auc ~ poly(r, rh, h, degree = 3), data = cd_data)
anova(mod3, mod2)

where mod2 is the model you summarize here. If significant, that'd be a similar finding to having a significant LOF test.