Strange error with using se = FALSE in geom_smooth()
ONKELINX, Thierry
I get a strange error when using se = FALSE in geom_smooth. Even
stranger is that it does not occur when using smaller subsamples of the
dataset. The dataset is attached.
Any ideas?
Best regards,
Thierry
> library(ggplot2)
> ggplot(test, aes(x = A, y = B, colour = C)) + geom_smooth()
> ggplot(test, aes(x = A, y = B, colour = C)) + geom_smooth(se = FALSE)
Error in class(output[[var]]) <- class(value) :
cannot set class to "array" unless the dimension attribute has length
> 0
> ggplot(tail(test, 3900), aes(x = A, y = B, colour = C)) +
geom_smooth(se = FALSE)
> ggplot(head(test, 3900), aes(x = A, y = B, colour = C)) +
geom_smooth(se = FALSE)
> ggplot(head(test, 4000), aes(x = A, y = B, colour = C)) +
geom_smooth(se = FALSE)
Error in class(output[[var]]) <- class(value) :
cannot set class to "array" unless the dimension attribute has length
> 0
> sessionInfo()
R version 2.11.1 (2010-05-31)
i386-pc-mingw32
locale:
[1] LC_COLLATE=Dutch_Belgium.1252 LC_CTYPE=Dutch_Belgium.1252
[3] LC_MONETARY=Dutch_Belgium.1252 LC_NUMERIC=C
[5] LC_TIME=Dutch_Belgium.1252
attached base packages:
[1] grid grDevices datasets splines graphics stats tcltk
[8] utils methods base
other attached packages:
[1] mgcv_1.6-2 ggplot2_0.8.8 proto_0.3-8 reshape_0.8.3
[5] plyr_1.0.3 survival_2.36-1
loaded via a namespace (and not attached):
[1] digest_0.4.2 lattice_0.19-9 Matrix_0.999375-43 nlme_3.1-96
[5] tools_2.11.1
> summary(test)
A B C
Min. : 1.00 Min. : 1.00 a:2723
1st Qu.: 13.00 1st Qu.:12.00 b:2723
Median : 19.00 Median :16.00 c:5446
Mean : 21.43 Mean :16.94
3rd Qu.: 27.00 3rd Qu.:21.00
Max. :130.00 Max. :50.00
------------------------------------------------------------------------
----
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium
Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
Thierry....@inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey
Druk dit bericht a.u.b. niet onnodig af.
Please do not print this message unnecessarily.
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in this message
and any annex are purely those of the writer and may not be regarded as stating
an official position of INBO, as long as the message is not confirmed by a duly
signed document.
Brian Diggs
On Jul 26, 3:37 am, "ONKELINX, Thierry" <Thierry.ONKEL...@inbo.be>
wrote:
>
> I get a strange error when using se = FALSE in geom_smooth. Even
> stranger is that it does not occur when using smaller subsamples of the
> dataset. The dataset is attached.
>
> Any ideas?
>
> Best regards,
>
> Thierry
down the conditions needed to trigger it, and I have a suggestion for
a temporary workaround until the real problem can be understood.
The conditions needed are:
* data with at least two group
* the first group with 1000 or more cases
* a geom_smooth with se=FALSE
If there is only a single group, it works with or without se=FALSE,
independent of the number of cases.
If the first group (where "first" is defined based on the order of the
levels of the grouping variable) has less than 1000 cases, it works
with or without se=FALSE. I assume the 1000 is significant because
that is the threshold which stat_smooth switches from loess to gam.
Without se=FALSE, any combination of size and order works.
As a workaround, if you add method="loess" to you geom_smooth call, it
should work.
=====Example code:
# create a data.frame that satisfies the criteria for the error to
occur
ex2 <- data.frame(A=c(1:1000,seq(1,1000,length=20)),
B=rnorm(1020),
C=rep(c("a","b"),c(1000,20)))
# has all three elements
ggplot(ex2, aes(x=A, y=B, group=C)) + geom_smooth(se=FALSE)
# Error in class(output[[var]]) <- class(value) :
# cannot set class to "array" unless the dimension attribute has
length > 0
# does not include the se=FALSE
ggplot(ex2, aes(x=A, y=B, group=C)) + geom_smooth()
# works
# rearrange the order of the levels on C so the smaller (<1000 cases)
# group is first;
ggplot(transform(ex2, C=factor(C,levels=c("b","a"))),
aes(x=A, y=B, group=C)) + geom_smooth(se=FALSE)
# works
# take away the grouping completely
ggplot(ex2, aes(x=A, y=B)) + geom_smooth(se=FALSE)
# works
# force the smoothing method to be "loess"
ggplot(ex2, aes(x=A, y=B, group=C)) + geom_smooth(se=FALSE,
method="loess")
# works
# Forced use of gam also triggers the error, even with a smaller data
set
ex3 <- data.frame(A=c(1:20,1:20),
B=rnorm(40),
C=rep(c("a","b"),each=20))
# request gam, using the same formula as default in StatSmooth
$calculate
# when more than 1000 cases
ggplot(ex3, aes(x=A, y=B, group=C)) +
geom_smooth(se=FALSE, method="gam", formula=y~s(x,bs="cs"))
#Error in class(output[[var]]) <- class(value) :
# cannot set class to "array" unless the dimension attribute has
length > 0
# take away the grouping
ggplot(ex3, aes(x=A, y=B)) +
geom_smooth(se=FALSE, method="gam", formula=y~s(x,bs="cs"))
# works, though I am not sure that the formula argument is doing
anything
# or take away the se=FALSE
ggplot(ex3, aes(x=A, y=B, group=C)) +
geom_smooth(method="gam", formula=y~s(x,bs="cs"))
# works, though I am not sure that the formula argument is doing
anything
=====End example code
Given these examples, the criteria could be restated as:
* data with at least two group
* a geom_smooth with se=FALSE
* the first group has its smoothing performed with "gam", either by
direct request or by virtue of having 1000 or more cases and using the
default auto selection of smoother
I'm going to let someone else with a better understanding of the guts
of the code figure out why these particular criteria lead to an error.
> sessionInfo()
R version 2.11.1 (2010-05-31)
i386-pc-mingw32
locale:
States.1252
[3] LC_MONETARY=English_United States.1252
LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] grid stats graphics grDevices utils datasets
other attached packages:
[1] mgcv_1.6-2 ggplot2_0.8.8 proto_0.3-8 reshape_0.8.3
loaded via a namespace (and not attached):
Matrix_0.999375-42 nlme_3.1-96
[6] tools_2.11.1
--
Brian Diggs
Senior Research Associate, Department of Surgery, Oregon Health &
Science University
> To call in the statistician after the experiment is done may be no more
> than asking him to perform a post-mortem examination: he may be able to
> say what the experiment died of.
> ~ Sir Ronald Aylmer Fisher
>
> The plural of anecdote is not data.
> ~ Roger Brinner
>
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of
> data.
> ~ John Tukey
>
> Druk dit bericht a.u.b. niet onnodig af.
> Please do not print this message unnecessarily.
>
> Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer
> en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
> door een geldig ondertekend document. The views expressed in this message
> and any annex are purely those of the writer and may not be regarded as stating
> an official position of INBO, as long as the message is not confirmed by a duly
> signed document.
>
> 51KViewDownload
takahashi kohske
probably this is a bug in ggplot2:::predictdf.glm
the function should be
function (model, xseq, se, level)
{
pred <- stats::predict(model, newdata = data.frame(x = xseq),
se = se, type = "link")
if (se) {
std <- qnorm(level/2 + 0.5)
data.frame(x = xseq, y = model$family$linkinv(as.vector(pred$fit)),
ymin = model$family$linkinv(as.vector(pred$fit -
std * pred$se)), ymax =
model$family$linkinv(as.vector(pred$fit +
std * pred$se)), se = as.vector(pred$se))
}
else {
# the class of pred is array, so it shoud be coerce to vector
# data.frame(x = xseq, y = model$family$linkinv(pred))
data.frame(x = xseq, y = model$family$linkinv(as.vector(pred)))
}
}
the strange behavior is because StatSmooth choose "loess" as fitting
method when the number of cases < 1000 and without specifying the
method.
possibly this bug will be fixed in the next version of ggplot2.
> --
> You received this message because you are subscribed to the ggplot2 mailing list.
> Please provide a reproducible example: http://gist.github.com/270442
>
> To post: email ggp...@googlegroups.com
> To unsubscribe: email ggplot2+u...@googlegroups.com
> More options: http://groups.google.com/group/ggplot2
>
Hadley Wickham
<takahash...@gmail.com> wrote:
> hi,
>
> probably this is a bug in ggplot2:::predictdf.glm
>
> the function should be
>
> function (model, xseq, se, level)
> {
> pred <- stats::predict(model, newdata = data.frame(x = xseq),
> se = se, type = "link")
> if (se) {
> std <- qnorm(level/2 + 0.5)
> data.frame(x = xseq, y = model$family$linkinv(as.vector(pred$fit)),
> ymin = model$family$linkinv(as.vector(pred$fit -
> std * pred$se)), ymax =
> model$family$linkinv(as.vector(pred$fit +
> std * pred$se)), se = as.vector(pred$se))
> }
> else {
> # the class of pred is array, so it shoud be coerce to vector
> # data.frame(x = xseq, y = model$family$linkinv(pred))
> data.frame(x = xseq, y = model$family$linkinv(as.vector(pred)))
> }
> }
>
> the strange behavior is because StatSmooth choose "loess" as fitting
> method when the number of cases < 1000 and without specifying the
> method.
>
> possibly this bug will be fixed in the next version of ggplot2.
I'm happy to make this change, but I would prefer to know why the
class of predict is an array - any ideas?
Hadley
--
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
http://had.co.nz/
takahashi kohske
> class of predict is an array - any ideas?
I don't know the reason. But help of predict.gam says they return
array, and actually if argument type="term", the return value has to
be array. Anyway, the author of gam would know the reason.
Hadley Wickham
>> class of predict is an array - any ideas?
>
> I don't know the reason. But help of predict.gam says they return
> array, and actually if argument type="term", the return value has to
> be array. Anyway, the author of gam would know the reason.
I see I've done it everywhere else, so I obviously just missed that line.
Thanks for the fix and it's now included in the devel version.