! SCRIPT NAME : lingrow.mx (jh)
! GOAL : To evaluate bivariate linear growth curve model for heritability estimation of latent variables: intercept and slope
! DATA : continuous
! INPUT : raw data
! UNI/BI/MULTI : uni (longitudinal observations)
! DATA-GROUPS : MZ, DZ-SS
! MEANS MODEL : grand mean, sex effect
! VARIANCE COVARIANCE MODEL(S) : 1.ACE , 2.ACE orthogonal, 3. AE,  4. ADE

! Downloading Mx software: <http://www.vcu.edu/mx>
! Mx script's library: <http://www.psy.vu.nl/mxbib>

! GenomEUtwin lingrow.mx
!
! Script for the analysis of (simulated) univariate
! continuous longitudinal twin data in wave representation.
! 
! Model: Linear growthcurve model. 
! Biometric model (ACE) for slopes and intercepts.

#define NT 4	  ! number of timepoints.
#define NTX2 8	  ! twice the number og timepoints.
#define NF 2      ! number of factors.
#define NFX2 4    ! twice the number of factors.
#define NINPUT 11 ! number of input variables (4 variables for each twin, plus family id, sex, zygosity).


GROUP #1 Initialization
CALCULATION NG=3
BEGIN MATRICES;
	F FULL   NT 2 FIXED ! matrix of factor loadings
        X LOWER   2 2 FREE  ! cholesky factor for genetic (co)variance
        Y LOWER   2 2 FREE  ! cholesky factor for variance due to shared covariance
        Z LOWER   2 2 FREE  ! cholesky factor for variance due to non-shared environment
        V LOWER   2 2 FIXED ! cholesky factor for genetic dominant (co)variance
 	L FULL    1 1  FREE ! cholesky factor for residual variance.
	B Iden  NT NT FIXED ! identity matrix
	J Iden  NF NF FIXED ! identity matrix
	W Zero  NT NT FIXED ! zero matrix
        H FULL    1 1 FIXED ! scalar = 1/2
        I FULL    1 1 FIXED ! scalar = 1/4
	Q FULL    2 1 FREE  ! mean intercept & slope 
	N FULL  1  NT FIXED ! Effect of covariate
END MATRICES;
MATRIX H .5
MATRIX I .25
MATRIX F
  1 0
  1 1
  1 2
  1 3


START -0.5 all
st 0.2766 X(1,1)
st -0.0033 X(2,1)
st -0.0014 X(2,2)

st 0.27152 Y(1,1)
st -0.0032231 Y(2,1)
st -0.000000007 Y(2,2)

st 0.4301 Z(1,1)
st -0.0052 Z(2,1)
st 0.0011 Z(2,2)

START 3.4923 Q(1,1) 
START -0.0028 Q(2,1)

START 0.0582  L(1,1) 

BOUND -100 100 ALL


!START  0.05 X(1,1) X(2,2) Y(1,1) Y(2,2) Z(1,1) Z(2,2)
!START  0.05 X(2,1) Y(2,1) Z(2,1)


BEGIN ALGEBRA;
	A = X*X';	
	C = Y*Y';
	E = Z*Z';
	D = V*V';
        K =((J.(A + D + C + E))~*A)|((J.(A + D + C + E))~*C)|((J.(A + D + C + E))~*E); ! scaled covariance matrix 
	M = A+D+C+E |  A+D+C _
	     A+D+C  | A+D+C+E;  	! covariance matrix of intercept and slope for mz pair
	S = A+D+C+E | H@A+I@D+C _
	    H@A+I@D+C | A+D+C+E; 	! covariance matrix of intercept and slope for dz pair
	G = (L*L')@B | W _
	     W | (L*L')@B; 		! matrix of residual variance    
	P = \stnd(A) | \stnd(E);
END ALGEBRA;
!INTERVAL
! K(1,1,1) K(1,2,2) K(1,1,3) K(1,2,4) K(1,1,5) K(1,2,6)
!INTERVAL
! Q(1,1,1) Q(1,2,1)
!INTERVAL
! N(1,1,1) N(1,1,2) N(1,1,3) N(1,1,4)
!INTERVAL
!  P(1,2,1) P(1,2,3)
END

GROUP #2 Mz twins
DATA NI=NINPUT                          ! number of inputvariables
REctangular file=example4.dat      	! input file containing raw data
LABELS
  pairid sex wave11 wave12 wave13 wave14 wave21 wave22 wave23 wave24 zygos
Select if  zygos = 1 ;                 ! mz's only for this group
Select sex wave11 wave12 wave13 wave14 wave21 wave22 wave23 wave24 ;
Definition sex;

BEGIN MATRICES;
	F FULL   NT NF 	=F1
	M COMPUTED 	=M1
	G COMPUTED 	=G1
	Q FULL    2 1 	=Q1 
	N FULL  1  NT 	=N1 		! Effect of covariate, female
	K FULL  1  1
	I IDEN  2  2
        U UNIT  1  2
END MATRICES;

Specify K
-1 


MEANS (F*Q)' + K*N | (F*Q)'+ K*N ;
COVARIANCE (I@F)&M + G ;
OPTIONS RS
END

GROUP #3 Dz twins
DATA NI=NINPUT
REctangular file=example4.dat
LABELS
  pairid sex wave11 wave12 wave13 wave14 wave21 wave22 wave23 wave24 zygos
Select if  zygos = 2 ;               ! dz's only for this group
Select sex wave11 wave12 wave13 wave14 wave21 wave22 wave23 wave24 ;
Definition sex;

BEGIN MATRICES;
	F FULL   NT 2 	=F1
	D COMPUTED 	=S1
	Q FULL    2 1 	=Q1
	G COMPUTED 	=G1
	N FULL  1  NT 	=N1 		! Effect of covariate, female
	K FULL  1  1
        H FULL  1  1 	=H1
	I IDEN  2  2
        U UNIT  1  2
END MATRICES;

Specify K
-1 


MEANS (F*Q)' + K*N | (F*Q)'+ K*N ;
COVARIANCE (I@F)&D + G ;
!Option  TH=-2
Option jiggle
!Option it=5000
OPTIONS MULTIPLE issat
END

SAVE lingrow.mxs

! FIT ORTHOGONAL 
GET lingrow.mxs
DROP @0 X(1,2,1) Y(1,2,1) Z(1,2,1) ! no covariance of intercept and slope.
END

!Significance of no heritability of intercept  (Under ACE model)
GET lingrow.mxs
DROP @0 X(1,1,1) X(1,2,1)
END

! Significance of no heritability of slope  (Under ACE model)
GET lingrow.mxs
DROP @0 X(1,2,1) X(1,2,2)
END

! Significance of no common environmental effect of intercept 
GET lingrow.mxs
DROP @0 Y(1,1,1) Y(1,2,1)
END

! Significance of no common environmental effect of slope 
GET lingrow.mxs
DROP @0 Y(1,2,1) Y(1,2,2)
END

! FIT AE MODEL
GET lingrow.mxs
DROP @0 Y(1,1,1) Y(1,2,1) Y(1,2,2)
END