On 2005-04-07, at 22:30, INAGAMI,SANAE wrote:
1) My sample was not a random sample. It oversampled certain
populations. Can I still use nlcom to calculate CI's on predicted
probabilities in multilevel ordered logistic regresssions using gllamm?
If you only include only fixed effects when you calculate the CIs then
it is based on an "average" subject in the population i.e. with all
random effects set to zero (the group mean). If you have a strong bias
in your sampling, this average subject might not be the one you think
it is. Also, if you have a strong bias in your sample, the random
effects might not be normally distributed which is assumed by the
model. However, -gllamm- takes sampling weights that might solve your
problem. See the -pweight()- option in the helpfile/manual.
You could also test if sampling predicts the estimated random effects,
after the fact, by regressing group dummys or continuos variables (e.g.
age if you oversampled older or younger subjects) on the latent
variables. see the -geqs()- option. If they did have an effect you
might have bias in your estimated variance components.
Another way is to use -gllapred u ,u- for subject specific empirical
Bayes predictions of the random effects. A significant ANOVA on these
predictions with your sampling groups as a factor would suggest that
you have group variation that might have biased your estimates. You
could also plot these estimates to assess the assumption of normality.
If you use this method, be sure you only use one observation/subject.
. gllapred u , u
. egen pick = tag(subject_id)
. anova u1 sampling_group if pick
. histogram u1 if pick, normal
2) Assuming the answer is yes...Having used the syntax that M.Ingre so
kindly provided...can I add the effect of another covariate into the
model.
The simple answer is yes. It is possible to modify the code for ANY
combination of the covariates (including random effects).
the code works...but is it ok methodology.
thanks.
This is a very important question. I think it is. Some might find it
odd to use a logistic model to calculate RRRs. But if I'm not mistaken,
this is exactly what a multinomial logit (mlogit) does. But for other
experts on this list to be able to comment, a replication of the code
might help:
clear
sysuse auto
// Create level 2 ID
gen make2 = substr(make, 1, index(make, " "))
egen pick =tag(make2)
gen id=_n if pick
sort make2 id
replace id = id[_n-1] if id == .
// Make turn>40 a dichotomous predictor variable
gen turn40 = turn > 40
// estimate gllamm model on rep78 with turn40 as predictor
gllamm rep78 turn40 , i(id) link(ologit) fam(bin) adapt
// Use local macros to temporary hold your input values
local pv1 = 1 // predictor value 1
local pv2 = 0 // predictor value 2
local cut = 11 // Choose cut point i.e. which level of the dependent
variable you want to predict
// calculate log of the ratio of the two probabilities predicted from
pv1/pv2
nlcom ln( (exp(`pv1' * [rep78]turn40 - [_cut`cut']_cons) / ///
(1 + exp(`pv1' * [rep78]turn40 - [_cut`cut']_cons))) / ///
(exp(`pv2' * [rep78]turn40 - [_cut`cut']_cons) / ///
(1 + exp(`pv2' * [rep78]turn40 - [_cut`cut']_cons))))
// Store the results from -nlcom-
matrix b = r(b)
matrix V = r(V)
// Estimates needs to be exponentiated to describe RRR
di "RRR: " exp(b[1,1])
di "Lower 95% CI: " exp(b[1,1] - sqrt(V[1,1])*1.96)
di "Upper 95% CI: " exp(b[1,1] + sqrt(V[1,1])*1.96)
*
* 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/
