[R] The argument 'eps.Pvalue' of `printCoefmat()`

2023-10-29 Thread Shu Fai Cheung
Hi all,

Just a minor issue that I am not sure whether this is considered a
"bug." It is about the help page.

In the help page of printCoefmat(), for the argument 'eps.Pvalue', the
description is as below:

number, ..

I have to read the source to figure out that this argument is to be
used by format.pval().

Maybe the description of 'eps.Pvalue' can be revised to refer users to
the help page of format.pval()?

Regards,
Shu Fai

__
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] Create a call but evaluate only some elements

2023-10-25 Thread Shu Fai Cheung
Dear Bert,

Many thanks for your suggestion! I am reading the section to
understand more about this topic. It is highly relevant to what I plan
to work on.

Regards,
Shu Fai

On Thu, Oct 26, 2023 at 5:38 AM Bert Gunter  wrote:
>
> As you seem to have a need for this sort of capability (e.g. bquote),
> see Section 6: "Computing on the Language" in the R Language
> Definition manual. Actually, if you are interested in a concise
> (albeit dense) overview of the R Language, you might consider going
> through the whole manual.
>
> Cheers,
> Bert
>
> On Wed, Oct 25, 2023 at 3:57 AM Shu Fai Cheung  
> wrote:
> >
> > Dear Iris,
> >
> > Many many thanks! This is exactly what I need! I have never heard
> > about bquote(). This function will also be useful to me on other
> > occasions.
> >
> > I still have a lot to learn about the R language ...
> >
> > Regards,
> > Shu Fai
> >
> >
> > On Wed, Oct 25, 2023 at 5:24 PM Iris Simmons  wrote:
> > >
> > > You can try either of these:
> > >
> > > expr <- bquote(lm(.(as.formula(mod)), dat))
> > > lm_out5 <- eval(expr)
> > >
> > > expr <- call("lm", as.formula(mod), as.symbol("dat"))
> > > lm_out6 <- eval(expr)
> > >
> > > but bquote is usually easier and good enough.
> > >
> > > On Wed, Oct 25, 2023, 05:10 Shu Fai Cheung  
> > > wrote:
> > >>
> > >> Hi All,
> > >>
> > >> I have a problem that may have a simple solution, but I am not
> > >> familiar with creating calls manually.
> > >>
> > >> This is example calling lm()
> > >>
> > >> ``` r
> > >> set.seed(1234)
> > >> n <- 10
> > >> dat <- data.frame(x1 = rnorm(n),
> > >>   x2 = rnorm(n),
> > >>   y = rnorm(n))
> > >>
> > >> lm_out <- lm(y ~ x1 + x2, dat)
> > >> lm_out
> > >> #>
> > >> #> Call:
> > >> #> lm(formula = y ~ x1 + x2, data = dat)
> > >> #>
> > >> #> Coefficients:
> > >> #> (Intercept)   x1   x2
> > >> #> -0.5755  -0.4151  -0.2411
> > >> lm_out$call
> > >> #> lm(formula = y ~ x1 + x2, data = dat)
> > >> ```
> > >>
> > >> The call is stored, "lm(formula = y ~ x1 + x2, data = dat)", and names
> > >> are not evaluated.
> > >>
> > >> I want to create a similar call, but only one of the elements is from a 
> > >> string.
> > >>
> > >> ```r
> > >> mod <- "y ~ x1 + x2"
> > >> ```
> > >>
> > >> This is what I tried but failed:
> > >>
> > >> ```r
> > >> lm_out2 <- do.call("lm",
> > >>list(formula = as.formula(mod),
> > >> data = dat))
> > >> lm_out2
> > >> #>
> > >> #> Call:
> > >> #> lm(formula = y ~ x1 + x2, data = structure(list(x1 = 
> > >> c(-1.20706574938542,
> > >> #> 0.27742924211066, 1.08444117668306, -2.34569770262935, 
> > >> 0.42912468881105,
> > >> #> 0.506055892157574, -0.574739960134649, -0.546631855784187,
> > >> -0.564451999093283,
> > >> #> -0.890037829044104), x2 = c(-0.477192699753547, -0.998386444859704,
> > >> #> -0.77625389463799, 0.0644588172762693, 0.959494058970771, 
> > >> -0.110285494390774,
> > >> #> -0.511009505806642, -0.911195416629811, -0.83717168026894, 
> > >> 2.41583517848934
> > >> #> ), y = c(0.134088220152031, -0.490685896690943, -0.440547872353227,
> > >> #> 0.459589441005854, -0.693720246937475, -1.44820491038647, 
> > >> 0.574755720900728,
> > >> #> -1.02365572296388, -0.0151383003641817, -0.935948601168394)), class
> > >> = "data.frame", row.names = c(NA,
> > >> #> -10L)))
> > >> #>
> > >> #> Coefficients:
> > >> #> (Intercept)   x1   x2
> > >> #> -0.5755  -0.4151  -0.2411
> > >> ```
> > >>
> > >> It does not have the formula, "as a formula": y ~ x1 + x2.
> > >> However, the name "dat" is evaluated. Therefore, the call stored does
> > >> not have the name 'dat',

Re: [R] Create a call but evaluate only some elements

2023-10-25 Thread Shu Fai Cheung
Dear Iris,

Many many thanks! This is exactly what I need! I have never heard
about bquote(). This function will also be useful to me on other
occasions.

I still have a lot to learn about the R language ...

Regards,
Shu Fai


