# program to simulate genotypes that differ in stability and in reaction norms

# model y = mu + G + (1+b)x + e

#genotype 1 - sensitive and variable

mu=6           # overall mean
b=0.5          # slope as deviation from 1
vare=0.5       # residual variance
g=1.0            # deviation from overall mean

n=100            # number of replicates per environment
numenv=5       # number of environments
x1<-array(0,c(n*numenv,2))

env=-3.0       # environmental scale is the mean performance scaled with a mean of 0.0 and a standard deviation of 1
r=0
for (j in 1:numenv) {
env=env+1.0
for (i in 1:n) {
r=r+1
x1[r,1]<-env
x1[r,2]<-mu+g+(1+b)*env+rnorm(1,mean=0,sd=sqrt(vare))

}
}

#genotype 2 - static stable

mu=6
b=-1
vare=0.1
g=-1.0

n=100
numenv=5
x2<-array(0,c(n*numenv,2))

env=-2.95  # set slightly different than for genotype 1 
r=0
for (j in 1:numenv) {
env=env+1.0
for (i in 1:n) {
r=r+1
x2[r,1]<-env
x2[r,2]<-mu+g+(1+b)*env+rnorm(1,mean=0,sd=sqrt(vare))

}
}

# fit linear model
RN1<-lm(x1[,2] ~ x1[,1])
RN2<-lm(x2[,2] ~ x2[,1])


#coefficient of determination
summary(RN1)$r.squared
summary(RN2)$r.squared


png('stability.png',units='in',width=3,height=3,res=600,pointsize=9)
plot(x1[,2] ~ x1[,1],xlab='environment',ylab='yield')
points(x2[,2] ~ x2[,1], col = 2)
abline(RN1)
abline(RN2,col=2)
dev.off()