Re: [R] Mapping 2D to 3D

[Jeff Newmiller]

No, but fortunately you are off in the weeds. Density has an 
internally-computed "z" coordinate... you should be looking at ?geom_contour.

On September 17, 2020 7:17:33 PM PDT, H  wrote:
>I am trying to understand how to map 2D to 3D using ggplot() and
>eventually plot_gg(). I am, however, stuck on understanding how to
>express the third variable to be mapped. This example:
>
>ggdiamonds = ggplot(diamonds, aes(x, depth)) +
>stat_density_2d(aes(fill = stat(nlevel)),
>geom = "polygon", n = 100, bins = 10,contour = TRUE) +
>facet_wrap(clarity~.) +
>scale_fill_viridis_c(option = "A")
>
>uses a variable nlevel that I now understand is calculated during the
>building of the ggplot but I have not figured out from where it is
>calculated or how to specify a variable of my choosing.
>
>Does anyone have a good reference for understanding how to specify this
>variable? Most examples on the 'net seem to use the same dataset but do
>not specify this particular aspect...
>
>__
>R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>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.

-- 
Sent from my phone. Please excuse my brevity.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Mapping 2D to 3D

[Abby Spurdle]

> But there's no reason for the user to do that when using the plotting 
> function.

I should amend the above.
There's no reason for the user to do that (compute a third "variable"
representing density), if using a high level plotting function, that's
designed to compute the density for you.

It is possible to do it in two or more steps.
(Compute the density matrix, then plot it).

Again, refer to the help file for kde2d.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Mapping 2D to 3D

[Abby Spurdle]

I'm not familiar with the gg graphics system.
However, I am familiar with density estimation, and density visualization.

There is *no* third variable, as such.
But rather, density estimates, which in this context, would usually be a matrix.
(And are computed inside the plotting or density estimation functions).

The documentation for the function you've used, says it uses  MASS::kde2d().
This does just that, returns an object, which contains a density matrix.
(Refer to the help file for kde2d).

Of course, there's no reason why one can't create a third variable,
from a mathematical perspective.
e.g. d, z, h, fv, or whatever you prefer...
And then set z = fh (x, y).

But there's no reason for the user to do that when using the plotting function.

Note that there are situations where one might want to set the limits
of the plot.
And set the breaks, colors, and color key.
e.g. Creating two or more plots, and putting them next to each other,
for comparison purposes.


On Fri, Sep 18, 2020 at 2:17 PM H  wrote:
>
> I am trying to understand how to map 2D to 3D using ggplot() and eventually 
> plot_gg(). I am, however, stuck on understanding how to express the third 
> variable to be mapped. This example:
>
> ggdiamonds = ggplot(diamonds, aes(x, depth)) +
> stat_density_2d(aes(fill = stat(nlevel)),
> geom = "polygon", n = 100, bins = 10,contour = TRUE) +
> facet_wrap(clarity~.) +
> scale_fill_viridis_c(option = "A")
>
> uses a variable nlevel that I now understand is calculated during the 
> building of the ggplot but I have not figured out from where it is calculated 
> or how to specify a variable of my choosing.
>
> Does anyone have a good reference for understanding how to specify this 
> variable? Most examples on the 'net seem to use the same dataset but do not 
> specify this particular aspect...
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Mapping 2D to 3D

[H]

I am trying to understand how to map 2D to 3D using ggplot() and eventually 
plot_gg(). I am, however, stuck on understanding how to express the third 
variable to be mapped. This example:

ggdiamonds = ggplot(diamonds, aes(x, depth)) +
stat_density_2d(aes(fill = stat(nlevel)),
geom = "polygon", n = 100, bins = 10,contour = TRUE) +
facet_wrap(clarity~.) +
scale_fill_viridis_c(option = "A")

uses a variable nlevel that I now understand is calculated during the building 
of the ggplot but I have not figured out from where it is calculated or how to 
specify a variable of my choosing.

Does anyone have a good reference for understanding how to specify this 
variable? Most examples on the 'net seem to use the same dataset but do not 
specify this particular aspect...

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] multiple time series plots

[Jim Lemon]

Hi Peri,
In that case, try this:

ylim<-range(c(df$Nation,df$East,df$Center),na.rm=TRUE)
plot(df$year, df$Nation, type="o", col="blue", lty=1, ylab="change",ylim=ylim)

Jim

On Fri, Sep 18, 2020 at 11:10 AM peri He  wrote:
>
> Hi Jim,
>
> Thank you for your reply.
>
> Sorry I am messaging you instead of asking publicly.
> Actually Nation, center and East are like medium and confidence intervals.
> So no, they are not in the limit of each other.
> I am trying to find a way to plot them in one chart. Most of the codes give 
> me error.
>
> 
> From: Jim Lemon 
> Sent: Friday, September 18, 2020 12:21 AM
> To: peri He 
> Cc: r-help@r-project.org 
> Subject: Re: [R] multiple time series plots
>
> Hi Peri,
> Without the data this is only a guess, but are the values of East and
> center within the limits of the plot generated by Nation?
>
> Jim
>
> On Fri, Sep 18, 2020 at 10:12 AM peri He  wrote:
> >
> > Dear Friends,
> >
> > I am trying to add two time series curves into one plot. I don't get any 
> > error but the lines are not added to my plot.
> > I am trying the following code;
> >
> > 1) df <- as.data.frame (read_excel("C:/Users/NO2.xlsx", "trend))
> >
> > 2) plot(df$year, df$Nation, type="o", col="blue", lty=1, ylab="change" )
> >
> > 3) lines(df$year, df$East, col="red",lty=2)
> >
> > 4) lines(df$year, df$center, col="dark red", lty=3)
> >
> > I need to mention that Nation, East and Center are some probability values.
> > And year, is a sequence of years. Lines 1 and 2 are run properly but lines 
> > 3 and 4 are not added.
> >
> > I would appreciate any idea.
> >
> > Regards,
> >
> > Peri
> >
> > [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] how to overlay two histograms

