Dear Imai,
Thank you for picking up on the syntax error. After correcting it, and making a couple
other minor changes, I think I'm running into the same limitation as I did attempting
to use "logit.mixed --- the response is 3-level ordinal, but here the logit function
is strictly defined for binary responses only.
I'm interpreting the output below "NA's introduced by coercion" to
mean that when the function encounters the 3rd response level, the response value is
replaced with NA.
I can try to the strategy suggested elsewhere: collapse the 2 lowest categories to make a
binary response.
Paul
z.out <- zelig(formula=Response ~ Product +
Application, id = "ID", data=Q3_soapfeel, model="logit.gee")
Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
running glm to get initial regression estimate
(Intercept) Product Application
31.26514 -0.02468 -0.70156
Warning in function (formula = formula(data), id = id, data = parent.frame(), :
NAs introduced by coercion
Error in function (formula = formula(data), id = id, data = parent.frame(), :
NA/NaN/Inf in foreign function call (arg 2)
Paul Prew | Statistician
651-795-5942 | fax 651-204-7504
Ecolab Research Center | Mail Stop ESC-F4412-A
655 Lone Oak Drive | Eagan, MN 55121-1560
-----Original Message-----
From: Kosuke Imai [mailto:kimai@Princeton.EDU]
Sent: Sunday, April 12, 2009 10:02 PM
To: Prew, Paul
Cc: zelig(a)lists.gking.harvard.edu
Subject: Re: [zelig] FW: ordinal response models with logit.mixed - plots
It seems:
id = "Subject"+ data=soap.data,
should be:
id = "Subject", data=soap.data,
Kosuke
--
Department of Politics
Princeton University
http://imai.princeton.edu
On Tue, 7 Apr 2009, Prew, Paul wrote:
I apologize if this was double-posted, I didn't see the original email
make it onto the list
Hello,
In short, I would like some advice on how to interpret and hopefully
plot the results from an ordinal mixed-effects model.
I am working on an analysis to predict the comfort of repeated
hand-washing with soap formulations. Ratings of Like, Indifferent or
Dislike (1, 2 or 3) were recorded after 10 consecutive hand-washings.
Dozens of people participated in the study. (This appears to be a good
candidate for logit.gee also, but this function appears to only accept
binomial responses, not multinomial. I got this response with logit.gee
>> "ERROR: Cgee" >>>
when trying the code: z.out <- zelig(formula=
Score ~ Product + Application,
id = "Subject"+ data=soap.data,
model="logit.gee")).
Subject Product Question Application Score Response
1 1 869 2 1 3 Like
2 1 869 2 2 3 Like
3 1 869 2 3 3 Like
4 1 869 2 4 3 Like
5 1 869 2 5 3 Like
6 1 869 2 6 3 Like
str(soap.data)
'data.frame':
560 obs. of 6 variables:
$ Subject : Factor w/ 56 levels
"1","2","3","4",..: 1 1 1 1 1 1 1 1
1 1 ...
$ Product : Factor w/ 2 levels "143","869": 2 2 2 2 2 2 2 2 2 2
...
$ Question : int 2 2 2 2 2 2 2 2 2 2 ...
$ Application: int 1 2 3 4 5 6 7 8 9 10 ...
$ Score : Factor w/ 3 levels "1","2","3": 3 3 3 3 3 3
3 3 3 3 ...
$ Response : Factor w/ 3 levels "Dislike","Indifferent",..: 3 3 3 3
3
3 3 3 3 3 ...
Logit.mixed
GENERICALLY: z.out <- zelig(formula= Response ~ x1 + x2 + tag(z1 + z2 |
g),
data=mydata, model="logit.mixed")
My interest is in knowing how the score changes as Applications
increase, comparing the 2 Products.
Using the random-intercepts model
z.out <- zelig(formula= Score ~ Product + Application + tag(1 |
Subject),
data=soap.data, model="logit.mixed")
I got the following output
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 9.0872 3.0488 2.981 0.0029 **
Product[T.869] -1.4118 3.5758 -0.395 0.6930
Application -0.0314 0.1147 -0.274 0.7839
Correlation of Fixed Effects:
(Intr) P[T.86
Prdc[T.869] -0.814
Application -0.216 0.004
The logit.mixed analysis worked fine, and I got what I wanted up to a
point. I now have 2 questions ---
1. Does anyone have any suggestions for using Zelig output to plot
the model predictions? I'd like to get something such as the confidence
interval plots that the Effects package produces for multinomial models.
But Effects doesn't recognize the "mer", "lme4" or
"gee" object classes.
Specifically, plots of confidence bands around the "Application" slope,
per each Product, are what would be helpful to the scientist I'm
assisting.
2. Could someone help interpret the simulation output? For the
last line of the output, it states Predicted Value Y|X. I assume this
is conditional on the average input values, ( not very interesting in
this case), but to me it seems to be predicting output values of 0 or 1.
The response scale was 1,2, or 3; zero was not possible. Clearly I
don't understand what Predicted Value Y|X output means.
x.out <- setx(z.out)
s.out <- sim(z.out, x = x.out)
summary(s.out)
Model: logit.mixed
Number of simulations: 1000
Values of X
(Intercept) Product[T.869] Application
1 1 0 5.481
Expected Values: E(Y|X)
mean sd 2.5% 97.5%
1 0.9962 0.01909 0.9694 1
Predicted Values: Y|X
0 1
1 0.066 0.934
CONFIDENTIALITY NOTICE:
This e-mail communication and any attachments may contain proprietary and privileged
information for the use of the designated recipients named above.
Any unauthorized review, use, disclosure or distribution is prohibited.
If you are not the intended recipient, please contact the sender by reply e-mail and
destroy all copies of the original message.
-
Zelig Mailing List, served by Harvard-MIT Data Center
Send messages: zelig(a)lists.gking.harvard.edu
[un]subscribe Options:
http://lists.gking.harvard.edu/?info=zelig
Zelig program information:
http://gking.harvard.edu/zelig/
CONFIDENTIALITY NOTICE:
This e-mail communication and any attachments may contain proprietary and privileged
information for the use of the designated recipients named above.
Any unauthorized review, use, disclosure or distribution is prohibited.
If you are not the intended recipient, please contact the sender by reply e-mail and
destroy all copies of the original message.
-
Zelig Mailing List, served by Harvard-MIT Data Center
Send messages: zelig(a)lists.gking.harvard.edu
[un]subscribe Options:
http://lists.gking.harvard.edu/?info=zelig
Zelig program information:
http://gking.harvard.edu/zelig/