####################################

## simulating the proband method

## for estimating segregation ratio

####################################

p.hat_rep(0,10)

p_rep(0,10)

# we will repeat the whole process 10 times to get

# and average and a sample sd

for(k in 1:10){

# create a "population" of 5 child families

# whose parents can have affected childeren with

# probability 1/4

mat <- matrix(rbinom(50000,1,.25),ncol=5)

# create a proband vector = # of probands in each family where

# the ascertainment prob for individual = .20

proband <- rep(0,10000)

for(i in 1:10000){

if(sum(mat[i,])>0){

proband[i] <- rbinom(1,sum(mat[i,]),.20)

}#if

}#for

# create matrix of selected families (at least one proband)

selected.fams <- mat[proband>=1,]

proband <- proband[proband>=1]

# vector of number affected in selected.fams

num.aff <- rep(0,dim(selected.fams)[1])

for(i in 1:dim(selected.fams)[1]){

num.aff[i] <- sum(selected.fams[i,])

}#for

# find out estimate for p

num <- sum(proband*(num.aff-1))

denom <- sum(proband*(5-1))

p.hat[k] <- (num/denom)

# what was actual p?  should be about .25

p[k] <- (sum(mat)/(dim(mat)[1]*dim(mat)[2]))

}#for

#### OUTPUT ################################################################

> p.hat

[1] 0.2431373 0.2445088 0.2500000 0.2598542 0.2574411 0.2503903 0.2426864

[8] 0.2469475 0.2401980 0.2487041

> p

[1] 0.24902 0.24880 0.24898 0.25124 0.25028 0.24848 0.24680 0.25112 0.24682

[10] 0.25106

> p.hat-p

[1] -0.005882745 -0.004291214  0.001020000  0.008614158  0.007161051

[6]  0.001910320 -0.004113556 -0.004172503 -0.006621980 -0.002355853

> stdev(p.hat-p)

[1] 0.005367616

> mean(p.hat-p)

[1] -0.0008732322

############################################################################