[R] How to create gridded data

2018-11-12 Thread lily li
Hi R users,

I have a question about manipulating data. For example, I have DF1 as the
following, how to transform it to a gridded dataset DF2? In DF2, each value
Precip is an attribute of the corresponding grid cell. So DF2 is like a
spatial surface, and can be imported to ArcGIS. Thanks for your help.

DF1
latitude   longitude   Precip
45.5   110.5 3.2
45.5   1115.0
45.5   111.5 1.8
45.5   1122.0
46  110.5 6.1
46  1114.5
46  111.5 7.8
46  1125.5
...


DF2
6.1   4.5   7.8   5.5
3.2   5.0   1.8   2.0
...

[[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] which element is duplicated?

2018-11-12 Thread Bert Gunter
"I'd like to see a cleverer solution that vectorizes..."

and Herve provided it.


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 Mon, Nov 12, 2018 at 9:43 PM Bert Gunter  wrote:

> It is not clear to what you want for the general case. Perhaps:
>
> > v <- letters[c(2,2,1,2,1,1)]
> > wh <- tapply(seq_along(v),factor(v), '[',1)
> > w <- wh[match(v,v[wh])]
> > w
> b b a b a a
> 1 1 3 1 3 3
> > ## and if you want NA's for the first occurences of unique values
> > ## of course:
> > w[wh] <- NA
> > w
>  b  b  a  b  a  a
> NA  1 NA  1  3  3
>
> I'd like to see a cleverer solution that vectorizes and avoids the
> tapply(), though.
>
> Cheers,
> Bert
>
>
>
>
> On Mon, Nov 12, 2018 at 8:33 PM Bert Gunter 
> wrote:
>
>> > match(v, unique(v))
>> [1] 1 2 2 1
>>
>> 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 Mon, Nov 12, 2018 at 5:08 PM Duncan Murdoch 
>> wrote:
>>
>>> The duplicated() function gives TRUE if an item in a vector (or row in a
>>> matrix, etc.) is a duplicate of an earlier item.  But what I would like
>>> to know is which item does it duplicate?
>>>
>>> For example,
>>>
>>> v <- c("a", "b", "b", "a")
>>> duplicated(v)
>>>
>>> returns
>>>
>>> [1] FALSE FALSE  TRUE  TRUE
>>>
>>> What I want is a fast way to calculate
>>>
>>>   [1] NA NA 2 1
>>>
>>> or (equally useful to me)
>>>
>>>   [1] 1 2 2 1
>>>
>>> The result should have the property that if result[i] == j, then v[i] ==
>>> v[j], at least for i != j.
>>>
>>> Does this already exist somewhere, or is it easy to write?
>>>
>>> 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.
>>>
>>

[[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] which element is duplicated?

2018-11-12 Thread Bert Gunter
It is not clear to what you want for the general case. Perhaps:

> v <- letters[c(2,2,1,2,1,1)]
> wh <- tapply(seq_along(v),factor(v), '[',1)
> w <- wh[match(v,v[wh])]
> w
b b a b a a
1 1 3 1 3 3
> ## and if you want NA's for the first occurences of unique values
> ## of course:
> w[wh] <- NA
> w
 b  b  a  b  a  a
NA  1 NA  1  3  3

I'd like to see a cleverer solution that vectorizes and avoids the
tapply(), though.

Cheers,
Bert




On Mon, Nov 12, 2018 at 8:33 PM Bert Gunter  wrote:

> > match(v, unique(v))
> [1] 1 2 2 1
>
> 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 Mon, Nov 12, 2018 at 5:08 PM Duncan Murdoch 
> wrote:
>
>> The duplicated() function gives TRUE if an item in a vector (or row in a
>> matrix, etc.) is a duplicate of an earlier item.  But what I would like
>> to know is which item does it duplicate?
>>
>> For example,
>>
>> v <- c("a", "b", "b", "a")
>> duplicated(v)
>>
>> returns
>>
>> [1] FALSE FALSE  TRUE  TRUE
>>
>> What I want is a fast way to calculate
>>
>>   [1] NA NA 2 1
>>
>> or (equally useful to me)
>>
>>   [1] 1 2 2 1
>>
>> The result should have the property that if result[i] == j, then v[i] ==
>> v[j], at least for i != j.
>>
>> Does this already exist somewhere, or is it easy to write?
>>
>> 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.
>>
>

[[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] which element is duplicated?

2018-11-12 Thread Pages, Herve
Hi,

On 11/12/18 17:08, Duncan Murdoch wrote:
> The duplicated() function gives TRUE if an item in a vector (or row in 
> a matrix, etc.) is a duplicate of an earlier item.  But what I would 
> like to know is which item does it duplicate?
>
> For example,
>
> v <- c("a", "b", "b", "a")
> duplicated(v)
>
> returns
>
> [1] FALSE FALSE  TRUE  TRUE
>
> What I want is a fast way to calculate
>
>  [1] NA NA 2 1
>
> or (equally useful to me)
>
>  [1] 1 2 2 1
>
> The result should have the property that if result[i] == j, then v[i] 
> == v[j], at least for i != j.
>
> Does this already exist somewhere, or is it easy to write?

I generally use match() for that:

 > v <- c("a", "b", "b", "a")

 > match(v, v)

[1] 1 2 2 1

H.

>
> Duncan Murdoch
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dhelp=DwICAg=eRAMFD45gAfqt84VtBcfhQ=BK7q3XeAvimeWdGbWY_wJYbW0WYiZvSXAJJKaaPhzWA=APEsp-OzJs6YdfshtiYe715BsAor8xTu26lpN4KGOrU=opxT_5og2YaWKdiXD-cRz0gWxGGMRG6kq20Jo8711qA=
>  
>
> PLEASE do read the posting guide 
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.R-2Dproject.org_posting-2Dguide.html=DwICAg=eRAMFD45gAfqt84VtBcfhQ=BK7q3XeAvimeWdGbWY_wJYbW0WYiZvSXAJJKaaPhzWA=APEsp-OzJs6YdfshtiYe715BsAor8xTu26lpN4KGOrU=ZaPnASTzuEmE8EHqFL6F5wYkPhhg_uv-CMrGjY2-_Q4=
> and provide commented, minimal, self-contained, reproducible code.

-- 
Hervé Pagès

Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024

E-mail: hpa...@fredhutch.org
Phone:  (206) 667-5791
Fax:(206) 667-1319

__
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] which element is duplicated?

2018-11-12 Thread Bert Gunter
> match(v, unique(v))
[1] 1 2 2 1

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 Mon, Nov 12, 2018 at 5:08 PM Duncan Murdoch 
wrote:

> The duplicated() function gives TRUE if an item in a vector (or row in a
> matrix, etc.) is a duplicate of an earlier item.  But what I would like
> to know is which item does it duplicate?
>
> For example,
>
> v <- c("a", "b", "b", "a")
> duplicated(v)
>
> returns
>
> [1] FALSE FALSE  TRUE  TRUE
>
> What I want is a fast way to calculate
>
>   [1] NA NA 2 1
>
> or (equally useful to me)
>
>   [1] 1 2 2 1
>
> The result should have the property that if result[i] == j, then v[i] ==
> v[j], at least for i != j.
>
> Does this already exist somewhere, or is it easy to write?
>
> 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.
>

[[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] semiparametric manova

2018-11-12 Thread David L Carlson
With two variables there are no combinations with less than 2 observations. 
Here's the part of the data you provided:

> df <- structure(list(V1 = c(200, 200, 200, 200, 200, 200, 200, 200, 
200, 200, 200, 200, 200, 200, 200, 500, 500, 500, 500, 500, 500, 
500, 500, 500, 500, 350, 350, 350, 350, 350, 200, 200, 200, 200, 
200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 500, 500, 
500, 500, 500, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 
350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 500, 500, 500, 
500, 500, 500, 500, 500, 500, 500, 350, 350, 350, 350, 350, 500, 
500, 500, 500, 500), V2 = c(16, 16, 16, 16, 16, 8, 8, 8, 8, 8, 
8, 8, 8, 8, 8, 16, 16, 16, 16, 16, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
8, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 16, 16, 16, 16, 16, 
8, 8, 8, 8, 8, 16, 16, 16, 16, 16, 8, 8, 8, 8, 8, 16, 16, 16, 
16, 16, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 16, 16, 16, 16, 
16, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24), V3 = c(16, 16, 16, 
16, 16, 16, 16, 16, 16, 16, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 23, 
23, 23, 23, 23, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 23, 23, 23, 23, 
23, 23, 23, 23, 23, 23, 16, 16, 16, 16, 16, 9, 9, 9, 9, 9, 23, 
23, 23, 23, 23, 16, 16, 16, 16, 16, 23, 23, 23, 23, 23, 16, 16, 
16, 16, 16, 23, 23, 23, 23, 23, 16, 16, 16, 16, 16, 9, 9, 9, 
9, 9)), class = "data.frame", row.names = c(NA, -90L))

> xtabs(~V1+V2+V3, df) # There are 9 cells with 0 entries. That is the problem.
, , V3 = 9

 V2
V18 16 24
  200 5  0  5
  350 5  5  0
  500 0  5  5

, , V3 = 16

 V2
V18 16 24
  200 5  5  0
  350 0  5  5
  500 5  0  5

, , V3 = 23

 V2
V18 16 24
  200 0  5  5
  350 5  0  5
  500 5  5  0

> xtabs(~V1+V2, df)  # No cells < 2 with V1, V2
 V2
V1 8 16 24
  200 10 10 10
  350 10 10 10
  500 10 10 10
> xtabs(~V1+V3, df)  # No cells < 2 with V1, V3
 V3
V1 9 16 23
  200 10 10 10
  350 10 10 10
  500 10 10 10
> xtabs(~V2+V3, df)  # No cells < 2 with V2, V3
V3
V29 16 23
  8  10 10 10
  16 10 10 10
  24 10 10 10


David L. Carlson
Department of Anthropology
Texas A University



-Original Message-
From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Yectli Huerta 
via R-help
Sent: Monday, November 12, 2018 6:25 PM
To: r-help@r-project.org
Subject: Re: [R] semiparametric manova

thanks for the replies.

i don't believe the data is the problem. here you see how i used 3 variables 
and it fails,
but when i use any combination of 2 variables, it does work

> head(df)
   V1 V2 V3 V4 V5   V6   V7
1 200 16 16  3 64 5.584092e+13 1.616745e+12
2 200 16 16  3 64 5.589262e+13 1.715906e+12
3 200 16 16  3 64 5.588578e+13 1.714084e+12
4 200 16 16  3 64 5.588061e+13 1.651920e+12
5 200 16 16  3 64 5.589810e+13 1.624824e+12
6 200  8 16  1 48 5.585124e+13 1.689478e+12
> library(MANOVA.RM)
> df$V1
[1] 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200 500 500 500 500
[20] 500 500 500 500 500 500 350 350 350 350 350 200 200 200 200 200 200 200 200
[39] 200 200 200 200 200 200 200 500 500 500 500 500 350 350 350 350 350 350 350
[58] 350 350 350 350 350 350 350 350 350 350 350 350 350 500 500 500 500 500 500
[77] 500 500 500 500 350 350 350 350 350 500 500 500 500 500
> df$V2
[1] 16 16 16 16 16  8  8  8  8  8  8  8  8  8  8 16 16 16 16 16  8  8  8  8  8
[26]  8  8  8  8  8 24 24 24 24 24 24 24 24 24 24 16 16 16 16 16  8  8  8  8  8
[51] 16 16 16 16 16  8  8  8  8  8 16 16 16 16 16 24 24 24 24 24 24 24 24 24 24
[76] 16 16 16 16 16 24 24 24 24 24 24 24 24 24 24
> df$V3
[1] 16 16 16 16 16 16 16 16 16 16  9  9  9  9  9  9  9  9  9  9 23 23 23 23 23
[26]  9  9  9  9  9  9  9  9  9  9 23 23 23 23 23 23 23 23 23 23 16 16 16 16 16
[51]  9  9  9  9  9 23 23 23 23 23 16 16 16 16 16 23 23 23 23 23 16 16 16 16 16
[76] 23 23 23 23 23 16 16 16 16 16  9  9  9  9  9
>
> MANOVA.wide(cbind(V6,V7)~V1*V2*V3,data=df,seed=1234)
Error in MANOVA.wide(cbind(V6, V7) ~ V1 * V2 * V3, data = df, seed = 1234) :
  There is at least one factor-level combination
   with less than 2 observations!

> MANOVA.wide(cbind(V6,V7)~V1*V2,data=df,seed=1234)
Call:
cbind(V6, V7) ~ V1 * V2

Wald-Type Statistic (WTS):
  Test statistic df p-value
V1    17.870  4   0.001
V2    20.392  4   0.000
V1:V2 24.127  8   0.002



> MANOVA.wide(cbind(V6,V7)~V1*V3,data=df,seed=1234)
Call:
cbind(V6, V7) ~ V1 * V3

Wald-Type Statistic (WTS):
  Test statistic df p-value
V1    18.566  4   0.001
V3    19.894  4   0.001
V1:V3 27.330  8   0.001

...
> MANOVA.wide(cbind(V6,V7)~V2*V3,data=df,seed=1234)
Call:
cbind(V6, V7) ~ V2 * V3

Wald-Type Statistic (WTS):
  Test statistic df p-value
V2    20.139  4   0.000
V3    19.947  4   0.001
V2:V3 32.088  8   0.000


__
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 

Re: [R] saveRDS() and readRDS() Why? [solved, pretty much anyway)

2018-11-12 Thread p_connolly

On 2018-11-13 12:55, William Dunlap wrote:

You wrote:
  ## On Windows 3.4.2


x <- airquality
saveRDS(x, file = "x.rds")
saveRDS(x, file = "y.rds")


  Files x.rds and y.rds are identical in size but utterly different in
content.

 Wow!  Can you show us the results of
  x <- datasets::airquality
  saveRDS(x, file="x.rds")
  saveRDS(x, file="y.rds")
  tools::md5sum(c("x.rds", "y.rds"))
  dput(readBin("x.rds", what="raw", n=file.size("x.rds")))
  dput(readBin("y.rds", what="raw", n=file.size("y.rds")))

(Copy and paste, as text, from the R session, so we can see the input
and the output in context.


If I do that on Linux or Windows, I get identical files.

If I copy x.rds and y.rds from Windows to Linux (which is what I wish to 
be able to do) using a shared folder between the Windows VirtualBox host 
to the Linux guest, I  get this:




x <- readRDS(file =  "x.rds")

Error in readRDS(file = "x.rds") : error reading from connection

x <- readRDS(file =  "y.rds")
  tools::md5sum(c("x.rds", "y.rds"))

 x.rds  y.rds
"5fef054848f39b4be02b7c54f1c71a20" "978a64d1dd342d16a381c9ca728d3665"

  dput(readBin("x.rds", what="raw", n=file.size("x.rds")))

as.raw(c(0x1f, 0x8b, 0x08, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x06, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 

Re: [R] which element is duplicated?

2018-11-12 Thread Michael Sumner
what about   as.integer(factor(v, levels = unique(v)))

I recall very clearly when I realized the power of this feature of
factor(), but I've not seen it discussed much.

Cheers, Mike.

On Tue, 13 Nov 2018 at 12:08 Duncan Murdoch 
wrote:

> The duplicated() function gives TRUE if an item in a vector (or row in a
> matrix, etc.) is a duplicate of an earlier item.  But what I would like
> to know is which item does it duplicate?
>
> For example,
>
> v <- c("a", "b", "b", "a")
> duplicated(v)
>
> returns
>
> [1] FALSE FALSE  TRUE  TRUE
>
> What I want is a fast way to calculate
>
>   [1] NA NA 2 1
>
> or (equally useful to me)
>
>   [1] 1 2 2 1
>
> The result should have the property that if result[i] == j, then v[i] ==
> v[j], at least for i != j.
>
> Does this already exist somewhere, or is it easy to write?
>
> 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.
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

[[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] which element is duplicated?

2018-11-12 Thread Duncan Murdoch
The duplicated() function gives TRUE if an item in a vector (or row in a 
matrix, etc.) is a duplicate of an earlier item.  But what I would like 
to know is which item does it duplicate?


For example,

v <- c("a", "b", "b", "a")
duplicated(v)

returns

[1] FALSE FALSE  TRUE  TRUE

What I want is a fast way to calculate

 [1] NA NA 2 1

or (equally useful to me)

 [1] 1 2 2 1

The result should have the property that if result[i] == j, then v[i] == 
v[j], at least for i != j.


Does this already exist somewhere, or is it easy to write?

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] semiparametric manova

2018-11-12 Thread Yectli Huerta via R-help
thanks for the replies.

i don't believe the data is the problem. here you see how i used 3 variables 
and it fails,
but when i use any combination of 2 variables, it does work

> head(df)
   V1 V2 V3 V4 V5   V6   V7
1 200 16 16  3 64 5.584092e+13 1.616745e+12
2 200 16 16  3 64 5.589262e+13 1.715906e+12
3 200 16 16  3 64 5.588578e+13 1.714084e+12
4 200 16 16  3 64 5.588061e+13 1.651920e+12
5 200 16 16  3 64 5.589810e+13 1.624824e+12
6 200  8 16  1 48 5.585124e+13 1.689478e+12
> library(MANOVA.RM)
> df$V1
[1] 200 200 200 200 200 200 200 200 200 200 200 200 200 200 200 500 500 500 500
[20] 500 500 500 500 500 500 350 350 350 350 350 200 200 200 200 200 200 200 200
[39] 200 200 200 200 200 200 200 500 500 500 500 500 350 350 350 350 350 350 350
[58] 350 350 350 350 350 350 350 350 350 350 350 350 350 500 500 500 500 500 500
[77] 500 500 500 500 350 350 350 350 350 500 500 500 500 500
> df$V2
[1] 16 16 16 16 16  8  8  8  8  8  8  8  8  8  8 16 16 16 16 16  8  8  8  8  8
[26]  8  8  8  8  8 24 24 24 24 24 24 24 24 24 24 16 16 16 16 16  8  8  8  8  8
[51] 16 16 16 16 16  8  8  8  8  8 16 16 16 16 16 24 24 24 24 24 24 24 24 24 24
[76] 16 16 16 16 16 24 24 24 24 24 24 24 24 24 24
> df$V3
[1] 16 16 16 16 16 16 16 16 16 16  9  9  9  9  9  9  9  9  9  9 23 23 23 23 23
[26]  9  9  9  9  9  9  9  9  9  9 23 23 23 23 23 23 23 23 23 23 16 16 16 16 16
[51]  9  9  9  9  9 23 23 23 23 23 16 16 16 16 16 23 23 23 23 23 16 16 16 16 16
[76] 23 23 23 23 23 16 16 16 16 16  9  9  9  9  9
>
> MANOVA.wide(cbind(V6,V7)~V1*V2*V3,data=df,seed=1234)
Error in MANOVA.wide(cbind(V6, V7) ~ V1 * V2 * V3, data = df, seed = 1234) :
  There is at least one factor-level combination
   with less than 2 observations!

> MANOVA.wide(cbind(V6,V7)~V1*V2,data=df,seed=1234)
Call:
cbind(V6, V7) ~ V1 * V2

Wald-Type Statistic (WTS):
  Test statistic df p-value
V1    17.870  4   0.001
V2    20.392  4   0.000
V1:V2 24.127  8   0.002



> MANOVA.wide(cbind(V6,V7)~V1*V3,data=df,seed=1234)
Call:
cbind(V6, V7) ~ V1 * V3

Wald-Type Statistic (WTS):
  Test statistic df p-value
V1    18.566  4   0.001
V3    19.894  4   0.001
V1:V3 27.330  8   0.001

...
> MANOVA.wide(cbind(V6,V7)~V2*V3,data=df,seed=1234)
Call:
cbind(V6, V7) ~ V2 * V3

Wald-Type Statistic (WTS):
  Test statistic df p-value
V2    20.139  4   0.000
V3    19.947  4   0.001
V2:V3 32.088  8   0.000



signature.asc
Description: OpenPGP digital 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] saveRDS() and readRDS() Why? [solved, pretty much anyway)

2018-11-12 Thread William Dunlap via R-help
You wrote:
  ## On Windows 3.4.2

>   x <- airquality
>   saveRDS(x, file = "x.rds")
>   saveRDS(x, file = "y.rds")
>

  Files x.rds and y.rds are identical in size but utterly different in
content.

Wow!  Can you show us the results of
  x <- datasets::airquality
  saveRDS(x, file="x.rds")
  saveRDS(x, file="y.rds")
  tools::md5sum(c("x.rds", "y.rds"))
  dput(readBin("x.rds", what="raw", n=file.size("x.rds")))
  dput(readBin("y.rds", what="raw", n=file.size("y.rds")))

(Copy and paste, as text, from the R session, so we can see the input and
the output in context.





Bill Dunlap
TIBCO Software
wdunlap tibco.com

On Mon, Nov 12, 2018 at 12:42 PM, p_connolly 
wrote:

> On 2018-11-12 22:49, peter dalgaard wrote:
>
>> Er, where, what, how? I can't reproduce that, at least not on 3.5.1 on
>> MacOS:
>>
>> x <- airquality
>>> saveRDS(x, file = "x.rds")
>>> x <- NULL
>>> x <- readRDS(file = "x.rds")
>>> x
>>>
>> Ozone Solar.R Wind Temp Month Day
>> 1  41 190  7.4   67 5   1
>> 2  36 118  8.0   72 5   2
>> 3  12 149 12.6   74 5   3
>> ...
>>
>> Looks fine to me.
>>
>>
> It seems to work fine using the same installation to read as used for the
> save.
> But it's a different story if the save was done on a Windows installation
> and read on a Linux installation.
>
> ## On Windows 3.4.2
>
>> x <- airquality
>> saveRDS(x, file = "x.rds")
>> saveRDS(x, file = "y.rds")
>>
>
> Files x.rds and y.rds are identical in size but utterly different in
> content.
>
> ## On Linux 3.5.1
>
> x <- readRDS(file =  "x.rds")
>>
> Error in readRDS(file = "x.rds") : error reading from connection
>
>> x <- readRDS(file =  "y.rds")
>> head(x)
>>
>   Ozone Solar.R Wind Temp Month Day
> 141 190  7.4   67 5   1
> 236 118  8.0   72 5   2
> 312 149 12.6   74 5   3
> 418 313 11.5   62 5   4
> 5NA  NA 14.3   56 5   5
> 628  NA 14.9   66 5   6
>
> It might just be the age of the Windows installation.  I don't have much
> use for Windows, so I haven't had much inclination to install a newer
> version.
>
> YMMV
>
>
>
> ?
>> -pd
>>
>>
>> On 12 Nov 2018, at 08:28 , Patrick Connolly 
>>> wrote:
>>>
>>> The solution was very simple.  Don't use the same name for the rds
>>> file  as used for the R object, viz a vie:
>>>
>>> saveRDS(x, file = "x.rds")
>>> and
>>> x <- readRDS(file = "x.rds")
>>>
>>> will not work; however
>>>
>>> saveRDS(x, file = "y.rds")
>>> and
>>> x <- readRDS(file = "y.rds")
>>> will work.
>>>
>>> An undocumented feature?
>>>
>>> Thanks to all who contributed.
>>>
>>
>> --
>> Peter Dalgaard, Professor,
>> Center for Statistics, Copenhagen Business School
>> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
>> Phone: (+45)38153501
>> Office: A 4.23
>> Email: pd@cbs.dk  Priv: pda...@gmail.com
>>
>
> __
> 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/posti
> ng-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] semiparametric manova

2018-11-12 Thread David Winsemius


On 11/12/18 12:37 PM, Yectli Huerta via R-help wrote:
> Hello,
>
> I was wondering if there are other packages like MANOVA.RM that could be used 
> to analysis non normal distributions. I have to analyze data with more than 2 
> predictor variables and a similar number of response variables. When I try 
> the function MANOVA.wide with more than 2 predictor variables, I get
>
> There is at least one factor-level combination
>     with less than 2 observations!
>
> Is there another package out there that can be used to analyze the 
> significance of more than 2 predictor variables?


That is a warning about a problem with your data in one or more of the 
"cells" defined by your factor levels. It's not going to be solved by 
choosing another package. You should be able to determine which variable 
is causing this warning by first looking at the data, a step should 
_always_ precede running any multivariate procedure.

-- 

David.

> thanks,
>
> yah
>
> __
> 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] semiparametric manova

2018-11-12 Thread David L Carlson
The error message does not say anything about having more than two predictor 
variables. It says that one of the combinations of the predictor variables has 
less than 2 observations (i.e. 1 or 0 observations). That is probably an issue 
of your sample size. You may need to consider combining some of the categories 
in your predictor variables or increasing your sample size. Are your predictor 
variables coded as factors? 


David L Carlson
Department of Anthropology
Texas A University
College Station, TX 77843-4352


-Original Message-
From: R-help  On Behalf Of Yectli Huerta via 
R-help
Sent: Monday, November 12, 2018 2:38 PM
To: r-help@r-project.org
Subject: [R] semiparametric manova

Hello,

I was wondering if there are other packages like MANOVA.RM that could be used 
to analysis non normal distributions. I have to analyze data with more than 2 
predictor variables and a similar number of response variables. When I try the 
function MANOVA.wide with more than 2 predictor variables, I get

There is at least one factor-level combination
   with less than 2 observations!

Is there another package out there that can be used to analyze the significance 
of more than 2 predictor variables?

thanks,

yah
__
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] Reporting binomial logistic regression from R results

2018-11-12 Thread Frodo Jedi
Dear Bert,
I understand and thanks for your recommendation. Unfortunately I do not
have any possibility to contact a statistical expert at the moment. So this
forum experts' recommendation would be crucial to me to understand how R
works in relation to my question.
I hope that someone could reply to my last questions.

Best regards

FJ

On Mon, Nov 12, 2018 at 7:48 PM Bert Gunter  wrote:

> Generally speaking, this list is about questions on R programming, not
> statistical issues. However, I grant you that your queries are in something
> of a gray area intersecting both.
>
> Nevertheless, based on your admitted confusion, I would recommend that you
> find a local statistical expert with whom you can consult 1-1 if at all
> possible. As others have already noted, you statistical understanding is
> muddy, and it can be quite difficult to resolve such confusion in online
> forums like this that cannot provide the close back and forth that may be
> required (as well as further appropriate study).
>
> Best,
> Bert
>
> On Mon, Nov 12, 2018 at 11:09 AM Frodo Jedi <
> frodojedi.mailingl...@gmail.com> wrote:
>
>> Dear Peter and Eik,
>> I am very grateful to you for your replies.
>> My current understanding is that from the GLM analysis I can indeed
>> conclude that the response predicted by System A is significantly
>> different
>> from that of System B, while the pairwise comparison A vs C leads to non
>> significance. Now the Wald test seems to be correct only for Systems B vs
>> C, indicating that the pairwise System B vs System C is significant. Am I
>> correct?
>>
>> However, my current understanding is also that I should use contrasts
>> instead of the wald test. So the default contrasts is with the System A,
>> now I should re-perform the GLM with another base. I tried to use the
>> option "contrasts" of the glm:
>>
>> > fit1 <- glm(Response ~ System, data = scrd, family = "binomial",
>> contrasts = contr.treatment(3, base=1,contrasts=TRUE))
>> > summary(fit1)
>>
>> > fit2 <- glm(Response ~ System, data = scrd, family = "binomial",
>> contrasts = contr.treatment(3, base=2,contrasts=TRUE))
>> > summary(fit2)
>>
>> > fit3 <- glm(Response ~ System, data = scrd, family = "binomial",
>> contrasts = contr.treatment(3, base=3,contrasts=TRUE))
>> > summary(fit3)
>>
>> However, the output of these three summary functions are identical. Why?
>> That option should have changed the base, but apparently this is not the
>> case.
>>
>>
>> Another analysis I found online (at this link
>>
>> https://stats.stackexchange.com/questions/60352/comparing-levels-of-factors-after-a-glm-in-r
>> )
>> to understand the differences between the 3 levels is to use glth with
>> Tuckey. I performed the following:
>>
>> > library(multcomp)
>> > summary(glht(fit, mcp(System="Tukey")))
>>
>> Simultaneous Tests for General Linear Hypotheses
>>
>> Multiple Comparisons of Means: Tukey Contrasts
>>
>>
>> Fit: glm(formula = Response ~ System, family = "binomial", data = scrd)
>>
>> Linear Hypotheses:
>>   Estimate Std. Error z value Pr(>|z|)
>> B - A == 0  -1.2715 0.3379  -3.763 0.000445 ***
>> C - A == 00.8588 0.4990   1.721 0.192472
>> C - B == 0 2.1303 0.4512   4.722  < 1e-04 ***
>> ---
>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>> (Adjusted p values reported -- single-step method)
>>
>>
>> Is this Tukey analysis correct?
>>
>>
>> I am a bit confused on what analysis I should do. I am doing my very best
>> to study all resources I can find, but I would really need some help from
>> experts, especially in using R.
>>
>>
>> Best wishes
>>
>> FJ
>>
>>
>>
>>
>>
>>
>> On Mon, Nov 12, 2018 at 1:46 PM peter dalgaard  wrote:
>>
>> > Yes, only one of the pairwise comparisons (B vs. C) is right. Also, the
>> > overall test has 3 degrees of freedom whereas a comparison of 3 groups
>> > should have 2. You (meaning Frodo) are testing that _all 3_ regression
>> > coefficients are zero, intercept included. That would imply that all
>> three
>> > systems have response probablilities og 0.5, which is not likely what
>> you
>> > want.
>> >
>> > This all suggests that you are struggling with the interpretation of the
>> > regression coefficients and their role in the linear predictor. This
>> should
>> > be covered by any good book on logistic regression.
>> >
>> > -pd
>> >
>> > > On 12 Nov 2018, at 14:15 , Eik Vettorazzi 
>> wrote:
>> > >
>> > > Dear Jedi,
>> > > please use the source carefully. A and C are not statistically
>> different
>> > at the 5% level, which can be inferred from glm output. Your last two
>> > wald.tests don't test what you want to, since your model contains an
>> > intercept term. You specified contrasts which tests A vs B-A, ie A-
>> > (B-A)==0 <-> 2*A-B==0 which is not intended I think. Have a look at
>> > ?contr.treatment and re-read your source doc to get an idea what dummy
>> > coding and indicatr variables are about.
>> > >
>> > > Cheers
>> > >
>> > >
>> > > Am 

Re: [R] missRanger package

2018-11-12 Thread MacQueen, Don via R-help
I could not find the word "censor" in the documentation for the missRanger 
package, so I think additional explanation is needed.

Also, I would expect information about censoring to be included in data 
provided to a function in a package -- inserting censoring into a package 
doesn't make sense.

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 11/12/18, 6:13 AM, "R-help on behalf of Rebecca Bingert" 
 wrote:

Hi,
does anybody know where I need to insert the censoring in the missRanger 
package?
Regards,
Rebecca

__
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] saveRDS() and readRDS() Why? [solved, pretty much anyway)

2018-11-12 Thread p_connolly

On 2018-11-12 22:49, peter dalgaard wrote:
Er, where, what, how? I can't reproduce that, at least not on 3.5.1 on 
MacOS:



x <- airquality
saveRDS(x, file = "x.rds")
x <- NULL
x <- readRDS(file = "x.rds")
x

Ozone Solar.R Wind Temp Month Day
1  41 190  7.4   67 5   1
2  36 118  8.0   72 5   2
3  12 149 12.6   74 5   3
...

Looks fine to me.



It seems to work fine using the same installation to read as used for 
the save.
But it's a different story if the save was done on a Windows 
installation and read on a Linux installation.


## On Windows 3.4.2

x <- airquality
saveRDS(x, file = "x.rds")
saveRDS(x, file = "y.rds")


Files x.rds and y.rds are identical in size but utterly different in 
content.


## On Linux 3.5.1


x <- readRDS(file =  "x.rds")

Error in readRDS(file = "x.rds") : error reading from connection

x <- readRDS(file =  "y.rds")
head(x)

  Ozone Solar.R Wind Temp Month Day
141 190  7.4   67 5   1
236 118  8.0   72 5   2
312 149 12.6   74 5   3
418 313 11.5   62 5   4
5NA  NA 14.3   56 5   5
628  NA 14.9   66 5   6

It might just be the age of the Windows installation.  I don't have much 
use for Windows, so I haven't had much inclination to install a newer 
version.


YMMV




?
-pd


On 12 Nov 2018, at 08:28 , Patrick Connolly 
 wrote:


The solution was very simple.  Don't use the same name for the rds
file  as used for the R object, viz a vie:

saveRDS(x, file = "x.rds")
and
x <- readRDS(file = "x.rds")

will not work; however

saveRDS(x, file = "y.rds")
and
x <- readRDS(file = "y.rds")
will work.

An undocumented feature?

Thanks to all who contributed.


--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd@cbs.dk  Priv: pda...@gmail.com


__
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] semiparametric manova

2018-11-12 Thread Yectli Huerta via R-help
Hello,

I was wondering if there are other packages like MANOVA.RM that could be used 
to analysis non normal distributions. I have to analyze data with more than 2 
predictor variables and a similar number of response variables. When I try the 
function MANOVA.wide with more than 2 predictor variables, I get

There is at least one factor-level combination
   with less than 2 observations!

Is there another package out there that can be used to analyze the significance 
of more than 2 predictor variables?

thanks,

yah

signature.asc
Description: OpenPGP digital 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] missRanger package

