For simplicity, first consider only one group. Your data can be modeled as $$Y_{ij} = \mu + \delta_i + \epsilon_{ij},$$ where $\delta_i \sim N(0, \tau^2), i=1,\ldots,k$ is the random day effect, and $\epsilon_{ij}\sim N(0, \sigma^2), j=1,\ldots, n_i$ is the random within-day replicate effect. Then you can obtain an unbiased estimate of the group average $$\hat\mu = \frac 1k \sum_{j=1}^j \bar{Y}_{i.} \sim N(\mu, \frac{\tau^2}{k} + \frac 1k \sum_{i=1}^k \frac{\sigma^2}{n_i}).$$ It is fairly simple to obtain estimates of the variance term with $\tau^2$ being the variance of the day means, and $\sigma^2$ the pooled within-day variance.
Start edit
The variance term in this formula gives you the square of the standard error for your plot. The point to show would be the mean of the within-day means. If you don't want to assume that the within-day variability is constant, just replace $\sigma$ with $\sigma_i$ - the within-day standard deviation on day 1, so the summation would be $\sum_iSE_i^2$, where $SE_i$ is the standard error of the within-day mean.
End edit
If for a moment you assume $n_i=n$, then the variance is $\frac{\tau^2}{k} + \frac{\sigma^2}{n}$. In your application the key observation is that $n >> k$, so the first term will dominate. In essence, you can ignore the within-day replicates, and just use them as an expensive way to obtain one (day-specific) observation.
I have run a "correct" mixed ANOVA analysis and a simple ANOVA on group means, and they give essentially the same result:
library(lme4)
m1 <- lmer(Meas ~ Treatment*Temp + (1|Rep), data=dat)
m1
Output:
Linear mixed model fit by REML
Formula: Meas ~ Treatment * Temp + (1 | Rep)
Data: dat
AIC BIC logLik deviance REMLdev
3351 3376 -1669 3343 3339
Random effects:
Groups Name Variance Std.Dev.
Rep (Intercept) 15.414 3.9260
Residual 42.000 6.4808
Number of obs: 504, groups: Rep, 12
Fixed effects:
Estimate Std. Error t value
(Intercept) 16.1579 5.2369 3.085
TreatmentT2 -9.3838 7.3903 -1.270
Temp -0.4723 0.3315 -1.425
TreatmentT2:Temp 1.0377 0.4683 2.216
And
m2 <- lm(mean ~ Treatment*Temp, data=dat2)
summary(m2)
with output
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.1037 5.2451 3.070 0.0153 *
TreatmentT2 -9.2846 7.4177 -1.252 0.2460
Temp -0.4694 0.3317 -1.415 0.1948
TreatmentT2:Temp 1.0302 0.4691 2.196 0.0594 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.063 on 8 degrees of freedom
Multiple R-squared: 0.5955, Adjusted R-squared: 0.4438
F-statistic: 3.926 on 3 and 8 DF, p-value: 0.05412