On Wed, Oct 25, 2023 at 5:24 PM Iris Simmons  wrote:
>
> You can try either of these:
>
> expr <- bquote(lm(.(as.formula(mod)), dat))
> lm_out5 <- eval(expr)
>
> expr <- call("lm", as.formula(mod), as.symbol("dat"))
> lm_out6 <- eval(expr)
>
> but bquote is usually easier and good enough.
>
> On Wed, Oct 25, 2023, 05:10 Shu Fai Cheung  wrote:
>>
>> Hi All,
>>
>> I have a problem that may have a simple solution, but I am not
>> familiar with creating calls manually.
>>
>> This is example calling lm()
>>
>> ``` r
>> set.seed(1234)
>> n <- 10
>> dat <- data.frame(x1 = rnorm(n),
>>   x2 = rnorm(n),
>>   y = rnorm(n))
>>
>> lm_out <- lm(y ~ x1 + x2, dat)
>> lm_out
>> #>
>> #> Call:
>> #> lm(formula = y ~ x1 + x2, data = dat)
>> #>
>> #> Coefficients:
>> #> (Intercept)   x1   x2
>> #> -0.5755  -0.4151  -0.2411
>> lm_out$call
>> #> lm(formula = y ~ x1 + x2, data = dat)
>> ```
>>
>> The call is stored, "lm(formula = y ~ x1 + x2, data = dat)", and names
>> are not evaluated.
>>
>> I want to create a similar call, but only one of the elements is from a 
>> string.
>>
>> ```r
>> mod <- "y ~ x1 + x2"
>> ```
>>
>> This is what I tried but failed:
>>
>> ```r
>> lm_out2 <- do.call("lm",
>>list(formula = as.formula(mod),
>> data = dat))
>> lm_out2
>> #>
>> #> Call:
>> #> lm(formula = y ~ x1 + x2, data = structure(list(x1 = c(-1.20706574938542,
>> #> 0.27742924211066, 1.08444117668306, -2.34569770262935, 0.42912468881105,
>> #> 0.506055892157574, -0.574739960134649, -0.546631855784187,
>> -0.564451999093283,
>> #> -0.890037829044104), x2 = c(-0.477192699753547, -0.998386444859704,
>> #> -0.77625389463799, 0.0644588172762693, 0.959494058970771, 
>> -0.110285494390774,
>> #> -0.511009505806642, -0.911195416629811, -0.83717168026894, 
>> 2.41583517848934
>> #> ), y = c(0.134088220152031, -0.490685896690943, -0.440547872353227,
>> #> 0.459589441005854, -0.693720246937475, -1.44820491038647, 
>> 0.574755720900728,
>> #> -1.02365572296388, -0.0151383003641817, -0.935948601168394)), class
>> = "data.frame", row.names = c(NA,
>> #> -10L)))
>> #>
>> #> Coefficients:
>> #> (Intercept)   x1   x2
>> #> -0.5755  -0.4151  -0.2411
>> ```
>>
>> It does not have the formula, "as a formula": y ~ x1 + x2.
>> However, the name "dat" is evaluated. Therefore, the call stored does
>> not have the name 'dat', but has the evaluated content.
>>
>> The following fits the same model. However, the call stores the name,
>> 'mod', not the evaluated result, y ~ x1 + x2.
>>
>> ```r
>> lm_out3 <- lm(mod, data = dat)
>> lm_out3
>> #>
>> #> Call:
>> #> lm(formula = mod, data = dat)
>> #>
>> #> Coefficients:
>> #> (Intercept)   x1   x2
>> #> -0.5755  -0.4151  -0.2411
>> ```
>>
>> The following method works. However, I have to do a dummy call,
>> extract the stored call, and set formula to the result of
>> as.formula(mod):
>>
>> ```r
>> lm_out3 <- lm(mod, data = dat)
>> lm_out3
>> #>
>> #> Call:
>> #> lm(formula = mod, data = dat)
>> #>
>> #> Coefficients:
>> #> (Intercept)   x1   x2
>> #> -0.5755  -0.4151  -0.2411
>>
>> call1 <- lm_out3$call
>> call1$formula <- as.formula(mod)
>> lm_out4 <- eval(call1)
>> lm_out4
>> #>
>> #> Call:
>> #> lm(formula = y ~ x1 + x2, data = dat)
>> #>
>> #> Coefficients:
>> #> (Intercept)   x1   x2
>> #> -0.5755  -0.4151  -0.2411
>> ```
>>
>> Is it possible to create the call directly, with only 'mod' evaluated,
>> and other arguments, e.g., 'dat', not evaluated?
>>
>> Regards,
>> Shu Fai
>>
>> __
>> 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] Create a call but evaluate only some elements

2023-10-25 Thread Shu Fai Cheung
Sorry for a typo, regarding the first attempt, lm_out2, using
do.call(), I meant:

'It does have the formula, "as a formula": y ~ x1 + x2.
However, the name "dat" is evaluated. ...'

Regards,
Shu Fai

On Wed, Oct 25, 2023 at 5:09 PM Shu Fai Cheung  wrote:
>
> Hi All,
>
> I have a problem that may have a simple solution, but I am not
> familiar with creating calls manually.
>
> This is example calling lm()
>
> ``` r
> set.seed(1234)
> n <- 10
> dat <- data.frame(x1 = rnorm(n),
>   x2 = rnorm(n),
>   y = rnorm(n))
>
> lm_out <- lm(y ~ x1 + x2, dat)
> lm_out
> #>
> #> Call:
> #> lm(formula = y ~ x1 + x2, data = dat)
> #>
> #> Coefficients:
> #> (Intercept)   x1   x2
> #> -0.5755  -0.4151  -0.2411
> lm_out$call
> #> lm(formula = y ~ x1 + x2, data = dat)
> ```
>
> The call is stored, "lm(formula = y ~ x1 + x2, data = dat)", and names
> are not evaluated.
>
> I want to create a similar call, but only one of the elements is from a 
> string.
>
> ```r
> mod <- "y ~ x1 + x2"
> ```
>
> This is what I tried but failed:
>
> ```r
> lm_out2 <- do.call("lm",
>list(formula = as.formula(mod),
> data = dat))
> lm_out2
> #>
> #> Call:
> #> lm(formula = y ~ x1 + x2, data = structure(list(x1 = c(-1.20706574938542,
> #> 0.27742924211066, 1.08444117668306, -2.34569770262935, 0.42912468881105,
> #> 0.506055892157574, -0.574739960134649, -0.546631855784187,
> -0.564451999093283,
> #> -0.890037829044104), x2 = c(-0.477192699753547, -0.998386444859704,
> #> -0.77625389463799, 0.0644588172762693, 0.959494058970771, 
> -0.110285494390774,
> #> -0.511009505806642, -0.911195416629811, -0.83717168026894, 2.41583517848934
> #> ), y = c(0.134088220152031, -0.490685896690943, -0.440547872353227,
> #> 0.459589441005854, -0.693720246937475, -1.44820491038647, 
> 0.574755720900728,
> #> -1.02365572296388, -0.0151383003641817, -0.935948601168394)), class
> = "data.frame", row.names = c(NA,
> #> -10L)))
> #>
> #> Coefficients:
> #> (Intercept)   x1   x2
> #> -0.5755  -0.4151  -0.2411
> ```
>
> It does not have the formula, "as a formula": y ~ x1 + x2.
> However, the name "dat" is evaluated. Therefore, the call stored does
> not have the name 'dat', but has the evaluated content.
>
> The following fits the same model. However, the call stores the name,
> 'mod', not the evaluated result, y ~ x1 + x2.
>
> ```r
> lm_out3 <- lm(mod, data = dat)
> lm_out3
> #>
> #> Call:
> #> lm(formula = mod, data = dat)
> #>
> #> Coefficients:
> #> (Intercept)   x1   x2
> #> -0.5755  -0.4151  -0.2411
> ```
>
> The following method works. However, I have to do a dummy call,
> extract the stored call, and set formula to the result of
> as.formula(mod):
>
> ```r
> lm_out3 <- lm(mod, data = dat)
> lm_out3
> #>
> #> Call:
> #> lm(formula = mod, data = dat)
> #>
> #> Coefficients:
> #> (Intercept)   x1   x2
> #> -0.5755  -0.4151  -0.2411
>
> call1 <- lm_out3$call
> call1$formula <- as.formula(mod)
> lm_out4 <- eval(call1)
> lm_out4
> #>
> #> Call:
> #> lm(formula = y ~ x1 + x2, data = dat)
> #>
> #> Coefficients:
> #> (Intercept)   x1   x2
> #> -0.5755  -0.4151  -0.2411
> ```
>
> Is it possible to create the call directly, with only 'mod' evaluated,
> and other arguments, e.g., 'dat', not evaluated?
>
> Regards,
> Shu Fai

