Wednesday, November 29, 2017

Object class & S3 methods for multiple Mantel test in phytools

The phytools package contains a small function, multi.mantel, for multiple matrix regression. The idea of this function is merely to fit a model in which distance or correlation matrix Y ~ X1 + X2 and so on, but then obtain our p-values for both the full model & the model coefficients via (Mantel) permutations of the rows & columns together in Y.

A colleague today contacted me with some questions about the function & so I decided it needed a small re-boot. None of the internal operation of the function has been changed; however, I have now updated it with a new object class, as well as print, residuals, and fitted S3 methods.

For fun, this is what the print method looks like:

print.multi.mantel<-function(x,...){
    if(hasArg(digits)) digits<-list(...)$digits
    else digits<-6
    star<-function(p){
        obj<-if(p>0.1) "" else if(p<=0.1&&p>0.05) "." else
            if(p<=0.05&&p>0.01) "*" else if(p<=0.01&&p>0.001) "**" else
            if(p<=0.001) "***"
        obj
    }
    cat("\nResults from a (multiple) Mantel regression using \"multi.mantel\":\n\n")
    cat("Coefficients:\n")
    object<-data.frame(x$coefficients,
        x$tstatistic,x$probt,
        sapply(x$probt,star))
    rownames(object)<-names(x$coefficients)
    colnames(object)<-c("Estimate","t value","Pr(>|t|)","")
    print(object)
    cat("---\n")
    cat(paste("Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01", 
        "\u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1\n"))
    cat(paste("Pr(>|t|) based on",x$nperm,
        "(Mantel) permutations of rows & columns together in Y.\n\n"))
    cat(paste("Multiple R-squared:",round(x$r.squared,digits),"\n"))
    cat(paste("F-statistic: ",round(x$fstatistic,digits),
        ", p-value (based on ",x$nperm," permutations): ",
        x$probF,"\n\n",sep=""))
}

Here's a small example. It takes as arguments a single dependent matrix, Y, and either a single independent matrix or a list of such matrices, X. Either matrices or class "dist" objects are acceptable. Here I use the latter.

