Function for Thorpe normalization following Bartels, P. J., Nelson, D. R., & Exline, R. P. (2011). Allometry and the removal of body size effects in the morphometric analysis of tardigrades. Journal of Zoological Systematics and Evolutionary Research, 49, 17-25.
BL is the column number of the trait to normalize for (Body length or Buccal Tube length)
traits is a vector with the column numbers of traits to normalize
thorpe.normalization = function(dataset,BL,traits){
BLs = dataset[,BL]
L_BLs = log10(BLs)
other_traits = dataset[,traits]
L_other_traits = log10(other_traits)
avg_BS = mean(BLs, na.rm=T)
linear_models = list()
for (i in 1:ncol(L_other_traits)){
linear_models[[i]] = lm(L_other_traits[,i]~L_BLs, na.action = na.omit)
names(linear_models)[i] = colnames(L_other_traits)[i]
}
results_table = data.frame(trait = names(linear_models),
b = unlist(lapply(linear_models, function(x){x$coefficients[2]})),
pval_b = unlist(lapply(linear_models, function(x){summary(x)$coefficients[2,4]})),
r_squared = unlist(lapply(linear_models, function(x){summary(x)$r.squared})),
signif = unlist(lapply(linear_models, function(x){summary(x)$coefficients[2,4]}))<=0.05
)
Y_star = t(other_traits) * matrix(rep(avg_BS/BLs,length(traits)),nrow=length(traits), byrow=T) * matrix(rep(results_table$b, nrow(other_traits)), nrow=length(traits), byrow=F)
colnames(Y_star) = paste0("Y_star_Ind", 1:ncol(Y_star))
Y_star_t = t(Y_star)
a_star_models = list()
for (i in 1:ncol(Y_star_t)){
a_star_models[[i]] = lm(Y_star_t[,i]~BLs, na.action = na.omit)
names(a_star_models)[i] = colnames(Y_star_t)[i]
}
results_table$a_star = round(unlist(lapply(a_star_models,function(x){x$coefficients[1]})),4)
results_table = cbind(results_table,Y_star)
return(results_table)
}
We validate the function by testing if the function works fine and can reproduce the same values obtained by bartels et al., 2011
# Run the function on the Bartels raw data
data.Bartels = read.table("Bartels_data.txt",header=T, sep="\t")
bartels_thorpe = thorpe.normalization(data.Bartels,2,c(3:9,11:13,15:17,20,21))
# Load the Bartels coefficnts and compare them to the ones obtained with the R function
original_Bartels = read.table("Bartels_original_thorpe.txt",header=T, sep="\t")
Plot the original vs. the newly obtained values
plot(bartels_thorpe$b~original_Bartels$b,main="Allometric exponent (b)",xlab="Values from Bartels et al. 2011", ylab="Values from thorpe.normalization function")
abline(lm(bartels_thorpe$b~original_Bartels$b))
# Get the r^2
summary(lm(bartels_thorpe$b~original_Bartels$b))$r.squared
## [1] 1
For the allometric exponent (b) we obtain an almost perfect correspondence. It has to be consider however that rounding can lead to slightly different results between the original Excel calculations and the R function.
plot(bartels_thorpe$a_star~original_Bartels$a_star,main="Y intercept (a*)",xlab="Values from Bartels et al. 2011", ylab="Values from thorpe.normalization function")
abline(lm(bartels_thorpe$a_star~original_Bartels$a_sta))
# Get the r^2
summary(lm(bartels_thorpe$a_star~original_Bartels$a_sta))$r.squared
## [1] 0.9977482
For the Y intercept (a*) we obtain an almost perfect correspondence. It has to be consider however that rounding can lead to slightly different results between the original Excel calculations and the R function.
data.annewintersae = read.table("Mac_annewintersae.txt",header=T, sep="\t")
annewintersae_thorpe = thorpe.normalization(data.annewintersae,5,6:30)
data.rybaki = read.table("Mac_rybaki.txt",header=T, sep="\t")
rybaki_thorpe = thorpe.normalization(data.rybaki,5,6:30)
It is possible to download all the table as Excel files by clicking on the “Excel” button.