[R] Weird SEs with effect()
John Fox
jfox at mcmaster.ca
Sun Feb 17 23:10:41 CET 2008
Dear Gustaf,
> -----Original Message-----
> From: Gustaf Granath [mailto:gustaf.granath at ebc.uu.se]
> Sent: February-17-08 4:18 PM
> To: John Fox
> Cc: 'Prof Brian Ripley'; r-help at r-project.org
> Subject: RE: [R] Weird SEs with effect()
>
> 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.
>
> I thought that the values produced by effect() were logs (e.g. in
> $fit). And then they were transformed (antilogged) with summary(). Was
> I wrong?
I'm sorry that you're continuing to have problems with this.
Yes, there is a point that you don't understand: The SEs are on the scale of
the log-counts, but you can't get correct SEs on the scale of the counts by
exponentiating the SEs on the scale of the log-counts. What summary(), etc.,
do (and you can do) to produce confidence intervals on the count scale is
first to compute the intervals on the log-count scale and then to transform
the end-points.
I'm afraid that I can't make the point more clearly than that.
I hope this helps,
John
>
> 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.
>
> #One of my GLM models (inf.level & def.level=factors, initial.size =
> covariate) #used as an example.
> #I was not able to make a reproducible example though. Sorry.
>
> model <-
> glm(tot.fruit~initial.size+inf.level+def.level,family=quasipoisson)
> summary(model)
> Coefficients:
> 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() which
> seems to be wrong.
>
> Regards,
>
> Gustaf
>
>
> > Quoting John Fox <jfox at 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.
>
> >
More information about the R-help
mailing list