procD.lm and advanced.procD.lm on geomorph dataframe returns: Error in cmd$points[, p] : subscript out of bounds

[CarlyMuletzWolz]

I am using procD.lm to determine if distances matrices differ among species, and then advanced.procD.lm.

These are community analyses using jaccard distances and bray-curtis distances. Jaccard distance is fine and same dimension as Bray distances. But, bray distances is returning the Error in cmd$points[, p] : subscript out of bounds.

These are all complex objects to share, so thought I would just share the code and output right now. I can spend more time getting you the data so that you could reproduce the error.

jacc_sp <- distance(GM_data_compare, "jaccard", binary = T)

bray_sp <- distance(GM_data_compare, "bray")

Species <- get_variable(GM_data_compare, "Species2")

Species
[1] Rhesus macaque Rhesus macaque Rhesus macaque Rhesus macaque Human
[6] Mantled howler monkey Rhesus macaque Rhesus macaque Rhesus macaque Rhesus macaque
[11] Mantled howler monkey Rhesus macaque Rhesus macaque Rhesus macaque Rhesus macaque
[16] Rhesus macaque Mantled howler monkey Rhesus macaque Rhesus macaque Rhesus macaque
[21] Rhesus macaque Rhesus macaque Rhesus macaque Rhesus macaque Mantled howler monkey
[26] Rhesus macaque Rhesus macaque Human Human Rhesus macaque
[31] Rhesus macaque Human Human Rhesus macaque Rhesus macaque
[36] Rhesus macaque Western lowland gorilla Owl monkey Rhesus macaque Human
[41] Human Rhesus macaque Sumatran orangutan Rhesus macaque Human
[46] Rhesus macaque Human Owl monkey Rhesus macaque Human
[51] Human Western lowland gorilla Human Human Sumatran orangutan
[56] Human Human Bornean orangutan Human Bornean orangutan
[61] Human Human Human Human Human
[66] Western lowland gorilla Human Western lowland gorilla Human Bornean orangutan
[71] Western lowland gorilla Human Human Bonobo Western lowland gorilla
[76] Chimpanzee Western lowland gorilla Human Human Bornean orangutan
[81] Bornean orangutan Bonobo Mantled howler monkey Human Western lowland gorilla
[86] Western lowland gorilla Chimpanzee Mantled howler monkey
9 Levels: Human Bonobo Chimpanzee Western lowland gorilla Bornean orangutan Sumatran orangutan ... Owl monkey

str(jacc_sp)
'dist' num [1:3828] 0.612 0.521 0.54 0.965 0.936 ...

