procrustes ANOVA

569 views
Skip to first unread message

Alex Sotola

unread,
Sep 20, 2016, 3:03:26 PM9/20/16
to geomorph R package
Hello,

I am trying to run a MANOVA using this package on some landmark/semi landmark shape data and classifiers of species and site. However, when I run the code I get an error (code below). Also, when I run the single factor ANOVA, they run just fine and give me results. Any help would be appreciated. 

Thanks!

> procD.lm(shape~species*site,data=gdf,iter=999,RRPP=FALSE,print.progress=FALSE)
Error in pfit$wXs[[k + 1]] : subscript out of bounds

Michael Collyer

unread,
Sep 21, 2016, 5:43:31 AM9/21/16
to Alex Sotola, geomorph R package
Alex,

This error could be caused by a number of issues, which are not easily recognizable without more information.  It looks like you are running the procD.lm example?  Did you set up the geomorph data frame, as in the example?

We certainly do not get an error when running the examples, and it would not be possible to publish the package through CRAN if there were errors.  To help you out we ned to see how you are setting up the objects you use as arguments in your procD.lm step that causes the error.

Cheers,
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 post to this group, send email to geomorph-...@googlegroups.com.
Visit this group at https://groups.google.com/group/geomorph-r-package.
To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/44a0b2b0-51e6-4eb2-b73b-56326557be3c%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Alex Sotola

unread,
Sep 21, 2016, 8:56:31 AM9/21/16
to geomorph R package, vaso...@gmail.com
Thanks for your response. My code is below. I've successfully run PCAs and one way ANOVAs but when I add the interaction of site*species I get the error. 

Thanks in advance!
Alex


d <- readland.tps("Chubs_Initial_16_SL_new.tps",readcurves=TRUE,specID="ID",warnmsg=TRUE)
dsliders <- read.csv("Sliders.csv",header=TRUE)
classifier <- read.csv("Classifiers2.csv",header=TRUE)
d

#define sliders
d1 <- define.sliders(d[,,1],nsliders=14)
d1
dsliders <- data.frame(d1)

#Generalized Procrustes Analysis (GPA)

gpa <- gpagen(d,curves=dsliders,ProcD=TRUE,print.progress=FALSE)
summary(gpa)
plot(gpa)

plotOutliers(gpa$coords)

gp <- as.factor(paste(classifier$Species))
gp2 <- as.factor(paste(classifier$Site))

##PCA
PCA <- plotTangentSpace(gpa$coords,groups=gp,legend=TRUE)
PCA2 <- plotTangentSpace(gpa$coords,groups=gp2,legend=TRUE)

pvar <- (PCA$sdev^2)/(sum(PCA$sdev^2))
names(pvar) <- seq(1:length(pvar))
barplot(pvar,xlab="Principal Components",ylab="% Variance",ylim=c(0,0.4))

pvar2 <- (PCA2$sdev^2)/(sum(PCA2$sdev^2))
names(pvar2) <- seq(1:length(pvar2))
barplot(pvar2,xlab="Principal Components",ylab="% Variance",ylim=c(0,0.4))

#ANOVA
gdf <- geomorph.data.frame(shape=gpa$coords,
                           site=classifier$Site,
                           species=classifier$Species)
gdf
attributes(gdf)

procD.lm(shape~site,data=gdf,iter=999,RRPP=FALSE,print.progress=FALSE)
procD.lm(shape~species,data=gdf,iter=999,RRPP=FALSE,print.progress=FALSE)

procD.lm(shape~species*site,data=gdf,iter=999,RRPP=FALSE,print.progress=FALSE) ##This is where the error occures


On Wednesday, September 21, 2016 at 4:43:31 AM UTC-5, Michael Collyer wrote:
Alex,

This error could be caused by a number of issues, which are not easily recognizable without more information.  It looks like you are running the procD.lm example?  Did you set up the geomorph data frame, as in the example?

We certainly do not get an error when running the examples, and it would not be possible to publish the package through CRAN if there were errors.  To help you out we ned to see how you are setting up the objects you use as arguments in your procD.lm step that causes the error.

