#delimit ;
log using twbiodem-analysis;
*! Hans-Peter Kohler
*! kohler@demogr.mpg.de
*! August 20, 2000


*! Demonstration of bivariate (ordered) probit models
*! for twin data;

*! load data
*! this is a wide dataset in which each twin pair is one line
*! in the dataset and the information for twin 1 and 2 in a
*! pair are indicated with the suffice 1 and 2 to each
*! individual-specific variable;
use twbiodem-wtwpairs, clear;
format sex zygos %12.0g;
format birthy nkids1 nkids2 oyelem1 oyelem2 %6.0g;


*! LIST PRIMARY VARIABLES IN DATA
*! Variables are:
*!    twpair:    id of twin pair
*!    sex        sex (same sex twin pairs only)
*!    zygos      zygosity
*!    birthy     birth year of twin pair
*!    nkids      number of children at age 35
*!    oyelem     years of elementary education
*!    oyseco     years of secondary education
*! ;

list sex zygos birthy nkids1 nkids2 oyelem1 oyelem2 in 1/10, nodisplay;

*! SUMMARY OF THESE VARIABLES;

tab1 sex zygos;

su birthy nkids1 nkids2 oyelem1 oyelem2;



*! CROSS TAB OF NUMBER OF CHILDREN AND
*! PEARSON CORRELATIONS BY SEX AND ZYGOSITY;

tab nkids*;

sort sex zygos;
by sex zygos: corr nkids*;


*! ESTIMATION OF BIVARIATE PROBIT MODEL FOR HAVING
*! AT LEAST ON CHILD
*!
*! because there may be a cohort trend in childlessness, we include
*! birth year among the covariates;

*! estimation only for females;
for num 1/2: g anykidX = (nkidsX > 0) if nkidsX != .;
g c2 = 1;
g h2 = Radd;

*! note: the option norhocons is included because we defined a
*! varable c2 that is a constant;
twbiopr anykid1 anykid2 if female == 1,
        xbeqc(birthy)
        rhoeq(c2 h2)
        norhocons;

*! the same estimation can also be combined for males and
*! females, where include a dummy for females to account for
*! differences in the level of childlessness across sex;
g male_c2 = female == 0;
g fem_c2 = female;
g male_h2 = Radd * (female == 0);
g fem_h2 = Radd * female;

twbiopr anykid1 anykid2,
        xbeqc(birthy female)
        rhoeq(male_c2 male_h2 fem_c2 fem_h2)
        norhocons;

*! since this model yields insignificant values for c2 for
*! both males and females, we may want to re-estimate as;

twbiopr anykid1 anykid2,
        xbeqc(birthy female)
        rhoeq(male_h2 fem_h2)
        norhocons;

*! and we can test of whether the male and female h2 are different
*! as;

test male_h2 = fem_h2;

*! finally, since the years of elementary education may be a
*! determant of having any children, we add the years of
*! education as an individual-specific covariate;

twbiopr anykid1 anykid2 if female == 1,
        xbeqc(birthy)
        xbeqi(oyelem)
        rhoeq(c2 h2)
        norhocons;



*! NUMBER OF CHILDREN IN CATEGORIES 0, 1, 2, 3 and more;

*! because only relatively few indiviuals have more than
*! three children (see above), we combine all of them into a
*! single category 3+ children
*! this can be done by adding the truncate option in twbiopr
*! ;

twbiopr nkids1 nkids2 if female == 1,
        xbeqc(birthy)
        rhoeq(c2 h2)
        maxcat(3) truncate
        norhocons;

*! again, this estimation can be combined for both sexes as
*! (given our estimates above, we also specify starting values);

twbiopr nkids1 nkids2,
        xbeqc(birthy female)
        rhoeq(male_c2 male_h2 fem_c2 fem_h2)
        maxcat(3) truncate
        rhogues(.2,.3,.15,.5)
        norhocons;

*! finally, we re-estimate the first model, adding education
*! (elementary and secondary) as an individual-specific covarate;

twbiopr nkids1 nkids2 if female == 1,
        xbeqc(birthy)
        xbeqi(oyelem oyseco)
        rhoeq(c2 h2)
        maxcat(3) truncate
        rhogues(0,.5)
        norhocons;

log close;