2018-11-12 Thread Bert Gunter
You have asked what I believe is an incoherent question, and thus are
unlikely to receive any useful replies (of course, I may be wrong about
this...).

Please read and follow the posting guide linked below to to ask a question
that can be answered.


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 Mon, Nov 12, 2018 at 12:03 PM Rebecca Bingert 
wrote:

> Hi,
> does anybody know where I need to insert the censoring in the missRanger
> package?
> Regards,
> Rebecca
>
> __
> 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] missRanger package

2018-11-12 Thread Rebecca Bingert

Hi,
does anybody know where I need to insert the censoring in the missRanger 
package?

Regards,
Rebecca

__
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] Reporting binomial logistic regression from R results

2018-11-12 Thread Bert Gunter
Generally speaking, this list is about questions on R programming, not
statistical issues. However, I grant you that your queries are in something
of a gray area intersecting both.

Nevertheless, based on your admitted confusion, I would recommend that you
find a local statistical expert with whom you can consult 1-1 if at all
possible. As others have already noted, you statistical understanding is
muddy, and it can be quite difficult to resolve such confusion in online
forums like this that cannot provide the close back and forth that may be
required (as well as further appropriate study).

Best,
Bert

On Mon, Nov 12, 2018 at 11:09 AM Frodo Jedi 
wrote:

