 Repeated measures Permanova and mean group differences for spectral data using PAVO and adonis function (vegan)

Hello everyone,

We are studying the iridescence (changes with viewing angle) of dorsal coloration in lizards. We took reflectance spectra at three different angles (0º, 60º, and 90º between illumination and viewing points).
My dataset takes the following form. The first column (wl) corresponds to the waveband visible to lizards (from 300 to 700 nm). Then, each column correponds to the reflectance of a lizard dorsum measured at three different angles (i.e. 3 measures per lizard). Each column except the first is hence labelled with the ID of the lizard (e.g. F04) and the viewing angle (0º, 60º, 90º), separated by a "_". This is an extract from my dataset.

#    wl      F04_90      F04_0     F04_60      F05_90      F05_0      F05_60      F22_90      F22_0     F22_60     F26_90       F26_0    F26_60      F32_0      F32_90
# 1 300 0.003454688 0.01715481 0.03353053 0.002307865 0.02507565 0.006666404 0.001536308 0.01859693 0.09439136 0.03034757 0.006715812 0.1117221 0.01111177 0.005094258
# 2 301 0.003454688 0.01715481 0.03353053 0.002307865 0.02507565 0.006666404 0.001536308 0.01859693 0.09439136 0.03034757 0.006715812 0.1117221 0.01111177 0.005094258
# 3 302 0.003454688 0.01715481 0.03353053 0.002307865 0.02507565 0.006666404 0.001536308 0.01859693 0.09439136 0.03034757 0.006715812 0.1117221 0.01111177 0.005094258
# 4 303 0.003454688 0.01715481 0.03353053 0.002307865 0.02507565 0.006666404 0.001536308 0.01859693 0.09439136 0.03034757 0.006715812 0.1117221 0.01111177 0.005094258
# 5 304 0.003454688 0.01715481 0.03353053 0.002307865 0.02507565 0.006666404 0.001536308 0.01859693 0.09439136 0.03034757 0.006715812 0.1117221 0.01111177 0.005094258
# 6 305 0.003454688 0.01715481 0.03353053 0.002307865 0.02507565 0.006666404 0.001536308 0.01859693 0.09439136 0.03034757 0.006715812 0.1117221 0.01111177 0.005094258

We want to estimate if the dorsal coloration differs with viewing angle. We tried to adapt the analysis suggested [here] (Electronic Supplementary Material for Comparing colours using visual models), but taking account of the repeated measures nature of our data.
First, we used the package PAVO to calculate the cone stimulation that each spectra elicits on each of the four types of cones present in the retina of lizards.

############ Visual model ##############

# Creating a visual model with the sensitivities of each of the four cones present in the retina of lizards

library(PAVO)
modelo<- sensmodel(peaksens=c(367, 457, 497, 562))
names(modelo)=c("wl", "u", "s", "m", "l")

modelo
# obtain the cone stimulation elicited by each spectra

vis<-vismodel(specs, visual=modelo, achromatic = c(562),illum = "D65", bkg = "ideal", relative=FALSE)
u           s           m          l      lum
F04_90 3.388424e-05 0.004203633 0.009196292 0.01474442 2098.286
F04_0  1.012247e-04 0.007537306 0.015045455 0.03407671 6530.361
F04_60 8.747567e-05 0.003736701 0.009886441 0.02395551 3664.618
F05_90 2.294536e-05 0.003674388 0.007913641 0.01324431 1973.242
F05_0  1.343674e-04 0.008426027 0.015407409 0.03125477 6076.853
F05_60 4.970419e-05 0.006446567 0.014428640 0.03139251 4967.756

# Obtain chromatic distances between each of the spectra

dis<-coldist(vis, achromatic=TRUE, n= c(1,1,1,4), weber=0.05)
head(dis) # Object with the chromatic (dS) and achromatic (dL) distances between each possible pairwise combination of spectra.

patch1 patch2       dS        dL
1 F04_90  F04_0 4.825718 11.353415
2 F04_90 F04_60 8.396470  5.576033
3 F04_90 F05_90 2.542436  0.614429
4 F04_90  F05_0 6.593240 10.633663
5 F04_90 F05_60 4.416959  8.618474
6 F04_90 F22_90 3.440865 10.924385

The unit of these chromatic distances are JND (just noticeable distances). Conventionally, a JND of 3 is considered large enough for an animal to discriminate between the colors even under dim light.

REPEATED MEASURES PERMANOVA

# Create dissimilarity matrix with the chromatic distances (dS) between every spectra

mat1 <- dist(coldist2mat(dis)[["dS"]])
mat1

F04_90       F04_0      F04_60      F05_90       F05_0      F05_60      F22_90       F22_0      F22_60      F26_90       F26_0      F26_60       F32_0
F04_0   24.5175568
F04_60  36.7366106  19.3193443
F05_90  14.3435559  35.5460223  45.0172783
F05_0   34.2073410  17.7556758  16.3868454  45.1874828
F05_60  18.6681426  33.0930678  41.5362896  11.4778527  44.3211641

# Create another dataframe (pred) with the factors "angle" and "ID" to fit in adonis and pairwise.adonis2 (vegan)

angle <- substring(rownames(as.matrix(mat1)),5,5) # Retain the 5º character of each column name in mat1.
angle

ID <- substring(rownames(as.matrix(mat1)),1,3) # Retain the first three characters of each columnname in mat1.
ID
pred<-cbind(ID,angle)

