**************************************************************** *! Version 9.0.1, 5 October 2006 *! Author: James Cui, Monash University *! Simulate three-generation families *! Original publication: Nov 2000 STB 58: 2-5 (dm82) **************************************************************** capture program drop simuped3 program simuped3 version 9.0 gettoken age1 0 : 0, parse(" ,") gettoken std1 0 : 0, parse(" ,") gettoken age2 0 : 0, parse(" ,") gettoken std2 0 : 0, parse(" ,") gettoken age3 0 : 0, parse(" ,") gettoken std3 0 : 0, parse(" ,") if (`age1'<0 | `age2'<0 | `age3'<0 | `std1'<0 | `std2'<0 | `std3'<0) { di in red "negative numbers invalid" exit 498 } syntax [, Saving(string) Reps(int 100) Alle(real 0.1) Sib(real 3) Si3(real 3)] if "`saving'" ~= "" { parse "`saving'", parse (" ") local output = "`1'" } else { local output = "temp.dta" } tempname nsimu mu mu1 p q simu tempvar degree id family nochild tlchild female age g x y i tempvar fchild mchild nogrand marry marriag 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 (`si3' < 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(5* (`sib' + `si3') + 4) qui set obs `nob' qui gen mu = `sib' qui gen mu1 = `si3' postfile simuped3 famid id degree family female marry g age using `output', replace qui gen `simu' = 0 while (`simu' < `nsimu') { /* begin simu */ qui replace `simu' = `simu' + 1 *--------------------------------------------------------------- * 1. 1st and 2nd generation *--------------------------------------------------------------- qui gen byte `degree' = cond(_n <= 4, 1, .) qui gen byte `id' = _n if `degree' == 1 qui gen byte `family' = cond(_n <= 2, 1, 2) if `degree' == 1 qui rndpoix mu qui gen byte `nochild' = xp qui replace `nochild'=. if `id'~=1 & `id'~=3 /* denote by 1st person in each family */ qui drop xp qui gen byte `tlchild' = `nochild'[1] + `nochild'[3] if `id' == 1 /* denote by 1st person */ qui replace `degree' = 2 if _n > 4 & _n <= 4 + `tlchild'[1] qui replace `id' = _n if `degree' == 2 qui replace `family' = cond(_n > 4 & _n <= 4 + `nochild'[1], 1, 2) if `degree' == 2 *--------------------------------------------------------------- * 2. SEX AND AGE *--------------------------------------------------------------- qui gen byte `female' = cond(_n == 2 | _n == 4, 1, 0) if `degree' == 1 /* the 2nd person female */ qui gen `x' = uniform() if `degree' == 2 qui replace `female'=cond(`x'<0.5, 0, 1) if `degree'==2 qui drop `x' 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 *--------------------------------------------------------------- * 3. GENERATION 1: GENOTYPE *--------------------------------------------------------------- qui gen `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 qui drop `x' *--------------------------------------------------------------- * 4. GENERATION 2: GENOTYPE *--------------------------------------------------------------- qui gen `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 qui drop `x' `y' *--------------------------------------------------------------- * 5. Number of Males and Females *--------------------------------------------------------------- qui sort `family' `degree' `female' qui by `family': gen `fchild' = sum(`female' == 1 & `degree' == 2) qui by `family': replace `fchild' = `fchild'[_N] if _n ~= _N qui replace `fchild' = . if `id' ~= 1 & `id' ~= 3 /* assign this number to grandpa */ qui gen `mchild' = `nochild' - `fchild' *--------------------------------------------------------------- * 5. Marriage * last F in family 1 marry first M in family 2 * first F in family 1 marry last M in family 2 * No marriag if either condition holds *--------------------------------------------------------------- qui sort `degree' `family' `female' qui gen `marry'=. /* check no. of females in family 1 > 1 */ if (1-`fchild'[1] <= 0 & 1-`mchild'[3] <= 0) { /* check no. of males in famiy 2 > 1 */ qui by `degree' `family': replace `marry' = 1 /* */ if `degree' == 2 & `family' == 1 & _n == _N & `female' == 1 qui by `degree' `family': replace `marry'=1 /* */ if `degree' == 2 & `family' == 2 & _n == 1 & `female' == 0 } else if (1-`mchild'[1] <= 0 & 1-`fchild'[3] <= 0) { qui by `degree' `family': replace `marry'=1 /* */ if `degree' == 2 & `family' == 1 & _n == 1 & `female' == 0 qui by `degree' `family': replace `marry'=1 /* */ if `degree' == 2 & `family' == 2 & _n == _N & `female' == 1 } *--------------------------------------------------------------- * 6. Generation 3, only if marriage exist *--------------------------------------------------------------- qui gen byte `nogrand' = . qui gen byte `marriag' = cond(1-`fchild'[1] <= 0 & 1-`mchild'[3] <= 0 | /* */ 1-`mchild'[1] <= 0 & 1-`fchild'[3] <= 0, 1, 0) if (`marriag'==1) { /* BEGIN marriag */ *--------------------------------------------------------------- * 7. Sibship size, generation 3 *--------------------------------------------------------------- qui rndpoix mu1 qui replace `nogrand' = xp if `id' == 1 /* assign to grandpa family 1 */ qui drop xp if (-`nogrand'[1] < 0) { qui replace `degree' = 3 if _n > 4 + `tlchild'[1] & /* */ _n <= 4 + `tlchild'[1] + `nogrand'[1] qui replace `id' = _n if `degree' == 3 qui replace `family' = 0 if `degree' == 3 *--------------------------------------------------------------- * 8. Sex *--------------------------------------------------------------- qui gen `x' = uniform() if `degree' == 3 qui replace `female' = cond(`x' < 0.5, 0, 1) if `degree' == 3 qui drop `x' *--------------------------------------------------------------- * 8. AGE and Genotype *--------------------------------------------------------------- qui replace `age' = max(0,int(`age3' + `std3' * invnorm(uniform()) + 0.5)) if `degree' == 3 qui sort `marry' `degree' `family' `female' /* Married couple listed first */ qui gen `x' = uniform() if `degree' == 3 qui gen `y' = uniform() if `degree' == 3 #delimit ; qui replace `g' = 11 if `degree' == 3 & (`g'[1] == 11 | `g'[1] == 12 & `x' < 0.5) & (`g'[2] == 11 | `g'[2] == 12 & `y' < 0.5) ; qui replace `g' = 22 if `degree' == 3 & (`g'[1] == 22 | `g'[1] == 12 & `x' >= 0.5) & (`g'[2] == 22 | `g'[2] == 12 & `y' >= 0.5) ; qui replace `g' = 12 if `degree' == 3 & (`g'[1] == 11 | `g'[1] == 12 & `x' < 0.5) & (`g'[2] == 22 | `g'[2] == 12 & `y' >= 0.5); qui replace `g' = 12 if `degree' == 3 & (`g'[1] == 22 | `g'[1] == 12 & `x' >= 0.5) & (`g'[2] == 11 | `g'[2] == 12 & `y' < 0.5) ; #delimit cr drop `x' `y' } } *--------------------------------------------------------------- * 9. OUTPUT TO STATA FILE *--------------------------------------------------------------- if (`marriag' == 1) { qui replace `marry' = 1 if `degree' == 1 qui replace `marry' = 0 if `marry' == . qui sort `id' qui gen `i' = 1 /* the ith observation */ while (`i' <= _N & `degree'[`i'] ~= .) { post simuped3 (`simu'[`i']) (`id'[`i']) (`degree'[`i']) (`family'[`i']) /* */ (`female'[`i']) (`marry'[`i']) (`g'[`i']) (`age'[`i']) qui replace `i'=`i' + 1 } qui drop `i' } if (`simu' / 10 == int(`simu' / 10)) { di in green `simu' _skip(1) _con /* display dots during simu */ } qui drop `degree' `id' `family' `nochild' `tlchild' `female' `age' `g' /* */ `fchild' `mchild' `nogrand' `marry' `marriag' } /* end SIMU */ postclose simuped3 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