__
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] Create a call but evaluate only some elements

2023-10-25 Thread Shu Fai Cheung
Hi All,

I have a problem that may have a simple solution, but I am not
familiar with creating calls manually.

This is example calling lm()

``` r
set.seed(1234)
n <- 10
dat <- data.frame(x1 = rnorm(n),
  x2 = rnorm(n),
  y = rnorm(n))

lm_out <- lm(y ~ x1 + x2, dat)
lm_out
#>
#> Call:
#> lm(formula = y ~ x1 + x2, data = dat)
#>
#> Coefficients:
#> (Intercept)   x1   x2
#> -0.5755  -0.4151  -0.2411
lm_out$call
#> lm(formula = y ~ x1 + x2, data = dat)
```

The call is stored, "lm(formula = y ~ x1 + x2, data = dat)", and names
are not evaluated.

I want to create a similar call, but only one of the elements is from a string.

```r
mod <- "y ~ x1 + x2"
```

This is what I tried but failed:

```r
lm_out2 <- do.call("lm",
   list(formula = as.formula(mod),
data = dat))
lm_out2
#>
#> Call:
#> lm(formula = y ~ x1 + x2, data = structure(list(x1 = c(-1.20706574938542,
#> 0.27742924211066, 1.08444117668306, -2.34569770262935, 0.42912468881105,
#> 0.506055892157574, -0.574739960134649, -0.546631855784187,
-0.564451999093283,
#> -0.890037829044104), x2 = c(-0.477192699753547, -0.998386444859704,
#> -0.77625389463799, 0.0644588172762693, 0.959494058970771, -0.110285494390774,
#> -0.511009505806642, -0.911195416629811, -0.83717168026894, 2.41583517848934
#> ), y = c(0.134088220152031, -0.490685896690943, -0.440547872353227,
#> 0.459589441005854, -0.693720246937475, -1.44820491038647, 0.574755720900728,
#> -1.02365572296388, -0.0151383003641817, -0.935948601168394)), class
= "data.frame", row.names = c(NA,
#> -10L)))
#>
#> Coefficients:
#> (Intercept)   x1   x2
#> -0.5755  -0.4151  -0.2411
```

It does not have the formula, "as a formula": y ~ x1 + x2.
However, the name "dat" is evaluated. Therefore, the call stored does
not have the name 'dat', but has the evaluated content.

The following fits the same model. However, the call stores the name,
'mod', not the evaluated result, y ~ x1 + x2.

```r
lm_out3 <- lm(mod, data = dat)
lm_out3
#>
#> Call:
#> lm(formula = mod, data = dat)
#>
#> Coefficients:
#> (Intercept)   x1   x2
#> -0.5755  -0.4151  -0.2411
```

The following method works. However, I have to do a dummy call,
extract the stored call, and set formula to the result of
as.formula(mod):

```r
lm_out3 <- lm(mod, data = dat)
lm_out3
#>
#> Call:
#> lm(formula = mod, data = dat)
#>
#> Coefficients:
#> (Intercept)   x1   x2
#> -0.5755  -0.4151  -0.2411

call1 <- lm_out3$call
call1$formula <- as.formula(mod)
lm_out4 <- eval(call1)
lm_out4
#>
#> Call:
#> lm(formula = y ~ x1 + x2, data = dat)
#>
#> Coefficients:
#> (Intercept)   x1   x2
#> -0.5755  -0.4151  -0.2411
```

Is it possible to create the call directly, with only 'mod' evaluated,
and other arguments, e.g., 'dat', not evaluated?

Regards,
Shu Fai

__
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] Plot to a device and examine the plot?

2023-10-15 Thread Shu Fai Cheung
Thanks for the suggestion. I believe it is exactly what I need. I will
try this function. Thanks!

Regards,
Shu Fai

On Mon, Oct 16, 2023 at 3:39 AM Paul Murrell  wrote:
>
> Hi
>
> You could also look at dev.capture(), depending on which screen device
> you are using.
>
> Paul
>
> On 16/10/23 05:24, Duncan Murdoch wrote:
> > On 15/10/2023 12:05 p.m., Shu Fai Cheung wrote:
> >  > Let me clarify my question:
> >  >
> >  > plot.new()
> >  > polygon(c(.5, .5, .75, .8), c(.25, .3, .4, .5))
> >  >
> >  > If the device is an on-screen device, can I check whether a
> > particular area
> >  > has anything drawn on it, or, to be precise, whether the color of a
> >  > particular area has all pixels equal to the background color. That is, if
> >  > the background is white, can I know whether a particular area is white?
> >  >
> >  > E.g., in the case above, the area bounded by rect(0, 0, .25, .25) is
> >  > completely white, while the area bounded by rect(0, 0, .75, .75) is not
> >  > because part of border of the polygon, black by default, has been
> > drawn in
> >  > this area.
> >  >
> >  > If the device is an image file, then I can check the color of a pixel. I
> >  > would like to know whether I can do the same with an on-screen device.
> >  >
> >
> > I think the answer is that you can't do that in general. However, in
> > general you can copy a plot to a different device using dev.copy() and
> > examine it there. It won't be pixel-by-pixel identical, but will
> > contain the same components, likely with slightly different scaling and
> > positioning if the new device isn't the same as the old one.
> >
> > You can also save the commands that drew the plot using recordPlot() and
> > redraw it using replayPlot() (which is essentially what dev.copy()
> > does), but the format of the object saved by recordPlot() is not
> > documented, and is subject to change with R version changes.
> >
> > Duncan Murdoch
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > <https://stat.ethz.ch/mailman/listinfo/r-help>
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > <http://www.R-project.org/posting-guide.html>
> > and provide commented, minimal, self-contained, reproducible code.
>
> --
> Dr Paul Murrell
> Te Kura Tatauranga | Department of Statistics
> Waipapa Taumata Rau | The University of Auckland
> Private Bag 92019, Auckland 1142, New Zealand
> 64 9 3737599 x85392
> p...@stat.auckland.ac.nz
> www.stat.auckland.ac.nz/~paul/
>

