The S-plus code below is for a sib-based permutation association test 
without any parental information, which is discussed in the following
reference:
 
Allison, D. B., Heo, M., Kaplan, N., & Martin, E. R.  (1999). Development
of Sibling-based Tests of Linkage in the Presence of Association for
Quantitative Traits That Do not Require Parental Information. American
Journal of Human Genetics, 64:1754-1764.
 
The three key parameters for the program are:
 
sibship: family id number -- must be the same for the sibling members in a
         same sibship;
geno:    the number of an allele ("A" or "a"), not the genotype -- must be
         0,1, and 2 for a bi-allelic marker;
pheno:   continuous phenotype values -- these could be residuals from some
         models for adjustments.
 
This program allows missing values for "geno" and "pheno". This program
is only for bi-allelic markers. 

We acknowledge that this program is not in a distributable form, and do not 
exclude any possibilities of coding error. 

The S-plus code
===============

sibbase.perm_function(sibship,geno,pheno){

ord.sibship_order(sibship)

sibship_sibship[ord.sibship]
geno_geno[ord.sibship]
pheno_pheno[ord.sibship]

NApheno_is.na(pheno)
NAgeno_is.na(geno)
NAyes_(NApheno|NAgeno)
sibship_sibship[!NAyes]
geno_geno[!NAyes]
pheno_pheno[!NAyes]

sibsize_tapply(sibship,list(sibship),length)
sib.number_rep(sibsize,sibsize)
mulsib_(sib.number > 1)

sibship_sibship[mulsib]
geno_geno[mulsib]
pheno_pheno[mulsib]
geno1_2-geno

sibsize_tapply(sibship,list(sibship),length)
pheno.sum1_tapply(pheno,list(sibship),sum)
pheno.sum2_tapply(pheno^2,list(sibship),sum)
geno.sum1_tapply(geno,list(sibship),sum)
geno.sum2_tapply(geno^2,list(sibship),sum)
geno1.sum1_tapply(geno1,list(sibship),sum)
geno1.sum2_tapply(geno1^2,list(sibship),sum)
phgn.sum_tapply(pheno*geno,list(sibship),sum)
phgn1.sum_tapply(pheno*geno1,list(sibship),sum)

muij_cbind(pheno.sum1*geno.sum1/sibsize,pheno.sum1*geno1.sum1/sibsize)
Vij1_cbind(pheno.sum2*geno.sum2/sibsize,pheno.sum2*geno1.sum2/sibsize)
Vij2_cbind((pheno.sum1^2-pheno.sum2)*(geno.sum1^2-geno.sum2)/sibsize/(sibsize-1),
     (pheno.sum1^2-pheno.sum2)*(geno1.sum1^2-geno1.sum2)/sibsize/(sibsize-1))
Vij_Vij1+Vij2-muij^2
pg_cbind(phgn.sum,phgn1.sum)
Snum_(apply(pg-muij,2,sum))^2
Sden_apply(Vij,2,sum)
Sperm_sum(Snum/Sden)/2
Pvalue_1-pchisq(Sperm,1)
list(chisq=Sperm,Pvalue=Pvalue,n=length(sibship),nsibship=length(unique(sibship)))

}