! Sensitivity to underreporting model ! drug abuse and dependence in male-male twins ! using self report and cotwin report #NGroups 6 ! number of groups in mx script #define $mzfreq 506 2 51 6 3 0 2 0 42 2 51 18 5 0 15 4 ! input mz freqs #define $dzfreq 292 6 56 9 3 0 1 1 56 4 26 10 13 0 7 4 ! input dz freqs ! -------------- Group ID 1---------------- G1: MZ Calc Calculation #define NumManifest 4 Begin Matrices; S Symm 16 16 A Full 16 16 F Full 4 16 I Iden 16 16 End Matrices; Matrix S 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 1 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0 0 1 Specify S 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 Label Row S T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Label Col S T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Matrix A 0 0 0 0 0 0 0 0 0 0 0.5 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.6 0.6 0.6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.6 0.6 0.6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Specify A 0 0 0 0 0 0 0 0 0 0 5 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 3 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Label Row A T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Label Col A T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Matrix F 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 Label Row F T1_T1 T2_T1 T2_T2 T1_T2 Label Col F T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Bound -0.99 0.99 1 ! bound correlation estimate between -1 and +1 Bound 0.001 0.999 2 3 4 5 6 7 8 ! bound parameter estimates Begin Algebra; V = F & ((I - A)~ & S); ! the covariance matrix of the reported variables End Algebra; Options RSidual End Group; ! -------------- Group ID 2---------------- G2: DZ Calc Calculation #define NumManifest 4 Begin Matrices; S Symm 16 16 A Full 16 16 F Full 4 16 I Iden 16 16 End Matrices; Matrix S 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 .5 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 1 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0 0 1 Specify S 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 Label Row S T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Label Col S T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Matrix A 0 0 0 0 0 0 0 0 0 0 0.5 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.6 0.6 0.6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.6 0.6 0.6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Specify A 0 0 0 0 0 0 0 0 0 0 5 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 3 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Label Row A T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Label Col A T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Matrix F ! used to peel off the cov matrix for the 4 reported variables from the larger cov matrix 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 Label Row F T1_T1 T2_T1 T2_T2 T1_T2 Label Col F T1_T1 T2_T1 T2_T2 T1_T2 A C E E C A LP1 LP2 RE1 RE2 RE1 RE2 Begin Algebra; V = F & ((I - A)~ & S); ! the covariance matrix of the reported variables End Algebra; Bound -0.99 0.99 1 ! bound correlation estimate between -1 and +1 Bound 0.001 0.999 2 3 4 5 6 7 8 ! bound parameter estimates Options RSidual End Group; G3: CONSTRAIN VARIANCES CONSTRAINT NI=1 BEGIN MATRICES = GROUP 1; K ZI 4 14 M Zero 4 2 H FULL 1 4 J IZ 2 16 G UNIT 1 3 END MATRICES; MATRIX H 1 0 0 0 BEGIN ALGEBRA; L = (K|M); W = (I - A)~ & S; ! the entire covariance matrix X = H & (L & W); ! Peel off the variance for latent phenotype for T1 from the covariance matrix Y = J & W; ! Peel off the variance for T1_on_T1 and T1_on_T2 from the covariance matrix B = \D2V (X) | \D2V (Y); END ALGEBRA; CONSTRAINT B = G; ! Set the variances for T1 latent phenotype, T1 on T1, and T1 on T2 equal to 1 OPTIONS RS END GROUP; G4: MZ Group Data NInput=1 Begin Matrices; C comp =V1 ! The correlation matrix for MZs I Iden 2 2 J Zero 2 2 K Zero 1 4 ! For means of normal (zeroes) B Unit 1 4 ! number of thresholds in each dimension (one) N Full 1 1 ! 2.0 for -2lnL O Full 16 1 ! observed data P Full 1 1 ! underreporting rate of Twin1 when Twin2 reports that Twin1 is negative for DAD Q Full 1 1 ! underreporting rate of Twin1 when Twin2 reports that Twin1 is positive for DAD R Full 1 1 ! underreporting rate of Twin2 when Twin1 reports that Twin2 is negative for DAD S Full 1 1 ! underreporting rate of Twin2 when Twin1 reports that Twin2 is positive for DAD T Full 1 4 Free ! Thresholds U Unit 1 1 Z Zero 1 1 End Matrices; Equate T 1 1 T 1 3 ! Because Tw1 and Tw2's self-reports should be equal Equate T 1 2 T 1 4 ! Because Tw1 and Tw2's co-twin reports should be equal Bound -3 3 T 1 1 - T 1 2 ! Put reasonable bounds on thresholds Matrix O $mzfreq Matrix N 2 ! P, Q, R, and S are the underreporting parameters Matrix P 0.1 ! P should equal R because parameters for T1 and T2 should be equal Matrix Q 0.0 ! Q should equal S because parameters for T1 and T2 should be equal Matrix R 0.1 Matrix S 0.0 Begin Algebra; X = ((Z|Z_Z|Z)|(Z|Z_Z|Z))_((Z|Z_Z|Z)|(Z|Z_Z|Z)); ! The transition matrix, Y, follows Y = (((U|Z_Z|U)|(Z|Z_Z|Z))_((((R.U)|Z)_(Z|(S.U)))|((((U-R).U)|Z)_(Z|((U-S).U)))))|X|X|X_ X|(((U|Z_Z|U)|(Z|Z_Z|Z))_((((R.U)|Z)_(Z|(S.U)))|((((U-R).U)|Z)_(Z|((U-S).U)))))|X|X_ (((P|Z_Z|P)|(Z|Z_Z|Z))_((((P.R)|Z)_(Z|(P.S)))|(((P.(U-R))|Z)_(Z|(P.(U-S))))))|X| ((((U-P)|Z_Z|(U-P))|(Z|Z_Z|Z))_((((R.(U-P))|Z)_(Z|(S.(U-P))))|((((U-P).(U-R))|Z)_(Z|((U-P).(U-S))))))|X_ X|(((Q|Z_Z|Q)|(Z|Z_Z|Z))_((((Q.R)|Z)_(Z|(Q.S)))|(((Q.(U-R))|Z)_(Z|(Q.(U-S))))))|X| ((((U-Q)|Z_Z|(U-Q))|(Z|Z_Z|Z))_((((R.(U-Q))|Z)_(Z|(S.(U-Q))))|((((U-Q).(U-R))|Z)_(Z|((U-Q).(U-S)))))); D = \allint(C_K_B_T) ; ! estimates the true cell proportions based on the correlations and thresholds E = D*Y; ! Computes the expected reported proportions in the cells F = \sum(E); ! just to check that this equals 1, else error G = \sum(O)@E; ! Given E, computes the expected reported frequencies L = \sum(O)@D; ! Given D, computes the expected true frequencies M = \sum(L); ! check that this sums to 1 H = O|L'|O-L'|G'|O-G'; ! Compare observed, expected true, and expected reported frequencies End Algebra; Compute -N.\sum(\ln(E).O'); ! Compute -2LL Option User-Defined End Group G5: DZ Group Data NInput=1 Begin Matrices; C comp =v2 ! The correlation matrix for DZs I Iden 2 2 J Zero 2 2 K Zero 1 4 ! For means of normal (zeroes) B Unit 1 4 ! number of thresholds in each dimension (one) N Full 1 1 ! 2.0 for -2lnL O Full 16 1 ! observed data ! P, Q, R, and S are the underreporting parameters P Full 1 1 =P4 ! underreporting rate of Twin1 when Twin2 reports that Twin1 is negative for DAD Q Full 1 1 =Q4 ! underreporting rate of Twin1 when Twin2 reports that Twin1 is positive for DAD R Full 1 1 =R4 ! underreporting rate of Twin2 when Twin1 reports that Twin2 is negative for DAD S Full 1 1 =S4 ! underreporting rate of Twin2 when Twin1 reports that Twin2 is positive for DAD T Full 1 4 Free ! Thresholds U Unit 1 1 Z Zero 1 1 End Matrices; Equate T 1 1 T 1 3 ! Because Tw1 and Tw2's self-reports should be equal Equate T 1 2 T 1 4 ! Because Tw1 and Tw2's co-twin reports should be equal Bound -3 3 T 1 1 - T 1 2 ! Put reasonable bounds on thresholds Matrix O $dzfreq Matrix N 2 Begin Algebra; X = ((Z|Z_Z|Z)|(Z|Z_Z|Z))_((Z|Z_Z|Z)|(Z|Z_Z|Z)); ! The transition matrix, Y, follows Y = (((U|Z_Z|U)|(Z|Z_Z|Z))_((((R.U)|Z)_(Z|(S.U)))|((((U-R).U)|Z)_(Z|((U-S).U)))))|X|X|X_ X|(((U|Z_Z|U)|(Z|Z_Z|Z))_((((R.U)|Z)_(Z|(S.U)))|((((U-R).U)|Z)_(Z|((U-S).U)))))|X|X_ (((P|Z_Z|P)|(Z|Z_Z|Z))_((((P.R)|Z)_(Z|(P.S)))|(((P.(U-R))|Z)_(Z|(P.(U-S))))))|X| ((((U-P)|Z_Z|(U-P))|(Z|Z_Z|Z))_((((R.(U-P))|Z)_(Z|(S.(U-P))))|((((U-P).(U-R))|Z)_(Z|((U-P).(U-S))))))|X_ X|(((Q|Z_Z|Q)|(Z|Z_Z|Z))_((((Q.R)|Z)_(Z|(Q.S)))|(((Q.(U-R))|Z)_(Z|(Q.(U-S))))))|X| ((((U-Q)|Z_Z|(U-Q))|(Z|Z_Z|Z))_((((R.(U-Q))|Z)_(Z|(S.(U-Q))))|((((U-Q).(U-R))|Z)_(Z|((U-Q).(U-S)))))); D = \allint(C_K_B_T) ; ! compute most likely true distribution E = D*Y; ! compute most likely reported distribution given underrept params F = \sum(E); ! just to check that this equals 1, else error G = \sum(O)@E; ! expected reported freqs L = \sum(O)@D; ! expected true freqs M = \sum(L); H = O|L'|O-L'|G'|O-G'; ! Compare observed, expected true, and expected reported frequencies End Algebra; Compute -N.\sum(\ln(E).O'); ! Compute -2LL Option User-Defined Option Multiple Issat jiggle Option func=1.e-11 Option feas=1.e-3 Option df=30 ! DF=30 because there are 16 freqs for MZs and DZs each, minus 2 for each group Option th=-20 ! sometimes the job is slow to converge End Group G6: Calc group calculation begin Matrices; v full 1 3 x full 1 1 = %f4 ! pull out the fit from group 4 y full 1 1 = %f5 ! pull out the fit from group 5 end matrices; specify v 2 3 4 ! put estimates for a, c and e into matrix v bound 0.001 0.999 v 1 1 - v 1 3 begin algebra; a = x+y; b = v.v; ! Square a, c, and e to calculate A, C, and E end algebra; Interval b 1 1 - b 1 3 ! calculate CIs on A, C and E End