[Ana Marija]

HI Jim,

fantastic solution!
Thank you so much!!!

Ana

On Thu, Sep 17, 2020 at 6:01 PM Jim Lemon  wrote:
>
> Hi Ana,
> Sorry it's not in ggplot, but it may help:
>
> d<-read.table(text="CHR counts name
>   1 193554  old
>   2 220816  old
>   3 174350  old
>   4 163112  old
>   5 168125  old
>   6 182366  old
>   7 143023  old
>   8 147410  old
>   9 122112  old
>  10 138394  old
>  11 130069  old
>  12 124850  old
>  13 104119  old
>  14  83931  old
>  15  72287  old
>  16  71550  old
>  17  58380  old
>  18  76812  old
>  19  37040  old
>  20  63407  old
>  21  33863  old
>  22  33812  old
>   1 202783  new
>   2 252124  new
>   3 213337  new
>   4 201001  new
>   5 207606  new
>   6 228133  new
>   7 147218  new
>   8 177518  new
>   9 121276  new
>  10 163447  new
>  11 158724  new
>  12 142183  new
>  13 89  new
>  14  83043  new
>  15  61063  new
>  16  55439  new
>  17  32883  new
>  18  69135  new
>  19  16624  new
>  20  48541  new
>  21  25479  new
>  22  19698  new",
> header=TRUE,stingsAsFactors=FALSE)
> barpos<-barplot(counts~name+CHR,data=d,beside=TRUE,names.arg=rep("",22))
> legend(40,22,c("new","old"),fill=c("gray20","gray80"))
> library(plotrix)
> staxlab(1,at=colMeans(barpos),labels=1:22)
>
> Jim
>
> On Fri, Sep 18, 2020 at 8:05 AM Ana Marija  
> wrote:
> >
> > Hello,
> >
> > I am trying to overlay two histograms with this:
> >
> > p <- ggplot(d, aes(CHR, counts, fill = name)) + geom_bar(position = "dodge")
> > p
> >
> > but I am getting this error:
> > Error: stat_count() can only have an x or y aesthetic.
> > Run `rlang::last_error()` to see where the error occurred.
> >
> > my data is this:
> >
> > > d
> >CHR counts name
> > 11 193554  old
> > 22 220816  old
> > 33 174350  old
> > 44 163112  old
> > 55 168125  old
> > 66 182366  old
> > 77 143023  old
> > 88 147410  old
> > 99 122112  old
> > 10  10 138394  old
> > 11  11 130069  old
> > 12  12 124850  old
> > 13  13 104119  old
> > 14  14  83931  old
> > 15  15  72287  old
> > 16  16  71550  old
> > 17  17  58380  old
> > 18  18  76812  old
> > 19  19  37040  old
> > 20  20  63407  old
> > 21  21  33863  old
> > 22  22  33812  old
> > 23   1 202783  new
> > 24   2 252124  new
> > 25   3 213337  new
> > 26   4 201001  new
> > 27   5 207606  new
> > 28   6 228133  new
> > 29   7 147218  new
> > 30   8 177518  new
> > 31   9 121276  new
> > 32  10 163447  new
> > 33  11 158724  new
> > 34  12 142183  new
> > 35  13 89  new
> > 36  14  83043  new
> > 37  15  61063  new
> > 38  16  55439  new
> > 39  17  32883  new
> > 40  18  69135  new
> > 41  19  16624  new
> > 42  20  48541  new
> > 43  21  25479  new
> > 44  22  19698  new
> >
> > Basically I need to show counts per CHR in "old" and "new" side by side.
> >
> > Please advise,
> > Ana
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] multiple time series plots

[Jim Lemon]

Hi Peri,
Without the data this is only a guess, but are the values of East and
center within the limits of the plot generated by Nation?

Jim

On Fri, Sep 18, 2020 at 10:12 AM peri He  wrote:
>
> Dear Friends,
>
> I am trying to add two time series curves into one plot. I don't get any 
> error but the lines are not added to my plot.
> I am trying the following code;
>
> 1) df <- as.data.frame (read_excel("C:/Users/NO2.xlsx", "trend))
>
> 2) plot(df$year, df$Nation, type="o", col="blue", lty=1, ylab="change" )
>
> 3) lines(df$year, df$East, col="red",lty=2)
>
> 4) lines(df$year, df$center, col="dark red", lty=3)
>
> I need to mention that Nation, East and Center are some probability values.
> And year, is a sequence of years. Lines 1 and 2 are run properly but lines 3 
> and 4 are not added.
>
> I would appreciate any idea.
>
> Regards,
>
> Peri
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] multiple time series plots

