data <- matrix(scan(),ncol=4,byrow=TRUE)
 1 1 1 41.2  1 1 1 42.6  1 1 2 41.2 1 1 2 41.4 
 2 1 1 135.7  2 1 1 136.8  2 1 2 143.0  2 1 2 143.3 
 3 1 1 163.2  3 1 1 163.3  3 1 2 181.4  3 1 2 180.3 
 1 2 3 39.8  1 2 3 40.3    1 2 4 41.5  1 2 4 43.0 
 2 2 3 132.4  2 2 3 130.3  2 2 4 134.4  2 2 4 130.0 
 3 2 3 173.6  3 2 3 173.9  3 2 4 174.9  3 2 4 175.6 
 1 3 5 41.9  1 3 5 42.7    1 3 6 45.5   1 3 6 44.7 
 2 3 5 137.4  2 3 5 135.2  2 3 6 141.1  2 3 6 139.1 
 3 3 5 166.6  3 3 5 165.5  3 3 6 175.0  3 3 6 172.0 

data

table712  <- data.frame(data)
colnames(table712)  <- c("conc","day","run","glucose")
rm(data)

table712$conc  <- as.factor(table712$conc)
table712$day  <- as.factor(table712$day)
table712$run  <- as.factor(table712$run)


aov(glucose ~ conc +day +conc:day + run%in%day +run:conc%in%day, data = table712)


 # here is an approach using the afex package

library(afex)

table712.afex1 <- mixed(glucose ~ conc +(1|day) +(1|conc:day) +(1|run:day) +(1|run:conc:day), data=table712)

summary(table712.afex1)

nice(table712.afex1)

library(emmeans)

emmeans(table712.afex1, ~conc)

# examining the conc x run(day) effect

boxplot(glucose ~ run*conc, data = table712)