> Dear Peter and Eik,
> I am very grateful to you for your replies.
> My current understanding is that from the GLM analysis I can indeed
> conclude that the response predicted by System A is significantly different
> from that of System B, while the pairwise comparison A vs C leads to non
> significance. Now the Wald test seems to be correct only for Systems B vs
> C, indicating that the pairwise System B vs System C is significant. Am I
> correct?
>
> However, my current understanding is also that I should use contrasts
> instead of the wald test. So the default contrasts is with the System A,
> now I should re-perform the GLM with another base. I tried to use the
> option "contrasts" of the glm:
>
> > fit1 <- glm(Response ~ System, data = scrd, family = "binomial",
> contrasts = contr.treatment(3, base=1,contrasts=TRUE))
> > summary(fit1)
>
> > fit2 <- glm(Response ~ System, data = scrd, family = "binomial",
> contrasts = contr.treatment(3, base=2,contrasts=TRUE))
> > summary(fit2)
>
> > fit3 <- glm(Response ~ System, data = scrd, family = "binomial",
> contrasts = contr.treatment(3, base=3,contrasts=TRUE))
> > summary(fit3)
>
> However, the output of these three summary functions are identical. Why?
> That option should have changed the base, but apparently this is not the
> case.
>
>
> Another analysis I found online (at this link
>
> https://stats.stackexchange.com/questions/60352/comparing-levels-of-factors-after-a-glm-in-r
> )
> to understand the differences between the 3 levels is to use glth with
> Tuckey. I performed the following:
>
> > library(multcomp)
> > summary(glht(fit, mcp(System="Tukey")))
>
> Simultaneous Tests for General Linear Hypotheses
>
> Multiple Comparisons of Means: Tukey Contrasts
>
>
> Fit: glm(formula = Response ~ System, family = "binomial", data = scrd)
>
> Linear Hypotheses:
>   Estimate Std. Error z value Pr(>|z|)
> B - A == 0  -1.2715 0.3379  -3.763 0.000445 ***
> C - A == 00.8588 0.4990   1.721 0.192472
> C - B == 0 2.1303 0.4512   4.722  < 1e-04 ***
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> (Adjusted p values reported -- single-step method)
>
>
> Is this Tukey analysis correct?
>
>
> I am a bit confused on what analysis I should do. I am doing my very best
> to study all resources I can find, but I would really need some help from
> experts, especially in using R.
>
>
> Best wishes
>
> FJ
>
>
>
>
>
>
> On Mon, Nov 12, 2018 at 1:46 PM peter dalgaard  wrote:
>
> > Yes, only one of the pairwise comparisons (B vs. C) is right. Also, the
> > overall test has 3 degrees of freedom whereas a comparison of 3 groups
> > should have 2. You (meaning Frodo) are testing that _all 3_ regression
> > coefficients are zero, intercept included. That would imply that all
> three
> > systems have response probablilities og 0.5, which is not likely what you
> > want.
> >
> > This all suggests that you are struggling with the interpretation of the
> > regression coefficients and their role in the linear predictor. This
> should
> > be covered by any good book on logistic regression.
> >
> > -pd
> >
> > > On 12 Nov 2018, at 14:15 , Eik Vettorazzi  wrote:
> > >
> > > Dear Jedi,
> > > please use the source carefully. A and C are not statistically
> different
> > at the 5% level, which can be inferred from glm output. Your last two
> > wald.tests don't test what you want to, since your model contains an
> > intercept term. You specified contrasts which tests A vs B-A, ie A-
> > (B-A)==0 <-> 2*A-B==0 which is not intended I think. Have a look at
> > ?contr.treatment and re-read your source doc to get an idea what dummy
> > coding and indicatr variables are about.
> > >
> > > Cheers
> > >
> > >
> > > Am 12.11.2018 um 02:07 schrieb Frodo Jedi:
> > >> Dear list members,
> > >> I need some help in understanding whether I am doing correctly a
> > binomial
> > >> logistic regression and whether I am interpreting the results in the
> > >> correct way. Also I would need an advice regarding the reporting of
> the
> > >> results from the R functions.
> > >> I want to report the results of a binomial logistic regression where I
> > want
> > >> to assess difference between the 3 levels of a factor (called System)
> on
> > >> the dependent variable (called Response) taking two values, 

