[R] Weird SEs with effect()
Gustaf Granath
gustaf.granath at ebc.uu.se
Sun Feb 17 22:17:57 CET 2008
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?
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