Predicting Shape Along a PC Axis for 3D Data

[Claudia]

Hi, there - I'm trying to predict shape at the minimum and maximum ends of a PC axis.  While I can do this using the 2D dataset "plethodon," I am having trouble making this work using a 3D dataset, such as "scallops."

The following code works "plethodon," but not for "scallops":

## Select PC1
PC <- PCA$x[, 1]

## Predict shape along PC 1
preds <- shape.predictor(Y.gpa$coords, x = PC, Intercept = FALSE,
                         pred1 = min(PC),  
                         pred2 = max(PC))

When I try to do this using the "scallops" data, I receive the following error message:

Error in model.frame.default(formula = form, data = dat, drop.unused.levels = TRUE) :
  variable lengths differ (found for 'X')
In addition: Warning message:
In data.frame(Y = Y, X = X) :
  row names were found from a short variable and have been discarded

Can you please tell me what I should do to make set predictions when using the "scallops" dataset?  I've read through the help file, but don't see any examples that run 3D data in doing so. 

Thank you for your help!

Best,

Claudia

Claudia M. Astorino 

Ph.D. Candidate

Biological Anthropology

NYCEP

CUNY Graduate Center

365 Fifth Avenue
New York, NY 10016

cmastorino.com

pronouns: she/her & they/them

[Mike Collyer]

Dear Claudia,

It worked fine for me using packages RRPP 2.1.2.999 and geomorph 4.0.9.999, plus the most current of their dependencies.  If I had to guess (since you did not say), you might be using older versions of these packages?  It is possible the issue was resolved before and updating to more recent versions of the packages would help.  Or, it could be something that you did to obtain the PCA object.  I used 

PCA <- gm.prcomp(Y.gpa$coords)

but maybe you did something different?

In other words, it is not possible to diagnose the issue you have because vital information is missing.  But these are the two places I would guess something could go wrong.

Best,

Mike

--
You received this message because you are subscribed to the Google Groups "geomorph R package" group.
To unsubscribe from this group and stop receiving emails from it, send an email to geomorph-r-pack...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/geomorph-r-package/958b5ded-fccc-4d96-99a3-76223fedaf3dn%40googlegroups.com.

[Claudia]

Hi - thank you for your quick reply!  For some reason, PCA was not saved as an object in my environment - tried it again after running PCA <- gm.prcomp(Y.gpa$coords) and it works fine.  Thank you!

Claudia

[mura...@gmail.com]

Mike,

I have a followup question about this. What's the difference between obtaining the PC extremas via shape.predict vs the coordinates returned by the gm.prcomp function? They are close, but not identical. The difference becomes much more noticeable if one puts the scale back in

library(geomorph)
data("plethodon")
Y.gpa <- gpagen(plethodon$land)    #GPA-alignment    
PCA <- gm.prcomp(Y.gpa$coords)
PC <- PCA$x[,1]

preds <- shape.predictor(Y.gpa$coords, x = PC, Intercept = FALSE,

                         pred1 = min(PC), pred2 = max(PC)) # PC 1 extremes, more technically
cbind(PCA$shapes$shapes.comp1$min, preds$pred1)

             X           Y           X           Y
1   0.16854339 -0.03004810  0.16924480 -0.03017315
2   0.21590881 -0.10565617  0.21680735 -0.10609588
3  -0.05589353 -0.01478701 -0.05612613 -0.01484855
4  -0.28239636 -0.09043202 -0.28357159 -0.09080836
5  -0.31217922 -0.05213652 -0.31347840 -0.05235350
6  -0.32823434 -0.02531190 -0.32960034 -0.02541724
7  -0.31326224  0.05433565 -0.31456593  0.05456178
8  -0.18070075  0.10199298 -0.18145276  0.10241744
9   0.03327212  0.08519855  0.03341059  0.08555312
10  0.20197495  0.07036741  0.20281550  0.07066026
11  0.27960272  0.06605489  0.28076633  0.06632978
12  0.57336445 -0.05957776  0.57575059 -0.05982570

#scale Y.gpa by centroid sizes
boas.gpa=Y.gpa$coords
for (i in 1:40) boas.gpa [,,i]=Y.gpa$coords[,,i] * Y.gpa$Csize[i]

PCA.boas <- gm.prcomp(boas.gpa)
PC.boas <- PCA.boas$x[,1]
preds.boas <- shape.predictor(boas.gpa, x = PC.boas, Intercept = FALSE,
                         pred1 = min(PC.boas), pred2 = max(PC.boas)) # PC 1 extremes, more technically

cbind(PCA.boas$shapes$shapes.comp1$min, preds.boas$pred1)
             X            Y          X          Y
1   0.14408750 -0.033671721  2.9157167 -0.6813721
2   0.18924267 -0.106759900  3.8294649 -2.1603652
3  -0.03423911 -0.006577069 -0.6928535 -0.1330918
4  -0.28383591 -0.085932779 -5.7436288 -1.7389131
5  -0.31181882 -0.054111868 -6.3098835 -1.0949935
6  -0.32806395 -0.028729714 -6.6386158 -0.5813669
7  -0.31597027  0.040436233 -6.3938912  0.8182570
8  -0.17965804  0.099201617 -3.6355129  2.0074178
9   0.02177795  0.096415451  0.4406929  1.9510377
10  0.19020766  0.070859408  3.8489921  1.4338923
11  0.35478953  0.060581591  7.1794276  1.2259131
12  0.55348079 -0.051711248 11.2000916 -1.0464152