Re: [R] Reporting binomial logistic regression from R results

2018-11-12 Thread Frodo Jedi
Dear Peter and Eik,
I am very grateful to you for your replies.
My current understanding is that from the GLM analysis I can indeed
conclude that the response predicted by System A is significantly different
from that of System B, while the pairwise comparison A vs C leads to non
significance. Now the Wald test seems to be correct only for Systems B vs
C, indicating that the pairwise System B vs System C is significant. Am I
correct?

However, my current understanding is also that I should use contrasts
instead of the wald test. So the default contrasts is with the System A,
now I should re-perform the GLM with another base. I tried to use the
option "contrasts" of the glm:

> fit1 <- glm(Response ~ System, data = scrd, family = "binomial",
contrasts = contr.treatment(3, base=1,contrasts=TRUE))
> summary(fit1)

> fit2 <- glm(Response ~ System, data = scrd, family = "binomial",
contrasts = contr.treatment(3, base=2,contrasts=TRUE))
> summary(fit2)

> fit3 <- glm(Response ~ System, data = scrd, family = "binomial",
contrasts = contr.treatment(3, base=3,contrasts=TRUE))
> summary(fit3)

However, the output of these three summary functions are identical. Why?
That option should have changed the base, but apparently this is not the
case.


Another analysis I found online (at this link
https://stats.stackexchange.com/questions/60352/comparing-levels-of-factors-after-a-glm-in-r
)
to understand the differences between the 3 levels is to use glth with
Tuckey. I performed the following:

> library(multcomp)
> summary(glht(fit, mcp(System="Tukey")))

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = Response ~ System, family = "binomial", data = scrd)

Linear Hypotheses:
  Estimate Std. Error z value Pr(>|z|)
B - A == 0  -1.2715 0.3379  -3.763 0.000445 ***
C - A == 00.8588 0.4990   1.721 0.192472
C - B == 0 2.1303 0.4512   4.722  < 1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


Is this Tukey analysis correct?


I am a bit confused on what analysis I should do. I am doing my very best
to study all resources I can find, but I would really need some help from
experts, especially in using R.


Best wishes

FJ






On Mon, Nov 12, 2018 at 1:46 PM peter dalgaard  wrote:

