[R] Weird SEs with effect()
Prof Brian Ripley
ripley at stats.ox.ac.uk
Mon Feb 18 10:28:23 CET 2008
On Mon, 18 Feb 2008, Gustaf Granath wrote:
> Dear John & Brian,
>
> Ok. Now I start to understand. So basically I cannot use the given SEs for
> my purposes. They only make sense on the scale of log-counts. However in a
> paper, the log-count scale is not very informative (the reader want to see
> the effect on the scale you measure). If I understand you right, the
> confidence intervals are fine though and maybe I will use them to illustrate
> the reliability of the estimate. The problem is that showing SEs of the
> adjusted means has become sort of the standard way to illustrate things in my
> field (too many SAS users ?) and I might run into trouble with the reviewers.
> I have several data sets with similar data and my plan was to use effect().
> That is why I want to figure this out. I hope I haven't been too annoying ;).
>
> Finally, is there a way to get correct SEs on the count scale (with adjusted
> means)??
> I guess not, judging by your answers.
Yes, there is, at least approximately. If X has mean mu and se s,
exp(X) has approximate mean exp(mu) and se s*exp(mu). As in
> X <- rnorm(1000, 10, 0.1)
> sd(X)
[1] 0.09811928
> sd(exp(X))
[1] 2172.124
> exp(10)
[1] 22026.47
You will find slightly more accurate formulae on the help page ?plnorm
(they are exact if you assume normality of the estimated effects).
>
> Thanks again for your help,
>
> Gustaf
>
>
> John Fox wrote:
>> 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
>>>>>>>>
>>>>>>>> ______________________________________________
>>>>>>>>
>>
>
> -
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-help
mailing list