library(phytools)
packageVersion("phytools")
## [1] '0.6.48'
dY
##             1          2          3          4          5          6
## 2  0.54603980                                                       
## 3  1.69212862 2.23816843                                            
## 4  0.16256209 0.70860190 1.52956653                                 
## 5  0.58192166 1.12796146 1.11020697 0.41935957                      
## 6  0.83933577 0.29329596 2.53146439 1.00189786 1.42125743           
## 7  0.58679271 0.04075291 2.27892134 0.74935480 1.16871437 0.25254306
## 8  1.75571440 2.30175420 0.06358577 1.59315231 1.17379274 2.59505017
## 9  0.93129804 1.47733784 0.76083058 0.76873595 0.34937638 1.77063381
## 10 0.97782473 1.52386454 0.71430389 0.81526264 0.39590307 1.81716050
##             7          8          9
## 2                                  
## 3                                  
## 4                                  
## 5                                  
## 6                                  
## 7                                  
## 8  2.34250711                      
## 9  1.51809075 0.82441636           
## 10 1.56461744 0.77788967 0.04652669
dX
## [[1]]
##             1          2          3          4          5          6
## 2  1.75416715                                                       
## 3  0.63607244 1.11809471                                            
## 4  3.05970185 1.30553469 2.42362940                                 
## 5  0.77190528 2.52607244 1.40797773 3.83160713                      
## 6  0.58102311 1.17314404 0.05504933 2.47867874 1.35292839           
## 7  2.09863486 0.34446771 1.46256242 0.96106698 2.87054015 1.51761176
## 8  1.45771945 0.29644770 0.82164701 1.60198239 2.22962474 0.87669634
## 9  0.46813108 1.28603607 0.16794136 2.59157077 1.24003636 0.11289203
## 10 1.23194515 0.52222200 0.59587271 1.82775670 2.00385043 0.65092204
##             7          8          9
## 2                                  
## 3                                  
## 4                                  
## 5                                  
## 6                                  
## 7                                  
## 8  0.64091541                      
## 9  1.63050378 0.98958837           
## 10 0.86668972 0.22577430 0.76381407
## 
## [[2]]
##             1          2          3          4          5          6
## 2  0.54790843                                                       
## 3  1.16512317 1.71303160                                            
## 4  0.26507653 0.28283189 1.43019970                                 
## 5  0.95864760 1.50655603 0.20647557 1.22372414                      
## 6  0.36213942 0.18576901 1.52726259 0.09706289 1.32078703           
## 7  0.03566574 0.51224269 1.20078891 0.22941080 0.99431334 0.32647369
## 8  1.29414610 1.84205453 0.12902294 1.55922264 0.33549850 1.65628553
## 9  1.73415785 2.28206628 0.56903468 1.99923439 0.77551025 2.09629728
## 10 0.30273318 0.24517525 1.46785635 0.03765665 1.26138078 0.05940624
##             7          8          9
## 2                                  
## 3                                  
## 4                                  
## 5                                  
## 6                                  
## 7                                  
## 8  1.32981184                      
## 9  1.76982359 0.44001175           
## 10 0.26706744 1.59687928 2.03689103
fit<-multi.mantel(dY,dX)
fit
## 
## Results from a (multiple) Mantel regression using "multi.mantel":
## 
## Coefficients:
##               Estimate   t value Pr(>|t|)  
## (intercept)  0.9667831  4.196181    0.618  
## X1          -0.1771598 -1.579220    0.103  
## X2           0.3950495  2.802433    0.023 *
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## Pr(>|t|) based on 1000 (Mantel) permutations of rows & columns together in Y.
## 
## Multiple R-squared: 0.213525 
## F-statistic: 5.701406, p-value (based on 1000 permutations): 0.028
residuals(fit)
##               1            2            3            4            5
## 2  -0.326426449                                                    
## 3   0.377750562  0.792734343                                       
## 4  -0.366883315 -0.138625616  0.427153292                          
## 5  -0.626824208  0.013532473  0.311292771 -0.352048579             
## 6  -0.167576456 -0.539041195  0.971089391  0.435892244  0.172382466
## 7  -0.022286474 -1.067365644  1.096874319 -0.137794548  0.317672448
## 8   0.535928683  0.659786867 -0.808605010  0.294205815  0.469470892
## 9  -0.637629371 -0.163140682 -0.400996988 -0.528721733 -0.704087161
## 10  0.109698119  0.552741765 -0.726790552  0.157408232 -0.714186278
##               6            7            8            9
## 2                                                     
## 3                                                     
## 4                                                     
## 5                                                     
## 6                                                     
## 7  -0.574353601                                       
## 8   1.129267512  0.963926842                          
## 9  -0.004290689  0.140999293 -0.140877962             
## 10  0.942226161  0.645871995 -0.779741779 -1.589612195
fitted(fit)
##            1         2         3         4         5         6         7
## 2  0.8724663                                                            
## 3  1.3143781 1.4454341                                                  
## 4  0.5294454 0.8472275 1.1024132                                        
## 5  1.2087459 1.1144290 0.7989142 0.7714081                              
## 6  1.0069122 0.8323372 1.5603750 0.5660056 1.2488750                    
## 7  0.6090792 1.1081186 1.1820470 0.8871494 0.8510419 0.8268967          
## 8  1.2197857 1.6419673 0.8721908 1.2989465 0.7043218 1.4657827 1.3785803
## 9  1.5689274 1.6404785 1.1618276 1.2974577 1.0534635 1.7749245 1.3770915
## 10 0.8681266 0.9711228 1.4410944 0.6578544 1.1100894 0.8749343 0.9187454
##            8         9
## 2                     
## 3                     
## 4                     
## 5                     
## 6                     
## 7                     
## 8                     
## 9  0.9652943          
## 10 1.5576314 1.6361389

You get the idea.

I simulated the distance matrices for the preceding as follows:

dY<-dist(y<-rnorm(10))
dX<-list(dist(rnorm(10)),dist(sqrt(0.5)*y+rnorm(10,sd=0.5)))

That's it.

1 comment:

  1. Liam, this is a really useful function. Thanks! Do you have any recommendations for how to use it for variation partitioning if you provide it with a list of independent X matrices?

    ReplyDelete