• attr(*, "Size")= int 88
• attr(*, "Labels")= chr [1:88] "29718" "32912" "30226" "29507" ...
• attr(*, "Diag")= logi FALSE
• attr(*, "Upper")= logi FALSE
• attr(*, "method")= chr "binary jaccard"
• attr(*, "call")= language vegdist(x = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.5842, 0, 0, 0, 0, 0, 0, 0, 0| truncated ...

str(bray_sp)
'dist' num [1:3828] 0.543 0.424 0.421 0.942 0.975 ...

• attr(*, "Size")= int 88
• attr(*, "Labels")= chr [1:88] "29718" "32912" "30226" "29507" ...
• attr(*, "Diag")= logi FALSE
• attr(*, "Upper")= logi FALSE
• attr(*, "method")= chr "bray"
• attr(*, "call")= language vegdist(x = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.5842, 0, 0, 0, 0, 0, 0, 0, 0| truncated ...

sp.gdf<-geomorph.data.frame( bray_sp=bray_sp, jacc_sp=jacc_sp,Species=Species)

attributes(sp.gdf)
$names
[1] "bray_sp" "jacc_sp" "Species"

$class
[1] "geomorph.data.frame"

Jacc_sp works with no problem

proc_sp <- procD.lm(jacc_sp~Species,data=sp.gdf)

Sums of Squares calculations: 1000 permutations.
|==================================================================================================================| 100%

proc_sp

Call:
procD.lm(f1 = jacc_sp ~ Species, data = sp.gdf)

Type I (Sequential) Sums of Squares and Cross-products
Randomized Residual Permutation Procedure Used
1000 Permutations
ANOVA effect sizes and P-values based on empirical F distributions

      Df     SS      MS     Rsq      F      Z Pr(>F)   

Species 8 10.719 1.33981 0.34705 5.2487 13.351 0.001 **
Residuals 79 20.166 0.25526 0.65295
Total 87 30.884

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

bray_sp gets an error

proc_sp_b <- procD.lm(bray_sp~Species,data=sp.gdf)
Error in cmd$points[, p] : subscript out of bounds

Again jacc_sp works without an error

advproc_sp <- advanced.procD.lm(jacc_sp~Species, ~1, groups = ~Species, data = sp.gdf)

Obtaining initial full model fit:
Preliminary Model Fit...
|==================================================================================================================| 100%

Coefficients estimation: 1000 permutations.
|==================================================================================================================| 100%

Sums of Squares calculations: 1000 permutations.
|==================================================================================================================| 100%

Obtaining initial reduced model fit:
Preliminary Model Fit...
|==================================================================================================================| 100%

No terms for ANOVA; only RSS calculated in each permutation

Coefficients estimation: 1000 permutations.
|==================================================================================================================| 100%

Sums of Squares calculations: 1000 permutations.
|==================================================================================================================| 100%

Performing ANOVA:
Sums of Squares calculations for 1 models: 1000 permutations.
|==================================================================================================================| 100%

Calculating LS means:
Coefficients estimation from RRPP on null model: 1000 permutations.
|==================================================================================================================| 100%

while bray_sp gets an error

advproc_sp_b <- advanced.procD.lm(bray_sp~Species, ~1, groups = ~Species, data = sp.gdf)
Error in cmd$points[, p] : subscript out of bounds

Any thoughts?

Thank you,

Carly

[mlcollyer]

Carly, While we try to adjust for negative eigen values, it is possible that the adjustment does not work. If you are sure that using bray distance produces real eigen vectors, you could use the PCoA scores instead, and I would like to maybe see your data to determine if there is a problem with our coding. If there are no real eigenvectors (no positive eigen values), this could be the problem; i.e., there is no mathematical solution. You might consider using adonis in vegan just to know if it computes a result. With no knowledge of data type and what appears to be trying all kinds of distances with your data, it is hard to say what is going on. The error suggests that the number of real eigenvectors is less than 1. Best, Mike

…

On Feb 4, 2019, at 1:23 PM, Carly Muletz Wolz ***@***.***> wrote: I am using procD.lm to determine if distances matrices differ among species, and then advanced.procD.lm. These are community analyses using jaccard distances and bray-curtis distances. Jaccard distance is fine and same dimension as Bray distances. But, bray distances is returning the Error in cmd$points[, p] : subscript out of bounds. These are all complex objects to share, so thought I would just share the code and output right now. I can spend more time getting you the data so that you could reproduce the error. jacc_sp <- distance(GM_data_compare, "jaccard", binary = T) bray_sp <- distance(GM_data_compare, "bray") Species <- get_variable(GM_data_compare, "Species2") Species [1] Rhesus macaque Rhesus macaque Rhesus macaque Rhesus macaque Human [6] Mantled howler monkey Rhesus macaque Rhesus macaque Rhesus macaque Rhesus macaque [11] Mantled howler monkey Rhesus macaque Rhesus macaque Rhesus macaque Rhesus macaque [16] Rhesus macaque Mantled howler monkey Rhesus macaque Rhesus macaque Rhesus macaque [21] Rhesus macaque Rhesus macaque Rhesus macaque Rhesus macaque Mantled howler monkey [26] Rhesus macaque Rhesus macaque Human Human Rhesus macaque [31] Rhesus macaque Human Human Rhesus macaque Rhesus macaque [36] Rhesus macaque Western lowland gorilla Owl monkey Rhesus macaque Human [41] Human Rhesus macaque Sumatran orangutan Rhesus macaque Human [46] Rhesus macaque Human Owl monkey Rhesus macaque Human [51] Human Western lowland gorilla Human Human Sumatran orangutan [56] Human Human Bornean orangutan Human Bornean orangutan [61] Human Human Human Human Human [66] Western lowland gorilla Human Western lowland gorilla Human Bornean orangutan [71] Western lowland gorilla Human Human Bonobo Western lowland gorilla [76] Chimpanzee Western lowland gorilla Human Human Bornean orangutan [81] Bornean orangutan Bonobo Mantled howler monkey Human Western lowland gorilla [86] Western lowland gorilla Chimpanzee Mantled howler monkey 9 Levels: Human Bonobo Chimpanzee Western lowland gorilla Bornean orangutan Sumatran orangutan ... Owl monkey str(jacc_sp) 'dist' num [1:3828] 0.612 0.521 0.54 0.965 0.936 ... attr(*, "Size")= int 88 attr(*, "Labels")= chr [1:88] "29718" "32912" "30226" "29507" ... attr(*, "Diag")= logi FALSE attr(*, "Upper")= logi FALSE attr(*, "method")= chr "binary jaccard" attr(*, "call")= language vegdist(x = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.5842, 0, 0, 0, 0, 0, 0, 0, 0| truncated ... str(bray_sp) 'dist' num [1:3828] 0.543 0.424 0.421 0.942 0.975 ... attr(*, "Size")= int 88 attr(*, "Labels")= chr [1:88] "29718" "32912" "30226" "29507" ... attr(*, "Diag")= logi FALSE attr(*, "Upper")= logi FALSE attr(*, "method")= chr "bray" attr(*, "call")= language vegdist(x = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.5842, 0, 0, 0, 0, 0, 0, 0, 0| truncated ... sp.gdf<-geomorph.data.frame( bray_sp=bray_sp, jacc_sp=jacc_sp,Species=Species) attributes(sp.gdf) $names [1] "bray_sp" "jacc_sp" "Species" $class [1] "geomorph.data.frame" Jacc_sp works with no problem proc_sp <- procD.lm(jacc_sp~Species,data=sp.gdf) Sums of Squares calculations: 1000 permutations. |==================================================================================================================| 100% proc_sp Call: procD.lm(f1 = jacc_sp ~ Species, data = sp.gdf) Type I (Sequential) Sums of Squares and Cross-products Randomized Residual Permutation Procedure Used 1000 Permutations ANOVA effect sizes and P-values based on empirical F distributions Df SS MS Rsq F Z Pr(>F) Species 8 10.719 1.33981 0.34705 5.2487 13.351 0.001 ** Residuals 79 20.166 0.25526 0.65295 Total 87 30.884 Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 bray_sp gets an error proc_sp_b <- procD.lm(bray_sp~Species,data=sp.gdf) Error in cmd$points[, p] : subscript out of bounds Again jacc_sp works without an error advproc_sp <- advanced.procD.lm(jacc_sp~Species, ~1, groups = ~Species, data = sp.gdf) Obtaining initial full model fit: Preliminary Model Fit... |==================================================================================================================| 100% Coefficients estimation: 1000 permutations. |==================================================================================================================| 100% Sums of Squares calculations: 1000 permutations. |==================================================================================================================| 100% Obtaining initial reduced model fit: Preliminary Model Fit... |==================================================================================================================| 100% No terms for ANOVA; only RSS calculated in each permutation Coefficients estimation: 1000 permutations. |==================================================================================================================| 100% Sums of Squares calculations: 1000 permutations. |==================================================================================================================| 100% Performing ANOVA: Sums of Squares calculations for 1 models: 1000 permutations. |==================================================================================================================| 100% Calculating LS means: Coefficients estimation from RRPP on null model: 1000 permutations. |==================================================================================================================| 100% while bray_sp gets an error advproc_sp_b <- advanced.procD.lm(bray_sp~Species, ~1, groups = ~Species, data = sp.gdf) Error in cmd$points[, p] : subscript out of bounds Any thoughts? Thank you, Carly — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#22>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AGjacqihYQ-tualJ4RLSZV9QNiSurjEIks5vKHqNgaJpZM4ahxAv>.

[CarlyMuletzWolz]

Mike,

I also forgot to mention that it was working fine a year ago. That exact same test and same data did not give this error message. I am guessing I updated the package since then.

I can do the adonis using vegan and I get the exact same answer for jaccard distances between adonis and procD.lm.

adonis(jacc_sp ~ Species)

Call:
adonis(formula = jacc_sp ~ Species)

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    

Species 8 10.719 1.33981 5.2487 0.34705 0.001 ***
Residuals 79 20.166 0.25526 0.65295
Total 87 30.884 1.00000

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

proc_sp <- procD.lm(jacc_sp~Species,data=sp.gdf2)

Sums of Squares calculations: 1000 permutations.
|==================================================================================================================| 100%

proc_sp

Call:
procD.lm(f1 = jacc_sp ~ Species, data = sp.gdf2)

Type I (Sequential) Sums of Squares and Cross-products
Randomized Residual Permutation Procedure Used
1000 Permutations
ANOVA effect sizes and P-values based on empirical F distributions

      Df     SS      MS     Rsq      F      Z Pr(>F)   

Species 8 10.719 1.33981 0.34705 5.2487 13.351 0.001 **
Residuals 79 20.166 0.25526 0.65295
Total 87 30.884

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

adonis for bray curtis distances does works

adonis(bray_sp ~ Species)

Call:
adonis(formula = bray_sp ~ Species)

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    

Species 8 13.407 1.67592 7.5931 0.43468 0.001 ***
Residuals 79 17.437 0.22072 0.56532
Total 87 30.844 1.00000

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

I can get bray to work if I convert the dist object to a matrix (as.matrix()), but it does not agree with the adonis results (and the same happens if I do that with jaccard - results are different than adonis). It seems that a dist object is the appropriate choice for the response matrix 'y'.

Attached are the files as .csvs. I exported as a matrix, so need to be reimported as as.dist() to reproduce error.
Archive.zip

Thank you,

Carly

[mlcollyer]

Carly, I found the issue and fixed it. The problem was with the R function, cmdscale, which performs the multidimensional scaling of the distance matrix. It returned a different number of eigenvectors and eigenvalues (which was triggered because of negative eigenvalues). I don’t know why it used to work and recently started to return a different number of eigenvectors and eigenvalues (I don’t believe we altered its use). Regardless, I updated our code to hopefully catch these instances and avoid issues. Should work now (after you update geomorph from Github). Cheers! Mike

…

On Feb 5, 2019, at 3:39 PM, Carly Muletz Wolz ***@***.***> wrote: Mike, I also forgot to mention that it was working fine a year ago. That exact same test and same data did not give this error message. I am guessing I updated the package since then. I can do the adonis using vegan and I get the exact same answer for jaccard distances between adonis and procD.lm. adonis(jacc_sp ~ Species) Call: adonis(formula = jacc_sp ~ Species) Permutation: free Number of permutations: 999 Terms added sequentially (first to last) Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) Species 8 10.719 1.33981 5.2487 0.34705 0.001 *** Residuals 79 20.166 0.25526 0.65295 Total 87 30.884 1.00000 Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 proc_sp <- procD.lm(jacc_sp~Species,data=sp.gdf2) Sums of Squares calculations: 1000 permutations. |==================================================================================================================| 100% proc_sp Call: procD.lm(f1 = jacc_sp ~ Species, data = sp.gdf2) Type I (Sequential) Sums of Squares and Cross-products Randomized Residual Permutation Procedure Used 1000 Permutations ANOVA effect sizes and P-values based on empirical F distributions Df SS MS Rsq F Z Pr(>F) Species 8 10.719 1.33981 0.34705 5.2487 13.351 0.001 ** Residuals 79 20.166 0.25526 0.65295 Total 87 30.884 Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 adonis for bray curtis distances works adonis(bray_sp ~ Species) Call: adonis(formula = bray_sp ~ Species) Permutation: free Number of permutations: 999 Terms added sequentially (first to last) Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) Species 8 13.407 1.67592 7.5931 0.43468 0.001 *** Residuals 79 17.437 0.22072 0.56532 Total 87 30.844 1.00000 Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1 I can get bray to work if I convert the dist object to a matrix (as.matrix()), but it does not agree with the adonis results (and the same happens if I do that with jaccard - results are different than adonis). It seems that a dist object is the appropriate choice for the response matrix 'y'. Attached are the files as .csvs. I exported as a matrix, so need to be reimported as as.dist() to reproduce error. Archive.zip <https://github.com/geomorphR/geomorph/files/2833765/Archive.zip> Thank you, Carly — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#22 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AGjacuURIcjEmNT7kyiolWUjM1qL7kgqks5vKev3gaJpZM4ahxAv>.