Statalist


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

Re: st: ttest or xtmelogit?


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/



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