From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>

Date: Mon, 18 Feb 2008 09:28:23 +0000 (GMT)

Date: Mon, 18 Feb 2008 09:28:23 +0000 (GMT)

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 ______________________________________________ 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 Mon 18 Feb 2008 - 09:35:27 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 Mon 18 Feb 2008 - 10: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.
*