# R program for additional multivariate exercise 14 (VHM 802)
# only includes code for simulation-based ANOVA

library(vegan)

# expanded prehistoric dogs data
prehdog2 <- read.csv("r:/prehist_dog2.csv", header=TRUE)
prehdog2std <- scale(prehdog2[,2:10])

# anosim
prehdog2.simF.eucl <- anosim( prehdog2std, prehdog2$species, distance="euclidean") 
summary(prehdog2.simF.eucl)
prehdog2.simF.manh <- anosim( prehdog2std, prehdog2$species, distance="manhattan") 
summary(prehdog2.simF.manh)

# permanova, using adonis
prehdog2.permF.eucl <- adonis (prehdog2std~as.factor(species), data=prehdog2, method="euclidean")
print(prehdog2.permF.eucl)
prehdog2.permF.manh <- adonis (prehdog2std~as.factor(species), data=prehdog2, method="manhattan")
print(prehdog2.permF.manh)

prehdog2sex <- prehdog2[prehdog2$sex>0,]
prehdog2sexstd <- scale(prehdog2sex[,2:10])
# permanova, using adonis
prehdog2sex.permF.eucl <- adonis (prehdog2sexstd~as.factor(species)*as.factor(sex), data=prehdog2sex, method="euclidean")
print(prehdog2sex.permF.eucl)
prehdog2sex.permF.manh <- adonis (prehdog2sexstd~as.factor(species)*as.factor(sex), data=prehdog2sex, method="manhattan")
print(prehdog2sex.permF.manh)

# Steneryd data
steneryd <- read.csv("r:/steneryd.csv", header=TRUE, row.names=1)
sten.spec <- steneryd[,1:25]
sten.specstd <- scale(sten.spec)

# permanova, using adonis and default (bray-curtis) distance
steneryd.permF.bc <- adonis( sten.spec ~ light+moisture+reaction+nitrogen, data=steneryd)
print(steneryd.permF.bc)

# permanova, using adonis2 and marginal tests, and default (bray-curtis) distance
steneryd.permF.bc2 <- adonis2( sten.spec ~ light+moisture+reaction+nitrogen, by="margin", data=steneryd)
print(steneryd.permF.bc2)

# permanova, using adonis and euclidean distances for correlations
steneryd.permF.eucl <- adonis( sten.specstd ~ light+moisture+reaction+nitrogen, data=steneryd, method="euclidean")
print(steneryd.permF.eucl)

# permanova, using adonis2 and marginal tests, and euclidean distances for correlations
steneryd.permF.eucl2 <- adonis2( sten.specstd ~ light+moisture+reaction+nitrogen, data=steneryd, by="margin", method="euclidean")
print(steneryd.permF.eucl2)