[peri He]

Dear Friends,

I am trying to add two time series curves into one plot. I don't get any error 
but the lines are not added to my plot.
I am trying the following code;

1) df <- as.data.frame (read_excel("C:/Users/NO2.xlsx", "trend))

2) plot(df$year, df$Nation, type="o", col="blue", lty=1, ylab="change" )

3) lines(df$year, df$East, col="red",lty=2)

4) lines(df$year, df$center, col="dark red", lty=3)

I need to mention that Nation, East and Center are some probability values.
And year, is a sequence of years. Lines 1 and 2 are run properly but lines 3 
and 4 are not added.

I would appreciate any idea.

Regards,

Peri

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] how to overlay two histograms

[Rasmus Liland]

On 2020-09-17 17:01 -0500, Ana Marija wrote:
> Hello,
> 
> I am trying to overlay two histograms with this:
> 
> p <- ggplot(d, aes(CHR, counts, fill = name)) + geom_bar(position = "dodge")
> p
> 
> but I am getting this error:
> Error: stat_count() can only have an x or y aesthetic.
> Run `rlang::last_error()` to see where the error occurred.

Dear Ana,

you need to specify stat="identity" [1], 
like so:

mapping <- ggplot2::aes(
  x=CHR,
  y=counts,
  fill=name)
p <-
  ggplot2::ggplot() +
  ggplot2::geom_bar(
data=d,
mapping = mapping,
position="dodge",
stat="identity")

Best,
Rasmus

[1] 
https://stackoverflow.com/questions/61068031/error-stat-count-can-only-have-an-x-or-y-aesthetic


signature.asc
Description: PGP signature
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] how to overlay two histograms

[Jim Lemon]

Hi Ana,
Sorry it's not in ggplot, but it may help:

d<-read.table(text="CHR counts name
  1 193554  old
  2 220816  old
  3 174350  old
  4 163112  old
  5 168125  old
  6 182366  old
  7 143023  old
  8 147410  old
  9 122112  old
 10 138394  old
 11 130069  old
 12 124850  old
 13 104119  old
 14  83931  old
 15  72287  old
 16  71550  old
 17  58380  old
 18  76812  old
 19  37040  old
 20  63407  old
 21  33863  old
 22  33812  old
  1 202783  new
  2 252124  new
  3 213337  new
  4 201001  new
  5 207606  new
  6 228133  new
  7 147218  new
  8 177518  new
  9 121276  new
 10 163447  new
 11 158724  new
 12 142183  new
 13 89  new
 14  83043  new
 15  61063  new
 16  55439  new
 17  32883  new
 18  69135  new
 19  16624  new
 20  48541  new
 21  25479  new
 22  19698  new",
header=TRUE,stingsAsFactors=FALSE)
barpos<-barplot(counts~name+CHR,data=d,beside=TRUE,names.arg=rep("",22))
legend(40,22,c("new","old"),fill=c("gray20","gray80"))
library(plotrix)
staxlab(1,at=colMeans(barpos),labels=1:22)

Jim

On Fri, Sep 18, 2020 at 8:05 AM Ana Marija  wrote:
>
> Hello,
>
> I am trying to overlay two histograms with this:
>
> p <- ggplot(d, aes(CHR, counts, fill = name)) + geom_bar(position = "dodge")
> p
>
> but I am getting this error:
> Error: stat_count() can only have an x or y aesthetic.
> Run `rlang::last_error()` to see where the error occurred.
>
> my data is this:
>
> > d
>CHR counts name
> 11 193554  old
> 22 220816  old
> 33 174350  old
> 44 163112  old
> 55 168125  old
> 66 182366  old
> 77 143023  old
> 88 147410  old
> 99 122112  old
> 10  10 138394  old
> 11  11 130069  old
> 12  12 124850  old
> 13  13 104119  old
> 14  14  83931  old
> 15  15  72287  old
> 16  16  71550  old
> 17  17  58380  old
> 18  18  76812  old
> 19  19  37040  old
> 20  20  63407  old
> 21  21  33863  old
> 22  22  33812  old
> 23   1 202783  new
> 24   2 252124  new
> 25   3 213337  new
> 26   4 201001  new
> 27   5 207606  new
> 28   6 228133  new
> 29   7 147218  new
> 30   8 177518  new
> 31   9 121276  new
> 32  10 163447  new
> 33  11 158724  new
> 34  12 142183  new
> 35  13 89  new
> 36  14  83043  new
> 37  15  61063  new
> 38  16  55439  new
> 39  17  32883  new
> 40  18  69135  new
> 41  19  16624  new
> 42  20  48541  new
> 43  21  25479  new
> 44  22  19698  new
>
> Basically I need to show counts per CHR in "old" and "new" side by side.
>
> Please advise,
> Ana
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] how to overlay two histograms

[Ana Marija]

Hello,

I am trying to overlay two histograms with this:

p <- ggplot(d, aes(CHR, counts, fill = name)) + geom_bar(position = "dodge")
p

