Statalist


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

Re: st: ttest or xtmelogit?


From   "Joseph Coveney" <jcoveney@bigplanet.com>
To   "Statalist" <statalist@hsphsun2.harvard.edu>
Subject   Re: st: ttest or xtmelogit?
Date   Tue, 11 Mar 2008 23:25:34 +0900

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/




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