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.
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.
library(mgcv) library(janitor) library(tidyverse)
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())
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)
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) library(compiler) 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)
How does one prove the coverage of the fitted intervals? Simply looking at the data they seem overly optimistic.