[Mike Collyer]

Hi Murat,

Good question!  Also an interesting one because you explored this with boas coordinates.  Your solution is missing one step that gm.prcomp does (that maybe is not strictly required).  Essentially, shape.predictor is nothing more than prediction from a linear model but it cannot be guaranteed that the predictions have unit size, like Procrustes coordinates do.  So by including one step after you obtain preds, as such:

preds <- lapply(preds, geomorph:::cs.scale)

the configurations are scaled to unit size.  (The cs.scale function is an internal function but essentially is a rescaling of a configuration by dividing by its centroid size.)

Whether using unit-size configurations or allowing the configurations to expand or contract by prediction, the deformation grids should look about the same except a degree of size-magnification.  Having that magnification would probably be important for inferences made using boas coordinates.  

Now, whether gm.prcomp should do this is a debatable topic.  The functional approach predates Bookstein’s advocacy for boas coordinates.  So if one does like you, uses gm.prcomp on boas coordinates, the step of scaling to centroid size produces erroneous results.  However, the gm.prcomp help file explicitly states that the argument should be Procrustes coordinates.  Using shape.predictor rather than the output provided by gm.prcomp will give you non-unit-size coordinates for which you could use PlotRefToTarget to obtain deformation grids.  Maybe this is an encumbrance but it would work.

We should consider adding an additional argument for whether to rescale shapes along PCs, after considering whether this is appropriate.  After all, shape.predictor is not predicting shapes if using boas coordinates; it is predicting form, perhaps.

Best,

Mike

To view this discussion, visit https://groups.google.com/d/msgid/geomorph-r-package/9fa2cdbf-b5d4-4136-9bf2-830adab5c561n%40googlegroups.com.

[mura...@gmail.com]

Thanks Mike,

To recap: predict.shape does not scale the coordinates after the prediction. If one want to enforce unit size constrain on prediction, they should run cs.scale on resultant coordinates, or directly use the min/max from gm.prcomp.

I think this makes total sense. I don't think you need to change anything, other than perhaps add an enter to the help of shape.predictor  the output will be in the scale of the input, and not necessarily has unit dimensions.

Best,

M

[Mike Collyer]

Hi Murat,

To recap: predict.shape does not scale the coordinates after the prediction.

Correct.

If one want to enforce unit size constrain on prediction, they should run cs.scale on resultant coordinates, or directly use the min/max from gm.prcomp.

Correct.  However, the consequence of not rescaling when using Procrustes coordinates is probably not great.  For example,

> preds <- shape.predictor(Y.gpa$coords, x = PC, Intercept = TRUE,

+                          pred1 = min(PC), pred2 = max(PC)) # PC 1 extremes, more technically

> preds2 <- lapply(preds, geomorph:::cs.scale)

> two.b.pls(preds[[1]], preds2[[1]])

Random PLS calculations: 1000 permutations.

  |==========================================================| 100%

Call:

two.b.pls(A1 = preds[[1]], A2 = preds2[[1]]) 

r-PLS: 1

Effect Size (Z): 3.1783

P-value: 0.001

Based on 1000 random permutations

> sapply(preds, geomorph:::csize)

   pred1    pred2 

1.004162 1.001530 

> sapply(preds2, geomorph:::csize)

pred1 pred2 

    1     1 

So by not using cs.scale, the resulting centroid size is nearly unit size and the two approaches — scaling or not scaling — produce perfectly correlated shapes, meaning an affine transformation of one to the other (and I am sure the affine transformation is just enlargement).  Although I have not verified, I believe using all PCs instead of just the first would result in predictions that have unit size.  The distortion caused by not rescaling is probably far less than would be made if someone made a figure in powerpoint and rescaled the deformation grids that are slapped onto a PC plot.  This is to say, this might be much ado about nothing, unless one wishes to predict shapes way beyond the realm of possibilities encompassed by the observed data, which one could do, for example, with picknplot.shape.  Not rescaling configurations would be weird, as some deformation grids would expand and others shrink, but it would just be a photocopy enlargement effect.

With boas coordinates, it’s not quite the same, and the photocopy enlargement effect would be important.

I think this makes total sense. I don't think you need to change anything, other than perhaps add an enter to the help of shape.predictor  the output will be in the scale of the input, and not necessarily has unit dimensions.

I’ll have to think about this.  This is a good idea, as a minimum, but it’s also low fruit to offer one additional argument that empowers users.  As you know, when empowering users, there is a delicate balance between offering additional tools and offering options to mess things up.  We can wrestle with that for a bit.  Regardless, happy that at least this makes sense for you!

Best,

Mike 

To view this discussion, visit https://groups.google.com/d/msgid/geomorph-r-package/d7213c93-1519-4c28-90ab-215875c51a8fn%40googlegroups.com.

[Mike Collyer]

Dear Colleagues,

After considering the discussion Murat and I had, yesterday, we decided to update geomorph to include the argument, unit.size, in both shape.predictor and gm.prcomp.  The default is FALSE in shape predictor but is TRUE in gm.prcomp, meaning nothing should change with using these functions.  However, now one can choose whether to scale configurations without having to resort to a post hoc manipulation of configurations.  This flexibility might be important if one is using Boas coordinates, for example.

These updates can be made by re-installing geomorph from Github.

Best,

Mike