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.
head(specs)
# 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)
head(vis)
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 from https://github.com/pmartinezarbizu/pairwiseAdonis/blob/master/pairwiseAdonis/R/pairwise.adonis2.R
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[[2]]
# extract factors on right hand side of formula
rhs <- x1[[3]]
# create model.frame matrix
x1[[2]] <- 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
#add parent call to res
res['parent_call'] <- list(paste(fostri[2],fostri[1],fostri[3],', 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[1])),
sep='~'))}
#pass new formula to adonis
if(is.null(strata)){
ad <- adonis2(xnew,data=mdat1, ... )
}else{
perm <- how(nperm = nperm)
setBlocks(perm) <- with(mdat1, mdat1[,ststri])
ad <- adonis2(xnew,data=mdat1,permutations = perm, ... )}
res[nameres[elem+1]] <- list(ad[1:5])
}
#names(res) <- names
class(res) <- c("pwadstrata", "list")
return(res)
}
### Method summary
summary.pwadstrata = function(object, ...) {
cat("Result of pairwise.adonis2:\n")
cat("\n")
print(object[1], ...)
cat("\n")
cat("Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n")
}
# Test on our data
pa1 <- pairwise.adonis2(mat1~angle, strata="ID", data=pred) # More info on https://www.researchgate.net/post/How_can_I_do_PerMANOVA_pairwise_contrasts_in_R
pa1
# Output:
$parent_call
[1] "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")
[1] "pwadstrata" "list"
My questions are:
-
Is this script correct? Are we effectively accounting for the repeated measures nature of our data?
-
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