|
|  |
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]
Re: st: nested design and logistic regression
Thanks Joseph for your advice. I have tried gllamm before asking the list
but was not sure how to setup the model and was a bit reluctant to use this
difficult feature of Stata (for me!). The gllamm model worked with my data
but, as a newby to gllamm, I don't really understand your last lines: the
purpose of setting up a matrix and re-run the gllamm command.
Denis
From: Joseph Coveney <[email protected]>
Reply-To: [email protected]
To: Statalist <[email protected]>
Subject: Re: st: nested design and logistic regression
Date: Thu, 06 Jul 2006 17:44:32 +0900
Denis Haine wrote:
I have a 2 factors nested design model, where treatment received by
siblings
(low dose vs. high dose) is nested in same treatment received by their
maternal parents (received or not). I could fit an anova for nested design
for a continuous dependent variable, e.g. anova y trt_mat / trt_sib|trt_mat
/, but how to fit a logistic regression taking into account this nesting of
the 2 treatments if I have a binary dependent variable?
--------------------------------------------------------------------------------
Wouldn't it be something like that below?
It's helpful to start out in familiar territory in order to get your
bearings, and progress in stages to terra incognita, with checkpoints along
the way.
So, the trail starts with -anova-,
and then passes through -xtmixed , reml- (comparing the variance
components)
to -xtmixed , ml-
to -gllamm , family(gaussian) link(identity)-
before arriving at -glamm , family(binomial) link(logit)-.
The ANOVA setup for nested factors (a.k.a. hierarchal design) can be found
in B. J. Winer, D. R. Brown & K. M. Michels, _Statistical Principles in
Experimental Design_ Third Edition. (New York: McGraw-Hill, 1991),
pp. 358--62. www.stata.com/bookstore/sped.html
I used -xtmixed- as the bridge, an illustrative analogy, to -gllamm-.
Joseph Coveney
clear
set more off
set memory 100M
set matsize 11000
set seed `=date("2006-07-06", "ymd")'
set obs 200
generate int mother = _n
generate float sigma_mother = 3 * invnorm(uniform())
generate treatment = mod(_n, 2)
forvalues dose = 0/1 {
generate float sigma_motherdose`dose' = 2 * ///
invnorm(uniform())
}
reshape long sigma_motherdose, i(mother) j(dose)
forvalues child = 1/2 {
generate float sigma_child`child' = ///
invnorm(uniform()) * _pi / sqrt(3)
}
reshape long sigma_child, i(mother dose) j(child)
generate float response = treatment + sigma_mother + ///
dose + sigma_motherdose + sigma_child
*
* Begin here
*
anova response treatment / mother|treatment ///
dose treatment*dose / mother|treatment*dose /
local n = e(N) / (e(df_m) + 1)
local sigma2_e = e(rss) / e(df_r)
local MSmotherdose = e(ssdenom_4) / e(dfdenom_4)
local sigma2_motherdose = (`MSmotherdose' - `sigma2_e' ) / `n'
local sigma2_mother = (e(ssdenom_1)/ e(dfdenom_1) - ///
`MSmotherdose') / `n' / (e(df_4) + 1)
display `sigma2_mother'
display `sigma2_motherdose'
display `sigma2_e'
*
* Bridge
*
quietly xi3 e.treatment*e.dose
quietly compress
egen int motherdose = group(mother dose)
xtmixed response _I* || mother: || motherdose:, reml ///
nolrtest variance nolog
xtmixed response _I* || mother: || motherdose:, mle ///
nolrtest variance nolog
gllamm response _I*, i(motherdose mother) nrf(1,1) ///
family(gaussian) link(identity) adapt nip(12) nolog
*
* Analogous logistic regression
*
summarize response, meanonly
generate byte binary_response = response > r(mean)
quietly gllamm binary_response _I*, i(motherdose mother) nrf(1,1) ///
family(binomial) link(logit)
matrix A = e(b)
gllamm binary_response _I*, i(motherdose mother) nrf(1,1) ///
family(binomial) link(logit) from(A) adapt nip(30)
exit
*
* 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/
*
* 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/
