###############################################################
#######   Genetica Quantitativa  - Equacao de Price    ########
####### Selecao Direcional - Disruptiva/Estabilizadora ########
###############################################################
#getwd() # informa que pasta que vc esta
#setwd("/Users/flaviamarquitti/Documents/NE441-F017_2019/Codigos") # altera a sua pasta

# Leitura de 3 dados de caracteristica (Phi) e fitness/aptidao (W)
d1=read.table("dados1.txt", header = TRUE)
d2=read.table("dados2.txt", header = TRUE)
d3=read.table("dados3.txt", header = TRUE)

### Neste caso, note que os 3 dos graficos de W x Phi sao positivamente relacionados,
### portanto temos que a media da caracteristica aumenta  com o tempo (direcionalidade). 
### Quanto ao tipo de selecao, apresentamos 3 cenario: disruptiva, estabiliazdora e neutro.
### Para isso, veja os graficos de W x (Phi-Phi_m)^2

par(mfcol=c(2,3))

### Dado 1:
plot(d1$W~d1$Phi, xlab=expression(paste(phi)), ylab="W", pch=16, col="lightblue")
model1=lm(d1$W~d1$Phi) # modelo de regressao linear de W com Pi
model1 # veja sinal do coeficiente angular
abline(model1, col="lightblue", lwd=2) # acrescentando a linha de regressao
abline(v=mean(d1$Phi), col="lightgray", lty=2) # acrescentando o Phi medio

d1.phi_phiM2=(d1$Phi-mean(d1$Phi))^2 # calculo de (Phi-Phi_m)^2 pelos dados de Phi
plot(d1$W~d1.phi_phiM2, xlab=expression(paste((phi-bar(phi))^2)), ylab="W", pch=16, col="lightblue")
model1.1=lm(d1$W~d1.phi_phiM2) # modelo de regressao linear de W com (Phi-Phi_m)^2
model1.1 # veja sinal do coeficiente angular
abline(model1.1, col="lightblue", lwd=2) # acrescentando a linha de regressao

### Dado 2:
plot(d2$W~d2$Phi, xlab=expression(paste(phi)), ylab="W", pch=16, col="lightcoral")
model2=lm(d2$W~d2$Phi)
model2
abline(model2, col="lightcoral", lwd=2)
abline(v=mean(d2$Phi), col="lightgray", lty=2)

d2.phi_phiM2=(d2$Phi-mean(d2$Phi))^2 
plot(d2$W~d2.phi_phiM2, xlab=expression(paste((phi-bar(phi))^2)), ylab="W", pch=16, col="lightcoral")
model2.2=lm(d2$W~d2.phi_phiM2)
model2.2
abline(model2.2, col="lightcoral", lwd=2)

### Dado 3:
plot(d3$W~d3$Phi, xlab=expression(paste(phi)), ylab="W", pch=16, col="sandybrown")
model3=lm(d3$W~d3$Phi)
model3
abline(model3, col="sandybrown", lwd=2)
abline(v=mean(d3$Phi), col="lightgray", lty=2)

d3.phi_phiM2=(d3$Phi-mean(d3$Phi))^2 
plot(d3$W~d3.phi_phiM2, xlab=expression(paste((phi-bar(phi))^2)), ylab="W", pch=16, col="sandybrown")
model3.3=lm(d3$W~d3.phi_phiM2)
model3.3
abline(model3.3, col="sandybrown", lwd=2)

#### So para mostrar que a caracteristica tbm pode diminuir e ter os diferentes efeitos quanto
#### a selecao disruptiva/estabilizadora

d4=read.table("dados4.txt", header = TRUE)
d5=read.table("dados5.txt", header = TRUE)

par(mfcol=c(2,2))

### Dado 4:
plot(d4$W~d4$Phi, xlab=expression(paste(phi)), ylab="W", pch=16, col="plum")
model4=lm(d4$W~d4$Phi)
model4
abline(model4, col="plum", lwd=2)
abline(v=mean(d4$Phi), col="lightgray", lty=2)

d4.phi_phiM2=(d4$Phi-mean(d4$Phi))^2 
plot(d4$W~d4.phi_phiM2, xlab=expression(paste((phi-bar(phi))^2)), ylab="W", pch=16, col="plum")
model4.4=lm(d4$W~d4.phi_phiM2)
model4.4
abline(model4.4, col="plum", lwd=2)

### Dado 5:
plot(d5$W~d5$Phi, xlab=expression(paste(phi)), ylab="W", pch=16, col="palegreen2")
model5=lm(d5$W~d5$Phi)
model5
abline(model5, col="palegreen2", lwd=2)
abline(v=mean(d5$Phi), col="lightgray", lty=2)

d5.phi_phiM2=(d5$Phi-mean(d5$Phi))^2 
plot(d5$W~d5.phi_phiM2, xlab=expression(paste((phi-bar(phi))^2)), ylab="W", pch=16, col="palegreen2")
model5.5=lm(d5$W~d5.phi_phiM2)
model5.5
abline(model5.5, col="palegreen2", lwd=2)