**************************************************************** *! Version 9.0.1, 5 October 2006 *! Author: James Cui, Monash University *! Simulate two-generation families *! Original publication: Nov 2000 STB 58: 2-5 (dm82) **************************************************************** capture program drop simuped2 program simuped2 version 9.0 gettoken age1 0 : 0, parse(" ,") gettoken std1 0 : 0, parse(" ,") gettoken age2 0 : 0, parse(" ,") gettoken std2 0 : 0, parse(" ,") if (`age1'<0 | `age2'<0 | `std1'<0 | `std2'<0) { di in red "negative numbers invalid" exit 498 } syntax [, Saving(string) Reps(int 100) Alle(real 0.1) Sib(real 3)] if "`saving'" ~= "" { parse "`saving'", parse (" ") local output = "`1'" } else { local output = "temp.dta" } tempname nsimu p q simu if (`reps' < 1) { di in red "reps() required" exit 198 } scalar `nsimu' = `reps' if (`sib' < 0) { di in red "Sibship size negative" exit 198 } if (`alle' < 0) { di in red "Allele frequency negative" exit 198 } scalar `p' = `alle' scalar `q' = 1-`p' local nob = int(10*`sib' + 2) qui set obs `nob' qui gen mu = `sib' tempvar simu degree id female nochild x g y age i postfile simuped2 famid id degree female g age using `output', replace qui gen `simu' = 0 while (`simu' < `nsimu') { /* begin simu */ qui replace `simu' = `simu' + 1 tempvar degree id female nochild x g y age i *--------------------------------------------------------------- * 1. SIZE OF SIBSHIP *--------------------------------------------------------------- qui rndpoix mu qui gen byte `nochild' = xp qui replace `nochild' = . if _n ~= 1 /* denoted by 1st person */ qui drop xp qui gen byte `degree' = 1 if _n <= 2 qui replace `degree' = 2 if _n > 2 & _n <= 2 + `nochild'[1] qui gen byte `id' = _n if _n <= 2 + `nochild'[1] *--------------------------------------------------------------- * 2. SEX AND AGE *--------------------------------------------------------------- qui gen byte `female' = 0 if `degree' == 1 qui replace `female' = 1 if _n == 2 & `degree' == 1 /* father first */ qui gen `x' = uniform() if `degree' == 2 qui replace `female' = 0 if `x' < 0.5 & `degree' == 2 qui replace `female' = 1 if `x' >= 0.5 & `degree' == 2 #delimit ; qui gen int `age' = max(0,int(`age1' + `std1' * invnorm(uniform()) + 0.5)) if `degree' == 1; qui replace `age' = max(0,int(`age2' + `std2' * invnorm(uniform()) + 0.5)) if `degree' == 2; #delimit cr *--------------------------------------------------------------- * 3. GENERATION 1: GENOTYPE *--------------------------------------------------------------- qui replace `x' = uniform() if `degree' == 1 qui gen byte `g' = 11 if `x' <= `p' * `p' & `degree' == 1 qui replace `g' = 12 if `x' > `p' * `p' & `x' <= `p' * `p' /* */ + 2 * `p' * `q' & `degree' == 1 qui replace `g' = 22 if `x' > `p' * `p' + 2 * `p' * `q' /* */ & `degree' == 1 *--------------------------------------------------------------- * 4. GENERATION 2: GENOTYPE *--------------------------------------------------------------- qui replace `x' = uniform() if `degree' == 2 qui gen `y' = uniform() if `degree' == 2 #delimit ; qui replace `g' = 11 if `degree' == 2 & (`g'[1] == 11 | `g'[1] == 12 & `x' < 0.5) & (`g'[2] == 11 | `g'[2] == 12 & `y' < 0.5); qui replace `g' = 22 if `degree' == 2 & (`g'[1] == 22 | `g'[1] == 12 & `x' >= 0.5) & (`g'[2] == 22 | `g'[2] == 12 & `y' >= 0.5); qui replace `g' = 12 if `degree' == 2 & (`g'[1] == 11 | `g'[1] == 12 & `x' < 0.5) & (`g'[2] == 22 | `g'[2] == 12 & `y' >= 0.5); qui replace `g' = 12 if `degree' == 2 & (`g'[1] == 22 | `g'[1] == 12 & `x' >= 0.5) & (`g'[2] == 11 | `g'[2] == 12 & `y' < 0.5); #delimit cr *--------------------------------------------------------------- * 5. OUTPUT TO STATA FILE *--------------------------------------------------------------- qui gen `i'=1 /* the ith observation */ while (`i' <= 2 + `nochild'[1]) { post simuped2 (`simu'[`i']) (`id'[`i']) (`degree'[`i']) /* */ (`female'[`i']) (`g'[`i']) (`age'[`i']) qui replace `i' = `i' + 1 } if (`simu' / 10 == int(`simu' / 10)) { di in green `simu' _skip(1) _con /* display dots during simu */ } qui drop `degree' `id' `female' `nochild' `x' `g' `y' `age' `i' } /* end SIMU */ postclose simuped2 use `output', clear qui gen str2 genotype = "AA" if g == 11 qui replace genotype = "Aa" if g == 12 qui replace genotype = "aa" if g == 22 qui drop g qui save `output', replace di _n in ye `"A file named "`output'" has been created"' drop _all end