[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