but I am getting this error:
Error: stat_count() can only have an x or y aesthetic.
Run `rlang::last_error()` to see where the error occurred.

my data is this:

> d
   CHR counts name
11 193554  old
22 220816  old
33 174350  old
44 163112  old
55 168125  old
66 182366  old
77 143023  old
88 147410  old
99 122112  old
10  10 138394  old
11  11 130069  old
12  12 124850  old
13  13 104119  old
14  14  83931  old
15  15  72287  old
16  16  71550  old
17  17  58380  old
18  18  76812  old
19  19  37040  old
20  20  63407  old
21  21  33863  old
22  22  33812  old
23   1 202783  new
24   2 252124  new
25   3 213337  new
26   4 201001  new
27   5 207606  new
28   6 228133  new
29   7 147218  new
30   8 177518  new
31   9 121276  new
32  10 163447  new
33  11 158724  new
34  12 142183  new
35  13 89  new
36  14  83043  new
37  15  61063  new
38  16  55439  new
39  17  32883  new
40  18  69135  new
41  19  16624  new
42  20  48541  new
43  21  25479  new
44  22  19698  new

Basically I need to show counts per CHR in "old" and "new" side by side.

Please advise,
Ana

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] linearHypothesis

[John Fox]

Dear Johan,

It's generally a good idea to keep the conversation on r-help to allow 
list members to follow it, and so I'm cc'ing this response to the list.


I hope that it's clear that car::linearHypothesis() computes the test as 
a Wald test of a linear hypothesis and not as a likelihood-ratio test by 
model comparison. As your example illustrates, however, the two tests 
are the same for a linear model, but this is not true more generally.


As I mentioned, you can find the details in many sources, including in 
Section 5.3.5 of Fox and Weisberg, An R Companion to Applied Regression, 
3rd Edition, the book with which the car package is associated.


Best,
 John

On 2020-09-17 4:03 p.m., Johan Lassen wrote:
Thank you John - highly appreciated! Yes, you are right, the less 
complex model may be seen as a restricted model of the starting model. 
Although the set of variables in the less complex model is not directly 
a subset of the variables of the starting model. What confused me at 
first was that I think of a subset model as a model having a direct 
subset of the set of variables of the starting model. Even though this 
is not the case in the example, the test still is on a restricted model 
of the starting model.

Thanks,
Johan

Den tor. 17. sep. 2020 kl. 15.55 skrev John Fox >:


Dear Johan,

On 2020-09-17 9:07 a.m., Johan Lassen wrote:
 > Dear R-users,
 >
 > I am using the R-function "linearHypothesis" to test if the sum
of all
 > parameters, but the intercept, in a multiple linear regression is
different
 > from zero.
 > I wonder if it is statistically valid to use the 
linearHypothesis-function

 > for this?

Yes, assuming of course that the hypothesis makes sense.


 > Below is a reproducible example in R. A multiple regression: y =
 > beta0*t0+beta1*t1+beta2*t2+beta3*t3+beta4*t4
 >
 > It seems to me that the linearHypothesis function does the
calculation as
 > an F-test on the extra residuals when going from the starting
model to a
 > 'subset' model, although all variables in the 'subset' model
differ from
 > the variables in the starting model.
 > I normally think of a subset model as a model built on the same
input data
 > as the starting model but one variable.
 >
 > Hence, is this a valid calculation?

First, linearHypothesis() doesn't literally fit alternative models, but
rather tests the linear hypothesis directly from the coefficient
estimates and their covariance matrix. The test is standard -- look at
the references in ?linearHypothesis or most texts on linear models.

Second, formulating the hypothesis using alternative models is also
legitimate, since the second model is a restricted version of the first.

 >
 > Thanks in advance,Johan
 >
 > # R-code:
 > y <-
 >

c(101133190,96663050,106866486,97678429,83212348,75719714,77861937,74018478,82181104,68667176,64599495,62414401,63534709,58571865,65222727,60139788,
 >

63355011,57790610,55214971,55535484,55759192,49450719,48834699,51383864,51250871,50629835,52154608,54636478,54942637)
 >
 > data <-
 >

data.frame(y,"t0"=1,"t1"=1990:2018,"t2"=c(rep(0,12),1:17),"t3"=c(rep(0,17),1:12),"t4"=c(rep(0,23),1:6))
 >
 > model <- lm(y~t0+t1+t2+t3+t4+0,data=data)

You need not supply the constant regressor t0 explicitly and suppress
the intercept -- you'd get the same test from linearHypothesis() for
lm(y~t1+t2+t3+t4,data=data).

 >
 > linearHypothesis(model,"t1+t2+t3+t4=0",test=c("F"))

test = "F" is the default.

 >
 > # Reproduce the result from linearHypothesis:
 > # beta1+beta2+beta3+beta4=0 -> beta4=-(beta1+beta2+beta3) ->
 > # y=beta0+beta1*t1+beta2*t2+beta3*t3-(beta1+beta2+beta3)*t4
 > # y = beta0'+beta1'*(t1-t4)+beta2'*(t2-t4)+beta3'*(t3-t4)
 >
 > data$t1 <- data$t1-data$t4
 > data$t2 <- data$t2-data$t4
 > data$t3 <- data$t3-data$t4
 >
 > model_reduced <- lm(y~t0+t1+t2+t3+0,data=data)
 >
 > anova(model_reduced,model)

