On Dec 25, 2009, at 9:38 AM, James Rome wrote:

Thanks for the help.

I tried making the pdf file as suggested. Acrobat said it was damaged
and could not be opened. Is this an R bug?

Hard to say. Graphics devices vary from OS to OS and I am on a Mac using a 64 bit bit version of R 2.10.1. I get no error with opening the pdf file so created using either the native Mac graphic viewer, Preview, or the Adobe Acrobat Reader. To decide whether you have identified a bug you would need to provide full sessionInfo and code, and then someone using your OS could try to reproduce the problem.

It did make a PostScript file that I was able to distill into PDF, but
it was gray scales. How do I get the color back?
And yes, I did do the layout I wanted so I could see how the days
compared for each hour.

On 12/24/09 4:56 PM, David Winsemius wrote:


pdf(test.pdf")
xyplot(Sepal.Length + Sepal.Width ~ Petal.Length + Petal.Width |
Species, data = iris, scales = "free", layout = c(2, 1, 2), auto.key =
list(x = .6, y = .7, corner = c(0, 0)))
dev.off()
You may not be getting what you expect, but it may be that your plots
are all being created, but too quickly to be seen. Try printing to a
more durable "canvas".

And I would like to add a Poisson Distribution fit to each of these
plots (see below), but am clueless as to how to go about it.

I would like to fit a distribution to the count data for each
combination of day and hour, and I am unable to see how to do this in a
vector manner.  For example, I tried
density((Arrival.Val | DAY*Hour), na.rm=TRUE)
which does not work.

I should think the this would be informative:

glm(Arrival.Val ~ DAY*Hour, family="poisson")

Since DAY and Hour are factors you will get a large number of
estimates. You can use the typical regression functions, such as
predict() and summary() to get the fitted values.

I tried glm:
---------
glm(Arrival.Val ~ DAY*as.factor(Hour), family="poisson")

Call:  glm(formula = Arrival.Val ~ DAY * as.factor(Hour), family =
"poisson")


This output came across rather mangled.

Coefficients:
                          (Intercept)
DAY[T.Monday]
                              3.15396
-0.61348
                      DAY[T.Saturday]
DAY[T.Sunday]
                             -0.43853
-0.93475


 snipped


0.39860
 DAY[T.Tuesday]:as.factor(Hour)[T.23]
DAY[T.Wednesday]:as.factor(Hour)[T.23]
                              0.43209
0.49274

Degrees of Freedom: 8124 Total (i.e. Null);  7963 Residual
 (18 observations deleted due to missingness)
Null Deviance:        40120
Residual Deviance: 17030     AIC: 59170
----------------
I am not sure what to make of this.

Those are estimated Poisson means (on a log-scale) for each of your factors, DAY and Hour.


So how do I get the fit plotted on top of my histograms?

See if this helps understand how a model relates to a data situation:
> testsim <- data.frame(Arrivals = rpois(24*3, lambda=c(rep(10,5),rep(40,4),rep(20,6), rep(40,4), rep(10, 24-5-4-6-4))), Hour= factor(rep(1:24, 3)), DAY=Sys.Date()+1:3 )

>testsim

> testsim
   Arrivals Hour        DAY
1         9    1 2009-12-26
2         6    2 2009-12-27
3        10    3 2009-12-28
4        10    4 2009-12-26
5        13    5 2009-12-27
6        34    6 2009-12-28
7        34    7 2009-12-26
8        35    8 2009-12-27
# output elided (72 lines)

> glm(Arrivals ~ Hour, data=testsim)

Call:  glm(formula = Arrivals ~ Hour, data = testsim)

Coefficients:
(Intercept) Hour2 Hour3 Hour4 Hour5 Hour6 Hour7 8.3333 -1.3333 0.3333 0.6667 4.0000 29.6667 28.0000 Hour8 Hour9 Hour10 Hour11 Hour12 Hour13 Hour14 30.3333 32.3333 11.3333 11.3333 10.6667 11.3333 12.6667 Hour15 Hour16 Hour17 Hour18 Hour19 Hour20 Hour21 8.6667 31.0000 26.6667 33.0000 28.6667 4.0000 -1.0000
     Hour22       Hour23       Hour24
     0.3333       4.6667       0.3333

Degrees of Freedom: 71 Total (i.e. Null);  48 Residual
Null Deviance:      12400
Residual Deviance: 1004         AIC: 444.1

The estimates are Intercept + factor_coefficient.


Is there a way to save the bin data from the histogram command?

I don't know if the lattice function supports that, but the base graphics function hist lets you get the breaks and counts.

?hist

You might need to wrap it in a call to tapply, since the help page does not say that a formula method is available, and I do not see a formula method with methods(hist). (There might be easier ways to get counts, such as xtabs() which supports a formula method.

?xtabs



David Winsemius, MD
Heritage Laboratories
West Hartford, CT

Again Thanks for the prompt holiday response.
Jim Rome

David Winsemius, MD
Heritage Laboratories
West Hartford, CT

______________________________________________
R-help@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.

Reply via email to