> Yes, only one of the pairwise comparisons (B vs. C) is right. Also, the
> overall test has 3 degrees of freedom whereas a comparison of 3 groups
> should have 2. You (meaning Frodo) are testing that _all 3_ regression
> coefficients are zero, intercept included. That would imply that all three
> systems have response probablilities og 0.5, which is not likely what you
> want.
>
> This all suggests that you are struggling with the interpretation of the
> regression coefficients and their role in the linear predictor. This should
> be covered by any good book on logistic regression.
>
> -pd
>
> > On 12 Nov 2018, at 14:15 , Eik Vettorazzi  wrote:
> >
> > Dear Jedi,
> > please use the source carefully. A and C are not statistically different
> at the 5% level, which can be inferred from glm output. Your last two
> wald.tests don't test what you want to, since your model contains an
> intercept term. You specified contrasts which tests A vs B-A, ie A-
> (B-A)==0 <-> 2*A-B==0 which is not intended I think. Have a look at
> ?contr.treatment and re-read your source doc to get an idea what dummy
> coding and indicatr variables are about.
> >
> > Cheers
> >
> >
> > Am 12.11.2018 um 02:07 schrieb Frodo Jedi:
> >> Dear list members,
> >> I need some help in understanding whether I am doing correctly a
> binomial
> >> logistic regression and whether I am interpreting the results in the
> >> correct way. Also I would need an advice regarding the reporting of the
> >> results from the R functions.
> >> I want to report the results of a binomial logistic regression where I
> want
> >> to assess difference between the 3 levels of a factor (called System) on
> >> the dependent variable (called Response) taking two values, 0 and 1. My
> >> goal is to understand if the effect of the 3 systems (A,B,C) in System
> >> affect differently Response in a significant way. I am basing my
> analysis
> >> on this URL: https://stats.idre.ucla.edu/r/dae/logit-regression/
> >> This is the result of my analysis:
> >>> fit <- glm(Response ~ System, data = scrd, family = "binomial")
> >>> summary(fit)
> >> Call:
> >> glm(formula = Response ~ System, family = "binomial", data = scrd)
> >> Deviance Residuals:
> >> Min   1Q   Median   3Q  Max
> >> -2.8840   0.1775   0.2712   0.2712   0.5008
> >> Coefficients:
> >>  Estimate Std. Error z value Pr(>|z|)
> >> (Intercept)3.2844 0.2825  11.626  < 2e-16 ***
> >> SystemB  -1.2715 0.3379  -3.763 0.000168 ***
> >> SystemC0.8588 0.4990   1.721 0.085266 .
> >> ---
> >> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 

Re: [R-es] Materiales - V Xornada de Usuarios de R en Galicia

2018-11-12 Thread Carlos Ortega
Gracias!.

Parece que hemos cogido un nuevo impulso si cabe al cambiar de sitio a uno
más céntrico..a ver si se confirma el buen momento... Tuvimos entre 50-60
personas en la última reunión, veremos este miércoles...

No creo que pueda pasarme a las Jornadas en Murcia, pero si el año que
viene son en Madrid ya nos veremos...!.

Seguiremos compartiendo...

Gracias por todo!
Carlos.

El lun., 12 nov. 2018 a las 14:28, 
escribió:

> Gracias a vosotros, Carlos.
>
> Que los del grupo de Madrid sois hiperactivos y, gracias a eso, tenemos un
> montón de materiales.
>
> :-)
>
>
> Buen trabajo
>
> (y seguid así, porfa!)
>
>
>
> --
> *De:* Carlos Ortega 
> *Enviado:* lunes, 12 de noviembre de 2018 14:25
> *Para:* Rodríguez Muíños, Miguel Ángel
> *Cc:* Lista R
> *Asunto:* Re: [R-es] Materiales - V Xornada de Usuarios de R en Galicia
>
> Gracias Miguel Ángel!.
>
> El lun., 12 nov. 2018 a las 11:58, <
> miguel.angel.rodriguez.mui...@sergas.es> escribió:
>
>> Hola.
>>
>> Por si es de vuestro interés, os informo de que ya están disponibles las
>> presentaciones y los videos de la Jornada de usuarios de Galicia, que tuvo
>> lugar el 25 de Octubre en Santiago de Compostela.
>>
>> https://www.r-users.gal/pagina/programa-2018
>>
>>
>>
>> Un Saludo,
>> --
>> Miguel Ángel Rodríguez Muíños
>> Dirección Xeral de Saúde Pública
>> Consellería de Sanidade
>> Xunta de Galicia
>> http://dxsp.sergas.es
>>
>>
>>
> --
> Saludos,
> Carlos Ortega
> www.qualityexcellence.es
>
> --
>
> Nota: A información contida nesta mensaxe e os seus posibles documentos
> adxuntos é privada e confidencial e está dirixida únicamente ó seu
> destinatario/a. Se vostede non é o/a destinatario/a orixinal desta mensaxe,
> por favor elimínea. A distribución ou copia desta mensaxe non está
> autorizada.
>
> Nota: La información contenida en este mensaje y sus posibles documentos
> adjuntos es privada y confidencial y está dirigida únicamente a su
> destinatario/a. Si usted no es el/la destinatario/a original de este
> mensaje, por favor elimínelo. La distribución o copia de este mensaje no
> está autorizada.
>
> See more languages: http://www.sergas.es/aviso-confidencialidad
>


-- 
Saludos,
Carlos Ortega
www.qualityexcellence.es

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


Re: [R] Reporting binomial logistic regression from R results

2018-11-12 Thread peter dalgaard
Yes, only one of the pairwise comparisons (B vs. C) is right. Also, the overall 
test has 3 degrees of freedom whereas a comparison of 3 groups should have 2. 
You (meaning Frodo) are testing that _all 3_ regression coefficients are zero, 
intercept included. That would imply that all three systems have response 
probablilities og 0.5, which is not likely what you want.

This all suggests that you are struggling with the interpretation of the 
regression coefficients and their role in the linear predictor. This should be 
covered by any good book on logistic regression.

-pd  

> On 12 Nov 2018, at 14:15 , Eik Vettorazzi  wrote:
> 
> Dear Jedi,
> please use the source carefully. A and C are not statistically different at 
> the 5% level, which can be inferred from glm output. Your last two wald.tests 
> don't test what you want to, since your model contains an intercept term. You 
> specified contrasts which tests A vs B-A, ie A- (B-A)==0 <-> 2*A-B==0 which 
> is not intended I think. Have a look at ?contr.treatment and re-read your 
> source doc to get an idea what dummy coding and indicatr variables are about.
> 
> Cheers
> 
> 
> Am 12.11.2018 um 02:07 schrieb Frodo Jedi:
>> Dear list members,
>> I need some help in understanding whether I am doing correctly a binomial
>> logistic regression and whether I am interpreting the results in the
>> correct way. Also I would need an advice regarding the reporting of the
>> results from the R functions.
>> I want to report the results of a binomial logistic regression where I want
>> to assess difference between the 3 levels of a factor (called System) on
>> the dependent variable (called Response) taking two values, 0 and 1. My
>> goal is to understand if the effect of the 3 systems (A,B,C) in System
>> affect differently Response in a significant way. I am basing my analysis
>> on this URL: https://stats.idre.ucla.edu/r/dae/logit-regression/
>> This is the result of my analysis:
>>> fit <- glm(Response ~ System, data = scrd, family = "binomial")
>>> summary(fit)
>> Call:
>> glm(formula = Response ~ System, family = "binomial", data = scrd)
>> Deviance Residuals:
>> Min   1Q   Median   3Q  Max
>> -2.8840   0.1775   0.2712   0.2712   0.5008
>> Coefficients:
>>  Estimate Std. Error z value Pr(>|z|)
>> (Intercept)3.2844 0.2825  11.626  < 2e-16 ***
>> SystemB  -1.2715 0.3379  -3.763 0.000168 ***
>> SystemC0.8588 0.4990   1.721 0.085266 .
>> ---
>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>> (Dispersion parameter for binomial family taken to be 1)
>> Null deviance: 411.26  on 1023  degrees of freedom
>> Residual deviance: 376.76  on 1021  degrees of freedom
>> AIC: 382.76
>> Number of Fisher Scoring iterations: 6
>> Following this analysis I perform the wald test in order to understand
>> whether there is an overall effect of System:
>> library(aod)
>>> wald.test(b = coef(fit), Sigma = vcov(fit), Terms = 1:3)
>> Wald test:
>> --
>> Chi-squared test:
>> X2 = 354.6, df = 3, P(> X2) = 0.0
>> The chi-squared test statistic of 354.6, with 3 degrees of freedom is
>> associated with a p-value < 0.001 indicating that the overall effect of
>> System is statistically significant.
>> Now I check whether there are differences between the coefficients using
>> again the wald test:
>> # Here difference between system B and C:
>>> l <- cbind(0, 1, -1)
>>> wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
>> Wald test:
>> --
>> Chi-squared test:
>> X2 = 22.3, df = 1, P(> X2) = 2.3e-06
>> # Here difference between system A and C:
>>> l <- cbind(1, 0, -1)
>>> wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
>> Wald test:
>> --
>> Chi-squared test:
>> X2 = 12.0, df = 1, P(> X2) = 0.00052
>> # Here difference between system A and B:
>>> l <- cbind(1, -1, 0)
>>> wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
>> Wald test:
>> --
>> Chi-squared test:
>> X2 = 58.7, df = 1, P(> X2) = 1.8e-14
>> My understanding is that from this analysis I can state that the three
>> systems lead to a significantly different Response. Am I right? If so, how
>> should I report the results of this analysis? What is the correct way?
>> Thanks in advance
>> Best wishes
>> FJ
>>  [[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.
> 
> -- 
> Eik Vettorazzi
> 
> Department of Medical Biometry and Epidemiology
> University Medical Center Hamburg-Eppendorf
> 
> Martinistrasse 52
> building W 34
> 20246 Hamburg
> 
> Phone: +49 (0) 40 7410 - 58243
> Fax:   +49 (0) 40 7410 - 57790
> Web: www.uke.de/imbe
> --
> 
> _
> 
> 

Re: [R-es] Materiales - V Xornada de Usuarios de R en Galicia

2018-11-12 Thread miguel.angel.rodriguez.muinos
Gracias a vosotros, Carlos.

Que los del grupo de Madrid sois hiperactivos y, gracias a eso, tenemos un 
mont�n de materiales.

:-)


Buen trabajo

(y seguid as�, porfa!)




De: Carlos Ortega 
Enviado: lunes, 12 de noviembre de 2018 14:25
Para: Rodr�guez Mu��os, Miguel �ngel
Cc: Lista R
Asunto: Re: [R-es] Materiales - V Xornada de Usuarios de R en Galicia

Gracias Miguel �ngel!.

El lun., 12 nov. 2018 a las 11:58, 
mailto:miguel.angel.rodriguez.mui...@sergas.es>>
 escribi�:
Hola.

Por si es de vuestro inter�s, os informo de que ya est�n disponibles las 
presentaciones y los videos de la Jornada de usuarios de Galicia, que tuvo 
lugar el 25 de Octubre en Santiago de Compostela.

https://www.r-users.gal/pagina/programa-2018



Un Saludo,
--
Miguel �ngel Rodr�guez Mu��os
Direcci�n Xeral de Sa�de P�blica
Conseller�a de Sanidade
Xunta de Galicia
http://dxsp.sergas.es



--
Saludos,
Carlos Ortega
www.qualityexcellence.es



Nota: A informaci�n contida nesta mensaxe e os seus posibles documentos 
adxuntos � privada e confidencial e est� dirixida �nicamente � seu 
destinatario/a. Se vostede non � o/a destinatario/a orixinal desta mensaxe, por 
favor elim�nea. A distribuci�n ou copia desta mensaxe non est� autorizada.

Nota: La informaci�n contenida en este mensaje y sus posibles documentos 
adjuntos es privada y confidencial y est� dirigida �nicamente a su 
destinatario/a. Si usted no es el/la destinatario/a original de este mensaje, 
por favor elim�nelo. La distribuci�n o copia de este mensaje no est� autorizada.

See more languages: http://www.sergas.es/aviso-confidencialidad

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


Re: [R] Reporting binomial logistic regression from R results

2018-11-12 Thread Eik Vettorazzi

Dear Jedi,
please use the source carefully. A and C are not statistically different 
at the 5% level, which can be inferred from glm output. Your last two 
wald.tests don't test what you want to, since your model contains an 
intercept term. You specified contrasts which tests A vs B-A, ie A- 
(B-A)==0 <-> 2*A-B==0 which is not intended I think. Have a look at 
?contr.treatment and re-read your source doc to get an idea what dummy 
coding and indicatr variables are about.


Cheers


Am 12.11.2018 um 02:07 schrieb Frodo Jedi:

Dear list members,
I need some help in understanding whether I am doing correctly a binomial
logistic regression and whether I am interpreting the results in the
correct way. Also I would need an advice regarding the reporting of the
results from the R functions.

I want to report the results of a binomial logistic regression where I want
to assess difference between the 3 levels of a factor (called System) on
the dependent variable (called Response) taking two values, 0 and 1. My
goal is to understand if the effect of the 3 systems (A,B,C) in System
affect differently Response in a significant way. I am basing my analysis
on this URL: https://stats.idre.ucla.edu/r/dae/logit-regression/

This is the result of my analysis:


fit <- glm(Response ~ System, data = scrd, family = "binomial")
summary(fit)


Call:
glm(formula = Response ~ System, family = "binomial", data = scrd)

Deviance Residuals:
 Min   1Q   Median   3Q  Max
-2.8840   0.1775   0.2712   0.2712   0.5008

Coefficients:
  Estimate Std. Error z value Pr(>|z|)
(Intercept)3.2844 0.2825  11.626  < 2e-16 ***
SystemB  -1.2715 0.3379  -3.763 0.000168 ***
SystemC0.8588 0.4990   1.721 0.085266 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

 Null deviance: 411.26  on 1023  degrees of freedom
Residual deviance: 376.76  on 1021  degrees of freedom
AIC: 382.76

Number of Fisher Scoring iterations: 6
Following this analysis I perform the wald test in order to understand
whether there is an overall effect of System:

library(aod)


wald.test(b = coef(fit), Sigma = vcov(fit), Terms = 1:3)

Wald test:
--

Chi-squared test:
X2 = 354.6, df = 3, P(> X2) = 0.0
The chi-squared test statistic of 354.6, with 3 degrees of freedom is
associated with a p-value < 0.001 indicating that the overall effect of
System is statistically significant.

Now I check whether there are differences between the coefficients using
again the wald test:

# Here difference between system B and C:


l <- cbind(0, 1, -1)
wald.test(b = coef(fit), Sigma = vcov(fit), L = l)

Wald test:
--

Chi-squared test:
X2 = 22.3, df = 1, P(> X2) = 2.3e-06



# Here difference between system A and C:


l <- cbind(1, 0, -1)
wald.test(b = coef(fit), Sigma = vcov(fit), L = l)

Wald test:
--

Chi-squared test:
X2 = 12.0, df = 1, P(> X2) = 0.00052



# Here difference between system A and B:


l <- cbind(1, -1, 0)
wald.test(b = coef(fit), Sigma = vcov(fit), L = l)

Wald test:
--

Chi-squared test:
X2 = 58.7, df = 1, P(> X2) = 1.8e-14

My understanding is that from this analysis I can state that the three
systems lead to a significantly different Response. Am I right? If so, how
should I report the results of this analysis? What is the correct way?


Thanks in advance

Best wishes

FJ

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



--
Eik Vettorazzi

Department of Medical Biometry and Epidemiology
University Medical Center Hamburg-Eppendorf

Martinistrasse 52
building W 34
20246 Hamburg

Phone: +49 (0) 40 7410 - 58243
Fax:   +49 (0) 40 7410 - 57790
Web: www.uke.de/imbe
--

_

Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; 
Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard Göke (Vorsitzender), Prof. Dr. Dr. Uwe 
Koch-Gromus, Joachim Prölß, Marya Verdel
_

SAVE PAPER - THINK BEFORE PRINTING
__
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-es] Materiales - V Xornada de Usuarios de R en Galicia

