PPFPA Data analysis
Faith Ndanu
2023-07-06
Goal: We are looking to see if there are any main effects of each individual variable as well as interaction effects between variables
1. load packages
2. LOAD DATA
PPFPA_Conjugation<- read_csv("PPC_full_data.csv")check if columns: protein, polymer, percent func and ratio are factors and percent conjugation is numbers
str(PPFPA_Conjugation)## spc_tbl_ [144 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Protein : chr [1:144] "Fn3" "Fn3" "Fn3" "Fn3" ...
## $ Polymer : chr [1:144] "mTEG" "mTEG" "mTEG" "mTEG" ...
## $ Percent_func : num [1:144] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ...
## $ Ratio(Pro:poly) : num [1:144] 0 0 0 1 1 1 2 2 2 0 ...
## $ percent_conjugation: num [1:144] 0.2157 0.0876 0.1753 0.7591 0.7658 ...
## - attr(*, "spec")=
## .. cols(
## .. Protein = col_character(),
## .. Polymer = col_character(),
## .. Percent_func = col_double(),
## .. `Ratio(Pro:poly)` = col_double(),
## .. percent_conjugation = col_double()
## .. )
## - attr(*, "problems")=<externalptr>3. DATA WRANGLING
Rename the “Ratio(Pro:poly)” column. Convert the protein, polymer, percent func and ratio to factors in new columns just to ensure that you have the original data
{
PPFPA_Conjugation_factor_combined <- PPFPA_Conjugation %>%
rename(Ratio = `Ratio(Pro:poly)`)%>%
mutate(Protein_f = as.factor(Protein),
Sidechain_f = as.factor(Polymer),
Percent_func_f = as.factor(Percent_func),
Ratio_f = as.factor(Ratio)) %>%
mutate (Protein_sidechain_f = as.factor(paste (Protein_f,
Sidechain_f, sep = "_"))) %>%
mutate (Protein_sidechain_func_f = as.factor(paste (Protein_sidechain_f, Percent_func_f,
sep = "_") )) %>%
mutate ( Protein_sidechain_func_ratio = as.factor(paste(Protein_sidechain_func_f, Ratio_f,
sep = "_"))) %>%
select(-Polymer,-Protein, -Ratio,-Percent_func, Protein_sidechain_func_ratio)
}confirm that the columns have been successfully converted to factors.
str(PPFPA_Conjugation_factor_combined)## tibble [144 × 8] (S3: tbl_df/tbl/data.frame)
## $ percent_conjugation : num [1:144] 0.2157 0.0876 0.1753 0.7591 0.7658 ...
## $ Protein_f : Factor w/ 2 levels "Fn3","Tf": 1 1 1 1 1 1 1 1 1 1 ...
## $ Sidechain_f : Factor w/ 3 levels "GLU","HPA","mTEG": 3 3 3 3 3 3 3 3 3 2 ...
## $ Percent_func_f : Factor w/ 4 levels "0.2","0.4","0.6",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Ratio_f : Factor w/ 3 levels "0","1","2": 1 1 1 2 2 2 3 3 3 1 ...
## $ Protein_sidechain_f : Factor w/ 6 levels "Fn3_GLU","Fn3_HPA",..: 3 3 3 3 3 3 3 3 3 2 ...
## $ Protein_sidechain_func_f : Factor w/ 16 levels "Fn3_GLU_0.2",..: 8 8 8 8 8 8 8 8 8 6 ...
## $ Protein_sidechain_func_ratio: Factor w/ 48 levels "Fn3_GLU_0.2_0",..: 22 22 22 23 23 23 24 24 24 16 ...Before doing the factorial anova we need to formulate some hypotheses for us to test. We can begin by first plotting means fo different variables.We have four different independent(2 x 3 x 2 x 3 comparison groups = 36 total) variables: -Protein: Fn3; Tf (2 categories) -> protein_f -Polymer functional group(side chains): mTEG; HPA; Glucamine(This was is not included in the dataframe yet) (3 categories) ->polymer_f -Percent functionalized: 0.8; 0.6 (2 categories) -Ratio of protein to polymer: 1:1; 1:10; 1: high (3 categories)
And one dependent variable: percent_conjugation( numeral)
comparing just the two proteins
boxplot( percent_conjugation ~ Protein_f,
data = PPFPA_Conjugation_factor_combined,
col = wes_palette(n =2, name = "GrandBudapest1"))
Protein and the side chains
boxplot( percent_conjugation ~ Protein_f * Sidechain_f, data = PPFPA_Conjugation_factor_combined,
col = c("white", "Grey"))
Protein, side chains and percent functionalization
boxplot( percent_conjugation ~ Protein_f * Sidechain_f * Percent_func_f,
data = PPFPA_Conjugation_factor_combined,
col = c("white", "Grey"))
Protein and the ratio
boxplot(percent_conjugation ~ Protein_f *Ratio_f,
data = PPFPA_Conjugation_factor_combined,
n.label = FALSE)
Clearly, the means for these groups are different. Now we need to do an anova to see if these differences are statistically significant. so we run factorial anova
FACTORIAL ANOVA
Model 1: Comparing the protein means( Fn3 Vs Tf)
model_1 <- aov(percent_conjugation ~ Protein_f , PPFPA_Conjugation_factor_combined)
summary( model_1)## Df Sum Sq Mean Sq F value Pr(>F)
## Protein_f 1 0.638 0.6383 4.728 0.0313 *
## Residuals 142 19.173 0.1350
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model 2: Comparing the Protein and sidechains percent conjugation while allowing interaction between the two variables
model_2 <- aov(percent_conjugation ~ Protein_f * Sidechain_f, PPFPA_Conjugation_factor_combined)
summary( model_2)## Df Sum Sq Mean Sq F value Pr(>F)
## Protein_f 1 0.638 0.6383 5.321 0.0226 *
## Sidechain_f 2 2.593 1.2967 10.809 4.36e-05 ***
## Protein_f:Sidechain_f 2 0.023 0.0115 0.096 0.9088
## Residuals 138 16.556 0.1200
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model 3: Comparng the Protein, sidechain, percent functionalization percent conjugation while allowing interaction between the three variables
model_3 <- aov(percent_conjugation ~ Protein_f* Sidechain_f * Percent_func_f , PPFPA_Conjugation_factor_combined)
summary( model_3)## Df Sum Sq Mean Sq F value Pr(>F)
## Protein_f 1 0.638 0.6383 6.242 0.01374 *
## Sidechain_f 2 2.593 1.2967 12.680 9.46e-06 ***
## Percent_func_f 3 0.787 0.2622 2.564 0.05761 .
## Protein_f:Sidechain_f 2 0.023 0.0115 0.112 0.89391
## Protein_f:Percent_func_f 3 0.914 0.3046 2.979 0.03399 *
## Sidechain_f:Percent_func_f 2 1.761 0.8806 8.610 0.00031 ***
## Protein_f:Sidechain_f:Percent_func_f 2 0.004 0.0022 0.022 0.97847
## Residuals 128 13.090 0.1023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For Model1,2 and 3, the p values are so big which means that the dfference in means is not that significant and we might need to in cooporate the protein polymer ratios to see if there is any statistical significance.
Model 4: Comparng the Protein, sidechain, and protein- polymer ration percent conjugation while allowing interaction between the three variables
model_4 <- aov(percent_conjugation ~ Protein_f* Sidechain_f * Ratio_f , PPFPA_Conjugation_factor_combined)
summary( model_4)## Df Sum Sq Mean Sq F value Pr(>F)
## Protein_f 1 0.638 0.638 12.428 0.00059 ***
## Sidechain_f 2 2.593 1.297 25.247 6.02e-10 ***
## Ratio_f 2 8.275 4.137 80.550 < 2e-16 ***
## Protein_f:Sidechain_f 2 0.023 0.011 0.224 0.80003
## Protein_f:Ratio_f 2 0.240 0.120 2.341 0.10044
## Sidechain_f:Ratio_f 4 1.391 0.348 6.769 5.73e-05 ***
## Protein_f:Sidechain_f:Ratio_f 4 0.179 0.045 0.870 0.48408
## Residuals 126 6.472 0.051
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model 5: Comparng the Protein, sidechain,percent_functionalization, and protein- polymer ratio percent conjugation while allowing interaction between the four variables
model_5 <- aov(percent_conjugation ~ Protein_f * Sidechain_f * Ratio_f * Percent_func_f , PPFPA_Conjugation_factor_combined)
summary( model_5)## Df Sum Sq Mean Sq F value Pr(>F)
## Protein_f 1 0.638 0.638 134.220 < 2e-16
## Sidechain_f 2 2.593 1.297 272.662 < 2e-16
## Ratio_f 2 8.275 4.137 869.934 < 2e-16
## Percent_func_f 3 0.787 0.262 55.130 < 2e-16
## Protein_f:Sidechain_f 2 0.023 0.011 2.414 0.0949
## Protein_f:Ratio_f 2 0.240 0.120 25.279 1.52e-09
## Sidechain_f:Ratio_f 4 1.391 0.348 73.108 < 2e-16
## Protein_f:Percent_func_f 3 0.914 0.305 64.056 < 2e-16
## Sidechain_f:Percent_func_f 2 1.761 0.881 185.153 < 2e-16
## Ratio_f:Percent_func_f 6 0.514 0.086 18.012 6.78e-14
## Protein_f:Sidechain_f:Ratio_f 4 0.179 0.045 9.395 1.88e-06
## Protein_f:Sidechain_f:Percent_func_f 2 0.004 0.002 0.468 0.6276
## Protein_f:Ratio_f:Percent_func_f 6 0.351 0.058 12.299 3.16e-10
## Sidechain_f:Ratio_f:Percent_func_f 4 1.536 0.384 80.767 < 2e-16
## Protein_f:Sidechain_f:Ratio_f:Percent_func_f 4 0.148 0.037 7.764 1.83e-05
## Residuals 96 0.457 0.005
##
## Protein_f ***
## Sidechain_f ***
## Ratio_f ***
## Percent_func_f ***
## Protein_f:Sidechain_f .
## Protein_f:Ratio_f ***
## Sidechain_f:Ratio_f ***
## Protein_f:Percent_func_f ***
## Sidechain_f:Percent_func_f ***
## Ratio_f:Percent_func_f ***
## Protein_f:Sidechain_f:Ratio_f ***
## Protein_f:Sidechain_f:Percent_func_f
## Protein_f:Ratio_f:Percent_func_f ***
## Sidechain_f:Ratio_f:Percent_func_f ***
## Protein_f:Sidechain_f:Ratio_f:Percent_func_f ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For both models 4 and 5, the p values for the individual variables are significantly small, which shows that the difference in means upon adding the ration variable is statistically significant. However, the p values for the interaction terms are very small, which could indicate that, the difference in means while allowing interaction between variables is not that significant( I’m not so sure if this is the correct interpretation of the p values for the interaction terms—-> need some help)
Let’s try adding the variables without the interactions for models 4 and 5
For model 4
model_6 <- aov(percent_conjugation ~ Protein_f + Sidechain_f + Ratio_f , PPFPA_Conjugation_factor_combined)
summary( model_6)## Df Sum Sq Mean Sq F value Pr(>F)
## Protein_f 1 0.638 0.638 10.61 0.00142 **
## Sidechain_f 2 2.593 1.297 21.55 7.18e-09 ***
## Ratio_f 2 8.275 4.137 68.75 < 2e-16 ***
## Residuals 138 8.305 0.060
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_7 <- aov(percent_conjugation ~ Protein_f + Sidechain_f + Ratio_f + Percent_func_f , PPFPA_Conjugation_factor_combined)
summary( model_7)## Df Sum Sq Mean Sq F value Pr(>F)
## Protein_f 1 0.638 0.638 11.462 0.00093 ***
## Sidechain_f 2 2.593 1.297 23.285 2.05e-09 ***
## Ratio_f 2 8.275 4.137 74.292 < 2e-16 ***
## Percent_func_f 3 0.787 0.262 4.708 0.00370 **
## Residuals 135 7.518 0.056
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In both models, the p values are very small. Therefore the difference in means is statistically significant.
RUN POST HOC TESTS
JUST EXPLORING THE DATA
PPFPA_Conjugation_Fn3_0.8 <- PPFPA_Conjugation_factor_combined %>%
group_by(Protein_sidechain_func_f, Ratio_f, Protein_f) %>%
filter(Protein_f =="Fn3") %>%
summarize( mean_percent_conjugation = mean(percent_conjugation))
ggplot(data = PPFPA_Conjugation_Fn3_0.8, mapping = aes(x = Protein_sidechain_func_f, y = mean_percent_conjugation, fill = Ratio_f)) +
geom_col(position = "dodge") +
scale_fill_manual(values = wes_palette(n =3, name = "GrandBudapest1")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
# facet_wrap(~Protein_f) +
# theme( panel.spacing = unit(5,"cm",data = NULL))PPFPA_Conjugation_Tf_0.8 <- PPFPA_Conjugation_factor_combined %>%
group_by(Protein_sidechain_func_f, Ratio_f) %>%
filter(Protein_f =="Tf") %>%
summarize( mean_percent_conjugation = mean(percent_conjugation))
ggplot(data = PPFPA_Conjugation_Tf_0.8, mapping = aes(x = Protein_sidechain_func_f, y = mean_percent_conjugation, fill = Ratio_f)) +
geom_col(position = "dodge") +
scale_fill_manual(values = wes_palette(n =3, name = "GrandBudapest1")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
Model_9 <- aov( mean_percent_conjugation ~ Protein_sidechain_func_f + Ratio_f,
data = PPFPA_Conjugation_Tf_0.8)
summary(Model_9 )## Df Sum Sq Mean Sq F value Pr(>F)
## Protein_sidechain_func_f 7 0.9464 0.1352 2.413 0.07618 .
## Ratio_f 2 0.9709 0.4854 8.663 0.00356 **
## Residuals 14 0.7845 0.0560
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TURKEY HOC
Turkey for the no interaction( parallel model)
TukeyHSD( model_7)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = percent_conjugation ~ Protein_f + Sidechain_f + Ratio_f + Percent_func_f, data = PPFPA_Conjugation_factor_combined)
##
## $Protein_f
## diff lwr upr p adj
## Tf-Fn3 -0.1331596 -0.2109442 -0.05537493 0.0009304
##
## $Sidechain_f
## diff lwr upr p adj
## HPA-GLU 0.1896453 0.075488715 0.3038018 0.0003850
## mTEG-GLU 0.3163669 0.202210315 0.4305234 0.0000000
## mTEG-HPA 0.1267216 -0.005095042 0.2585382 0.0623364
##
## $Ratio_f
## diff lwr upr p adj
## 1-0 0.3740308 0.25987424 0.4881874 0.0000000
## 2-0 0.5790047 0.46484812 0.6931612 0.0000000
## 2-1 0.2049739 0.09081732 0.3191304 0.0001145
##
## $Percent_func_f
## diff lwr upr p adj
## 0.4-0.2 -0.2750834 -0.479708535 -0.070458352 0.0035269
## 0.6-0.2 -0.1008254 -0.267901126 0.066250249 0.3993034
## 0.8-0.2 -0.1629620 -0.330037724 0.004113651 0.0587601
## 0.6-0.4 0.1742580 0.007182317 0.341333692 0.0373609
## 0.8-0.4 0.1121214 -0.054954281 0.279197094 0.3043373
## 0.8-0.6 -0.0621366 -0.180276950 0.056003754 0.5214624From the factorial anova test that we did earlier, it showed that there were differences in groups. The above Turkey hoc test shows the differences indicated by the ” diff” column for each group. The p-values are quite small too. For instance, Fn 3 is much better than Tf, a high protein polymer ratio is high and so on.
Next we need to check and see if there are any interactions at all and if so are they significant. The factorial anova that we did on model 5( shows interaction for all the groups) showed that there were differences between interatcing groups( as revealed by smaller p value). Now let’s run the Turkey Hoc test to see where the difference is.
Turkey Hoc for the interaction model
turkey.mod5<-TukeyHSD(model_5)
#turkey.mod5Only significant values
sig.turkey <-turkey.mod5 %>%
tidy() %>% #library(broom)
as.data.frame() %>%
filter(adj.p.value < 0.05)break down into individuals
1:10 80%
visualize to see the 3d effect