editet statistics to dunnett test to have only time effect control group not all verses all

code/plotting_quantitative_dti_values.py

@@ -7,6 +7,7 @@ from statsmodels.stats.multicomp import MultiComparison
 from itertools import combinations
 import statsmodels
 import scipy
+from scipy.stats import dunnett
 
 # Set up warning suppression for convergence
 import warnings
@@ -58,7 +59,7 @@ code_dir = os.path.dirname(os.path.abspath(__file__))
 parent_dir = os.path.dirname(code_dir)
 
 # Create a new folder for mixed model analysis
-mixed_model_dir = os.path.join(parent_dir, 'output', "Final_Quantitative_output", 'mixed_model_analysis')
+mixed_model_dir = os.path.join(parent_dir, 'output', "Final_Quantitative_output", 'mixed_model_analysis_dunnett')
 os.makedirs(mixed_model_dir, exist_ok=True)
 
 # Step 4: Save the resulting dataframe to a CSV file
@@ -115,21 +116,26 @@ for dd in df["dialation_amount"].unique():
                 # Plot and save boxplots for significant models
                 PlotBoxplot(filtered_df, qq, mm,dd ,mixed_model_dir,sig)
                 
-                # Perform Tukey's HSD test for post hoc analysis
+                # Perform Dunnet test for post hoc analysis
+                log.write("Dunnett Multiple Comparisons for time effect:\n")
+                log.write("Group\tTP1\tTP2\tPValue\tSignificant\n")
                 for gg in filtered_df["Group"].unique():
-
-                    filtered_df2 = filtered_df[filtered_df["Group"]==gg]
-                    mc = MultiComparison(filtered_df2["Value"], filtered_df2["merged_timepoint"])
-                    tukey_result = mc.tukeyhsd()
-                    
-                    # Print the summary of the test to log file
-                    log.write(f"Posthoc Analysis for {gg}:\n")
-                    log.write(str(tukey_result.summary()) + '\n')
-                
+                    control_data = filtered_df[(filtered_df["Group"]==gg) & (filtered_df["merged_timepoint"]==0)]["Value"]
+                    for tt in filtered_df["merged_timepoint"].unique():
+                        treatment_data = filtered_df[(filtered_df["Group"]==gg) & (filtered_df["merged_timepoint"]==tt)]["Value"]
+                        result = dunnett(treatment_data, control=control_data)
+                        if result.pvalue < 0.05:  # You can adjust the significance level as needed
+                            significant = "Yes"
+                        else:
+                            significant = "No"
+                        
+                        # Write the results to the log file
+                        log.write(f"{gg}\t0\t{tt}\t{result.pvalue}\t{significant}\n")
+                                    
                 # Perform the Sidak multiple comparisons for the group effect
                 group_combinations = combinations(filtered_df["Group"].unique(), 2)
-                log.write("Sidak Multiple Comparisons:\n")
-                log.write("Group1\tGroup2\tTimePoint\tPValue\tSignificantDifference\n")
+                log.write("Sidak Multiple Comparisons for group effect:\n")
+                log.write("Group1\tGroup2\tTimePoint\tPValue\tSignificant\n")
                 for group1, group2 in group_combinations:
                     for time_point in filtered_df["merged_timepoint"].unique():
                         group1_data = filtered_df[(filtered_df["Group"] == group1) & (filtered_df["merged_timepoint"] == time_point)]["Value"]