#ngroups 10
!
! Ascertainment correction, estimating pi
!
G1: Model parameters
 Calc
 Begin Matrices;
  H Full 1 1
  Q Full 1 1
  U unit 1 1
!  M Stan 2 2 Free
!  L Stan 2 2 Free
  X Lower 1 1 Free
  Y Lower 1 1 Free
  Z Lower 1 1 Free
  W Lower 1 1 
  O Zero 1 1
 End Matrices;
 Begin Algebra;
  A= X*X'; 
  C= Y*Y';
  E= Z*Z';
  D= W*W';
  L=  A+C+D+E   | H@A+C+Q@D _
      H@A+C+Q@D | A+C+D+E ;  ! DZ cov matrix
  M=  A+C+D+E | A+C+D _ 
      A+C+D   | A+C+D+E ;    ! MZ cov matrix
 End Algebra;
 Matrix H .5
 Matrix Q .25
 Start .6 all
 Bound 0 1 x 1 1 y 1 1
 Bound 0.1 1 z 1 1
 Intervals A 1 1 1 C 1 1 1 E 1 1 1 
End

#define repeat_number 8
G2: MZ twin pairs
 Data Ninput=9
 Labels studyid studyyr zygosity twin1 twin2 asc freq pihat threshold
 Ordinal File=scz_mx.dat
 Select if studyid = $repeat_number
 Select if zygosity = 1;
 Select twin1 twin2 asc freq pihat threshold ;
 Definition Asc Freq Pihat Threshold ;
 Begin Matrices= Group 1;
  F full 1 1 ! Frequency
  G full 1 4 ! to select ascertainment prob
  P full 1 1 ! ascertainment probability p
  T full 1 1 ! threshold
 End Matrices;
 Begin Algebra;
  S= (P) _ (P) _ P+P-P*P ; ! Vector of ascertainment probabilities
  J= (P+P)*\mnor(U_O_T_T_U)-P*P*\mnor(M_O|O_T|T_T|T_U|U);
 End Algebra:

 Specify F Freq
 Specify P Pihat
 Specify T threshold
 Matrix G 1 1 1 1
 Specify G asc 0 asc 0

 Covariances M;
 Thresholds T|T ;
 Frequency F ;
#if repeat_number < 4
! pihat is not 1 so weight is needed
 Weight \part(S,G)%J ;
#elseif repeat_number = 9
! pihat is not 1 so weight is needed
 Weight \part(S,G)%J ;
#elseif repeat_number = 7
! population sample, so no weight required
! Weight \part(S,G)%J ;
#elseif repeat_number = 10
! population sample, so no weight required
! Weight \part(S,G)%J ;
#else
! pihat is 1 but -- not sampled
 Weight \part(S,G)%J ;
#endif
 Options RSidual
End

G3: DZ twin pairs
 Data Ninput=9
 Labels studyid studyyr zygosity twin1 twin2 asc freq pihat threshold
 Ordinal File=scz_mx.dat
 Select if studyid = $repeat_number
 Select if zygosity = 2
 Select twin1 twin2 asc freq pihat threshold ;
 Definition Asc Freq Pihat Threshold ;
 Begin Matrices= Group 1 ;
  F full 1 1     ! Frequency
  G full 1 4 =G2 ! To select ascertainment probability
  P full 1 1 ! ascertainment probability p
  T full 1 1 ! threshold
 End Matrices;
 Specify F Freq
 Matrix G 1 1 1 1
 Specify G asc 0 asc 0
 Specify P Pihat
 Specify T threshold
 Begin Algebra;
  S= (P) _ (P) _ P+P-P*P ; ! Vector of ascertainment probabilities
  K= (P+P)*\mnor(U_O_T_T_U)-P*P*\mnor(L_O|O_T|T_T|T_U|U);
 End Algebra:
 Covariances L;
 Thresholds T|T ;
 Frequency F ;
#if repeat_number = 9
! pihat is not 1 so weight is needed
 Weight \part(S,G)%K ;
#elseif repeat_number = 10
! population sample, so no weight required (could estimate threshold)
! Weight \part(S,G)%K ;
#elseif repeat_number = 11
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#elseif repeat_number = 12
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#elseif repeat_number = 8
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#endif
 Options RSidual NDecimals=4 
 Option jiggle func=1.e-12
End