pred<-as.data.frame(pred)

# Test homogeneity of dispersion among angles

var1 <- betadisper(mat1, group=pred\$angle, type="centroid") #Esto es para comprobar la homogeneidad de varianzas.
permutest(var1)
var1 # Results show homogeonous dispersion

# Perform Permanova to test (within each individual) if perceived dorsal coloration varies chromatically with viewing angles (iridescence).

perm<-how(nperm=999)
setBlocks(perm) <- with(pred, ID) # Restrict permutations within the repeated measures of each individual (we are not interested in how the coloration of lizard A at 0º differs from coloration of lizard B at 60º).

manova1<- adonis2(mat1~angle, data=pred, contr.sum=TRUE, permutations = perm)
manova1

# Permutation test for adonis under reduced model
# Terms added sequentially (first to last)
# Blocks:  with(pred, ID)
# Permutation: free
# Number of permutations: 999

# adonis2(formula = mat1 ~ angle, data = pred, permutations = perm, contr.sum = TRUE)
#          Df SumOfSqs      R2      F Pr(>F)
# angle     2     7543 0.16212 4.6438  0.003 **
# Residual 48    38983 0.83788
# Total    50    46526 1.00000

# Now for each pairwise contrast between the different angles, but taking account of the repeated measures (within-individual) nature of our data.

pairwise.adonis2 <- function(x, data, strata = NULL, nperm=999, ... ) {

#describe parent call function
ststri <- ifelse(is.null(strata),'Null',strata)
fostri <- as.character(x)
#list to store results

#copy model formula
x1 <- x
# extract left hand side of formula
lhs <- x1[]
# extract factors on right hand side of formula
rhs <- x1[]
# create model.frame matrix
x1[] <- NULL
rhs.frame <- model.frame(x1, data, drop.unused.levels = TRUE)

# create unique pairwise combination of factors
co <- combn(unique(as.character(rhs.frame[,1])),2)

# create names vector
nameres <- c('parent_call')
for (elem in 1:ncol(co)){
nameres <- c(nameres,paste(co[1,elem],co[2,elem],sep='_vs_'))
}
#create results list
res <- vector(mode="list", length=length(nameres))
names(res) <- nameres

res['parent_call'] <- list(paste(fostri,fostri,fostri,', strata =',ststri, ', permutations',nperm ))

#start iteration trough pairwise combination of factors
for(elem in 1:ncol(co)){

#reduce model elements
if(inherits(eval(lhs),'dist')){
xred <- as.dist(as.matrix(eval(lhs))[rhs.frame[,1] %in% c(co[1,elem],co[2,elem]),
rhs.frame[,1] %in% c(co[1,elem],co[2,elem])])
}else{
xred <- eval(lhs)[rhs.frame[,1] %in% c(co[1,elem],co[2,elem]),]
}

mdat1 <-  data[rhs.frame[,1] %in% c(co[1,elem],co[2,elem]),]

# redefine formula
if(length(rhs) == 1){
xnew <- as.formula(paste('xred',as.character(rhs),sep='~'))
}else{
xnew <- as.formula(paste('xred' ,
paste(rhs[-1],collapse= as.character(rhs)),
sep='~'))}

if(is.null(strata)){
}else{
perm <- how(nperm = nperm)
setBlocks(perm) <- with(mdat1, mdat1[,ststri])

}
#names(res) <- names
return(res)
}

### Method summary
cat("\n")
print(object, ...)
cat("\n")

cat("Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n")
}

# Test on our data

pa1

# Output:

\$parent_call
 "mat1 ~ angle , strata = ID , permutations 999"

\$`9_vs_0`
Df SumOfSqs      R2    F Pr(>F)
angle     1   4557.9 0.22255 9.16  0.002 **
Residual 32  15922.7 0.77745
Total    33  20480.6 1.00000
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

\$`9_vs_6`
Df SumOfSqs      R2      F Pr(>F)
angle     1     1888 0.05276 1.7823  0.047 *
Residual 32    33892 0.94724
Total    33    35779 1.00000
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

\$`0_vs_6`
Df SumOfSqs      R2      F Pr(>F)
angle     1     4869 0.14745 5.5345  0.022 *
Residual 32    28151 0.85255
Total    33    33020 1.00000
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

attr(,"class")

My questions are:

1. Is this script correct? Are we effectively accounting for the repeated measures nature of our data?

2. We would like to obtain mean differences (with confidence intervals) between every pair of angles, although based exclusively on within-individual distances. Something based on the following:

boot <- bootcoldist(vis, by=pred\$angle, n=c(1,1,1,4), weber=0.05, weber.achro=0.05)

#     dS.mean   dS.lwr   dS.upr   dL.mean    dL.lwr    dL.upr
# 0-6 2.647124 1.163678 6.509701  6.183066  3.203439  9.275883
# 0-9 5.401745 4.072397 7.621412 17.700675 13.223363 22.452985
# 6-9 4.767297 3.767327 7.384998 11.517609  7.023582 16.149197

Where 0-6 stands for the mean difference between all spectra measured at 0º and all spectra measured at 60º, 0-9 = 0º-90º, 6-9 = 60º-90º.

Is there a way to code the argument in "by=" so that I can obtain a similar output although based exclusively on within-individual distances?

Thank you very much,

Javier

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.