__
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] Plot to a device and examine the plot?

2023-10-15 Thread Shu Fai Cheung
Thanks a lot for introducing these functions! I am not aware of them
but it seems that they can help me to do what I want to do.

Regards,
Shu Fai

Regards,
Shu Fai Cheung (張樹輝)


On Mon, Oct 16, 2023 at 12:24 AM Duncan Murdoch
 wrote:
>
> On 15/10/2023 12:05 p.m., Shu Fai Cheung wrote:
> > Let me clarify my question:
> >
> > plot.new()
> > polygon(c(.5, .5, .75, .8), c(.25, .3, .4, .5))
> >
> > If the device is an on-screen device, can I check whether a particular area
> > has anything drawn on it, or, to be precise, whether the color of a
> > particular area has all pixels equal to the background color. That is, if
> > the background is white, can I know whether a particular area is white?
> >
> > E.g.,  in the case above, the area bounded by rect(0, 0, .25, .25) is
> > completely white, while the area bounded by rect(0, 0, .75, .75) is not
> > because part of border of the polygon, black by default, has been drawn in
> > this area.
> >
> > If the device is an image file, then I can check the color of a pixel. I
> > would like to know whether I can do the same with an on-screen device.
> >
>
> I think the answer is that you can't do that in general.  However, in
> general you can copy a plot to a different device using dev.copy() and
> examine it there.  It won't be pixel-by-pixel identical, but will
> contain the same components, likely with slightly different scaling and
> positioning if the new device isn't the same as the old one.
>
> You can also save the commands that drew the plot using recordPlot() and
> redraw it using replayPlot() (which is essentially what dev.copy()
> does), but the format of the object saved by recordPlot() is not
> documented, and is subject to change with R version changes.
>
> Duncan Murdoch
>

__
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] Plot to a device and examine the plot?

2023-10-15 Thread Shu Fai Cheung
Let me clarify my question:

plot.new()
polygon(c(.5, .5, .75, .8), c(.25, .3, .4, .5))

If the device is an on-screen device, can I check whether a particular area
has anything drawn on it, or, to be precise, whether the color of a
particular area has all pixels equal to the background color. That is, if
the background is white, can I know whether a particular area is white?

E.g.,  in the case above, the area bounded by rect(0, 0, .25, .25) is
completely white, while the area bounded by rect(0, 0, .75, .75) is not
because part of border of the polygon, black by default, has been drawn in
this area.

If the device is an image file, then I can check the color of a pixel. I
would like to know whether I can do the same with an on-screen device.

Thanks.

Regards,
Shu Fai


On Sun, Oct 15, 2023 at 11:44 PM Shu Fai Cheung 
wrote:

> Sorry that I did not make my question clear enough.
>
> If the device is file based (e.g., a PNG file), then I believe I can do
> that, by using functions that can inspect an image. This is the solution I
> originally wanted to try. However, it requires creating a file in the
> process.
>
> I would like to see whether it is possible to inspect the canvas if the
> device is an on-screen one, like the pop-up window in R for Windows.
>
> Regards,
> Shu Fai
>
>
> On Sun, Oct 15, 2023 at 11:31 PM Jeff Newmiller 
> wrote:
>
>> This question is not clear to me. What is it you hope to retrieve from
>> the device?
>>
>> Note that the type of device in your example is system-dependent. The
>> content in a png() would be different than the content in a win.graph()
>> device.
>>
>> On October 15, 2023 8:04:00 AM PDT, Shu Fai Cheung <
>> shufai.che...@gmail.com> wrote:
>> >Hi All,
>> >
>> >I want to inspect the content of a plot generated by another function.
>> >
>> >For example:
>> >
>> >plot.new()
>> >polygon(c(.5, .5, .75, .8), c(.25, .3, .4, .5))
>> >
>> >A polygon will be drawn. If I do not know what has been done to generate
>> >the plot, is it possible to query the content in the active device?
>> >
>> >Regards,
>> >Shu Fai
>> >
>> >   [[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.
>>
>> --
>> Sent from my phone. Please excuse my brevity.
>>
>

[[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] Plot to a device and examine the plot?

2023-10-15 Thread Shu Fai Cheung
Sorry that I did not make my question clear enough.

If the device is file based (e.g., a PNG file), then I believe I can do
that, by using functions that can inspect an image. This is the solution I
originally wanted to try. However, it requires creating a file in the
process.

I would like to see whether it is possible to inspect the canvas if the
device is an on-screen one, like the pop-up window in R for Windows.

Regards,
Shu Fai


On Sun, Oct 15, 2023 at 11:31 PM Jeff Newmiller 
wrote:

> This question is not clear to me. What is it you hope to retrieve from the
> device?
>
> Note that the type of device in your example is system-dependent. The
> content in a png() would be different than the content in a win.graph()
> device.
>
> On October 15, 2023 8:04:00 AM PDT, Shu Fai Cheung <
> shufai.che...@gmail.com> wrote:
> >Hi All,
> >
> >I want to inspect the content of a plot generated by another function.
> >
> >For example:
> >
> >plot.new()
> >polygon(c(.5, .5, .75, .8), c(.25, .3, .4, .5))
> >
> >A polygon will be drawn. If I do not know what has been done to generate
> >the plot, is it possible to query the content in the active device?
> >
> >Regards,
> >Shu Fai
> >
> >   [[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.
>
> --
> Sent from my phone. Please excuse my brevity.
>

[[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] Plot to a device and examine the plot?

2023-10-15 Thread Shu Fai Cheung
Hi All,

I want to inspect the content of a plot generated by another function.

For example:

plot.new()
polygon(c(.5, .5, .75, .8), c(.25, .3, .4, .5))

A polygon will be drawn. If I do not know what has been done to generate
the plot, is it possible to query the content in the active device?

Regards,
Shu Fai

[[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] save(), load(), saveRDS(), and readRDS()

2023-09-28 Thread Shu Fai Cheung
Hi All,

There is a thread about the use of save(), load(), saveRDS(), and
loadRDS(). It led me to think about a question regarding them.

In my personal work, I prefer using saveRDS() and loadRDS() as I don't like
the risk of overwriting anything in the global environment. I also like the
freedom to name an object when reading it from a file.

However, for teaching, I have to teach save() and load() because, in my
discipline, it is common for researchers to share their datasets on the
internet using the format saved by save(), and so students need to know how
to use load() and what will happen when using it. Actually, I can't recall
encountering datasets shared by the .rds format. I have been wondering why
save() was usually used in that case.

That discussion led me to read the help pages again and I noticed the
following warning, from the help page of saveRDS():

"Files produced by saveRDS (or serialize to a file connection) are not
suitable as an interchange format between machines, for example to download
from a website. The files produced by save
 have a header identifying
the file type and so are better protected against erroneous use."

When will the problem mentioned in the warning occur? That is, when will a
file saved by saveRDS() not be read correctly? Saved in Linux and then read
in Windows? Is it possible to create a reproducible error?

Regards,
Shu Fai

[[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] save() and load()

2023-09-25 Thread Shu Fai Cheung
Hi,

You can try this:

head(irisdata)

Objects loaded by load() keep their names when being saved. In your
case, it is 'irisdata'.

You can also use verbose = TRUE to show the names of objects loaded:

load(file = "irisdataTest.RData", verbose = TRUE)

Hope this helps.

Regards,
Shu Fai

On Tue, Sep 26, 2023 at 9:24 AM AbouEl-Makarim Aboueissa
 wrote:
>
> Dear ALL:
>
> I am teaching statistical packages class this semester, in R programing I
> am trying to explain the use of save() and load() with an example using the
> iris data. It seems that the save() function works, BUT when I tried to
> load the data back to R, it seems that there is a problem(s), I could not
> figure out what went wrong.
>
> Any help would be highly appreciated.
>
>
> I saved the iris data in my computer in the text format, "iris.with.head.txt
> ".
>
> Here are my R codes:
>
> > irisdata<-read.table("G:/iris.with.head.txt", header=T)
> >
> > head(irisdata)
>   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
> 1  5.1 3.5  1.4 0.2  setosa
> 2  4.9 3.0  1.4 0.2  setosa
> 3  4.7 3.2  1.3 0.2  setosa
> 4  4.6 3.1  1.5 0.2  setosa
> 5  5.0 3.6  1.4 0.2  setosa
> 6  5.4 3.9  1.7 0.4  setosa
>
>
>
> *# saving the data as an .rda*
>
> save(irisdata,file="G:/irisdataTest.rda")
>
> *# load the data back to R*
>
> load(file="G:/irisdataTest.rda")
>
>
> >head(irisdataTest)
> Error in head(irisdataTest) : object 'irisdataTest' not found
>
> > irisdataTest
> Error: object 'irisdataTest' not found
>
>
>
> with many thanks
> abou
> __
>
>
> *AbouEl-Makarim Aboueissa, PhD*
>
> *Professor, Mathematics and Statistics*
> *Graduate Coordinator*
>
> *Department of Mathematics and Statistics*
> *University of Southern Maine*
>
> [[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] MASS::mvrnorm() on MKL may produce different numbers even when the seed is the same?

2023-08-17 Thread Shu Fai Cheung
Thanks a lot! I am not aware of this behavior! I will take this into
account in my work.

Regards,
Shu Fai

On Thu, Aug 17, 2023 at 10:31 PM Bill Dunlap  wrote:
>
> MKL's results can depend on the number of threads running and perhaps other 
> things.  They blame it on the non-associativity of floating point arithmetic. 
>  This article gives a way to make results repeatable:
>
> https://www.intel.com/content/www/us/en/developer/articles/technical/introduction-to-the-conditional-numerical-reproducibility-cnr.html
>
> -Bill
>
> On Wed, Aug 16, 2023 at 8:11 PM Shu Fai Cheung  
> wrote:
>>
>> Hi All,
>>
>> When addressing an error in one of my packages in the CRAN check at CRAN,
>> I found something strange which, I believe, is unrelated to my package.
>> I am not sure whether it is a bug or a known issue. Therefore, I would like
>> to have advice from experts here.
>>
>> The error at CRAN check occurred in a test, and only on R-devel with MKL. No
>> issues on all other platforms. The test compares two sets of random numbers,
>> supposed to be identical because generated from the same seed. However, when
>> I tried that in R-devel (2023-08-15 r84957) on Ubuntu 22.04.3 LTS, linked to
>> MKL (2023.2.0), I found that this is not the case in one situation.
>>
>> I don't know why but somehow only one particular set of means and
>> covariance matrix, among many others in the code, led to that problem.
>> Please find
>> below the code to reproduce the means and the covariance matrix (pardon me
>> for the long code):
>>
>> mu0 <- structure(c(0.52252416853188655, 0.39883382931927774, 
>> 1.6296642535174521,
>> 2.9045763671121816, -0.19816874840500939, 0.42610841566522556,
>> 0.30155498531316366, 1.0339601619394503, 3.4125587827873192,
>> 13.125481598475405, 19.275480386183503, 658.94225353462195, 
>> 1.0997193726987393,
>> 9.9980286642877214, 6.4947188998602092, -12.952617780813679,
>> 8.4991882784024799, 106.82209520278377, 0.19623712449219838), class =
>> c("numeric"))
>>
>> sigma0 <- structure(c(0.0047010182215945009, 1.4688424622163808e-17,
>> -2.5269697696885822e-17,
>> -1.2451516929358655e-18, 3.6553475725253408e-19, 3.2092363914356227e-19,
>> -2.6341938382508503e-18, 8.6439582878287556e-19, -1.295240602054123e-17,
>> 1.1154057497258702e-18, -8.2407621704807367e-33, -1.0891686096304908e-15,
>> -2.6800671450262819e-18, -0.0009225142979911, -1.4833018148344697e-16,
>> 2.292242132203e-16, -3.5128615476876035e-18, 2.8004770623734002e-33,
>> 2.0818099301348832e-34, 1.294907062174661e-17, 0.012788608283099183,
>> 1.5603789410838518e-15, 1.0425182851251724e-16, 2.758318745465066e-17,
>> -9.7047963858148572e-18, -7.4685438142667792e-17, -1.9426723638193857e-17,
>> -7.8775307262071738e-17, -4.0773523140359949e-17, 4.1212975037936416e-18,
>> -1.199934828824835e-14, 2.301740948367742e-17, -2.5410883836579325e-18,
>> -0.129172198695584, -1.9922957470809647e-14, 1.7077690739402761e-16,
>> 1.663702957392817e-30, 5.5362608501514495e-32, -2.5575000656645995e-17,
>> 7.3879960338026333e-16, 0.052246980157830372, -0.007670030643844231,
>> -1.5468402204507243e-17, 1.0518757786432683e-17, 2.9992131812421253e-16,
>> 1.6588830885081527e-17, -8.3276045185697992e-16, -3.3713008849128695e-15,
>> 1.5823980389961169e-17, 2.2243227896948194e-15, -3.2071920299461956e-17,
>> 5.0187645877462975e-18, -7.4622951185527265e-15, -0.44701112744184102,
>> -1.854066605407503e-16, -5.2511356793649759e-30, -6.7993385117653912e-32,
>> -1.5343740968067595e-18, -3.6785921616583949e-17, -0.0076700306438433133,
>> 0.019231143599164325, 5.7818784939170702e-18, 1.7819897158350214e-17,
>> -2.3927039320969442e-17, -3.9630951242056644e-18, 1.9779107704573469e-15,
>> 7.8103367010164485e-15, -1.2351275675875924e-17, -8.6959160183013082e-15,
>> -1.601002787863e-18, 3.0110116065267407e-19, 3.7155867714997651e-16,
>> -0.12490087179942209, -2.2029631919898945e-16, -8.0869824211018481e-31,
>> -5.058673986695e-33, 2.0667225433033359e-19, 1.0141735013231455e-18,
>> 3.2845977347050714e-17, -2.9662734037779067e-19, 9.9595082331331027e-05,
>> -1.9982466114987023e-18, 2.6420113469249318e-18, 5.2415327077690342e-19,
>> 3.007533948014542e-18, -5.1128158961673558e-17, 6.7791246559883483e-19,
>> -3.042327665547861e-15, 1.9523738271941284e-19, -4.0556768902105205e-20,
>> -1.024372770865438e-17, -0.010638955366526863, 2.5786757042756338e-17,
>> 1.3546130005558746e-32, 1.7947249868721723e-33, 1.9228405363682382e-19,
>> -4.1209014001548407e-18, 1.5228291529275058e-17, 2.0443162438349059e-17,
>> -1.9572693345249383e-18, 0.001270958102

Re: [R] MASS::mvrnorm() on MKL may produce different numbers even when the seed is the same?

2023-08-17 Thread Shu Fai Cheung
Sorry for the confusion caused by that sentence. I did
not expect using the same seed generates the same
sequence of numbers unconditionally. What puzzled me
is not having the same sequence when running the
same code repeatedly in the same computer in the same
session of R. The following four sets of results were generated
from a script with four copies of the two lines, ran in the
same session:

> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1]  0.4851605  0.5704446  1.6873036  2.7645014 -0.2020908
> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1]  0.4851605  0.5704446  1.6872818  2.7644787 -0.2023346
> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1]  0.4851605  0.5704446  1.6872818  2.7644787 -0.2023346
> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1]  0.4851605  0.5704446  1.6873036  2.7645014 -0.2020908

I could not figure out the pattern. I could not even
reproduce the order above. Sometimes, the order can be
different. Sometimes, the four sets of results can be identical.

I did a few more tests, on another computer (Windows 10),
with several copies of R installed. For R 4.3.1 and
R-devel (2023-08-12 r84939 ucrt), for which I did not touch
anything about BLAS and LAPACK, they consistently gave
the following results:

> set.seed(1234)
> MASS::mvrnorm(n = 5, mu = mu0, Sigma = sigma0)[1, 1:5]
[1]  0.4851605  0.5704446  1.6872818  2.7644787 -0.2023346

Therefore, it is possible that in my test with MKL, somehow
BLAS and LAPACK came with R were sometimes used, even
though the sessionInfo() output shows that MKL is used.

I am not sure whether it is because I did anything wrong in
setting up R to use MKL.

More on the machine on which I found that issue. It was a virtual
machine built using VirtualBox. The host system is Windows 10
(though I guess it should not matter). The OS is Ubuntu 22.04.3 LTS.
I installed MKL (2023.2.0) following these steps:

https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-download.html?operatingsystem=linux=aptpackagemanager

I compiled R-devel (2023-08-15 r84957) and then link MKL to it
following these steps (with some modifications, as the description
seems to be outdated):

https://www.intel.com/content/www/us/en/developer/articles/technical/quick-linking-intel-mkl-blas-lapack-to-r.html

Ideally, I should compile R with MKL, instead of using the quick
linking method. However, I am occupied with something else and
could not find an updated guide on how to do this quickly. The
guides I located were several years old.

Regards,
Shu Fai


On Thu, Aug 17, 2023 at 8:57 PM Ben Bolker  wrote:

>
>  > However, should the numbers
>  > generated identical if the same seed is used?
>
>I don't see how using the same seed can overcome floating-point
> differences across platforms (compilers etc.) stemming from differences
> in an eigen() computation (based on arcane details like use of
> registers, compiler optimizations, etc.).  set.seed() guarantees that
> the normal deviates generated by rnorm() will be identical, but those
> random numbers are then multiplied by the eigenvectors derived from
> eigen(), at which point the differences will crop up.
>
>This has been discussed on Twitter:
> https://twitter.com/wviechtb/status/1230078883317387264
>
> Wolfgang Viechtbauer, 2020-02-18:
>
> Another interesting reproducibility issue that came up. MASS::mvrnorm()
> can give different values despite setting the same seed. The problem:
> Results of eigen() (which is used by mvrnorm) can be machine dependent
> (help(eigen) does warn about this).
>
> Interestingly, mvtnorm::rmvnorm() with method="eigen" gives the same
> values across all machines I tested this on (and method="svd" give the
> same values). With method="chol" I get different values, but again
> consistent across machines.
>
> Ah, mvtnorm::rmvnorm() applies the results from eigen() in a different
> way that appears to be less (not?) affected by the indeterminacy in the
> eigenvectors. Quite clever.
>
>
> On 2023-08-16 11:10 p.m., Shu Fai Cheung wrote:
> > Hi All,
> >
> > When addressing an error in one of my packages in the CRAN check at CRAN,
> > I found something strange which, I believe, is unrelated to my package.
> > I am not sure whether it is a bug or a known issue. Therefore, I would
> like
> > to have advice from experts here.
> >
> > The error at CRAN check occurred in a test, and only on R-devel with
> MKL. No
> > issues on all other platforms. The test compares two sets of random
> numbers,
> > supposed to be identical because generated from the same seed. However,
> when
> > I tried that in R-devel (2023-08-15 r84957) on Ubuntu 22.04.3 LTS,
> 

[R] MASS::mvrnorm() on MKL may produce different numbers even when the seed is the same?

2023-08-16 Thread Shu Fai Cheung
Hi All,

When addressing an error in one of my packages in the CRAN check at CRAN,
I found something strange which, I believe, is unrelated to my package.
I am not sure whether it is a bug or a known issue. Therefore, I would like
to have advice from experts here.

The error at CRAN check occurred in a test, and only on R-devel with MKL. No
issues on all other platforms. The test compares two sets of random numbers,
supposed to be identical because generated from the same seed. However, when
I tried that in R-devel (2023-08-15 r84957) on Ubuntu 22.04.3 LTS, linked to
MKL (2023.2.0), I found that this is not the case in one situation.

I don't know why but somehow only one particular set of means and
covariance matrix, among many others in the code, led to that problem.
Please find
below the code to reproduce the means and the covariance matrix (pardon me
for the long code):

mu0 <- structure(c(0.52252416853188655, 0.39883382931927774, 1.6296642535174521,
2.9045763671121816, -0.19816874840500939, 0.42610841566522556,
0.30155498531316366, 1.0339601619394503, 3.4125587827873192,
13.125481598475405, 19.275480386183503, 658.94225353462195, 1.0997193726987393,
9.9980286642877214, 6.4947188998602092, -12.952617780813679,
8.4991882784024799, 106.82209520278377, 0.19623712449219838), class =
c("numeric"))

sigma0 <- structure(c(0.0047010182215945009, 1.4688424622163808e-17,
-2.5269697696885822e-17,
-1.2451516929358655e-18, 3.6553475725253408e-19, 3.2092363914356227e-19,
-2.6341938382508503e-18, 8.6439582878287556e-19, -1.295240602054123e-17,
1.1154057497258702e-18, -8.2407621704807367e-33, -1.0891686096304908e-15,
-2.6800671450262819e-18, -0.0009225142979911, -1.4833018148344697e-16,
2.292242132203e-16, -3.5128615476876035e-18, 2.8004770623734002e-33,
2.0818099301348832e-34, 1.294907062174661e-17, 0.012788608283099183,
1.5603789410838518e-15, 1.0425182851251724e-16, 2.758318745465066e-17,
-9.7047963858148572e-18, -7.4685438142667792e-17, -1.9426723638193857e-17,
-7.8775307262071738e-17, -4.0773523140359949e-17, 4.1212975037936416e-18,
-1.199934828824835e-14, 2.301740948367742e-17, -2.5410883836579325e-18,
-0.129172198695584, -1.9922957470809647e-14, 1.7077690739402761e-16,
1.663702957392817e-30, 5.5362608501514495e-32, -2.5575000656645995e-17,
7.3879960338026333e-16, 0.052246980157830372, -0.007670030643844231,
-1.5468402204507243e-17, 1.0518757786432683e-17, 2.9992131812421253e-16,
1.6588830885081527e-17, -8.3276045185697992e-16, -3.3713008849128695e-15,
1.5823980389961169e-17, 2.2243227896948194e-15, -3.2071920299461956e-17,
5.0187645877462975e-18, -7.4622951185527265e-15, -0.44701112744184102,
-1.854066605407503e-16, -5.2511356793649759e-30, -6.7993385117653912e-32,
-1.5343740968067595e-18, -3.6785921616583949e-17, -0.0076700306438433133,
0.019231143599164325, 5.7818784939170702e-18, 1.7819897158350214e-17,
-2.3927039320969442e-17, -3.9630951242056644e-18, 1.9779107704573469e-15,
7.8103367010164485e-15, -1.2351275675875924e-17, -8.6959160183013082e-15,
-1.601002787863e-18, 3.0110116065267407e-19, 3.7155867714997651e-16,
-0.12490087179942209, -2.2029631919898945e-16, -8.0869824211018481e-31,
-5.058673986695e-33, 2.0667225433033359e-19, 1.0141735013231455e-18,
3.2845977347050714e-17, -2.9662734037779067e-19, 9.9595082331331027e-05,
-1.9982466114987023e-18, 2.6420113469249318e-18, 5.2415327077690342e-19,
3.007533948014542e-18, -5.1128158961673558e-17, 6.7791246559883483e-19,
-3.042327665547861e-15, 1.9523738271941284e-19, -4.0556768902105205e-20,
-1.024372770865438e-17, -0.010638955366526863, 2.5786757042756338e-17,
1.3546130005558746e-32, 1.7947249868721723e-33, 1.9228405363682382e-19,
-4.1209014001548407e-18, 1.5228291529275058e-17, 2.0443162438349059e-17,
-1.9572693345249383e-18, 0.0012709581028842473, -0.0018515998737074948,
9.3128024867175073e-20, -5.1895788671618993e-17, 2.7373981615289174e-16,
1.2812977711597223e-17, -2.2792319486263267e-15, 4.1599503721278813e-19,
-3.7733269771394201e-20, 4.16234419478765e-17, -1.5986158133468129e-16,
-0.016037670900735834, -5.0763064708173244e-33, -1.0176066166236902e-50,
-5.9296704797665404e-18, -8.0698801402861772e-17, 2.9646619457173492e-16,
-2.8879431119718227e-17, 2.9393253663483117e-18, -0.0018515998737074959,
0.090335687736246548, 5.6849412731758135e-19, 8.8934458836007799e-17,
-4.1390858756690365e-16, 4.120323677410211e-16, 2.8000915545933503e-15,
2.8094462743052983e-17, 1.1636214841356605e-18, 8.151036749771e-16,
-3.004558574117108e-15, 0.0061666744014873013, -2.5381532354086619e-33,
-2.7110175619881585e-34, -3.9720857469692968e-18, -2.2722246209861298e-17,
6.5087631702810744e-18, -5.8833898464959171e-18, 6.0598974659958357e-19,
-6.4324667436451943e-19, -4.9865458549184915e-19, 0.010690736164778532,
-2.5717220776886191e-17, -2.9932234433250055e-17, -1.7586566364625091e-32,
-2.0572361612722435e-15, -4.1554892682395136e-18, 7.7947068522170098e-19,
2.2950757715435305e-16, -6.8563401956907788e-17, 8.3986032618135276e-18,
-3.0929447793865389e-45, 

Re: [R] printCoefmat() and zap.ind

2023-07-08 Thread Shu Fai Cheung
(Sorry for sending it twice. I forgot to reply
to the mailing list.)

Many many thanks for the comments and examples!

I could write my own function to achieve what
I want to do. However, I would like to find a method that
uses built-in functions only and prints the output in a format
identical to that of the default output of print.summary.lm(),
which uses printCoefmat() internally.

It seems that this cannot be done easily for now. This
is a workaround.

```r

set.seed(5689417)
n <- 1
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- .5 * x1 + .6 * x2 + rnorm(n, -0.0002366, .2)
dat <- data.frame(x1, x2, y)
out <- lm(y ~ x1 + x2, dat)
out_summary <- summary(out)
out_summary$coefficients[, "Estimate"] <-
  round(out_summary$coefficients[, "Estimate"], 4)
out_summary$coefficients[, "Std. Error"] <-
  round(out_summary$coefficients[, "Std. Error"], 4)

printCoefmat(out_summary$coefficients)
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept)   0. 0.00200.001
#> x10.5021 0.0020  254.70   <2e-16 ***
#> x20.6002 0.0020  301.23   <2e-16 ***

#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

I have to round the two columns first before calling
printCoefmat(). Not nice but works for now.

Regards,
Shu Fai Cheung

在 2023年7月8日週六 00:41,Martin Maechler  寫道:

> >>>>> Martin Maechler
> >>>>> on Fri, 7 Jul 2023 18:12:24 +0200 writes:
>
> >>>>> Shu Fai Cheung
> >>>>> on Thu, 6 Jul 2023 17:14:27 +0800 writes:
>
> >> Hi All,
>
> >> I would like to ask two questions about printCoefmat().
>
> > Good... this function, originally named print.coefmat(),
> > is 25 years old (in R) now:
>
> > 
> > r1902 | maechler | 1998-08-14 19:19:05 +0200 (Fri, 14 Aug 1998) |
> > Changed paths:
> > M R-0-62-patches/CHANGES
> > M R-0-62-patches/src/library/base/R/anova.R
> > M R-0-62-patches/src/library/base/R/glm.R
> > M R-0-62-patches/src/library/base/R/lm.R
> > M R-0-62-patches/src/library/base/R/print.R
>
> > print.coefmat(.) about ok
> > 
>
> > (yes, at the time, the 'stats' package did not exist yet ..)
>
> > so it may be a good time to look at it.
>
>
> >> First, I found a behavior of printCoefmat() that looks strange to
> me,
> >> but I am not sure whether this is an intended behavior:
>
> >> ``` r
> >> set.seed(5689417)
> >> n <- 1
> >> x1 <- rnorm(n)
> >> x2 <- rnorm(n)
> >> y <- .5 * x1 + .6 * x2 + rnorm(n, -0.0002366, .2)
> >> dat <- data.frame(x1, x2, y)
> >> out <- lm(y ~ x1 + x2, dat)
> >> out_summary <- summary(out)
> >> printCoefmat(out_summary$coefficients)
> >> #>   Estimate Std. Error t value Pr(>|t|)
> >> #> (Intercept) 1.7228e-08 1.9908e-030.001
> >> #> x1  5.0212e-01 1.9715e-03  254.70   <2e-16 ***
> >> #> x2  6.0016e-01 1.9924e-03  301.23   <2e-16 ***
> >> #> ---
> >> #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> >> printCoefmat(out_summary$coefficients,
> >> zap.ind = 1,
> >> digits = 4)
> >> #> Estimate Std. Error t value Pr(>|t|)
> >> #> (Intercept) 0.00   0.001991 0.01
> >> #> x1  0.502100   0.001971   254.7   <2e-16 ***
> >> #> x2  0.600200   0.001992   301.2   <2e-16 ***
> >> #> ---
> >> #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >> ```
>
> >> With zap.ind = 1, the values in "Estimate" were correctly
> >> zapped using digits = 4. However, by default, "Estimate"
> >> and "Std. Error" are formatted together. Because the
> >> standard errors are small, with digits = 4, zero's were added
> >> to values in "Estimate", resulting in "0.502100" and
> >> "0.600200", which are misleading because, if rounded to
> >> the 6th decimal place, the values to be displayed should
> >> be "0.502122" and "0.600162".
>
> >> Is this behavior of printCoefmat() intended/normal?
>
> 

[R] printCoefmat() and zap.ind

2023-07-06 Thread Shu Fai Cheung
Hi All,

I would like to ask two questions about printCoefmat().

First, I found a behavior of printCoefmat() that looks strange to me,
but I am not sure whether this is an intended behavior:

``` r
set.seed(5689417)
n <- 1
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- .5 * x1 + .6 * x2 + rnorm(n, -0.0002366, .2)
dat <- data.frame(x1, x2, y)
out <- lm(y ~ x1 + x2, dat)
out_summary <- summary(out)
printCoefmat(out_summary$coefficients)
#>   Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.7228e-08 1.9908e-030.001
#> x1  5.0212e-01 1.9715e-03  254.70   <2e-16 ***
#> x2  6.0016e-01 1.9924e-03  301.23   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
printCoefmat(out_summary$coefficients,
 zap.ind = 1,
 digits = 4)
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.00   0.001991 0.01
#> x1  0.502100   0.001971   254.7   <2e-16 ***
#> x2  0.600200   0.001992   301.2   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

With zap.ind = 1, the values in "Estimate" were correctly
zapped using digits = 4. However, by default, "Estimate"
and "Std. Error" are formatted together. Because the
standard errors are small, with digits = 4, zero's were added
to values in "Estimate", resulting in "0.502100" and
"0.600200", which are misleading because, if rounded to
the 6th decimal place, the values to be displayed should
be "0.502122" and "0.600162".

Is this behavior of printCoefmat() intended/normal?

Second, how can I use zap without this behavior?
In cases like the one above, I need to use zap such that
the intercept will not be displayed in scientific notation.
Disabling scientific notation cannot achieve the desired
goal.

I tried adding cs.ind = 1:

```r
printCoefmat(out_summary$coefficients,
 zap.ind = 1,
 digits = 4,
 cs.ind = 1)
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept)   0.   0.001991 0.01
#> x10.5021   0.001971   254.7   <2e-16 ***
#> x20.6002   0.001992   301.2   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

However, this solution is not ideal because the numbers
of decimal places of "Estimate" and "Std. Error" are
different. How can I get the output like this one?

```r
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept)   0.   0.0020 0.01
#> x10.5021   0.0020   254.7   <2e-16 ***
#> x20.6002   0.0020   301.2   <2e-16 ***
```

Thanks for your attention.

Regards,
Shu Fai Cheung

__
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.