#define repeat_number 9
G2: MZ twin pairs
 Data Ninput=9
 Labels studyid studyyr zygosity twin1 twin2 asc freq pihat threshold
 Ordinal File=scz_mx.dat
 Select if studyid = $repeat_number
 Select if zygosity = 1;
 Select twin1 twin2 asc freq pihat threshold ;
 Definition Asc Freq Pihat Threshold ;
 Begin Matrices= Group 1;
  F full 1 1 ! Frequency
  G full 1 4 ! to select ascertainment prob
  P full 1 1 ! ascertainment probability p
  T full 1 1 ! threshold
 End Matrices;
 Begin Algebra;
  S= (P) _ (P) _ P+P-P*P ; ! Vector of ascertainment probabilities
  J= (P+P)*\mnor(U_O_T_T_U)-P*P*\mnor(M_O|O_T|T_T|T_U|U);
 End Algebra:

 Specify F Freq
 Specify P Pihat
 Specify T threshold
 Matrix G 1 1 1 1
 Specify G asc 0 asc 0

 Covariances M;
 Thresholds T|T ;
 Frequency F ;
#if repeat_number < 4
! pihat is not 1 so weight is needed
 Weight \part(S,G)%J ;
#elseif repeat_number = 9
! pihat is not 1 so weight is needed
 Weight \part(S,G)%J ;
#elseif repeat_number = 7
! population sample, so no weight required
! Weight \part(S,G)%J ;
#elseif repeat_number = 10
! population sample, so no weight required
! Weight \part(S,G)%J ;
#else
! pihat is 1 but -- not sampled
 Weight \part(S,G)%J ;
#endif
 Options RSidual
End

G3: DZ twin pairs
 Data Ninput=9
 Labels studyid studyyr zygosity twin1 twin2 asc freq pihat threshold
 Ordinal File=scz_mx.dat
 Select if studyid = $repeat_number
 Select if zygosity = 2
 Select twin1 twin2 asc freq pihat threshold ;
 Definition Asc Freq Pihat Threshold ;
 Begin Matrices= Group 1 ;
  F full 1 1     ! Frequency
  G full 1 4 =G2 ! To select ascertainment probability
  P full 1 1 ! ascertainment probability p
  T full 1 1 ! threshold
 End Matrices;
 Specify F Freq
 Matrix G 1 1 1 1
 Specify G asc 0 asc 0
 Specify P Pihat
 Specify T threshold
 Begin Algebra;
  S= (P) _ (P) _ P+P-P*P ; ! Vector of ascertainment probabilities
  K= (P+P)*\mnor(U_O_T_T_U)-P*P*\mnor(L_O|O_T|T_T|T_U|U);
 End Algebra:
 Covariances L;
 Thresholds T|T ;
 Frequency F ;
#if repeat_number = 9
! pihat is not 1 so weight is needed
 Weight \part(S,G)%K ;
#elseif repeat_number = 10
! population sample, so no weight required (could estimate threshold)
! Weight \part(S,G)%K ;
#elseif repeat_number = 11
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#elseif repeat_number = 12
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#elseif repeat_number = 8
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#endif
 Options RSidual NDecimals=4 
 Option jiggle func=1.e-12
End


#define repeat_number 11
G2: MZ twin pairs
 Data Ninput=9
 Labels studyid studyyr zygosity twin1 twin2 asc freq pihat threshold
 Ordinal File=scz_mx.dat
 Select if studyid = $repeat_number
 Select if zygosity = 1;
 Select twin1 twin2 asc freq pihat threshold ;
 Definition Asc Freq Pihat Threshold ;
 Begin Matrices= Group 1;
  F full 1 1 ! Frequency
  G full 1 4 ! to select ascertainment prob
  P full 1 1 ! ascertainment probability p
  T full 1 1 ! threshold
 End Matrices;
 Begin Algebra;
  S= (P) _ (P) _ P+P-P*P ; ! Vector of ascertainment probabilities
  J= (P+P)*\mnor(U_O_T_T_U)-P*P*\mnor(M_O|O_T|T_T|T_U|U);
 End Algebra:

 Specify F Freq
 Specify P Pihat
 Specify T threshold
 Matrix G 1 1 1 1
 Specify G asc 0 asc 0

 Covariances M;
 Thresholds T|T ;
 Frequency F ;
#if repeat_number < 4
! pihat is not 1 so weight is needed
 Weight \part(S,G)%J ;
#elseif repeat_number = 9
! pihat is not 1 so weight is needed
 Weight \part(S,G)%J ;
#elseif repeat_number = 7
! population sample, so no weight required
! Weight \part(S,G)%J ;
#elseif repeat_number = 10
! population sample, so no weight required
! Weight \part(S,G)%J ;
#else
! pihat is 1 but -- not sampled
 Weight \part(S,G)%J ;
#endif
 Options RSidual
End

