[banner]

Summary and Analysis of Extension Program Evaluation in R

Salvatore S. Mangiafico

Permutation Tests for Linear Models

Packages used in this chapter

 

The packages used in this chapter include:

•  lmPerm

•  vegan

 

The following commands will install these packages if they are not already installed:


if(!require(lmPerm)){install.packages("lmPerm")}
if(!require(vegan)){install.packages("vegan")}


Permutation tests for linear models with the lmPerm package

 

The lmPerm package can be used to test coefficients from a linear model with a permutation test.  By default, non-sequential sums of squares (that is, type 2 or type 3 sums of squares) are returned by the anova() call. Model p-value and r-squared can be reported with the summary() call.  At the time of writing, I don’t know if there are methods for post-hoc tests.  Residuals and predicted values can be asked for from the model object.

 

This example uses the Lung Capacity data set, but we’ll truncate the data set to the first 200 observations, just to get more interesting results.


LungCapData = read.csv(header=TRUE, stringsAsFactors = TRUE,
                       file="https://rcompanion.org/documents/LungCapData.csv")

LungCapData = LungCapData[1:200,]

library(lmPerm)

model = lmp(LungCap ~ Age + Smoke + Gender + Smoke:Gender, data=LungCapData)

anova(model)


Analysis of Variance Table

Response: LungCap
              Df R Sum Sq R Mean Sq Iter Pr(Prob)   
Age            1   968.59    968.59 5000   <2e-16 ***
Smoke          1     3.11      3.11  357   0.2213   
Gender         1    17.51     17.51 5000   <2e-16 ***
Smoke:Gender   1     0.10      0.10   51   0.7059   
Residuals    195   364.53      1.87


summary(model)


Residual standard error: 1.367 on 195 degrees of freedom
Multiple R-Squared: 0.7377,   Adjusted R-squared: 0.7323
F-statistic: 137.1 on 4 and 195 DF,  p-value: < 2.2e-16


Permutation tests for linear models with the vegan package

 

The vegan package is tailor-made for applications in ecology. The adonis2() function can test model parameters with permutation tests for models with either a univariate or multivariate left-hand side (dependent variable).  The vegdist() function is used to create a distance matrix for the dependent variable(s). This is then used as the left-hand side of the model formula.

 

This example uses the Lung Capacity data set, but again we’ll truncate the data set to the first 200 observations to compare results to the lmPerm results.

 

Here, model2 uses the by="margin" option to ask for non-sequential sums of squares (that is, type 2 or type 3 sums of squares) for the model terms.  Note that the p-value for Smoke:Gender is returned, but not for Smoke or Gender individually.


LungCapData = read.csv(header=TRUE, stringsAsFactors = TRUE,
                       file="https://rcompanion.org/documents/LungCapData.csv")

LungCapData = LungCapData[1:200,]

library(vegan)

Yvar   = LungCapData$LungCap

LHS = vegdist(Yvar)

model2 = adonis2(LHS ~ Age + Smoke + Gender + Smoke:Gender, data=LungCapData,
                 permutations=1e4,
                 by="margin") 

model2


Permutation test for adonis under reduced model
Marginal effects of terms

              Df SumOfSqs      R2        F    Pr(>F)   
Age            1   3.5762 0.51068 216.2090 9.999e-05 ***
Smoke:Gender   1   0.0030 0.00043   0.1817    0.7763   
Residual     195   3.2254 0.46058                      
Total        199   7.0029 1.00000   00


The overall p-value and r-squared value for the overall model can be obtained by fitting a model with the by=NULL option.


model3 = adonis2(LHS ~ Age + Smoke + Gender + Smoke:Gender, data=LungCapData,
                 permutations=1e4,
                 by=NULL)

model3


Permutation test for adonis under reduced model

          Df SumOfSqs      R2      F    Pr(>F)   
Model      4   3.7775 0.53942 57.094 9.999e-05 ***
Residual 195   3.2254 0.46058                    
Total    199   7.0029 1.00000


Tests for sequentially-added terms will be reported with the by="term" option. Note that the p-value for all of Smoke, Gender, and Smoke:Gender is returned.


model4 = adonis2(LHS ~ Age + Smoke + Gender + Smoke:Gender, data=LungCapData,
                 permutations=1e4,
                 by="term")

model4


Permutation test for adonis under reduced model
Terms added sequentially (first to last)

Df SumOfSqs      R2        F    Pr(>F)   
Age            1   3.5988 0.51390 217.5745 9.999e-05 ***
Smoke          1   0.0207 0.00296   1.2531    0.2629   
Gender         1   0.1549 0.02212   9.3664    0.0004 ***
Smoke:Gender   1   0.0030 0.00043   0.1817    0.7797   
Residual     195   3.2254 0.46058                      
Total        199   7.0029 1.00000