! Script for estimating polychoric correlations : correlation between two 3-category variables
! Variable: Age of smoking onset.  
! Categories = non-smoker | smoker after 18 years | smoker before/at 18 years
! The age at onset variable has 3 smoking categories, 2 thresholds per variable, 2 for twin1 and 2 for twin2
! Equality in Thresholds are tested 
! THRESHOLD MODEL 2x2 : (thresholds twin1, thresholds twin2)

! Exercise: run, interpret results and then fill-in table
#define nvarx2   2      	! number of variables x number of twins
#define nthresh 2 		! number of thresholds
#ngroups 2

G1: Data and model for MZM 
DAta NInput_vars=3  
Missing=.
Ordinal File=smk_prac.ord
Labels  
  zyg ageon_t1 ageon_t2
SELECT IF 	zyg = 2
SELECT 		ageon_t1 ageon_t2  /

Begin Matrices;  		
R STAN nvarx2 nvarx2 FREE		! Correlation matrix
T FULL nthresh nvarx2 FREE		! thresh tw1, thresh tw2
L Lower nthresh nthresh
End matrices;

Value 1 L 1 1 to L nthresh nthresh 	! initialize L

COV 		R / 		! Predicted Correlation matrix for MZM pairs

Thresholds 	L*T /		! threshold model to ensure t1>t2>t3 etc.......

Start .1 T 1 1 T 1 2 			! starting values for the thresholds
start .2 T 2 1 T 2 2
Bound 0.01 1 T 1 1 T 1 2
Bound 0.1 5 T 2 1 T 2 2 
Start .7 R 2 1 	 		! starting value for the correlation
Option RS
END
 
G2: Data and model for DZM 
DAta NInput_vars=3
Missing=.
Ordinal File=smk_prac.ord
Labels  
  zyg ageon_t1 ageon_t2
SELECT IF 	zyg = 4
SELECT 		ageon_t1 ageon_t2  /

Begin Matrices;  		
R STAN nvarx2 nvarx2 FREE		! Correlation matrix
T FULL nthresh nvarx2 FREE		! thresh tw1, thresh tw2
L Lower nthresh nthresh =L1
End matrices;

COV 		R / 		! Predicted Correlation matrix for MZM pairs

Thresholds 	L*T /		! threshold model to ensure t1>t2>t3 etc.......

Start .05 T 1 1 T 1 2 			! starting values for the thresholds
start .2 T 2 1 T 2 2
Bound 0.01 1 T 1 1 T 1 2
Bound 0.1 5 T 2 1 T 2 2 
Start .3 R 2 1 	 		! starting value for the correlation
Option func=1.E-10  		! function precision is less than usual
Option RS issat multiple
End

Save cor.mxs        		!this saves the full model so you can test reduced models

! Test birthorder
EQ T 1 1 1 T 1 1 2 		!constrain TH1 to be equal across Tw1 and Tw2 MZM 
EQ T 1 2 1 T 1 2 2 		!constrain TH2 to be equal across Tw1 and Tw2 MZM 
EQ T 2 1 1 T 2 1 2 		!constrain TH1 to be equal across Tw1 and Tw2 DZM 
EQ T 2 2 1 T 2 2 2 		!constrain TH2 to be equal across Tw1 and Tw2 DZM 
End

Get cor.mxs

! Test equality of thresholds between MZM & DZM
EQ T 1 1 1 T 1 1 2 T 2 1 1 T 2 1 2 	!constrain TH1 to be equal across all Males
EQ T 1 2 1 T 1 2 2 T 2 2 1 T 2 2 2 	!constrain TH2 to be equal across all Males
End