G3: DZ twin pairs
 Data Ninput=9
 Labels studyid studyyr zygosity twin1 twin2 asc freq pihat threshold
 Ordinal File=scz_mx.dat
 Select if studyid = $repeat_number
 Select if zygosity = 2
 Select twin1 twin2 asc freq pihat threshold ;
 Definition Asc Freq Pihat Threshold ;
 Begin Matrices= Group 1 ;
  F full 1 1     ! Frequency
  G full 1 4 =G2 ! To select ascertainment probability
  P full 1 1 ! ascertainment probability p
  T full 1 1 ! threshold
 End Matrices;
 Specify F Freq
 Matrix G 1 1 1 1
 Specify G asc 0 asc 0
 Specify P Pihat
 Specify T threshold
 Begin Algebra;
  S= (P) _ (P) _ P+P-P*P ; ! Vector of ascertainment probabilities
  K= (P+P)*\mnor(U_O_T_T_U)-P*P*\mnor(L_O|O_T|T_T|T_U|U);
 End Algebra:
 Covariances L;
 Thresholds T|T ;
 Frequency F ;
#if repeat_number = 9
! pihat is not 1 so weight is needed
 Weight \part(S,G)%K ;
#elseif repeat_number = 10
! population sample, so no weight required (could estimate threshold)
! Weight \part(S,G)%K ;
#elseif repeat_number = 11
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#elseif repeat_number = 12
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#elseif repeat_number = 8
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#endif
 Options RSidual NDecimals=4 
 Option jiggle func=1.e-12
End


#define repeat_number 12
G2: MZ twin pairs
 Data Ninput=9
 Labels studyid studyyr zygosity twin1 twin2 asc freq pihat threshold
 Ordinal File=scz_mx.dat
 Select if studyid = $repeat_number
 Select if zygosity = 1;
 Select twin1 twin2 asc freq pihat threshold ;
 Definition Asc Freq Pihat Threshold ;
 Begin Matrices= Group 1;
  F full 1 1 ! Frequency
  G full 1 4 ! to select ascertainment prob
  P full 1 1 ! ascertainment probability p
  T full 1 1 ! threshold
 End Matrices;
 Begin Algebra;
  S= (P) _ (P) _ P+P-P*P ; ! Vector of ascertainment probabilities
  J= (P+P)*\mnor(U_O_T_T_U)-P*P*\mnor(M_O|O_T|T_T|T_U|U);
 End Algebra:

 Specify F Freq
 Specify P Pihat
 Specify T threshold
 Matrix G 1 1 1 1
 Specify G asc 0 asc 0

 Covariances M;
 Thresholds T|T ;
 Frequency F ;
#if repeat_number < 4
! pihat is not 1 so weight is needed
 Weight \part(S,G)%J ;
#elseif repeat_number = 9
! pihat is not 1 so weight is needed
 Weight \part(S,G)%J ;
#elseif repeat_number = 7
! population sample, so no weight required
! Weight \part(S,G)%J ;
#elseif repeat_number = 10
! population sample, so no weight required
! Weight \part(S,G)%J ;
#else
! pihat is 1 but -- not sampled
 Weight \part(S,G)%J ;
#endif
 Options RSidual
End

G3: DZ twin pairs
 Data Ninput=9
 Labels studyid studyyr zygosity twin1 twin2 asc freq pihat threshold
 Ordinal File=scz_mx.dat
 Select if studyid = $repeat_number
 Select if zygosity = 2
 Select twin1 twin2 asc freq pihat threshold ;
 Definition Asc Freq Pihat Threshold ;
 Begin Matrices= Group 1 ;
  F full 1 1     ! Frequency
  G full 1 4 =G2 ! To select ascertainment probability
  P full 1 1 ! ascertainment probability p
  T full 1 1 ! threshold
 End Matrices;
 Specify F Freq
 Matrix G 1 1 1 1
 Specify G asc 0 asc 0
 Specify P Pihat
 Specify T threshold
 Begin Algebra;
  S= (P) _ (P) _ P+P-P*P ; ! Vector of ascertainment probabilities
  K= (P+P)*\mnor(U_O_T_T_U)-P*P*\mnor(L_O|O_T|T_T|T_U|U);
 End Algebra:
 Covariances L;
 Thresholds T|T ;
 Frequency F ;
#if repeat_number = 9
! pihat is not 1 so weight is needed
 Weight \part(S,G)%K ;
#elseif repeat_number = 10
! population sample, so no weight required (could estimate threshold)
! Weight \part(S,G)%K ;
#elseif repeat_number = 11
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#elseif repeat_number = 12
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#elseif repeat_number = 8
! pihat is 1 but -- not sampled
 Weight \part(S,G)%K ;
#endif
 Options RSidual NDecimals=4 
 Option jiggle func=1.e-12
End


Group 24 Constraint
 Constraint
 Begin Matrices = Group 2
 U Unit 1 1
 End matrices;

 Constraint U = A+C+E+D;
 Option multiple issat
End