Cheers,
Mike

On Sep 20, 2016, at 3:03 PM, Alex Sotola <vaso...@gmail.com> wrote:

Hello,

I am trying to run a MANOVA using this package on some landmark/semi landmark shape data and classifiers of species and site. However, when I run the code I get an error (code below). Also, when I run the single factor ANOVA, they run just fine and give me results. Any help would be appreciated. 

Thanks!

> procD.lm(shape~species*site,data=gdf,iter=999,RRPP=FALSE,print.progress=FALSE)
Error in pfit$wXs[[k + 1]] : subscript out of bounds


--
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-package+unsub...@googlegroups.com.

Michael Collyer

unread,
Sep 21, 2016, 9:10:43 AM9/21/16
to Alex Sotola, geomorph R package
It looks like the problem could be that you failed to indicate species and sites are factors, as you did when you called them gp and gp2.  This is just a guess, without results.  But once you get to the part where you have a factor interaction, it might be that the underlying model design matrix is not properly defined and the function searches for but does not find the appropriate matrix.  The error statement, "pfit$wXs[[k + 1]] “ is a result of the function trying to use the last design matrix for the various procD.lm fits (there would actually be four: intercept, species, species + site, species + site + species:site).  The k + 1 tells me that it is the last model as k = number of model terms (species, site, species:site).  Something is wrong with this matrix.  What is wrong, I do not know.  But if I had to guess, it’s that the variables you entered did not make sense, and that it did not happen with the single-factor MANOVAs might be because you actually did multivariate regression, and only once you had an interaction, did your model matrix blow up.

Try making these factors and let me know if that does not fix the problem.

Cheers!
Mike


To unsubscribe from this group and stop receiving emails from it, send an email to geomorph-r-pack...@googlegroups.com.

To post to this group, send email to geomorph-...@googlegroups.com.
Visit this group at https://groups.google.com/group/geomorph-r-package.

Alex Sotola

unread,
Sep 21, 2016, 9:36:05 AM9/21/16
to geomorph R package, vaso...@gmail.com
Thank you again for the response. Unfortunately, that had occurred to me, and I did check to make sure that site and species were both factors, which they were, but I neglected to put that in the code in the previous reply, sorry for that. Any other ideas? 
Mike


Michael Collyer

unread,
Sep 21, 2016, 11:14:38 AM9/21/16
to Alex Sotola, geomorph R package
Alex,

Thank you for sending your data and script.  Although your issue is specific to your data, this is something everyone can learn from.

The reason you received an error - one which is hard to detect and alert - is that you do not have any factor interactions.  The levels of your site variable are unique to levels in your species variable.  The most complex model you can have is one that includes species and site, but not these plus their interaction, as there are no coefficients for interactions to estimate.  What R and geomorph does in this situation is just pump out possible model design matrices, but one is missing because you asked for an interaction model but it could not be generated.  The “out of bounds” error is a common one, when R is asked to do something say, k + 1 times, when there are only k objects for such operations.

That was technical; here is something more practical.  You do not want to have a factorial model; you want a nested factor model.

This does not work
procD.lm(shape~species*site,data=gdf,iter=999,RRPP=FALSE,print.progress=FALSE)

But this does
procD.lm(shape~species/site,data=gdf,iter=999,RRPP=FALSE,print.progress=FALSE)

I also recommend using nested update, e.g,

myANOVA <- procD.lm(shape~species/site,data=gdf,iter=999,RRPP=FALSE,print.progress=FALSE)
nested.update(myANOVA, ~species/site)

This will adjust F values for fixed effects that have nested terms.

Hope that helps!

Mike


To unsubscribe from this group and stop receiving emails from it, send an email to geomorph-r-pack...@googlegroups.com.

To post to this group, send email to geomorph-...@googlegroups.com.
Visit this group at https://groups.google.com/group/geomorph-r-package.
Reply all
Reply to author
Forward
0 new messages