From: Gustaf Granath <gustaf.granath_at_ebc.uu.se>

Date: Sun, 17 Feb 2008 22:17:57 +0100

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 17 Feb 2008 - 21:21:59 GMT

Date: Sun, 17 Feb 2008 22:17:57 +0100

Dear John and Brian,

Thank you for your help. I get the feeling that it is something
fundamental that I do not understand here. Furthermore, a day of
reading did not really help so maybe we have reached a dead end here.
Nevertheless, here comes one last try.

What I want:

I am trying to make a barplot with adjusted means with SEs (error
bars), with the y axis labeled on the response scale.

Estimate Std. Error t value Pr(>|t|) (Intercept) 1.9368528 0.1057948 18.308 < 2e-16 *** initial.size 0.0015245 0.0001134 13.443 < 2e-16 *** inf.level50 -0.3142688 0.0908063 -3.461 0.000612 *** def.level12.5 -0.2329221 0.1236992 -1.883 0.060620 .def.level25 -0.1722354 0.1181993 -1.457 0.146062 def.level50 -0.3543826 0.1212906 -2.922 0.003731 **

(Dispersion parameter for quasipoisson family taken to be 6.431139)

Null deviance: 2951.5 on 322 degrees of freedom Residual deviance: 1917.2 on 317 degrees of freedom

library(effects)

def <- effect("def.level",model,se=TRUE)
summary(def)

$effect

def.level

0 12.5 25 50
11.145382 8.829541 9.381970 7.819672

$lower

def.level

0 12.5 25 50

9.495220 7.334297 7.867209 6.467627

$upper

def.level

0 12.5 25 50

13.08232 10.62962 11.18838 9.45436

#Confidence intervals makes sense and are in line with the glm model
result. Now #lets look at the standard errors. Btw, why aren't they
given with summary?

def$se

324 325 326 327 0.08144281 0.09430438 0.08949864 0.09648573

# As you can see, the SEs are very very very small. #In a graph it would look weird in combination with the glm result. #I thought that these values were logs. Thats why I used exp() whichseems to be wrong.

Regards,

Gustaf

> Quoting John Fox <jfox@mcmaster.ca>:

*> Dear Brian and Gustaf,
**>
**> I too have a bit of trouble following what Gustaf is doing, but I think that
**> Brian's interpretation -- that Gustaf is trying to transform the standard
**> errors via the inverse link rather than transforming the ends of the
**> confidence intervals -- is probably correct. If this is the case, then what
**> Gustaf has done doesn't make sense.
**>
**> It is possible to get standard errors on the scale of the response (using,
**> e.g., the delta method), but it's probably better to work on the scale of
**> the linear predictor anyway. This is what the summary, print, and plot
**> methods in the effects package do (as is documented in the help files for
**> the package -- see the transformation argument under ?effect and the type
**> argument under ?summary.eff).
**>
**> Regards,
**> John
**>
**> --------------------------------
**> John Fox, Professor
**> Department of Sociology
**> McMaster University
**> Hamilton, Ontario, Canada L8S 4M4
**> 905-525-9140x23604
**> http://socserv.mcmaster.ca/jfox
**>
**>
**>> -----Original Message-----
**>> From: Prof Brian Ripley [mailto:ripley_at_stats.ox.ac.uk]
**>> Sent: February-17-08 6:42 AM
**>> To: Gustaf Granath
**>> Cc: John Fox; r-help_at_r-project.org
**>> Subject: Re: [R] Weird SEs with effect()
**>>
**>> On Sun, 17 Feb 2008, Gustaf Granath wrote:
**>>
**>> > Hi John,
**>> >
**>> > In fact I am still a little bit confused because I had read the
**>> > ?effect help and the archives.
**>> >
**>> > ?effect says that the confidence intervals are on the linear
**>> predictor
**>> > scale as well. Using exp() on the untransformed confidence intervals
**>> > gives me the same values as summary(eff). My confidence intervals
**>> > seems to be correct and reflects the results from my glm models.
**>> >
**>> > But when I use exp() to get the correct SEs on the response scale I
**>> > get SEs that sometimes do not make sense at all. Interestingly I have
**>>
**>> What exactly are you doing here? I suspect you are not using the
**>> correct
**>> formula to transform the SEs (you do not just exponeniate them), but
**>> without the reproducible example asked for we cannot tell.
**>>
**>> > found a trend. For my model with adjusted means ~ 0.5-1.5 I get huge
**>> > SEs (SEs > 1, but my glm model shows significant differences between
**>> > level 1 = 0.55 and level 2 = 1.15). Models with means around 10-20 my
**>> > SEs are fine with exp(). Models with means around 75-125 my SEs get
**>> > way too small with exp().
**>> >
**>> > Something is not right here (or maybe they are but I don not
**>> > understand it) so I think my best option will be to use the
**>> confidence
**>> > intervals instead of SEs in my plot.
**>>
**>> If you want confidence intervals, you are better off computing those on
**>> a
**>> reasonable scale and transforming then. Or using a profile likelihood
**>> to
**>> compute them (which will be equivariant under monotone scale
**>> transformations).
**>>
**>> > Regards,
**>> >
**>> > Gustaf
**>> >
**>> >
**>> >> Quoting John Fox <jfox_at_mcmaster.ca>:
**>> >>
**>> >> Dear Gustaf,
**>> >>
**>> >> From ?effect, "se: a vector of standard errors for the effect, on
**>> the scale
**>> >> of the linear predictor." Does that help?
**>> >>
**>> >> Regards,
**>> >> John
**>> >>
**>> >> --------------------------------
**>> >> John Fox, Professor
**>> >> Department of Sociology
**>> >> McMaster University
**>> >> Hamilton, Ontario, Canada L8S 4M4
**>> >> 905-525-9140x23604
**>> >> http://socserv.mcmaster.ca/jfox
**>> >>
**>> >>
**>> >>> -----Original Message-----
**>> >>> From: r-help-bounces_at_r-project.org [mailto:r-help-bounces_at_r-
**>> >>> project.org] On Behalf Of Gustaf Granath
**>> >>> Sent: February-16-08 11:43 AM
**>> >>> To: r-help_at_r-project.org
**>> >>> Subject: [R] Weird SEs with effect()
**>> >>>
**>> >>> Hi all,
**>> >>>
**>> >>> Im a little bit confused concerning the effect() command, effects
**>> >>> package.
**>> >>> I have done several glm models with family=quasipoisson:
**>> >>>
**>> >>> model <-glm(Y~X+Q+Z,family=quasipoisson)
**>> >>>
**>> >>> and then used
**>> >>>
**>> >>> results.effects <-effect("X",model,se=TRUE)
**>> >>>
**>> >>> to get the "adjusted means". I am aware about the debate concerning
**>> >>> adjusted means, but you guys just have to trust me - it makes sense
**>> >>> for me.
**>> >>> Now I want standard error for these means.
**>> >>>
**>> >>> results.effects$se
**>> >>>
**>> >>> gives me standard error, but it is now it starts to get confusing.
**>> The
**>> >>> given standard errors are very very very small - not realistic. I
**>> >>> thought that maybe these standard errors are not back transformed
**>> so I
**>> >>> used exp() and then the standard errors became realistic. However,
**>> for
**>> >>> one of my glm models with quasipoisson the standard errors make
**>> kind
**>> >>> of sense without using exp() and gets way to big if I use exp(). To
**>> be
**>> >>> honest, I get the feeling that Im on the wrong track here.
**>> >>>
**>> >>> Basically, I want to know how SE is calculated in effect() (all I
**>> know
**>> >>> is that the reported standard errors are for the fitted values) and
**>> if
**>> >>> anyone knows what is going on here.
**>> >>>
**>> >>> Regards,
**>> >>>
**>> >>> Gustaf Granath
**>> >>>
**>> >>> ______________________________________________
**>> >>> R-help_at_r-project.org mailing list
**>> >>> https://stat.ethz.ch/mailman/listinfo/r-help
**>> >>> PLEASE do read the posting guide http://www.R-project.org/posting-
**>> >>> guide.html
**>> >>> and provide commented, minimal, self-contained, reproducible code.
*

*>
*

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 17 Feb 2008 - 21:21:59 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sun 17 Feb 2008 - 23:30:14 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*