I'm grappling with a question that I think I know the answer to, but not able to prove it mathematically.

I have two oils, which are each run once on an engine across ~800 seconds while the engine speeds are increasing over time. At each time point, a measurement of the aeration is taken, which is the response variable. The relationship between engine speed and the aeration is quite non-linear, and I want to fit a spline to model the relationship of engine speed and aeration, and look for differences across speeds in the samples.

As noted, each sample is only run once. From a traditional ANOVA point of view, one will argue that there is only "one replication" and we can't make any inferences about it (that's what my colleagues argue). However, if we are modeling the relationship between the response and the engine speed, aren't the measurments at each time point conditionally independent on the predictor variable?
Do we really need more than one run in this case to be able to make any inferences?

Also, to make it clear, this is not about power. My colleagues argue that since there is only one replication, we cannot make ANY inference.

Just to make it concrete what the data/experiment looks like, I've attached a plot of the data.

enter image description here

Follow up question

Confidence intervals vs Prediction Intervals

I'm modeling the above data like so: A global shared smooth + group level smooths + group level intercepts and in addition, modeling the variance as a function of group level smooths and group level intercepts.
Based on AIC, a model assuming constant variance has a much larger AIC than the one modeling variance as a function of engine speed and sample.

mod_gam4 <- gam(list(aeration_scale ~ s(speed_scale, bs = "cs") + s(speed_scale, bs = "cs", by = sample) +
                 s(sample, bs = "re"),~ s(speed_scale, bs = "cs", by = sample) + s(sample, bs = "re")), 
                data = dta_cut, method = "REML", family = gaulss())

Confidence Intervals

I use the predict.gam to get fitted values and the associated intervals. Looking at the data, there's quite a lot of uncertainty but the CI below seems too optimistic - this sounds unscientific, but I would have expected wider intervals given the amount of noise.

grid <- expand.grid(speed_scale = seq(-2,2,by = 0.005), 
                    sample = levels(dta_cut$sample))
pred_list <- predict.gam(mod_gam5, newdata = grid, unconditional = TRUE, 
                     se.fit = TRUE, type = "response")
preds <- pred_list$fit %>%
  data.frame() %>%
  janitor::clean_names() %>%
  rename(fit = x)

cis <- pred_list$se.fit %>%
  data.frame() %>%
  janitor::clean_names() %>%
  rename(se_fit = x)

bind_cols(preds, cis, grid) %>%
  mutate(lower = fit - 1.96*se_fit, upper = fit + 1.96*se_fit) %>%
  ggplot(.,aes(speed_scale, fit, group = sample)) + geom_line(aes(color = sample)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2)

95 % Confidence Intervals

Prediction intervals

Instead, I wanted to look at prediction intervals which incorpoarate the error of the responses. This seems to capture the wide uncertainty in responses.

Designmat = predict.gam(mod_gam4,newdata=grid,type="lpmatrix", unconditional = TRUE)
tmp <- Designmat %*% vcov(mod_gam4)

diagMult <- cmpfun(function(m1, m2) sapply(seq_len(nrow(m1)), 
                                            function(i) m1[i,] %*% m2[,i]))
predvar <-  diagMult(tmp, t(Designmat))
SE <- sqrt(predvar) 
SE2 <- sqrt(predvar+mod_gam4$sig2) 
tfrac <- qt(.975, mod_gam4$df.residual)
interval = tfrac*SE2 %>% data.frame() %>%clean_names()

bind_cols(preds, cis, grid, interval) %>%
  mutate(lower = fit - x, upper = fit + x) %>%
  ggplot(.,aes(speed_scale, fit, group = sample)) + geom_line(aes(color = sample)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2)

enter image description here

Question -

How does one prove the coverage of the fitted intervals? Simply looking at the data they seem overly optimistic.

  • $\begingroup$ This is quite like a time series, except that engine_speed is a surrogate for time. So no, I don't see any problem with you only having the data shown in the plot. $\endgroup$ – Reinstate Monica - G. Simpson Oct 2 '19 at 15:50
  • $\begingroup$ Thank you for the response! Exactly my thoughts, except, my colleagues are trained in the traditional DOE and ANOVA which harps on replications, replications, replications. $\endgroup$ – Dilsher Singh Dhillon Oct 2 '19 at 16:02
  • $\begingroup$ Your follow-up is probably best as a new question. But, you can't prove it as you don't know the truth. When people normally do this they start with some truth, sample data from said truth, then form intervals, then look at the coverage of those intervals. Repeat many times and you have some idea as to the behaviour of the intervals. You can simplify what you wrote to use gratia::predicted_samples() or gratia::simulate() for GAMs, which will generate simulated data (responses) from the model & then do quantile() on the simulated values for each observation in newdata to get intervals. $\endgroup$ – Reinstate Monica - G. Simpson Oct 2 '19 at 17:35

Your Answer

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

Browse other questions tagged or ask your own question.