Tests in the new transducer and rubber caps

---

1 - Database

---

data=read.csv("C:/Users/Samsung/OneDrive/Documentos/1 - PósDoc/2 - FAPESP - Silvipastoril/Responsabilidades/Teste do Transducer novo/data.csv")
data

#factors
data$Rubber_cap=as.factor(data$Rubber_cap)
data$Transducer=as.factor(data$Transducer)
str(data)

'data.frame':   140 obs. of  7 variables:
 $ ID              : chr  "PadraoLANA" "PadraoLANA" "PadraoLANA" "PadraoLANA" ...
 $ Bottle          : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Rubber_cap      : Factor w/ 2 levels "NR","OR": 1 1 1 1 1 2 2 2 2 2 ...
 $ Time            : num  2 2 2 2 2 2 2 2 2 2 ...
 $ Transducer      : Factor w/ 2 levels "New","Old": 2 2 2 2 2 2 2 2 2 2 ...
 $ psi             : num  NA NA NA NA NA NA NA NA NA NA ...
 $ Actual_volume_mL: num  NA NA NA NA NA NA NA NA NA NA ...

---

2 - Old vs New transducer and rubber caps - PSI

---

data1 = subset(data, !Time == 2)
data1

#model
mod = lmer(psi~Transducer*Rubber_cap + (1|Time), data = data1)
#P value e Médias
Pvalue = anova(mod)
Pvalue

Type III Analysis of Variance Table with Satterthwaite's method
                       Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)    
Transducer            0.56709 0.56709     1   108 18.0396 4.6e-05 ***
Rubber_cap            0.00580 0.00580     1   108  0.1846  0.6683    
Transducer:Rubber_cap 0.01695 0.01695     1   108  0.5391  0.4644    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias1 = emmeans(mod, ~ Transducer)

NOTE: Results may be misleading due to involvement in interactions

medias1

 Transducer emmean    SE   df lower.CL upper.CL
 New         0.938 0.324 5.03    0.106     1.77
 Old         1.077 0.324 5.02    0.245     1.91

Results are averaged over the levels of: Rubber_cap 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

medias2 = emmeans(mod, ~ Rubber_cap)

NOTE: Results may be misleading due to involvement in interactions

medias2

 Rubber_cap emmean    SE   df lower.CL upper.CL
 NR           1.00 0.324 5.03    0.168     1.83
 OR           1.01 0.324 5.02    0.182     1.85

Results are averaged over the levels of: Transducer 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

#Old vs New transducer - PSI

library(ggplot2)
plot1=ggplot(data1, aes(x=Time, y=psi, shape=Transducer))+
  geom_line(stat="summary",fun="mean")+
  geom_point(aes(shape = Transducer), size = 3, stat="summary",fun="mean")+
  scale_shape_manual(values = c(1, 15), name= "Transducer",labels = c("New", "Old"))+
  scale_y_continuous(name="Pressure (psi)", breaks = seq(0,3,0.1))+
  scale_x_continuous(name="Fermentation time (h)", breaks=seq(0,24,2))
plot1

plot1.1=plot1+theme(axis.line = element_line(colour = "black", size = 0.7, linetype = "solid"),
 panel.background = element_rect(fill = "transparent"),
 legend.background = element_rect(fill = "transparent", size=0.5, linetype="solid",colour ="black"),legend.position = c(0.1, 0.75),legend.key.size = unit(0.42, 'cm'))+
  annotate(geom="text", y=2.6, x=4,label="P value = <0.001",size=4,color="black")
plot1.1
ggsave("Transducer.png", plot = plot1.1, width = 16, height = 11, units = "cm", dpi = 300)

#Old vs New rubber caps - PSI

plot2=ggplot(data1, aes(x=Time, y=psi, shape=Rubber_cap))+
  geom_line(stat="summary",fun="mean")+
  geom_point(aes(shape = Rubber_cap), size = 3, stat="summary",fun="mean")+
  scale_shape_manual(values = c(1, 15), name= "Rubber cap",labels = c("New", "Old"))+
  scale_y_continuous(name="Pressure (psi)", breaks = seq(0,3,0.1))+
  scale_x_continuous(name="Fermentation time (h)", breaks=seq(0,24,2))
plot2

plot2.1=plot2+theme(axis.line = element_line(colour = "black", size = 0.7, linetype = "solid"),
 panel.background = element_rect(fill = "transparent"),
 legend.background = element_rect(fill = "transparent", size=0.5, linetype="solid",colour ="black"),legend.position = c(0.1, 0.75),legend.key.size = unit(0.42, 'cm'))+
  annotate(geom="text", y=2.6, x=4,label="P value = 0.668",size=4,color="black")
plot2.1
ggsave("Rubber caps.png", plot = plot2.1, width = 16, height = 11, units = "cm", dpi = 300)

---

3 - Calibration curve for the new transducer

---

data2 = subset(data, !Transducer == "Old") 
data2

mod1 = lm(Actual_volume_mL ~ psi, data = data2)
summary(mod1)

Call:
lm(formula = Actual_volume_mL ~ psi, data = data2)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4377 -0.3393  0.0277  0.4118  3.6288 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.6399     0.2838   2.255   0.0275 *  
psi           5.3324     0.1937  27.534   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.401 on 65 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.921, Adjusted R-squared:  0.9198 
F-statistic: 758.1 on 1 and 65 DF,  p-value: < 2.2e-16

colSums(!is.na(data2[, c("psi", "Actual_volume_mL")]))

             psi Actual_volume_mL 
              67               67 

plot3 = ggplot(data2, aes(x = psi, y = Actual_volume_mL))+
  geom_point(shape=19, color='darkblue', size=1.3)+
  geom_smooth(method = "lm", se = TRUE, color='black')+
  scale_x_continuous(name ="Pressure (psi)", breaks=seq(0,4,0.5))+
  scale_y_continuous(name = "Actual volume (mL)",breaks=seq(0,20,2))
plot3 

  
plot3.1=plot3+ 
  theme(panel.background = element_rect(fill = "transparent", color = "black",size = 0.3),
  axis.line = element_line(size = 0.3, linetype = "solid",color = "black"),
  axis.ticks = element_line(size = 0.3,color = "black"),
  axis.title.y = element_text(size = 12,color = "black"),
  axis.title.x = element_text(size = 12,color = "black"),
  axis.text.x = element_text(size = 10,color = "black"),
  axis.text.y = element_text(size = 10,color = "black"))+
annotate(geom="text", y=17, x=.7,label=expression(paste(Transducer, " = Newlab 2025")),size=5,color="red4")+
annotate(geom="text", y=15, x=.7,label=expression(paste(P-value, " = 0.0275")),size=5,color="black")+
annotate(geom="text", y=13, x=.7,label=expression(paste(y, " = 5.33x + 0.6399")),size=5,color="black")+
annotate(geom="text", y=11, x=.7,label=expression(paste(R^2, " = 0.92")),size=5,color="black")+
annotate(geom="text", y=9, x=.7,label=expression(paste(n, " = 67")),size=5,color="black")
plot3.1

ggsave("curve.png", plot = plot3.1, width = 16, height = 11, units = "cm", dpi = 300)