2018-11-12 Thread Carlos Ortega
Gracias Miguel Ángel!.

El lun., 12 nov. 2018 a las 11:58, 
escribió:

> Hola.
>
> Por si es de vuestro interés, os informo de que ya están disponibles las
> presentaciones y los videos de la Jornada de usuarios de Galicia, que tuvo
> lugar el 25 de Octubre en Santiago de Compostela.
>
> https://www.r-users.gal/pagina/programa-2018
>
>
>
> Un Saludo,
> --
> Miguel Ángel Rodríguez Muíños
> Dirección Xeral de Saúde Pública
> Consellería de Sanidade
> Xunta de Galicia
> http://dxsp.sergas.es
>
>
>
>
>
>
>
>
>
> 
>
> Nota: A información contida nesta mensaxe e os seus posibles documentos
> adxuntos é privada e confidencial e está dirixida únicamente ó seu
> destinatario/a. Se vostede non é o/a destinatario/a orixinal desta mensaxe,
> por favor elimínea. A distribución ou copia desta mensaxe non está
> autorizada.
>
> Nota: La información contenida en este mensaje y sus posibles documentos
> adjuntos es privada y confidencial y está dirigida únicamente a su
> destinatario/a. Si usted no es el/la destinatario/a original de este
> mensaje, por favor elimínelo. La distribución o copia de este mensaje no
> está autorizada.
>
> See more languages: http://www.sergas.es/aviso-confidencialidad
>
> ___
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>


-- 
Saludos,
Carlos Ortega
www.qualityexcellence.es

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


Re: [R-es] curso

2018-11-12 Thread miguel.angel.rodriguez.muinos
Hola David.

Creo que te confundes de destinatario.
Esta es la lista R-help-es (R en castellano) y no organiza cursos.
Seguramente te has anotado a algún curso (Coursera, quizás?) y han referenciado 
esta lista pero no tiene nada que ver con la organización del curso.

Un saludo.



De: R-help-es  en nombre de David Montes 

Enviado: lunes, 12 de noviembre de 2018 14:09
Para: r-help-es@r-project.org
Asunto: [R-es] curso

Buenas tardes:
No puedo pagar curso, he hecho ejercicios y he visto todos los videos, anularme 
si quereis si sacais otra vez curso cuando trabaje si me gustaria pagar y 
sacarme el titulo, gracias.

Atentamente David Montes Navarro.

Good day:

I can not pay for the course, but you can also watch all the videos, cancel and 
want if you take one more course when I work and I would like to pay and get 
the title, thanks.

Sincerely David Montes Navarro.


