# -------------
# Plot results
# -------------

# Vector with likelihood results
likelihood <- c(0.1555633, 0.1345785, 0.1239500, 0.1177795, 0.1114660,
                0.1114660, 0.1177795, 0.1239500, 0.1345785, 0.1555633)

# Matrix to store all results, 5 columns: GROUP X_twin1 X_twin2 cov(X_tw1,X_tw2) Likelihood
group <- 1:10
X_tw1 <- rep(0.5,t=10)
X_tw2 <- rep(c(-0.5,0.5),each=5)
covTw1Tw2 <- c(seq(-0.5,0.5,by=0.25),seq(-0.5,0.5,by=0.25))
results <- print(matrix(c(group,X_tw1,X_tw2,covTw1Tw2,likelihood),ncol=5))

# Plot Likelihood results for twin pair X=[0.5 -0.5]
plot(results[1:5,4], results[1:5,5],col='red',ylim=c(0.10,0.16),type='b',
     ylab="Twin-pair Likelihood",xlab="cov(Tw1,Tw2)")

# Add to the same plot Likelihood results for twin pair X=[0.5 0.5]
lines(results[1:5,4], results[6:10,5],type='b',col='blue')

# Add legend
text(-0.25,0.16,pos=1,labels="X=[0.5 -0.5]",col='red')
text(-0.25,0.115,pos=1,labels="X=[0.5 0.5]",col='blue')

stop()

# If you want to save to a JPEG 
jpeg("likelihood.jpeg",width=400,height=400,quality=100)
plot(results[1:5,4], results[1:5,5],col='red',ylim=c(0.10,0.16),type='b',
     ylab="Twin-pair Likelihood",xlab="cov(Tw1,Tw2)")
lines(results[1:5,4], results[6:10,5],type='b',col='blue')
text(-0.25,0.16,pos=1,labels="X=[0.5 -0.5]",col='red')
text(-0.25,0.115,pos=1,labels="X=[0.5 0.5]",col='blue')
dev.off()