This script generates data with familly per site interactions. This data is then analysed with breedR and ecovalences of the different familes estimated & compared to each other.

• Simulate data w/ given parameters
• Perform Two-Way analysis
• Store BLUPS and BLUES
• Aggregate results per family
• Calculate ecovalences
• Rank families given ecovalences
• Do graphic outputs

library(breedR)

## Loading required package: sp

# family effects assumed to be independent
num_fam <- 35
mean_eff_fam <- 0
sd_eff_fam <- 0.11
#generate some random normal deviates
fam_eff <- rnorm(num_fam,mean_eff_fam,sd_eff_fam)

# family size and residuals
fam_size <- 20
mean_residuals <- 0
sd_residuals <- 0.55

# site effects are sampled from a normal distribution
# but resulting deviates are kept as fixed
num_sites <- 5
mean_eff_site <- 0
site_eff <- rnorm(num_sites,mean_eff_site,0.5)

# global_mean
global_mean <- 0

# family by site interaction
mean_eff_inter <- 0
sd_eff_inter <- 1-sd_eff_fam-sd_residuals

total_obs <- num_fam*num_sites*fam_size

self <- matrix(data = NA, nrow = total_obs, ncol = 1,
               dimnames = NULL)
fam <- matrix(data = NA, nrow = total_obs, ncol = 1,
              dimnames = NULL)
site <- matrix(data = NA, nrow = total_obs, ncol = 1,
               dimnames = NULL)
inter <- matrix(data = NA, nrow = total_obs, ncol = 1,
                dimnames = NULL)
rep <- matrix(data = NA, nrow = total_obs, ncol = 1,
              dimnames = NULL)
phenot <- matrix(data = NA, nrow = total_obs, ncol = 1,
                 dimnames = NULL)
counter <- 0

# generate data
for (i in 1:num_fam){
  # generate interaction terms for family across sites
  # in the form of some random normal deviates
  inter_eff <- rnorm(num_sites,mean_eff_inter,sd_eff_inter)
  for (j in 1:num_sites){
    # generates residuals for family at given site
    # in the form of some random normal deviates
    res_eff <- rnorm(fam_size,mean_residuals,sd_residuals)
    for (k in 1:fam_size){
      counter <- counter+1
      self[counter] <- counter
      fam[counter] <- i
      site[counter] <- j
      inter[counter] <- paste(i, j, sep = "_")
      rep[counter] <- k
      phenot[counter] <- global_mean+
        site_eff[j]+
        fam_eff[i]+
        inter_eff[j]+
        res_eff[k]
    }
  }
}

data_toy <- data.frame(cbind(self,fam,site,inter,rep,as.numeric(phenot)))
names(data_toy) <- c("self","fam","site","inter","rep","phenot")
write.table(data_toy, file="data_toy.txt", quote = F,
            sep = " ", row.names = F)
rm(data_toy)
data_toy <- read.table(file="data_toy.txt",header = T)

# TWO-WAY anova
two_way <- remlf90(fixed = phenot ~ factor(site),
                   random = ~ factor(fam) + factor(inter),
                   data = data_toy)

## No specification of initial variances.
##      Using default value of 1 for all variance components.
##      See ?breedR.getOption.

summary(two_way)

## Linear Mixed Model with pedigree and spatial effects fit by AI-REMLF90 ver. 1.110 
##    Data: data_toy 
##   AIC  BIC logLik
##  6273 6292  -3134
## 
## Parameters of special components:
## 
## 
## Variance components:
##               Estimated variances     S.E.
## factor(fam)               0.01963 0.011341
## factor(inter)             0.11125 0.015395
## Residual                  0.31310 0.007679
## 
## Fixed effects:
##                     value   s.e.
## factor(site).1  0.0047806 0.0647
## factor(site).2  0.5641400 0.0647
## factor(site).3  0.0339751 0.0647
## factor(site).4 -0.1478764 0.0647
## factor(site).5 -0.0590882 0.0647

# str(two_way)

# adds blups to data.frame
PBV.full <- ranef(two_way)$`factor(fam)`
PBV <- model.matrix(two_way)$`factor(fam)` %*% PBV.full
data_toy$blup_fam <- PBV
PBV.full <- ranef(two_way)$`factor(inter)`
PBV <- model.matrix(two_way)$`factor(inter)` %*% PBV.full
data_toy$blup_inter <- PBV
data_toy$blup_gi <- data_toy$blup_fam+data_toy$blup_inter

# adds blues to data.frame
fixed_eff <- data.frame(fixef(two_way))
data_toy$ef_site <- fixed_eff[data_toy$site,1]
data_toy$fitted <- fitted(two_way)

heading="family by site NoRs"
x_labels <- as.character(c(1:num_sites))
for (i in 1:num_sites){
  x_labels[i] <- paste("site",i,sep = "")
}
colores <- rainbow(n = num_fam)
linetype <- c(1:num_fam) 

aggdata1 <-aggregate(data_toy$site, by=list(data_toy$fam,data_toy$site)
                     ,FUN=mean)
aggdata2 <-aggregate(data_toy$fam, by=list(data_toy$fam,data_toy$site)
                     ,FUN=mean)
aggdata3 <-aggregate(data_toy$fitted, by=list(data_toy$fam,data_toy$site)
                     ,FUN=mean)
aggdata4 <-aggregate(data_toy$blup_inter, by=list(data_toy$fam,data_toy$site)
                     ,FUN=mean)
aggdata <- data.frame(cbind(aggdata1$x,aggdata2$x,aggdata3$x,aggdata4$x))
names(aggdata) <- c("site","fam","fitted","inter")

max_val <- max(aggdata$fitted)
min_val <- min(aggdata$fitted)

# calculates ecovalences from interaction blups
aggdata$sq_inter <- aggdata$inter^2
aggdata_inter <- data.frame(aggregate(aggdata$sq_inter,
                                      by=list(aggdata$fam),FUN=sum))
names(aggdata_inter) <- c("fam","sq_inter")
total_sq_inter <- sum(aggdata_inter$sq_inter)
aggdata_inter$ecoval <- aggdata_inter$sq_inter/total_sq_inter

# rank in ascending order families by their ecovalence
rank <- aggdata_inter$fam[order(aggdata_inter$ecoval)]
upper_rank <- 4
upper_threshold <- aggdata_inter$ecoval[rank[upper_rank]]
lower_rank <- num_fam-upper_rank-1
lower_threshold <- aggdata_inter$ecoval[rank[lower_rank]]