I don't understand the problem that these differences would pose for analyzing the treatments over time. You are interested in comparing the effects of the treatment over time across samples, right? That can be done as a trajectory analysis, comparing both the rates and directions of the change from each sample mean. So, they are different at the outset, but the analysis would compare the changes from the different outsets from time to time. Given that you have multiple time points, that would also allow you to compare the "shape" of the trajectories across samples over time. I don't see how a PCA would help you with the design of this analysis, although it might provide a useful low-dimensional visual summary of your data. The analysis can be done using the trajectory.analysis function in geomorph/RRPP.

If the PCA scores are a good representation of the variation in shape, they could be used as a covariate. However, if initial shape is to have an effect on facial deformity progression (which I assume it will given the question to use a covariate), then the experiment may be confounded because the different treatments also differ in the initial shape. Better to restart the experiment and randomly assign individudals to the different treatments. This will make sure that there is no potential confounding factors in the analysis from the onset of the experiment. My two cents.

I think that it would help to see your statistical model's design. It is not clear from your question why the difference at 1 week would matter (or complicate the analysis) unless you think that there would be an interaction between shape at that stage and the subsequent treatment effects. If those differences at 1 week are already evidence of treatment effects, then why would you want to control for them? Using PC1 as a covariate wouldn't control for the differences in mean unless all of the difference between the means at 1 week is along that one axis, and nothing else is. But once it's clear how you've designed your model, it may be more clear what the issue is and how to deal with it.

Hello,

Restarting the project is not possible. Essentially, there were potential issues in the family groups involved in the study and how they were sorted. We were sorting out 80,000 fish by hand into 4 different treatments with 3 replicates (12 tanks total). We divided 8 source tanks (multiple family groups in each) into these 12 treatment tanks, placing half a source tank in a low-density tank and half into a high-density tank. There was no overlap in source tanks between each tank (ie. source tank combination A, B, C could not be found in any other tank which could be said for all 12 tanks). We understand there was a flaw in our methods, but we are hoping to salvage what we can so I can get my thesis completed. We have been looking into the trajectory.analysis function in geomorph/RRPP (thank you, Miriam). However, if there are any other potential fixes, we are open to ideas.

Thanks for all your feedback by the way! I greatly appreciate it.

You should be able to complete the analysis using RRPP (which can handle mixed models). The pairwise function in RRPP is more flexible than trajectory analysis because you use it after running the model to set up the pairwise comparisons that you want to do. You run the model, see if you have interesting interactions, and design the pairwise analyses to show you whether those are because of different magnitudes of directions of the response. I thought of the trajectory analysis because of the "shape" parameter which compares the shapes of the curves (of the response, over time), but there are many possibilities for designing analyses. You know that there are differences at the outset, but you can compare population means over the course of the treatment to find out if they increase/decrease of if the differences, over time, are in the same direction. So even if your design isn't optimal, there is a lot that you can do.