##18-10-2017
##statistical analysis note for paper in titled
##"Genetic variation and genotype by environment interaction in susceptibility 
##of Abies nordmanniana to the fungus Neonectria neomacrospora"

##read data
total = read.table("//a00143.science.domain/prw187/Documents/JINGXU/IGN/5.paper/p2.Neonectria genetics/8.revision/R statistical analysis/data storage upload/total.txt",
                  sep="\t", header=TRUE)
str(total)

##set charactor of variables
total<-transform(total, blok=as.factor(blok), clone=as.factor(clone),
                row=as.factor(row), col=as.factor(col))

##########Single site spatial analysis##############################

####FP.251###

##select data only for FP.251
tuse <- subset(total, site=="FP251")

##single univariate spatial analysis for FP.251

tuse.asr<- asreml(NNecScr ~ Edge , random = ~ blok+clone+units,
                  rcov = ~ ar1(row):ar1(col),
                  data = tuse, na.method.X='include',maxit=50)

##check model validation
plot(tuse.asr)

##test the significance of fixed effect
wald(tuse.asr)

##view variance components
summary(tuse.asr)$varcomp

##extract BLUP Etimates
BLUP<-tuse_spa$coefficients$random

##plot variogram of spatial analysis
plot(variogram(tuse.asr))

####FP.259###

###select data only for FP.259
silk<- subset(total, site=="FP259")

####single univariate spatial analysis for FP.251

silk.asr<- asreml(NNecScr ~ 1, random = ~ blok+clone+units,
                  rcov = ~ ar1(row):ar1(col),
                  data = silk, na.method.X='include',maxit=50)

##check model validation
plot(silk.asr)

##view variance components
summary(silk.asr)$varcomp

##extract BLUP Etimates
BLUP<-silk.asr$coefficients$random

##plot variogram of spatial analysis
plot(variogram(silk.asr))

####FP.272###
###select data only for FP.272
skel<- subset(total, site=="FP272")

####single univariate spatial analysis for FP.272

skel.asr<- asreml(NNecScr ~ 1, random = ~ blok+clone+units,
                  rcov = ~ ar1(row):ar1(col),
                  data = skel, na.method.X='include',maxit=50)

##check model validation
plot(skel.asr)

##view variance components
summary(skel.asr)$varcomp

##extract BLUP Etimates
BLUP<-skel.asr$coefficients$random

##plot variogram of spatial analysis
plot(variogram(skel.asr))

#################multi-site analysis (G BY E interaction)#######################################

##multisite analysis
Multi.asr <- asreml(NNecScr ~ site+ at(site, 1):Edge,
                  random = ~ diag(site):blok+ corgh(site):clone+
                    at(site,1):units+at(site,2):units+at(site,3):units, 
                  rcov = ~ at(site):ar1(row):ar1(col),
                  data = total, na.method.X='include', maxit=100)
##view variance components
summary(Multi.asr)$varcomp

plot(variogram(Multi.asr))