# Install package mvtnorm if not already installed
# Then...
require(mvtnorm)

# -----------
# Parameters
# -----------

# Consider a phenotype X measured in a large sample of twin pairs. We'll assume for this exercise that:
# (i)   mean(X_tw1) = 0, var(X_tw1) = 1
# (ii)  mean(X_tw2) = 0, var(X_tw2) = 1
# (iii) cov(X_tw1,X_tw2) = 0.5

pop_means <- c(0,0)	  				# mean(X_tw1)=0; mean(X_tw2)=0
pop_varcov <- matrix(c(1,0.5,0.5,1),ncol=2)	# var(X_tw1)=0; cov(X_tw1,X_tw2)=0.5;
                                                # cov(X_tw2,X_tw1)=0.5; var(X_tw2)


# -------------------------------------------------------------------------------------
# Exercise 1: estimate the likelihood of observing in the population a twin pair with 
#              traits 0.5 and -0.3 under the above assumptions and using standard R functions
# -------------------------------------------------------------------------------------

X_tw1 <- 0.5
X_tw2 <- -0.3

# The likelihood of an individual bivariate observation (ie. a twin pair with values X=[tw1 tw2]) is:
# L = inv( sqrt( det(2 * pi * pop_varcov) ) ) * e^( -0.5 * (X-pop_means) %*% inv(pop_varcov) %*% t(X-pop_means) ) 

X <- matrix(c(X_tw1,X_tw2),ncol=2)	# Twin pair values - in matrix format
X_Likelihood <- solve( sqrt( det(2 * pi * pop_varcov) ) ) * exp( -0.5 * (X-pop_means) %*% solve(pop_varcov) %*% t(X-pop_means) )

# Print likelihood and -2LogL to screen:
c( X_Likelihood, -2*log(X_Likelihood) )


# -------------------------------------------------------------------------------------
# Exercise 2: estimate the likelihood of observing in the population a twin pair with 
#              traits 0.5 and -0.3 under the above assumptions and using openMx functions
# -------------------------------------------------------------------------------------
# The formula for the bivariate likelihood is the same as above:
# L = inv( sqrt( det(2 * pi * pop_varcov) ) ) * e^( -0.5 * (X-pop_means) %*% inv(pop_varcov) %*% t(X-pop_means) ) 

X <- c(X_tw1,X_tw2)				# Twin pair values - in vector format
X_Likelihood <- mxModel("X_L",
        mxMatrix( type="Full", nrow=1, ncol=2, values=X, name="X" ), 
        mxMatrix( type="Full", nrow=1, ncol=2, values=pop_means, name="M" ), 
        mxMatrix( type="Full", nrow=2, ncol=2, values=as.numeric(pop_varcov), name="S" ),
        # Compute likelihood for data vector X in two steps (c1,c2) to make it easier:
        # L = inv( 2*pi*sqrt( det(S) ) ) * e^( -0.5 * (X-M) %*% inv(S) %*% t(X-M) ) 
        #     <------------------------>   <-------------------------------------->
        #                  c1                                    c2
        mxAlgebra( expression= solve( 2 * pi * sqrt( det(S) ) ), name="c1" ),
        mxAlgebra( expression= exp( -0.5 * (X-M) %*% solve(S) %*% t(X-M) ), name="c2" ),
        mxAlgebra( expression= c1 * c2, name="X_Likelihood" ),
        mxAlgebraObjective("X_Likelihood")
)

# Run the mxModel that will compute the likelihood of twin pair X
output <- mxRun(X_Likelihood)

# Print likelihood to screen:
output$X_Likelihood


# -------------------------------------------------------------------------------------
# Exercise 3: estimate the likelihood of observing in the population a twin pair with 
#             traits 0.5 and -0.3 under the above assumptions using the dmvnorm() function 
#             from the mvnorm package
# -------------------------------------------------------------------------------------

X <- c(X_tw1,X_tw2)				# Twin pair values - in vector format
X_Likelihood <- dmvnorm(X, mean=pop_means,sigma=pop_varcov)

# Print likelihood to screen:
X_Likelihood



# ----------------------
# Exercise for the class
# ----------------------

# Divide class into 10 groups (~5 pairs of students for each, ie 1 half-row)
# Each group to compute to estimate Likelihood for a twin pair considering:
# twin pair values of [0.5 -0.5] or [0.5 0.5] AND
# cov(tw1,tw2) of -0.5 -0.25 0 0.25 or 0.5

L=rep(NA,t=10)
groups=matrix(c(1:10,rep(0.5,t=10),rep(c(-0.5,0.5),each=5),
       seq(-0.5,0.5,by=0.25),seq(-0.5,0.5,by=0.25),L),ncol=5)
colnames(groups)=c("group","X_tw1","X_tw2","cov(tw1,tw2)","Likelihood")