Re: st: permutations

[David Airey <david.airey@vanderbilt.edu>]

.

I calculated a test of a difference between any two SNP genotypes. You  
are saying you want to know which pairwise comparisons are different.  
If you just want those two probabilities (genotype 2 = genotype 0,  
genotype 1 = genotype 0), you get these by adding a couple lines to  
our code, and modifying the permute command as below. To get the other  
comparison (genotype 2 = genotype 1), you would need to use a  
different base comparison, or you would need to use the command lincom.

***do_begin***

describe
table age caco snp

program plogistic, rclass
	version 10
	args caco snp age
	xi: logistic `caco' i.`snp' i.`age'
	test _Isnp_1 _Isnp_2
	return scalar chi2 = r(chi2)
	test _Isnp_1
	return scalar chi2_1 = r(chi2)
	test _Isnp_2
	return scalar chi2_2 = r(chi2)
end

set seed 3132008

permute caco chi2=r(chi2) chi2_1=r(chi2_1) chi2_2=r(chi2_2), ///
	reps(1000): plogistic caco snp age

***do_end***


Monte Carlo permutation results                   Number of obs    
=       2394

      command:  plogistic caco snp age
         chi2:  r(chi2)
       chi2_1:  r(chi2_1)
       chi2_2:  r(chi2_2)
  permute var:  caco

------------------------------------------------------------------------------
T            |     T(obs)       c       n   p=c/n   SE(p) [95% Conf.  
Interval]
------------- 
+----------------------------------------------------------------
        chi2 |   .1560796     920    1000  0.9200  0.0086  . 
9014202    .936058
      chi2_1 |   .1184664     724    1000  0.7240  0.0141  . 
6951618   .7515127
      chi2_2 |   .1448991     719    1000  0.7190  0.0142   . 
690024   .7466803
------------------------------------------------------------------------------
Note:  confidence intervals are with respect to p=c/n.
Note:  c = #{|T| >= |T(obs)|}

.


On Mar 17, 2008, at 5:23 AM, Christopher Intemann wrote:


> Hello David,
> thank you for the example, it was very helpful.
> However, still the result is not what I'm actually up to.
> In your example, you're calculating:
> > xi: logistic caco i.snp i.age
> > i.snp             _Isnp_0-2           (naturally coded; _Isnp_0  
> > omitted)
> > i.age             _Iage_35-60         (naturally coded; _Iage_35  
> > omitted)
> >
> > Logistic regression                               Number of obs    
> > =       2116
> >                                                 LR chi2(7)       
> > =      21.17
> >                                                 Prob > chi2      
> > =     0.0035
> > Log likelihood = -1455.9021                       Pseudo R2        
> > =     0.0072
> >
> > ------------------------------------------------------------------------------
> >       caco | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf.  
> > Interval]
> > ------------- 
> > +----------------------------------------------------------------
> >    _Isnp_1 |   .9580097   .1193996    -0.34   0.731     .7503816     
> > 1.223088
> >    _Isnp_2 |   .9524632   .1218649    -0.38   0.703     .7412068     
> > 1.223931
> >   _Iage_40 |    1.63372   .2251451     3.56   0.000      
> > 1.247017     2.14034
> >   _Iage_45 |   1.545441    .216941     3.10   0.002     1.173722     
> > 2.034885
> >   _Iage_50 |   1.634428   .2582141     3.11   0.002     1.199194     
> > 2.227627
> >   _Iage_55 |    1.09121    .184666     0.52   0.606     .7831755     
> > 1.520399
> >   _Iage_60 |   1.202247   .1956237     1.13   0.258     .8739581     
> > 1.653854
> > ------------------------------------------------------------------------------
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/