Yes, this is equivalent to the test performed by linearHypothesis()
using the coefficients and their covariances from the original model.

I hope this helps,
   John

-- 
John Fox, Professor Emeritus

McMaster University
Hamilton, Ontario, Canada
web: https://socialsciences.mcmaster.ca/jfox/
 >



--
Johan Lassen

"In the cities people live in time -
in the mountains people live in space" (Budistisk munk).


--
John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
web: https://socialsciences.mcmaster.ca/jfox/

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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

Re: [R] Seasonal adjustment with seas - Error: X-13 has returned a non-zero exist status

[Jeff Newmiller]

I don't know definitively, but if I were in your shoes I would put my personal 
package library on a local drive. UNC paths are not supported by all of the 
Windows file-related system calls.

On September 17, 2020 6:59:33 AM PDT, "Weber Sylvain (HES)" 
 wrote:
>Dear all,
>
>I am using function "seas" from package "seasonal" for seasonal
>adjustment of time series. I am using the latest version of R (4.0.2,
>64-bit).
>The following lines provide a small working example:
>## EXAMPLE ##
>library(seasonal)
>m <- seas(AirPassengers)
>m
>
>On my private laptop, everything is going smoothly, and I obtain the
>following output:
>## OUTPUT ##
>Call:
>seas(x = AirPassengers)
>
>Coefficients:
> Weekday  Easter[1] AO1951.May  MA-Nonseasonal-01  
>-0.002950.017770.100160.11562  
>   MA-Seasonal-12  
>  0.49736  
>
>But on my office computer, I receive the following error message after
>issuing " m <- seas(AirPassengers)":
>## ERROR ##
>Error: X-13 has returned a non-zero exist status, which means that the
>current spec file cannot be processed for an unknown reason.
>In addition: Warning message:
>In system(cmd, intern = intern, wait = wait | intern,
>show.output.on.console = wait,  :
>running command 'C:\Windows\system32\cmd.exe /c
>"\\hes-nas-prairie.hes.adhes.hesge.ch/Home_S/sylvain.weber/Documents/R/win-library/4.0/x13binary/bin/x13ashtml.exe"
>C:\Users\SYLVAI~1.WEB\AppData\Local\Temp\Rtmpe0MDE7\x132d0cb8579b9/iofile
>-n -s' had status 1
>
>A similar issue has been already mentioned in various blogs, for
>instance:
>https://rdrr.io/github/christophsax/seasonal/src/travis/test-x13messages.R
>http://freerangestats.info/blog/2015/12/21/m3-and-x13
>https://github.com/christophsax/seasonal/issues/154
>but I couldn't find any solution yet.
>
>Does anyone have any idea why this issue is occurring?
>
>Thanks very much.
>
>Sylvain
>
>__
>R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>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.

-- 
Sent from my phone. Please excuse my brevity.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Seasonal adjustment with seas - Error: X-13 has returned a non-zero exist status

[Jeff Newmiller]

I don't know definitively, but if I were in your shoes I would put my personal 
package library on a local drive. UNC paths are not supported by all of the 
Windows file-related system calls.

On September 17, 2020 6:59:33 AM PDT, "Weber Sylvain (HES)" 
 wrote:
>Dear all,
>
>I am using function "seas" from package "seasonal" for seasonal
>adjustment of time series. I am using the latest version of R (4.0.2,
>64-bit).
>The following lines provide a small working example:
>## EXAMPLE ##
>library(seasonal)
>m <- seas(AirPassengers)
>m
>
>On my private laptop, everything is going smoothly, and I obtain the
>following output:
>## OUTPUT ##
>Call:
>seas(x = AirPassengers)
>
>Coefficients:
> Weekday  Easter[1] AO1951.May  MA-Nonseasonal-01  
>-0.002950.017770.100160.11562  
>   MA-Seasonal-12  
>  0.49736  
>
>But on my office computer, I receive the following error message after
>issuing " m <- seas(AirPassengers)":
>## ERROR ##
>Error: X-13 has returned a non-zero exist status, which means that the
>current spec file cannot be processed for an unknown reason.
>In addition: Warning message:
>In system(cmd, intern = intern, wait = wait | intern,
>show.output.on.console = wait,  :
>running command 'C:\Windows\system32\cmd.exe /c
>"\\hes-nas-prairie.hes.adhes.hesge.ch/Home_S/sylvain.weber/Documents/R/win-library/4.0/x13binary/bin/x13ashtml.exe"
>C:\Users\SYLVAI~1.WEB\AppData\Local\Temp\Rtmpe0MDE7\x132d0cb8579b9/iofile
>-n -s' had status 1
>
>A similar issue has been already mentioned in various blogs, for
>instance:
>https://rdrr.io/github/christophsax/seasonal/src/travis/test-x13messages.R
>http://freerangestats.info/blog/2015/12/21/m3-and-x13
>https://github.com/christophsax/seasonal/issues/154
>but I couldn't find any solution yet.
>
>Does anyone have any idea why this issue is occurring?
>
>Thanks very much.
>
>Sylvain
>
>__
>R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>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.

-- 
Sent from my phone. Please excuse my brevity.
-- 
Sent from my phone. Please excuse my brevity.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Seasonal adjustment with seas - Error: X-13 has returned a non-zero exist status

[Bill Dunlap]

The problem might to due to using an UNC path (//machine//blah/...) instead
of the traditional DOS path ("H:/blah/...").

E.g., if my working directory has a UNC path, cmd.exe will not work as
expected:
> getwd()
[1] "server/dept/devel/bill-sandbox"
> system("C:\\WINDOWS\\SYSTEM32\\cmd.exe /c echo foo bar", intern=TRUE)
[1] "'server\\dept\\devel\\bill-sandbox'"
[2] "CMD.EXE was started with the above path as the current directory."
[3] "UNC paths are not supported.  Defaulting to Windows directory."
[4] "foo bar"

You can test that by mapping "\\hes-nas-prairie.hes.adhes.hesge.ch/Home_S
"
to a driver letter and using letter-colon everywhere instead of "\\...".

-Bill

On Thu, Sep 17, 2020 at 8:46 AM Weber Sylvain (HES) 
wrote:

> Dear all,
>
> I am using function "seas" from package "seasonal" for seasonal adjustment
> of time series. I am using the latest version of R (4.0.2, 64-bit).
> The following lines provide a small working example:
> ## EXAMPLE ##
> library(seasonal)
> m <- seas(AirPassengers)
> m
>
> On my private laptop, everything is going smoothly, and I obtain the
> following output:
> ## OUTPUT ##
> Call:
> seas(x = AirPassengers)
>
> Coefficients:
>   Weekday  Easter[1] AO1951.May
> MA-Nonseasonal-01
>  -0.002950.017770.10016
> 0.11562
>MA-Seasonal-12
>   0.49736
>
> But on my office computer, I receive the following error message after
> issuing " m <- seas(AirPassengers)":
> ## ERROR ##
> Error: X-13 has returned a non-zero exist status, which means that the
> current spec file cannot be processed for an unknown reason.
> In addition: Warning message:
> In system(cmd, intern = intern, wait = wait | intern,
> show.output.on.console = wait,  :
>   running command 'C:\Windows\system32\cmd.exe /c "\\
> hes-nas-prairie.hes.adhes.hesge.ch/Home_S/sylvain.weber/Documents/R/win-library/4.0/x13binary/bin/x13ashtml.exe"
> C:\Users\SYLVAI~1.WEB\AppData\Local\Temp\Rtmpe0MDE7\x132d0cb8579b9/iofile
> -n -s' had status 1
>
> A similar issue has been already mentioned in various blogs, for instance:
> https://rdrr.io/github/christophsax/seasonal/src/travis/test-x13messages.R
> http://freerangestats.info/blog/2015/12/21/m3-and-x13
> https://github.com/christophsax/seasonal/issues/154
> but I couldn't find any solution yet.
>
> Does anyone have any idea why this issue is occurring?
>
> Thanks very much.
>
> Sylvain
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Seasonal adjustment with seas - Error: X-13 has returned a non-zero exist status

[Bert Gunter]

Consider contacting the package maintainer, found by:
maintainer("seasonal") for such a package specific question.

Cheers,

Bert Gunter

"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Thu, Sep 17, 2020 at 8:46 AM Weber Sylvain (HES) 
wrote:

> Dear all,
>
> I am using function "seas" from package "seasonal" for seasonal adjustment
> of time series. I am using the latest version of R (4.0.2, 64-bit).
> The following lines provide a small working example:
> ## EXAMPLE ##
> library(seasonal)
> m <- seas(AirPassengers)
> m
>
> On my private laptop, everything is going smoothly, and I obtain the
> following output:
> ## OUTPUT ##
> Call:
> seas(x = AirPassengers)
>
> Coefficients:
>   Weekday  Easter[1] AO1951.May
> MA-Nonseasonal-01
>  -0.002950.017770.10016
> 0.11562
>MA-Seasonal-12
>   0.49736
>
> But on my office computer, I receive the following error message after
> issuing " m <- seas(AirPassengers)":
> ## ERROR ##
> Error: X-13 has returned a non-zero exist status, which means that the
> current spec file cannot be processed for an unknown reason.
> In addition: Warning message:
> In system(cmd, intern = intern, wait = wait | intern,
> show.output.on.console = wait,  :
>   running command 'C:\Windows\system32\cmd.exe /c "\\
> hes-nas-prairie.hes.adhes.hesge.ch/Home_S/sylvain.weber/Documents/R/win-library/4.0/x13binary/bin/x13ashtml.exe"
> C:\Users\SYLVAI~1.WEB\AppData\Local\Temp\Rtmpe0MDE7\x132d0cb8579b9/iofile
> -n -s' had status 1
>
> A similar issue has been already mentioned in various blogs, for instance:
> https://rdrr.io/github/christophsax/seasonal/src/travis/test-x13messages.R
> http://freerangestats.info/blog/2015/12/21/m3-and-x13
> https://github.com/christophsax/seasonal/issues/154
> but I couldn't find any solution yet.
>
> Does anyone have any idea why this issue is occurring?
>
> Thanks very much.
>
> Sylvain
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] Seasonal adjustment with seas - Error: X-13 has returned a non-zero exist status

[Weber Sylvain (HES)]

Dear all,

I am using function "seas" from package "seasonal" for seasonal adjustment of 
time series. I am using the latest version of R (4.0.2, 64-bit).
The following lines provide a small working example:
## EXAMPLE ##
library(seasonal)
m <- seas(AirPassengers)
m

On my private laptop, everything is going smoothly, and I obtain the following 
output:
## OUTPUT ##
Call:
seas(x = AirPassengers)

Coefficients:
  Weekday  Easter[1] AO1951.May  MA-Nonseasonal-01  
 -0.002950.017770.100160.11562  
   MA-Seasonal-12  
  0.49736  

But on my office computer, I receive the following error message after issuing 
" m <- seas(AirPassengers)":
## ERROR ##
Error: X-13 has returned a non-zero exist status, which means that the current 
spec file cannot be processed for an unknown reason.
In addition: Warning message:
In system(cmd, intern = intern, wait = wait | intern, show.output.on.console = 
wait,  :
  running command 'C:\Windows\system32\cmd.exe /c 
"\\hes-nas-prairie.hes.adhes.hesge.ch/Home_S/sylvain.weber/Documents/R/win-library/4.0/x13binary/bin/x13ashtml.exe"
 C:\Users\SYLVAI~1.WEB\AppData\Local\Temp\Rtmpe0MDE7\x132d0cb8579b9/iofile -n 
-s' had status 1

A similar issue has been already mentioned in various blogs, for instance:
https://rdrr.io/github/christophsax/seasonal/src/travis/test-x13messages.R
http://freerangestats.info/blog/2015/12/21/m3-and-x13
https://github.com/christophsax/seasonal/issues/154
but I couldn't find any solution yet.

Does anyone have any idea why this issue is occurring?

Thanks very much.

Sylvain

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] Unnesting JSON using R

[William Michels via R-help]

Hi Fred, I believe the preferred package is jsonlite:

https://cran.r-project.org/package=jsonlite
https://jeroen.cran.dev/jsonlite/index.html

HTH, Bill.

W. Michels, Ph.D.

On Tue, Sep 15, 2020 at 1:48 PM Fred Kwebiha  wrote:
>
> Source=https://jsonformatter.org/e038ec
>
> The above is nested json.
>
> I want the output to be as below
> dataElements.name,dataElements.id,categoryOptionCombos.name,categoryOptionCombos.id
>
> Any help in r?
> *Best Regards,*
>
> *FRED KWEBIHA*
> *+256-782-746-154*
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] linearHypothesis

[John Fox]

Dear Johan,

On 2020-09-17 9:07 a.m., Johan Lassen wrote:

Dear R-users,

I am using the R-function "linearHypothesis" to test if the sum of all
parameters, but the intercept, in a multiple linear regression is different
from zero.
I wonder if it is statistically valid to use the  linearHypothesis-function
for this?


Yes, assuming of course that the hypothesis makes sense.



Below is a reproducible example in R. A multiple regression: y =
beta0*t0+beta1*t1+beta2*t2+beta3*t3+beta4*t4

It seems to me that the linearHypothesis function does the calculation as
an F-test on the extra residuals when going from the starting model to a
'subset' model, although all variables in the 'subset' model differ from
the variables in the starting model.
I normally think of a subset model as a model built on the same input data
as the starting model but one variable.

Hence, is this a valid calculation?


First, linearHypothesis() doesn't literally fit alternative models, but 
rather tests the linear hypothesis directly from the coefficient 
estimates and their covariance matrix. The test is standard -- look at 
the references in ?linearHypothesis or most texts on linear models.


Second, formulating the hypothesis using alternative models is also 
legitimate, since the second model is a restricted version of the first.




Thanks in advance,Johan

# R-code:
y <-
c(101133190,96663050,106866486,97678429,83212348,75719714,77861937,74018478,82181104,68667176,64599495,62414401,63534709,58571865,65222727,60139788,
63355011,57790610,55214971,55535484,55759192,49450719,48834699,51383864,51250871,50629835,52154608,54636478,54942637)

data <-
data.frame(y,"t0"=1,"t1"=1990:2018,"t2"=c(rep(0,12),1:17),"t3"=c(rep(0,17),1:12),"t4"=c(rep(0,23),1:6))

model <- lm(y~t0+t1+t2+t3+t4+0,data=data)


You need not supply the constant regressor t0 explicitly and suppress 
the intercept -- you'd get the same test from linearHypothesis() for 
lm(y~t1+t2+t3+t4,data=data).




linearHypothesis(model,"t1+t2+t3+t4=0",test=c("F"))


test = "F" is the default.



# Reproduce the result from linearHypothesis:
# beta1+beta2+beta3+beta4=0 -> beta4=-(beta1+beta2+beta3) ->
# y=beta0+beta1*t1+beta2*t2+beta3*t3-(beta1+beta2+beta3)*t4
# y = beta0'+beta1'*(t1-t4)+beta2'*(t2-t4)+beta3'*(t3-t4)

data$t1 <- data$t1-data$t4
data$t2 <- data$t2-data$t4
data$t3 <- data$t3-data$t4

model_reduced <- lm(y~t0+t1+t2+t3+0,data=data)

anova(model_reduced,model)


Yes, this is equivalent to the test performed by linearHypothesis() 
using the coefficients and their covariances from the original model.


I hope this helps,
 John

--
John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
web: https://socialsciences.mcmaster.ca/jfox/




__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[R] linearHypothesis

[Johan Lassen]

Dear R-users,

I am using the R-function "linearHypothesis" to test if the sum of all
parameters, but the intercept, in a multiple linear regression is different
from zero.
I wonder if it is statistically valid to use the  linearHypothesis-function
for this?
Below is a reproducible example in R. A multiple regression: y =
beta0*t0+beta1*t1+beta2*t2+beta3*t3+beta4*t4

It seems to me that the linearHypothesis function does the calculation as
an F-test on the extra residuals when going from the starting model to a
'subset' model, although all variables in the 'subset' model differ from
the variables in the starting model.
I normally think of a subset model as a model built on the same input data
as the starting model but one variable.

Hence, is this a valid calculation?

Thanks in advance,Johan

# R-code:
y <-
c(101133190,96663050,106866486,97678429,83212348,75719714,77861937,74018478,82181104,68667176,64599495,62414401,63534709,58571865,65222727,60139788,
63355011,57790610,55214971,55535484,55759192,49450719,48834699,51383864,51250871,50629835,52154608,54636478,54942637)

data <-
data.frame(y,"t0"=1,"t1"=1990:2018,"t2"=c(rep(0,12),1:17),"t3"=c(rep(0,17),1:12),"t4"=c(rep(0,23),1:6))

model <- lm(y~t0+t1+t2+t3+t4+0,data=data)

linearHypothesis(model,"t1+t2+t3+t4=0",test=c("F"))

# Reproduce the result from linearHypothesis:
# beta1+beta2+beta3+beta4=0 -> beta4=-(beta1+beta2+beta3) ->
# y=beta0+beta1*t1+beta2*t2+beta3*t3-(beta1+beta2+beta3)*t4
# y = beta0'+beta1'*(t1-t4)+beta2'*(t2-t4)+beta3'*(t3-t4)

data$t1 <- data$t1-data$t4
data$t2 <- data$t2-data$t4
data$t3 <- data$t3-data$t4

model_reduced <- lm(y~t0+t1+t2+t3+0,data=data)

anova(model_reduced,model)

-- 
Johan Lassen

"In the cities people live in time -
in the mountains people live in space" (Budistisk munk).

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Re: [R] How to reduce the sparseness in a TDM to make a cluster plot readable?

[Abby Spurdle]

I'm not familiar with these subjects.
And hopefully, someone who is, will offer some better suggestions.

But to get things started, maybe...
(1) What packages are you using (re: tdm)?
(2) Where does the problem happen, in dist, hclust, the plot method
for hclust, or in the package(s) you are using?
(3) Do you think you could produce a small reproducible example,
showing what is wrong, and explaining you would like it to do instead?

Note that if the problem relates to hclust, or the plot method, then
you should be able to produce a much simpler example.
e.g.

mycount.matrix <- matrix (rpois (25000, 20),, 5)
head (mycount.matrix, 3)
tail (mycount.matrix, 3)

plot (hclust (dist (mycount.matrix) ) )

On Tue, Sep 15, 2020 at 6:54 AM Andrew  wrote:
>
> Hello all
>
> I am doing some text mining on a set of five plain text files and have
> run into a snag when I run hclust in that there are just too many leaves
> for anything to be read. It returns a solid black line.
>
> The texts have been converted into a TDM which has a dim of 5,292 and 5
> (as per 5 docs).
>
> My code for removing sparsity is as follows:
>
>  > tdm2 <- removeSparseTerms(tdm, sparse=0.9)
>
>  > inspect(tdm2)
>
> <>
> Non-/sparse entries: 10415/16045
> Sparsity   : 61%
> Maximal term length: 22
> Weighting  : term frequency (tf)
>
> While the tf-idf weighting returns this when 0.9 sparseness is removed:
>
>  > inspect(tdm.tfidf)
> <>
> Non-/sparse entries: 7915/18545
> Sparsity   : 70%
> Maximal term length: 22
> Weighting  : term frequency - inverse document frequency
> (normalized) (tf-idf)
>
> I have experimented by decreasing the value I use for decreasing
> sparseness, and that helps a bit, for example:
>
>  > tdm2 <- removeSparseTerms(tdm, sparse=0.215)
>  > inspect(tdm2)
> <>
> Non-/sparse entries: 3976/369
> Sparsity   : 8%
> Maximal term length: 14
> Weighting  : term frequency (tf)
>
> But, no matter what I do, the resulting plot is unreadable. The code for
> plotting the cluster is:
>
>  > hc <- hclust(dist(tdm2, method = "euclidean"), method = "complete")
>  > plot(hc, yaxt = 'n', main = "Hierarchical clustering")
>
> Can someone kindly either advise me what I am doing wrong and/ or
> signpost me to some detailed info on how to fix this.
>
> Many thanks in anticipation.
>
> Andy
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.