[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
David Airey <david.airey@Vanderbilt.Edu> |

To |
statalist@hsphsun2.harvard.edu |

Subject |
Re: st: ttest or xtmelogit? |

Date |
Tue, 11 Mar 2008 11:32:06 -0500 |

.

I'm running now at lower values of rho, which are more typical of my data (.05 < rho < .15).

I agree simulating at larger rho's is smart because estimating rho accurately is not easy with small samples, so I should look at the upper bound of the estimated CI of rho.

This is fascinating stuff! I wonder if reviewers will buy this kind of supporting documentation for analysis choice down the road.

Thanks Joseph.

-Dave

On Mar 11, 2008, at 9:25 AM, Joseph Coveney wrote:

I wrote:

generate double mu = 0.5 + trt * `delta' + ///

`s' * _pi / sqrt(3) * invnormal(uniform())

It would be better to have:

generate double mu = invlogit(trt * `delta' + ///

sqrt(`rho' / (1 - `rho') * _pi^2 / 3) * invnormal(uniform()))

Conclusions regarding the t test with sum scores versus transformed mean scores are the same.

Joseph Coveney

clear *

set more off

set seed `=date("2008-03-11", "YMD")'

set seed0 `=date("2008-03-11", "YMD")'

*

capture program drop simem

program define simem, rclass

version 10

syntax [, Delta(real 0) Rho(real 0)]

drop _all

tempname p_t p_l

set obs 12

generate int pid = _n

generate byte trt = mod(_n, 2)

generate double mu = invlogit(trt * `delta' + ///

sqrt(`rho' / (1 - `rho') * _pi^2 / 3) * invnormal(uniform()))

assert !missing(mu) // if rho is misspecified

generate byte den = 50

rndbinx mu den

*

generate double pos_prim = bnlx

replace pos_prim = 50 - 1 / 2 / 50 if pos_prim == 50

replace pos_prim = 1 / 2 / 50 if !pos_prim

generate double logit_pi = logit(pos_prim / den)

ttest bnlx, by(trt)

scalar define `p_t' = r(p)

ttest logit_pi, by(trt)

scalar define `p_l' = r(p)

*

rename bnlx count1

generate byte count0 = den - count1

reshape long count, i(pid) j(positive)

drop if !count

expand count

xtlogit positive trt, i(pid) intpoints(30)

return scalar rho = e(rho)

return scalar p_r = chi2tail(1, e(chi2))

return scalar delta = _b[trt]

return scalar p_t = `p_t'

return scalar p_l = `p_l'

end

*

forvalues delta = 0/2 {

forvalues rho = 0(0.35)0.7 {

display in smcl as text _newline(1) ///

"delta = " %03.1f `delta' " rho = " %03.1f `rho'

quietly simulate delta = r(delta) rho = r(rho) p_r = r(p_r) ///

p_t = r(p_t) p_l = r(p_l), reps(300): ///

simem , rho(`rho') delta(`delta')

summarize delta rho

foreach var of varlist p_* {

generate byte pos_`var' = `var' < 0.05

}

summarize pos_*

}

}

exit

It takes time to run the do-file, because of the addition of - xtlogit- in order to verify delta and rho. For anyone interested but time-pressed, results are listed below:

Null hypothesis; independence

delta (mean SD): 0.01 0.16

rho (mean SD): 0.00 0.01

Type I error rate (nominal 5%)

random effects logistic regresion: 3%

Student's t test on untransformed: 4%

t test on logistic-transformed mean: 4%

Null hypothesis; rho = 30%

delta (mean SD): -0.02 0.77

rho (mean SD): 0.30 0.10

Type I error rate (nominal 5%)

random effects logistic regresion: 9%

Student's t test on untransformed: 2%

t test on logistic-transformed mean: 2%

Null hypothesis; rho = 70%

delta (mean SD): 0.08 1.71

rho (mean SD): 0.64 0.12

Type I error rate (nominal 5%)

random effects logistic regresion: 10%

Student's t test on untransformed: 4%

t test on logistic-transformed mean: 5%

Delta = 1; independence

delta (mean SD): 1.0 0.18

rho (mean SD): 0.00 0.00

Power (12 mice)

random effects logistic regresion: 100% (300/300)

Student's t test on untransformed: 100%

t test on logistic-transformed mean: 100%

Delta = 1; rho = 30%

delta (mean SD): 0.99 0.83

rho (mean SD): 0.29 0.11

Power (12 mice)

Student's t test on untransformed: 24%

t test on logistic-transformed mean: 21%

Delta = 1; rho = 70%

delta (mean SD): 1.1 1.6

rho (mean SD): 0.64 0.12

Power (12 mice)

Student's t test on untransformed: 9%

t test on logistic-transformed mean: 6%

Delta = 2; independence

delta (mean SD): 2.0 0.21

rho (mean SD): 0.00 0.01

Power (12 mice)

As for delta = 1

Delta = 2; rho = 30%

delta (mean SD): 1.9 0.82

rho (mean SD): 0.28 0.11

Power (12 mice)

Student's t test on untransformed: 58%

t test on logistic-transformed mean: 51%

Delta = 2; rho = 70%

delta (mean SD): 2.1 1.7

rho (mean SD): 0.61 0.13

Power (12 mice)

Student's t test on untransformed: 23%

t test on logistic-transformed mean: 23%

*

* 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/

-- David C. Airey, Ph.D. Pharmacology Research Assistant Professor Center for Human Genetics Research Member Department of Pharmacology School of Medicine Vanderbilt University Rm 8158A Bldg MR3 465 21st Avenue South Nashville, TN 37232-8548 TEL (615) 936-1510 FAX (615) 936-3747 EMAIL david.airey@vanderbilt.edu URL http://people.vanderbilt.edu/~david.c.airey/dca_cv.pdf URL http://www.vanderbilt.edu/pharmacology * * 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/

**References**:**Re: st: ttest or xtmelogit?***From:*"Joseph Coveney" <jcoveney@bigplanet.com>

- Prev by Date:
**st: booststraping in two stage models** - Next by Date:
**Re: st: RE: Stata has disappointing support of contrasts/multiple comparisons in mixed ANOVA** - Previous by thread:
**Re: st: ttest or xtmelogit?** - Next by thread:
**Re: st: ttest or xtmelogit?** - Index(es):

© Copyright 1996–2014 StataCorp LP | Terms of use | Privacy | Contact us | What's new | Site index |