[https://ipmcdn.avast.com/images/icons/icon-envelope-tick-round-orange-animated-no-repeat-v1.gif]
  Libre de virus. 
www.avast.com

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es



Nota: A información contida nesta mensaxe e os seus posibles documentos 
adxuntos é privada e confidencial e está dirixida únicamente ó seu 
destinatario/a. Se vostede non é o/a destinatario/a orixinal desta mensaxe, por 
favor elimínea. A distribución ou copia desta mensaxe non está autorizada.

Nota: La información contenida en este mensaje y sus posibles documentos 
adjuntos es privada y confidencial y está dirigida únicamente a su 
destinatario/a. Si usted no es el/la destinatario/a original de este mensaje, 
por favor elimínelo. La distribución o copia de este mensaje no está autorizada.

See more languages: http://www.sergas.es/aviso-confidencialidad

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


[R-es] curso

2018-11-12 Thread David Montes
Buenas tardes:
No puedo pagar curso, he hecho ejercicios y he visto todos los videos, anularme 
si quereis si sacais otra vez curso cuando trabaje si me gustaria pagar y 
sacarme el titulo, gracias.

Atentamente David Montes Navarro.

Good day:

I can not pay for the course, but you can also watch all the videos, cancel and 
want if you take one more course when I work and I would like to pay and get 
the title, thanks.

Sincerely David Montes Navarro.


[https://ipmcdn.avast.com/images/icons/icon-envelope-tick-round-orange-animated-no-repeat-v1.gif]
  Libre de virus. 
www.avast.com

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


Re: [R] Reporting binomial logistic regression from R results

2018-11-12 Thread Frodo Jedi
Dear Petr,
thank you very much for your feedback.

Can anyone in the list advise me if the way I report the results is correct?

Kind regards

FJ


On Mon, Nov 12, 2018 at 1:02 PM PIKAL Petr  wrote:

> Hi Frodo
>
>
>
> I do not consider myself as an arbiter in statistical results and their
> presentation. Again your text seems to as good as any other.
>
>
>
> You should keep responses to mailing list as others could have another
> opinion.
>
>
>
> Cheers
>
> Petr
>
>
>
>
>
> *From:* Frodo Jedi 
> *Sent:* Monday, November 12, 2018 1:48 PM
> *To:* PIKAL Petr 
> *Subject:* Re: [R] Reporting binomial logistic regression from R results
>
>
>
> Dear Petr,
>
> many thanks for your reply. I was wondering whether in your opinion it is
> correct to report in a journal the following text:
>
>
>
>
>
> “A logistic regression was performed to ascertain the effects of the
> system type on the likelihood that participants report correct
> identifications. The logistic regression model was statistically
> significant, χ2(3) = 354.6, p < 0.001, indicating an overall effect of the
> system type on participants' identification performances. The Wald test was
> used to compare the model coefficients related to the three systems.
> Results showed that participants’ accuracy was significantly lower for the
> system B compared to both the system C (χ2(1) = 22.3, p < 0.001) and the
> system A (χ2(1) = 58.7, p < 0.001), as well as that the system C led to
> significantly higher identification accuracies than the system A (χ2(1) =
> 12, p < 0.001).”
>
>
>
>
>
> Best wishes
>
>
>
> FJ
>
>
>
>
>
>
>
>
>
>
>
> On Mon, Nov 12, 2018 at 10:05 AM PIKAL Petr 
> wrote:
>
> Dear Frodo (or Jedi)
>
> The results seems to confirm your assumption that 3 systems are different.
> How you should present results probably depends on how it is usual to
> report such results in your environment.
>
> BTW. It seems to me like homework and this list has no homework policy
> (Sorry, if I am mistaken).
>
> Cheers
> Petr
> > -Original Message-
> > From: R-help  On Behalf Of Frodo Jedi
> > Sent: Monday, November 12, 2018 2:08 AM
> > To: r-help@r-project.org
> > Subject: [R] Reporting binomial logistic regression from R results
> >
> > Dear list members,
> > I need some help in understanding whether I am doing correctly a binomial
> > logistic regression and whether I am interpreting the results in the
> correct way.
> > Also I would need an advice regarding the reporting of the results from
> the R
> > functions.
> >
> > I want to report the results of a binomial logistic regression where I
> want to
> > assess difference between the 3 levels of a factor (called System) on the
> > dependent variable (called Response) taking two values, 0 and 1. My goal
> is to
> > understand if the effect of the 3 systems (A,B,C) in System affect
> differently
> > Response in a significant way. I am basing my analysis on this URL:
> > https://stats.idre.ucla.edu/r/dae/logit-regression/
> >
> > This is the result of my analysis:
> >
> > > fit <- glm(Response ~ System, data = scrd, family = "binomial")
> > > summary(fit)
> >
> > Call:
> > glm(formula = Response ~ System, family = "binomial", data = scrd)
> >
> > Deviance Residuals:
> > Min   1Q   Median   3Q  Max
> > -2.8840   0.1775   0.2712   0.2712   0.5008
> >
> > Coefficients:
> >  Estimate Std. Error z value Pr(>|z|)
> > (Intercept)3.2844 0.2825  11.626  < 2e-16 ***
> > SystemB  -1.2715 0.3379  -3.763 0.000168 ***
> > SystemC0.8588 0.4990   1.721 0.085266 .
> > ---
> > Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> >
> > (Dispersion parameter for binomial family taken to be 1)
> >
> > Null deviance: 411.26  on 1023  degrees of freedom Residual deviance:
> > 376.76  on 1021  degrees of freedom
> > AIC: 382.76
> >
> > Number of Fisher Scoring iterations: 6
> > Following this analysis I perform the wald test in order to understand
> whether
> > there is an overall effect of System:
> >
> > library(aod)
> >
> > > wald.test(b = coef(fit), Sigma = vcov(fit), Terms = 1:3)
> > Wald test:
> > --
> >
> > Chi-squared test:
> > X2 = 354.6, df = 3, P(> X2) = 0.0
> > The chi-squared test statistic of 354.6, with 3 degrees of freedom is
> associated
> > with a p-value < 0.001 indicating that the overall effect of System is
> statistically
> > significant.
> >
> > Now I check whether there are differences between the coefficients using
> again
> > the wald test:
> >
> > # Here difference between system B and C:
> >
> > > l <- cbind(0, 1, -1)
> > > wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
> > Wald test:
> > --
> >
> > Chi-squared test:
> > X2 = 22.3, df = 1, P(> X2) = 2.3e-06
> >
> >
> >
> > # Here difference between system A and C:
> >
> > > l <- cbind(1, 0, -1)
> > > wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
> > Wald test:
> > --
> >
> > Chi-squared test:
> > X2 = 12.0, df = 1, P(> X2) = 0.00052
> >
> >
> >
> > # Here 

Re: [R] Reporting binomial logistic regression from R results

2018-11-12 Thread PIKAL Petr
Hi Frodo

I do not consider myself as an arbiter in statistical results and their 
presentation. Again your text seems to as good as any other.

You should keep responses to mailing list as others could have another opinion.

Cheers
Petr


From: Frodo Jedi 
Sent: Monday, November 12, 2018 1:48 PM
To: PIKAL Petr 
Subject: Re: [R] Reporting binomial logistic regression from R results

Dear Petr,
many thanks for your reply. I was wondering whether in your opinion it is 
correct to report in a journal the following text:


“A logistic regression was performed to ascertain the effects of the system 
type on the likelihood that participants report correct identifications. The 
logistic regression model was statistically significant, χ2(3) = 354.6, p < 
0.001, indicating an overall effect of the system type on participants' 
identification performances. The Wald test was used to compare the model 
coefficients related to the three systems. Results showed that participants’ 
accuracy was significantly lower for the system B compared to both the system C 
(χ2(1) = 22.3, p < 0.001) and the system A (χ2(1) = 58.7, p < 0.001), as well 
as that the system C led to significantly higher identification accuracies than 
the system A (χ2(1) = 12, p < 0.001).”


Best wishes

FJ





On Mon, Nov 12, 2018 at 10:05 AM PIKAL Petr 
mailto:petr.pi...@precheza.cz>> wrote:
Dear Frodo (or Jedi)

The results seems to confirm your assumption that 3 systems are different. How 
you should present results probably depends on how it is usual to report such 
results in your environment.

BTW. It seems to me like homework and this list has no homework policy (Sorry, 
if I am mistaken).

Cheers
Petr
> -Original Message-
> From: R-help 
> mailto:r-help-boun...@r-project.org>> On Behalf 
> Of Frodo Jedi
> Sent: Monday, November 12, 2018 2:08 AM
> To: r-help@r-project.org
> Subject: [R] Reporting binomial logistic regression from R results
>
> Dear list members,
> I need some help in understanding whether I am doing correctly a binomial
> logistic regression and whether I am interpreting the results in the correct 
> way.
> Also I would need an advice regarding the reporting of the results from the R
> functions.
>
> I want to report the results of a binomial logistic regression where I want to
> assess difference between the 3 levels of a factor (called System) on the
> dependent variable (called Response) taking two values, 0 and 1. My goal is to
> understand if the effect of the 3 systems (A,B,C) in System affect differently
> Response in a significant way. I am basing my analysis on this URL:
> https://stats.idre.ucla.edu/r/dae/logit-regression/
>
> This is the result of my analysis:
>
> > fit <- glm(Response ~ System, data = scrd, family = "binomial")
> > summary(fit)
>
> Call:
> glm(formula = Response ~ System, family = "binomial", data = scrd)
>
> Deviance Residuals:
> Min   1Q   Median   3Q  Max
> -2.8840   0.1775   0.2712   0.2712   0.5008
>
> Coefficients:
>  Estimate Std. Error z value Pr(>|z|)
> (Intercept)3.2844 0.2825  11.626  < 2e-16 ***
> SystemB  -1.2715 0.3379  -3.763 0.000168 ***
> SystemC0.8588 0.4990   1.721 0.085266 .
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> (Dispersion parameter for binomial family taken to be 1)
>
> Null deviance: 411.26  on 1023  degrees of freedom Residual deviance:
> 376.76  on 1021  degrees of freedom
> AIC: 382.76
>
> Number of Fisher Scoring iterations: 6
> Following this analysis I perform the wald test in order to understand whether
> there is an overall effect of System:
>
> library(aod)
>
> > wald.test(b = coef(fit), Sigma = vcov(fit), Terms = 1:3)
> Wald test:
> --
>
> Chi-squared test:
> X2 = 354.6, df = 3, P(> X2) = 0.0
> The chi-squared test statistic of 354.6, with 3 degrees of freedom is 
> associated
> with a p-value < 0.001 indicating that the overall effect of System is 
> statistically
> significant.
>
> Now I check whether there are differences between the coefficients using again
> the wald test:
>
> # Here difference between system B and C:
>
> > l <- cbind(0, 1, -1)
> > wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
> Wald test:
> --
>
> Chi-squared test:
> X2 = 22.3, df = 1, P(> X2) = 2.3e-06
>
>
>
> # Here difference between system A and C:
>
> > l <- cbind(1, 0, -1)
> > wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
> Wald test:
> --
>
> Chi-squared test:
> X2 = 12.0, df = 1, P(> X2) = 0.00052
>
>
>
> # Here difference between system A and B:
>
> > l <- cbind(1, -1, 0)
> > wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
> Wald test:
> --
>
> Chi-squared test:
> X2 = 58.7, df = 1, P(> X2) = 1.8e-14
>
> My understanding is that from this analysis I can state that the three systems
> lead to a significantly different Response. Am I right? If so, how should I 
> report
> the results of this analysis? What is the 

Re: [R-es] Materiales - V Xornada de Usuarios de R en Galicia

2018-11-12 Thread José Luis Cañadas
Gracias Miguel Angel

El lun., 12 nov. 2018 a las 12:08, ceveve () escribió:

> Mil gracias.
>
> El lun., 12 nov. 2018 a las 4:58,  >
> escribió:
>
> > Hola.
> >
> > Por si es de vuestro interés, os informo de que ya están disponibles las
> > presentaciones y los videos de la Jornada de usuarios de Galicia, que
> tuvo
> > lugar el 25 de Octubre en Santiago de Compostela.
> >
> > https://www.r-users.gal/pagina/programa-2018
> >
> >
> >
> > Un Saludo,
> > --
> > Miguel Ángel Rodríguez Muíños
> > Dirección Xeral de Saúde Pública
> > Consellería de Sanidade
> > Xunta de Galicia
> > http://dxsp.sergas.es
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > 
> >
> > Nota: A información contida nesta mensaxe e os seus posibles documentos
> > adxuntos é privada e confidencial e está dirixida únicamente ó seu
> > destinatario/a. Se vostede non é o/a destinatario/a orixinal desta
> mensaxe,
> > por favor elimínea. A distribución ou copia desta mensaxe non está
> > autorizada.
> >
> > Nota: La información contenida en este mensaje y sus posibles documentos
> > adjuntos es privada y confidencial y está dirigida únicamente a su
> > destinatario/a. Si usted no es el/la destinatario/a original de este
> > mensaje, por favor elimínelo. La distribución o copia de este mensaje no
> > está autorizada.
> >
> > See more languages: http://www.sergas.es/aviso-confidencialidad
> >
> > ___
> > R-help-es mailing list
> > R-help-es@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-help-es
> >
>
>
> --
> *César O. Velázquez Vega**.*
>
> [[alternative HTML version deleted]]
>
> ___
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


Re: [R-es] Materiales - V Xornada de Usuarios de R en Galicia

2018-11-12 Thread ceveve
Mil gracias.

El lun., 12 nov. 2018 a las 4:58, 
escribió:

> Hola.
>
> Por si es de vuestro interés, os informo de que ya están disponibles las
> presentaciones y los videos de la Jornada de usuarios de Galicia, que tuvo
> lugar el 25 de Octubre en Santiago de Compostela.
>
> https://www.r-users.gal/pagina/programa-2018
>
>
>
> Un Saludo,
> --
> Miguel Ángel Rodríguez Muíños
> Dirección Xeral de Saúde Pública
> Consellería de Sanidade
> Xunta de Galicia
> http://dxsp.sergas.es
>
>
>
>
>
>
>
>
>
> 
>
> Nota: A información contida nesta mensaxe e os seus posibles documentos
> adxuntos é privada e confidencial e está dirixida únicamente ó seu
> destinatario/a. Se vostede non é o/a destinatario/a orixinal desta mensaxe,
> por favor elimínea. A distribución ou copia desta mensaxe non está
> autorizada.
>
> Nota: La información contenida en este mensaje y sus posibles documentos
> adjuntos es privada y confidencial y está dirigida únicamente a su
> destinatario/a. Si usted no es el/la destinatario/a original de este
> mensaje, por favor elimínelo. La distribución o copia de este mensaje no
> está autorizada.
>
> See more languages: http://www.sergas.es/aviso-confidencialidad
>
> ___
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>


-- 
*César O. Velázquez Vega**.*

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


[R-es] Materiales - V Xornada de Usuarios de R en Galicia

2018-11-12 Thread miguel.angel.rodriguez.muinos
Hola.

Por si es de vuestro interés, os informo de que ya están disponibles las 
presentaciones y los videos de la Jornada de usuarios de Galicia, que tuvo 
lugar el 25 de Octubre en Santiago de Compostela.

https://www.r-users.gal/pagina/programa-2018



Un Saludo,
--
Miguel Ángel Rodríguez Muíños
Dirección Xeral de Saúde Pública
Consellería de Sanidade
Xunta de Galicia
http://dxsp.sergas.es











Nota: A información contida nesta mensaxe e os seus posibles documentos 
adxuntos é privada e confidencial e está dirixida únicamente ó seu 
destinatario/a. Se vostede non é o/a destinatario/a orixinal desta mensaxe, por 
favor elimínea. A distribución ou copia desta mensaxe non está autorizada.

Nota: La información contenida en este mensaje y sus posibles documentos 
adjuntos es privada y confidencial y está dirigida únicamente a su 
destinatario/a. Si usted no es el/la destinatario/a original de este mensaje, 
por favor elimínelo. La distribución o copia de este mensaje no está autorizada.

See more languages: http://www.sergas.es/aviso-confidencialidad

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


[R-es] 55 Reunión Grupo Usuarios de R de Madrid...

2018-11-12 Thread Carlos Ortega
Buenas a todos,

Para los que estéis en Madrid el próximo miércoles tendremos la siguiente
reunión del Grupo de R de Madrid.

https://www.meetup.com/es-ES/Grupo-de-Usuarios-de-R-de-Madrid/events/256121363/

Y para los que no puedan pasarse, el material: presentaciones y videos,
estarán disponibles a la máxima brevedad.

Gracias!
Carlos Ortega
www.qualityexcellence.es

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


Re: [R] Reporting binomial logistic regression from R results

2018-11-12 Thread PIKAL Petr
Dear Frodo (or Jedi)

The results seems to confirm your assumption that 3 systems are different. How 
you should present results probably depends on how it is usual to report such 
results in your environment.

BTW. It seems to me like homework and this list has no homework policy (Sorry, 
if I am mistaken).

Cheers
Petr
> -Original Message-
> From: R-help  On Behalf Of Frodo Jedi
> Sent: Monday, November 12, 2018 2:08 AM
> To: r-help@r-project.org
> Subject: [R] Reporting binomial logistic regression from R results
>
> Dear list members,
> I need some help in understanding whether I am doing correctly a binomial
> logistic regression and whether I am interpreting the results in the correct 
> way.
> Also I would need an advice regarding the reporting of the results from the R
> functions.
>
> I want to report the results of a binomial logistic regression where I want to
> assess difference between the 3 levels of a factor (called System) on the
> dependent variable (called Response) taking two values, 0 and 1. My goal is to
> understand if the effect of the 3 systems (A,B,C) in System affect differently
> Response in a significant way. I am basing my analysis on this URL:
> https://stats.idre.ucla.edu/r/dae/logit-regression/
>
> This is the result of my analysis:
>
> > fit <- glm(Response ~ System, data = scrd, family = "binomial")
> > summary(fit)
>
> Call:
> glm(formula = Response ~ System, family = "binomial", data = scrd)
>
> Deviance Residuals:
> Min   1Q   Median   3Q  Max
> -2.8840   0.1775   0.2712   0.2712   0.5008
>
> Coefficients:
>  Estimate Std. Error z value Pr(>|z|)
> (Intercept)3.2844 0.2825  11.626  < 2e-16 ***
> SystemB  -1.2715 0.3379  -3.763 0.000168 ***
> SystemC0.8588 0.4990   1.721 0.085266 .
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> (Dispersion parameter for binomial family taken to be 1)
>
> Null deviance: 411.26  on 1023  degrees of freedom Residual deviance:
> 376.76  on 1021  degrees of freedom
> AIC: 382.76
>
> Number of Fisher Scoring iterations: 6
> Following this analysis I perform the wald test in order to understand whether
> there is an overall effect of System:
>
> library(aod)
>
> > wald.test(b = coef(fit), Sigma = vcov(fit), Terms = 1:3)
> Wald test:
> --
>
> Chi-squared test:
> X2 = 354.6, df = 3, P(> X2) = 0.0
> The chi-squared test statistic of 354.6, with 3 degrees of freedom is 
> associated
> with a p-value < 0.001 indicating that the overall effect of System is 
> statistically
> significant.
>
> Now I check whether there are differences between the coefficients using again
> the wald test:
>
> # Here difference between system B and C:
>
> > l <- cbind(0, 1, -1)
> > wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
> Wald test:
> --
>
> Chi-squared test:
> X2 = 22.3, df = 1, P(> X2) = 2.3e-06
>
>
>
> # Here difference between system A and C:
>
> > l <- cbind(1, 0, -1)
> > wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
> Wald test:
> --
>
> Chi-squared test:
> X2 = 12.0, df = 1, P(> X2) = 0.00052
>
>
>
> # Here difference between system A and B:
>
> > l <- cbind(1, -1, 0)
> > wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
> Wald test:
> --
>
> Chi-squared test:
> X2 = 58.7, df = 1, P(> X2) = 1.8e-14
>
> My understanding is that from this analysis I can state that the three systems
> lead to a significantly different Response. Am I right? If so, how should I 
> report
> the results of this analysis? What is the correct way?
>
>
> Thanks in advance
>
> Best wishes
>
> FJ
>
> [[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.
Osobní údaje: Informace o zpracování a ochraně osobních údajů obchodních 
partnerů PRECHEZA a.s. jsou zveřejněny na: 
https://www.precheza.cz/zasady-ochrany-osobnich-udaju/ | Information about 
processing and protection of business partner’s personal data are available on 
website: https://www.precheza.cz/en/personal-data-protection-principles/
Důvěrnost: Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a 
podléhají tomuto právně závaznému prohláąení o vyloučení odpovědnosti: 
https://www.precheza.cz/01-dovetek/ | This email and any documents attached to 
it may be confidential and are subject to the legally binding disclaimer: 
https://www.precheza.cz/en/01-disclaimer/

__
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] saveRDS() and readRDS() Why? [solved, pretty much anyway)

2018-11-12 Thread Patrick Connolly
The solution was very simple.  Don't use the same name for the rds
file  as used for the R object, viz a vie:

saveRDS(x, file = "x.rds")
and
x <- readRDS(file = "x.rds")

will not work; however

saveRDS(x, file = "y.rds")
and
x <- readRDS(file = "y.rds")
will work.

An undocumented feature?

Thanks to all who contributed.



On Sat, 10-Nov-2018 at 08:48PM +1300, Patrick Connolly wrote:

|> On Thu, 08-Nov-2018 at 11:06AM +0100, Martin Maechler wrote:
|> 
|> |> > Patrick Connolly 
|> |> > on Thu, 8 Nov 2018 20:27:24 +1300 writes:
|> 
|> [...]
|> 
|> |> > 
|> |> > I still don't know Why, but I know How.
|> |> 
|> |> Hmm.. and nobody has been able to reproduce your problem, right?
|> |> 
|> |> IIUC, currently you are suggesting that [on Windows], if you do
|> |> 
|> |>   saveRDS(rawdata, file="rawdata.rds")
|> |> 
|> |> the resulting file is does not work withreadRDS()  on Linux.
|> |> What again are your R versions on the two platforms?
|> 
|> 
|> It's an old version on Windows.  I haven't used Windows R since then.
|> 
|> major  3  
|> minor  2.4
|> year   2016   
|> month  03 
|> day16 
|> 
|> 
|> I've tried R-3.5.0 and R-3.5.1 Linux versions.  The problem might be
|> entirely because of the ancient Windows version. 
|> 
|> 
|> |> 
|> |> Could you  dput() -- provide a (short if possible) version of rawdata 
where
|> |> that problem occurs ?
|> 
|> I can't make a smaller version of rawdata which comes from scraping a
|> non-public web address, but next week when I'm back where those
|> machines are, I'll try it with a small data frame which is
|> reproducible.
|> 
|> 
|> 
|> |> 
|> |> Best,
|> |> Martin
|> |> 
|> |> 
|> |> > 
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.   
|> |> >___Patrick Connolly   
|> |> >  {~._.~}   Great minds discuss ideas
|> |> >  _( Y )_ Average minds discuss events 
|> |> > (:_~*~_:)  Small minds discuss people  
|> |> >  (_)-(_)  . Eleanor Roosevelt
|> |> >   
|> |> > ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
|> |> > 
|> |> > __
|> |> > 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.
|> 
|> -- 
|> ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.   
|>___Patrick Connolly   
|>  {~._.~}   Great minds discuss ideas
|>  _( Y )_  Average minds discuss events 
|> (:_~*~_:)  Small minds discuss people  
|>  (_)-(_)   . Eleanor Roosevelt
|>
|> ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
|> 
|> __
|> 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.

-- 
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.   
   ___Patrick Connolly   
 {~._.~}   Great minds discuss ideas
 _( Y )_ Average minds discuss events 
(:_~*~_:)  Small minds discuss people  
 (_)-(_)  . Eleanor Roosevelt
  
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.

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