date:20151214

[R] Hexbin: Counting Bins That Meet Certain Criteria

2015-12-14 Thread Sidoti, Salvatore A.

Greetings!

Is there a way to count the bins in a hexbin plot that meet certain criteria? 
For instance, what if I wanted to count the bins (hexes) that have a datapoint 
density of some number x?

Thank you!

__
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] stopifnot with logical(0)

2015-12-14 Thread Martin Maechler

> "DM" == Duncan Murdoch 
> on Sat, 12 Dec 2015 09:05:04 -0500 writes:

DM> On 12/12/2015 8:44 AM, peter dalgaard wrote:
>> 
>>> On 12 Dec 2015, at 10:54 , Martin Maechler
>>>  wrote:
>>> 
>>> My conclusion: Breaking such a fundamental lemma of
>>> logic as "the empty set is always true"
>> 
>> Umm, that doesn't make sense to me. Surely you mean that
>> "an AND-operation over an empty index set is TRUE"? A
>> similar OR operation is FALSE, i.e. they behave like
>> empty products and sums, respectively.
>> 

DM> How about "the empty set is all true, and all false."

or, what the I *meant* with the above:

  "All statements about elements of the empty set are true"

  ((and I still like the short form, even though it is not
correct strictly logically/mathematically))

Of course, Peter is correct,  and that
   any(logical(0))   is  FALSE
is really the only sensical way.

__
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] matrix which results singular but at the same time positive definite

2015-12-14 Thread Stefano Sofia

Dear John and dear Peter,
I needed time to understand better some practical implications derived from 
your hints.

Thank you
Stefano


Da: Fox, John [j...@mcmaster.ca]
Inviato: giovedì 10 dicembre 2015 17.24
A: peter dalgaard; Stefano Sofia
Cc: r-help@r-project.org
Oggetto: RE: [R] matrix which results singular but at the same time positive 
definite

Dear Peter,

> -Original Message-
> From: peter dalgaard [mailto:pda...@gmail.com]
> Sent: Thursday, December 10, 2015 11:09 AM
> To: Stefano Sofia
> Cc: Fox, John; r-help@r-project.org
> Subject: Re: [R] matrix which results singular but at the same time
> positive definite
>
> Looks like the ill-conditioning is almost entirely due to scaling, e.g.

Yes, that's my point. Sorry I didn't make it clearer.

Best,
 John

>
> > eigen(cov2cor(m))
> $values
> [1] 1.7234899 1.3380701 0.6619299 0.2765101
> ...
>
> This is an annoyance in several parts of numerical linear algebra:
> Routines assume that R^n has all coordinates on a similar scale and
> therefore think that anything on the order of 1e-7 or so is effectively
> zero.
>
> Condition numbers do this too:
>
> > kappa(m)
> [1] 1.066582e+13
> > kappa(cov2cor(m))
> [1] 5.489243
>
>
> -pd
>
> On 10 Dec 2015, at 16:41 , Stefano Sofia
>  wrote:
>
> > Dear John,
> > thank you for your considerations.
> > This matrix (which is a variance matrix) is part of an algorithm for
> forward-filtering and backward-sampling of Dynamic Linear Models (West
> and Harrison, 1997), applied to DLM representation of ARIMA processes
> (Petris, Petrone, Campagnoli).  It is therefore very difficult to
> explain why this variance matrix becomes so ill conditioned. This
> already happens at the first iteration of the algorithm. I will try to
> work on initial conditions and some fixed parameters.
> >
> > Thank you again
> > Stefano
> >
> >
> > 
> > Da: Fox, John [j...@mcmaster.ca]
> > Inviato: giovedì 10 dicembre 2015 14.41
> > A: Stefano Sofia; r-help@r-project.org
> > Oggetto: RE: matrix which results singular but at the same time
> positivedefinite
> >
> > Dear Stefano,
> >
> > You've already had a couple of informative responses directly
> addressing your question, but are you aware how ill-conditioned the
> matrix is (one of the responses alluded to this)?
> >
> >> kappa(X, exact=TRUE)
> > [1] 7.313338e+12
> >
> >> eigen(X)$values
> > [1] 4.964711e+00 9.356881e-01 4.863392e-12 6.788344e-13
> >
> > Two of the variables have variances around 10^0 and the other two
> around 10^-12. Of course, you haven't said anything about the context,
> but does it really make sense to analyze the data on these scales?
> >
> > Best,
> > John
> >
> > -
> > John Fox, Professor
> > McMaster University
> > Hamilton, Ontario
> > Canada L8S 4M4
> > Web: socserv.mcmaster.ca/jfox
> >
> >
> >
> >
> >> -Original Message-
> >> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of
> Stefano Sofia
> >> Sent: December 10, 2015 5:08 AM
> >> To: r-help@r-project.org
> >> Subject: [R] matrix which results singular but at the same time
> positive definite
> >>
> >> Dear list users,
> >> through the "matrixcalc" package I am performing some checks of
> variance
> >> matrices (which must be positive definite).
> >> In this example, it happens that the matrix A here reported is
> singular but
> >> positive definite. Is it possible?
> >>
> >>  [,1]  [,2]  [,3]  [,4]
> >> [1,]  1.904255e-12 -1.904255e-12 -8.238960e-13 -1.240294e-12 [2,] -
> >> 1.904255e-12  3.637979e-12  1.364242e-12  1.818989e-12 [3,] -
> 8.238960e-13
> >> 1.364242e-12  4.809988e+00  7.742369e-01 [4,] -1.240294e-12
> 1.818989e-12
> >> 7.742369e-01  1.090411e+00
> >>
> >> print(is.non.singular.matrix(A, tol = 1e-18)) FALSE
> print(is.positive.definite(A,
> >> tol=1e-18)) TRUE
> >>
> >> Is there something wrong with this matrix?
> >> Any comment will be appreciated.
> >> Stefano
> >>
> >>
> >> 
> >>
> >> AVVISO IMPORTANTE: Questo messaggio di posta elettronica può
> contenere
> >> informazioni confidenziali, pertanto è destinato solo a persone
> autorizzate alla
> >> ricezione. I messaggi di posta elettronica per i client di Regione
> Marche
> >> possono contenere informazioni confidenziali e con privilegi legali.
> Se non si è il
> >> destinatario specificato, non leggere, copiare, inoltrare o
> archiviare questo
> >> messaggio. Se si è ricevuto questo messaggio per errore, inoltrarlo
> al mittente
> >> ed eliminarlo completamente dal sistema del proprio computer. Ai
> sensi
> >> dell'art. 6 della DGR n. 1394/2008 si segnala che, in caso di
> necessità ed
> >> urgenza, la risposta al presente messaggio di posta elettronica può
> essere
> >> visionata da persone estranee al destinatario.
> >> IMPORTANT NOTICE: This e-mail message is intended to be received only
> by
>

[R] R: forest plot metafor

2015-12-14 Thread Mario Petretta

Many tanks to Wolfgang Viechtbauer for his time and help.

The suggested code works very well

Mario Petretta.

__
Message: 11
Date: Sat, 12 Dec 2015 17:22:20 +0100
From: "Viechtbauer Wolfgang (STAT)"

To: "r-help@r-project.org" 
Subject: Re: [R] R:  forest plot metafor
Message-ID:
<077e31a57da26e46ab0d493c9966ac730f2460c...@um-mail4112.unimaas.nl>
Content-Type: text/plain; charset="us-ascii"

No, this is not possible. But you can just add the weights yourself. One 
difficulty here is that you want the weights to the right of the estimates and 
corresponding CIs and those are always placed at the right of the figure. One 
possibility would be to make the right margin large and then use mtext() to 
place the weights in that margin. An example:

library(metafor)
data(dat.bcg)
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
res <- rma(yi, vi, data=dat, method="FE")
wi <- formatC(weights(res), format="f", digits=1, width=4)
par(mar=c(5,4,4,6))
forest(res, xlim=c(-7,8), digits=c(2,0))
mtext(wi, side=4, at=13:1, line=4, las=2, adj=1)
par(xpd=TRUE)
abline(h=c(0,14))
par(xpd=FALSE)

Adjust as needed for your data.

Best,
Wolfgang

From: R-help [r-help-boun...@r-project.org] On Behalf Of Mario Petretta 
[petre...@unina.it]
Sent: Saturday, December 12, 2015 3:57 PM
To: r-help@r-project.org
Subject: [R] R:  forest plot metafor

Many thanks to Professor Michael Dewey for his time.

I apologize for the error about the claim to obtain weights from escalc and I 
realise that the weights are a function of the fit, not the data.

weights(res) is OK to extract the weights from the fitted object and to put 
them in the data  frame.

However, I again ask:

Using:

 forest(res, showweights=TRUE)

it is possible to directly order the columns, simply placing the weights after 
effect size and CI?

Sincerely

Mario Petretta

---
Questa e-mail è stata controllata per individuare virus con Avast antivirus.
https://www.avast.com/antivirus

__
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] [R-pkgs] New Package "radiomics": First and Second Order Matrix Statistics

2015-12-14 Thread joel carlson

Dear R Users,

I am pleased to inform you that my package, 'radiomics' is now
available on CRAN:
https://cran.r-project.org/web/packages/radiomics/


The radiomics package offers 4 classes of texture matrices, and
associated feature sets. These matrices can be used for image
classification, or for extraction of meaningful information from the
input image (as in the case of the research field of radiomics). The
texture matrices are:

 - Gray level co-occurrence matrix (GLCM)

 - Gray level run-length matrix (GLRLM)

 - Gray level size-zone matrix (GLSZM)

 - Multiple gray level size-zone matrix (MGLSZM)

Each of these matrices are calculated from an input grayscale matrix,
and return objects of class 'glcm', 'glrlm', 'glszm' or 'mglszm'.
Features of the matrices can be calculated using the calc_features()
function. calc_features() can also be called on an image matrix, and
will return first-order summary statistics.

The github page can be viewed
here:https://github.com/joelcarlson/radiomics


Lists of features, and an introduction to the package can be found in
the package vignette.

In-depth explanation of texture analysis matrices can be found here:

http://joelcarlson.me/2015/07/10/radiomics-package/

Any comments, contributions, and feedback are always welcome.

I hope you find this helpful!
Regards
---
Joel Carlson

[[alternative HTML version deleted]]

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

__
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] Problems when using estimable() (gmodels) on lme()-objects

2015-12-14 Thread CG Pettersson

R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree", W7-32

Dear all,
I just upgraded R and downloaded a new version of "gmodels", to use estimable() 
for extracting fixed effects from lme()-objects.
Now it seems like estimable() have changed its view of dimensions of my 
lme-objects and contrast matrix dimensions, as I get this error message:

> estimable(y_15,Mat15)
Error in estimable.default(y_15, Mat15) :
 Dimension of structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, : 
25x26, not compatible with no of parameters in y_15: 25
Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, : 
25x26, not compatible with no of parameters in y_15: 25
Dimension of 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, : 
25x26, not compatible with no of parameters in y_15: 25
Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, : 
25x26, not compatible with no of parameters in y_15: 25
Dimension of 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, : 
25x26, not compatible with no of parameters in y_15: 25
Dimension of 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, : 
25x26, not compatible with no of parameters in y_15: 25
Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, : 
25x26, not compatible with no of parameters in y_15: 25
Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>

Mat15 is made as a 25x25 contrast matrix, but the updated estimable() sees it 
as a 25x26.
Why, and what should I do about it?

I do think the change has occurred in estimable() as I tested one year old 
lme() objects and contrast matrix who worked perfectly then, giving the same 
error message when run with updated R.
Does estimable() see the row-names as the first column nowdays if I don�t 
specify anything else?

Best regards
CG Pettersson
Scientific Project Manager, PhD
__
Lantm�nnen Corporate R
Phone:  +46 10 556 19 85
Mobile: + 46 70 330 66 85
Email: cg.petters...@lantmannen.com
Visiting Address: S:t G�ransgatan 160 A
Address: Box 30192, SE-104 25 Stockholm
Webb: http://www.lantmannen.com
Registered Office: Stockholm
Before printing, think about the environment


[[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] Problems using estimable() (gmodels) on lme()-objects

2015-12-14 Thread CG Pettersson

R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree", W7-32



Dear all,

I just upgraded R and downloaded a new version of “gmodels”, to use
estimable() for extracting fixed effects from lme()-objects.

Now it seems like estimable() have changed its view of dimensions of my
lme-objects and contrast matrix dimensions, as I get this error message:



> estimable(y_15,Mat15)

Error in estimable.default(y_15, Mat15) :

 Dimension of structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, : 25x26, not compatible with no of parameters in y_15: 25

Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
: 25x26, not compatible with no of parameters in y_15: 25

Dimension of 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
: 25x26, not compatible with no of parameters in y_15: 25

Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
: 25x26, not compatible with no of parameters in y_15: 25

Dimension of 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
: 25x26, not compatible with no of parameters in y_15: 25

Dimension of 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
: 25x26, not compatible with no of parameters in y_15: 25

Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
: 25x26, not compatible with no of parameters in y_15: 25

Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

>



Mat15 is made as a 25x25 contrast matrix, but the updated estimable() sees
it as a 25x26.

Why, and what should I do about it?



I do think the change has occurred in estimable() as I tested one year old
lme() objects and contrast matrix who worked perfectly then, giving the
same error message when run with updated R.
Does estimable() see the row-names as the first column nowdays if I don´t
specify anything else?

All the best
/CG Pettersson

[[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] stopifnot with logical(0)

2015-12-14 Thread Hadley Wickham

On Sat, Dec 12, 2015 at 1:51 PM, Martin Maechler
 wrote:
>> Hadley Wickham 
>> on Sat, 12 Dec 2015 08:08:54 -0600 writes:
>
> > On Sat, Dec 12, 2015 at 3:54 AM, Martin Maechler
> >  wrote:
> >>> Henrik Bengtsson  on
> >>> Fri, 11 Dec 2015 08:20:55 -0800 writes:
> >>
> >> > On Fri, Dec 11, 2015 at 8:10 AM, David Winsemius
> >>  wrote:
> >> >>
> >> >>> On Dec 11, 2015, at 5:38 AM, Dario Beraldi
> >>  wrote:
> >> >>>
> >> >>> Hi All,
> >> >>>
> >> >>> I'd like to understand the reason why
> >> stopifnot(logical(0) == x) doesn't >>> (never?) throw an
> >> exception, at least in these cases:
> >> >>
> >> >> The usual way to test for a length-0 logical object is
> >> to use length():
> >> >>
> >> >> x <- logical(0)
> >> >>
> >> >> stopifnot( !length(x) & mode(x)=="logical" )
> >>
> >> > I found
> >>
> >> > stopifnot(!length(x), mode(x) == "logical")
> >>
> >> > more helpful when troubleshooting, because it will tell
> >> you whether > it's !length(x) or mode(x) == "logical"
> >> that is FALSE.  It's as if you > wrote:
> >>
> >> > stopifnot(!length(x)) > stopifnot(mode(x) == "logical")
> >>
> >> > /Henrik
> >>
> >> Yes, indeed, thank you Henrik --- and Jeff Newmiller
> >> who's nice humorous reply added other relevant points.
> >>
> >> As author stopifnot(), I do agree with Dario's "gut
> >> feeling" that stopifnot() "somehow ought to do the right
> >> thing" in cases such as
> >>
> >> stopifnot(dim(x) == c(3,4))
> >>
> >> which is really subtle version of his cases {But the gut
> >> feeling is wrong, as I argue from now on}.
>
> > Personally, I think the problem there is that people
> > forget that == is vectorised, and for a non-vectorised
> > equality check you really should use identical:
>
> > stopifnot(identical(dim(x), c(3,4)))
>
> You are right "in theory"  but practice is less easy:
> identical() tends to be  too subtle for many users ... even
> yourself (;-), not really of course!),  Hadley, in the above case:
>
> Your stopifnot() would *always* stop, i.e., signal an error
> because typically all dim() methods return integer, and c(3,4)
> is double.
> So, if even Hadley gets it wrong so easily, I wonder if its good
> to advertize to always use  identical() in such cases.
> I indeed would quite often use identical() in such tests, and
> you'd too and would quickly find and fix the "trap" of course..
> So you are mostly right also in my opinion...

Ooops, yes - but you would discover this pretty quickly if you weren't
coding in a email client ;)

I wonder if R is missing an equality operator for this case. Currently:

* == is suboptimal because it's vectorised
* all.equal is suboptimal because it returns TRUE or a text string
* identical is suboptimal because it doesn't do common coercions

Do we need another function (equals()?) that uses the same coercion
rules as == but isn't vectorised? (Like == it would only work with
vectors, so you'd still need identical() for (e.g.) comparing
environments)

Hadley

-- 
http://had.co.nz/

__
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] repeated measure with quantitative independent variable

2015-12-14 Thread Cristiano Alessandro


Hi all,

I am new to R, and I am trying to set up a repeated measure analysis 
with a quantitative (as opposed to factorized/categorical) 
within-subjects variable. For a variety of reasons I am not using 
linear-mixed models, rather I am trying to fit a General Linear Model (I 
am aware of assumptions and limitations) to assess whether the value of 
the within-subjects variable affects statistically significantly the 
response variable. I have two questions. To make myself clear I propose 
the following exemplary dataset (where myfactor_nc is the quantitative 
within-subjects variable; i.e. each subject performs the experiment 
three times -- nc_factor=1,2,3 -- and produces the response in variable dv).


dv <- c(1,3,4,2,2,3,2,5,6,3,4,4,3,5,6);
subject <- 
factor(c("s1","s1","s1","s2","s2","s2","s3","s3","s3","s4","s4","s4","s5","s5","s5"));

myfactor_nc <- c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3)
mydata_nc <- data.frame(dv, subject, myfactor_nc)

*Question 1 (using function aov)*

Easily done...

am1_nc <- aov(dv ~ myfactor_nc + Error(subject/myfactor_nc), data=mydata_nc)
summary(am1_nc)

Unlike the case when myfactor_nc is categorical, this produces three 
error strata: Error: subject, Error: subject:myfactor_nc, Error: Within. 
I cannot understand the meaning of the latter. How is that computed?


*Question 2 (using function lm)*

Now I would like to do the same with the functions lm() and Anova() 
(from the car package). What I have done so far (please correct me if I 
am mistaking) is the following:


# Unstack the dataset
dvm <- with(mydata_nc, cbind(dv[myfactor_nc==1],dv[myfactor_nc==2], 
dv[myfactor_nc==3]))


#Fit the linear model
mlm1 <- lm(dvm ~ 1)

(is that model above correct for my design?)

Now I should use the Anova function, but it seems that it only accepts 
factors, and not quantitative within-subject variable.


Any help is highly appreciated!

Thanks
Cristiano

__
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] Extract survival probabilities from an Extended Cox model

2015-12-14 Thread Justine Nasejje

Dear All,
   How can I extract survival probabilities from a fitted extended
Cox model in R? Thank you in advance.

[[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] R-Help

2015-12-14 Thread Bello A rasheed via R-help

Dear R users 

Kindly help me on how run this code because there is error when ever i run it.
It give me an error at Error: lent >= 1 is not TRUE 
Error: object 'hubPsi' not found
Error in cbind(b = b, A = A, B = B, V = A/B^2, e = B^2/A, gamma. = b/B,  :   
dims [product 3] do not match the length of object [21] Thank you in 
advance
R-CODE
Hf0 <- function(x, c=1.35) pmax2(-c,pmin2(x,c))

Hf <- new("functionX",Hf0)

stopifnot(validObject(Hf)) # ok !


 

 
### psiFunc() examples## classical {trivial, notinteresting}:

F1 <- function(x) rep.int(1,length(x))

cPsi <- psiFunc(rho = function(x)x^2 / 2, psi = function(x) x,

    wgt = F1, Dpsi = F1,

    Erho = function(x) rep.int(1/2,length(x)),

    Epsi2 = F1, EDpsi = F1)


 
### MASS  -- ?rlm --- has

##

##- psi.huber   (u, k = 1.345, deriv =0)

##- psi.hampel  (u, a = 2, b = 4, c =8, deriv = 0)

##- psi.bisquare(u, c = 4.685, deriv = 0)

##  where deriv = 0 :  psi(x)/xi.e.  'wgt'


 
## MM has more in psi-funs.R (seeabove)

## Reproduce Table 1, p.138 of Hampel,Rousseeuw, Ronchetti, Stahel (1986):
b <- c(seq(0, 3, by = 0.1), 4, 5,Inf)

A <- hubPsi@Epsi2(b)

B <- hubPsi@EDpsi(b)


 
options(width = 100, digits = 7)


 
cbind(b = b, A = A, B = B, V = A/B^2, e= B^2/A,

 gamma. = b/B, k. = 1 + b^2/A, lambda. = 1/B)


 
(hubPsi2 <- chgDefaults(hubPsi, k =2)) # works too!

 Bello Abdulkadir RasheedPh.D candidate ( mathematics/statistics)
Department of mathematical Science/ faculty of Science
University Technology Malaysia
81310 UTM, Johor Bahru, Johor, Malaysia
IC. NO. 201204M10026
Matric No. PS113124
HP. No.  0060149144837E-mail: arasheedbe...@yahoo.com
supervisor- Assoc. Prof. Robiah Binti Adnan
[[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] Problems using estimable() (gmodels) on lme()-objects

2015-12-14 Thread Bert Gunter

If you do not receive a satisfactory reply here, this is the sort of
question that you should email the maintainer for, found by:

> maintainer("gmodels")
[1] "Gregory R. Warnes "

Cheers,
Bert
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, Dec 14, 2015 at 2:23 AM, CG Pettersson  wrote:
> R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree", W7-32
>
>
>
> Dear all,
>
> I just upgraded R and downloaded a new version of “gmodels”, to use
> estimable() for extracting fixed effects from lme()-objects.
>
> Now it seems like estimable() have changed its view of dimensions of my
> lme-objects and contrast matrix dimensions, as I get this error message:
>
>
>
>> estimable(y_15,Mat15)
>
> Error in estimable.default(y_15, Mat15) :
>
>  Dimension of structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
> 0, : 25x26, not compatible with no of parameters in y_15: 25
>
> Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> : 25x26, not compatible with no of parameters in y_15: 25
>
> Dimension of 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
> : 25x26, not compatible with no of parameters in y_15: 25
>
> Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> : 25x26, not compatible with no of parameters in y_15: 25
>
> Dimension of 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> : 25x26, not compatible with no of parameters in y_15: 25
>
> Dimension of 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> : 25x26, not compatible with no of parameters in y_15: 25
>
> Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> : 25x26, not compatible with no of parameters in y_15: 25
>
> Dimension of 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>
>>
>
>
>
> Mat15 is made as a 25x25 contrast matrix, but the updated estimable() sees
> it as a 25x26.
>
> Why, and what should I do about it?
>
>
>
> I do think the change has occurred in estimable() as I tested one year old
> lme() objects and contrast matrix who worked perfectly then, giving the
> same error message when run with updated R.
> Does estimable() see the row-names as the first column nowdays if I don´t
> specify anything else?
>
> All the best
> /CG Pettersson
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Draw a dendrogram with ROCK clustering (cba package)

2015-12-14 Thread Ahreum Lee

 Dear all. 
 
Hi, 
Currently, i do clustering analysis with several algorithms in R.  
one of them is ROCK clustering.  
 
Thanks to "cba" package in R, i could easily analyze with ROCK.
(https://cran.r-project.org/web/packages/cba/cba.pdf)
 
but what I want to see is a dendrogram which could be intuitive to define the 
number of clusters. 
 
I tried to find possible functions which could draw a dendrogram within cba 
packages.
but i couldn't...
 
Is there any way to draw a dendrogram with ROCK clustering outcome ?  
 
Any comments would be helpful for me to keep on this research.


Thanks !




All the best, 


Ahreum Lee. 

[[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] [R-pkgs] randomForestSRC 2.0.0

2015-12-14 Thread Udaya B. Kogalur

randomForestSRC 2.0.0 is now available on CRAN.   The package is able
to handle multivariate regression (data sets in which the response is
one or more continuous variables), multivariate classification (data
sets in which the response is one or more factor variables),
multivariate mixed (data sets in which the response is a combination
of continuous and factor variables), survival, competing risk, and
unsupervised forests (data sets in which no response is chosen).

The memory and CPU profiles for the package have improved.  The
package is OpenMP compliant on all popular platforms.  Instructions
for enabling this functionality is available in Rd files, and
pre-compiled OpenMP binaries are also available at:

http://www.ccs.miami.edu/~hishwaran/rfsrc.html.

Linear gains in computation times have been observed on multi-core
desktops and cluster environments using OpenMP.

Finally, the ability to implement user defined split rules is now
available. Instructions and examples for all families are provided in
the package's src directory in splitCustom.c.

Udaya B. Kogalur, Ph.D.
Contract Staff, Dept. of Quantitative Health Sciences, Cleveland
Clinic Foundation

ubkoga...@gmail.com

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

__
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] repeated measure with quantitative independent variable

2015-12-14 Thread Fox, John

Dear Cristiano,

If I understand correctly what you want to do, you should be able to use 
Anova() in the car package (your second question) by treating your numeric 
repeated-measures predictor as a factor and defining a single linear contrast 
for it.

Continuing with your toy example:

> myfactor_nc <- factor(1:3)
> contrasts(myfactor_nc) <- matrix(-1:1, ncol=1)
> idata <- data.frame(myfactor_nc)
> Anova(mlm1, idata=idata, idesign=~myfactor_nc)
Note: model has only an intercept; equivalent type-III tests substituted.

Type III Repeated Measures MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df   Pr(>F)   
(Intercept)  1   0.93790   60.409  1  4 0.001477 **
myfactor_nc  1   0.834787.579  2  3 0.067156 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

With just 3 distinct levels, however, you could just make myfactor_nc an 
ordered factor, not defining the contrasts explicitly, and then you'd get both 
linear and quadratic contrasts.

I hope this helps,
 John

---
John Fox, Professor
McMaster University
Hamilton, Ontario, Canada
http://socserv.socsci.mcmaster.ca/jfox/



> -Original Message-
> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of
> Cristiano Alessandro
> Sent: Monday, December 14, 2015 8:43 AM
> To: r-help@r-project.org
> Subject: [R] repeated measure with quantitative independent variable
> 
> Hi all,
> 
> I am new to R, and I am trying to set up a repeated measure analysis
> with a quantitative (as opposed to factorized/categorical)
> within-subjects variable. For a variety of reasons I am not using
> linear-mixed models, rather I am trying to fit a General Linear Model (I
> am aware of assumptions and limitations) to assess whether the value of
> the within-subjects variable affects statistically significantly the
> response variable. I have two questions. To make myself clear I propose
> the following exemplary dataset (where myfactor_nc is the quantitative
> within-subjects variable; i.e. each subject performs the experiment
> three times -- nc_factor=1,2,3 -- and produces the response in variable
> dv).
> 
> dv <- c(1,3,4,2,2,3,2,5,6,3,4,4,3,5,6);
> subject <-
> factor(c("s1","s1","s1","s2","s2","s2","s3","s3","s3","s4","s4","s4","s5
> ","s5","s5"));
> myfactor_nc <- c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3)
> mydata_nc <- data.frame(dv, subject, myfactor_nc)
> 
> *Question 1 (using function aov)*
> 
> Easily done...
> 
> am1_nc <- aov(dv ~ myfactor_nc + Error(subject/myfactor_nc),
> data=mydata_nc)
> summary(am1_nc)
> 
> Unlike the case when myfactor_nc is categorical, this produces three
> error strata: Error: subject, Error: subject:myfactor_nc, Error: Within.
> I cannot understand the meaning of the latter. How is that computed?
> 
> *Question 2 (using function lm)*
> 
> Now I would like to do the same with the functions lm() and Anova()
> (from the car package). What I have done so far (please correct me if I
> am mistaking) is the following:
> 
> # Unstack the dataset
> dvm <- with(mydata_nc, cbind(dv[myfactor_nc==1],dv[myfactor_nc==2],
> dv[myfactor_nc==3]))
> 
> #Fit the linear model
> mlm1 <- lm(dvm ~ 1)
> 
> (is that model above correct for my design?)
> 
> Now I should use the Anova function, but it seems that it only accepts
> factors, and not quantitative within-subject variable.
> 
> Any help is highly appreciated!
> 
> Thanks
> Cristiano
> 
> __
> 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] stopifnot with logical(0)

2015-12-14 Thread Duncan Murdoch


On 14/12/2015 11:10 AM, Hadley Wickham wrote:

On Sat, Dec 12, 2015 at 1:51 PM, Martin Maechler
 wrote:
>> Hadley Wickham 
>> on Sat, 12 Dec 2015 08:08:54 -0600 writes:
>
> > On Sat, Dec 12, 2015 at 3:54 AM, Martin Maechler
> >  wrote:
> >>> Henrik Bengtsson  on
> >>> Fri, 11 Dec 2015 08:20:55 -0800 writes:
> >>
> >> > On Fri, Dec 11, 2015 at 8:10 AM, David Winsemius
> >>  wrote:
> >> >>
> >> >>> On Dec 11, 2015, at 5:38 AM, Dario Beraldi
> >>  wrote:
> >> >>>
> >> >>> Hi All,
> >> >>>
> >> >>> I'd like to understand the reason why
> >> stopifnot(logical(0) == x) doesn't >>> (never?) throw an
> >> exception, at least in these cases:
> >> >>
> >> >> The usual way to test for a length-0 logical object is
> >> to use length():
> >> >>
> >> >> x <- logical(0)
> >> >>
> >> >> stopifnot( !length(x) & mode(x)=="logical" )
> >>
> >> > I found
> >>
> >> > stopifnot(!length(x), mode(x) == "logical")
> >>
> >> > more helpful when troubleshooting, because it will tell
> >> you whether > it's !length(x) or mode(x) == "logical"
> >> that is FALSE.  It's as if you > wrote:
> >>
> >> > stopifnot(!length(x)) > stopifnot(mode(x) == "logical")
> >>
> >> > /Henrik
> >>
> >> Yes, indeed, thank you Henrik --- and Jeff Newmiller
> >> who's nice humorous reply added other relevant points.
> >>
> >> As author stopifnot(), I do agree with Dario's "gut
> >> feeling" that stopifnot() "somehow ought to do the right
> >> thing" in cases such as
> >>
> >> stopifnot(dim(x) == c(3,4))
> >>
> >> which is really subtle version of his cases {But the gut
> >> feeling is wrong, as I argue from now on}.
>
> > Personally, I think the problem there is that people
> > forget that == is vectorised, and for a non-vectorised
> > equality check you really should use identical:
>
> > stopifnot(identical(dim(x), c(3,4)))
>
> You are right "in theory"  but practice is less easy:
> identical() tends to be  too subtle for many users ... even
> yourself (;-), not really of course!),  Hadley, in the above case:
>
> Your stopifnot() would *always* stop, i.e., signal an error
> because typically all dim() methods return integer, and c(3,4)
> is double.
> So, if even Hadley gets it wrong so easily, I wonder if its good
> to advertize to always use  identical() in such cases.
> I indeed would quite often use identical() in such tests, and
> you'd too and would quickly find and fix the "trap" of course..
> So you are mostly right also in my opinion...

Ooops, yes - but you would discover this pretty quickly if you weren't
coding in a email client ;)

I wonder if R is missing an equality operator for this case. Currently:

* == is suboptimal because it's vectorised
* all.equal is suboptimal because it returns TRUE or a text string
* identical is suboptimal because it doesn't do common coercions

Do we need another function (equals()?) that uses the same coercion
rules as == but isn't vectorised? (Like == it would only work with
vectors, so you'd still need identical() for (e.g.) comparing
environments)
I don't think so.  We already have all(), so all(x == y) would do what 
you want.


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] stopifnot with logical(0)

2015-12-14 Thread William Dunlap

Hadley wrote
> * all.equal is suboptimal because it returns TRUE or a text string

That feature works ok with stopifnot():
  > stopifnot(all.equal("one", 1))
  Error: all.equal("one", 1) are not all TRUE
and I suppose stopifnot could be enhanced to print the text strings that
all.equal() returns so the user has a better idea of what went wront.

I find that all.equal's comparing of names gets in the way of the intent
of stopifnot
  > stopifnot(all.equal(1:2, c(X=1,Y=2)))
  Error: all.equal(1:2, c(X = 1, Y = 2)) is not TRUE

When we use things like all.equal() or new operators that do
non-recycycling binary operations we may make things simpler
for the programmer, but harder for the user.  I think the user would
rather see a specific message about what is wrong rather than a
note that an obscure function did not return TRUE.
> f <- function(x) {
 stopifnot(length(x)==1, x>0)
 seq_len(x)
  }
> f(1:3)
Error: length(x) == 1 is not TRUE
> f(-4)
Error: x > 0 is not TRUE
> f(4)
[1] 1 2 3 4

Another step in this direction is to change the call to stop() in stopifnot
from what I assume is
stop(message, call.=FALSE)
to
stop(simpleError(message, sys.call(-1)))
so the error message included where error was:
> f(-1)
Error in f(-1) : x > 0 is not TRUE
(This is Bug 16188.)

Bill Dunlap
TIBCO Software
wdunlap tibco.com


On Mon, Dec 14, 2015 at 8:10 AM, Hadley Wickham  wrote:
> On Sat, Dec 12, 2015 at 1:51 PM, Martin Maechler
>  wrote:
>>> Hadley Wickham 
>>> on Sat, 12 Dec 2015 08:08:54 -0600 writes:
>>
>> > On Sat, Dec 12, 2015 at 3:54 AM, Martin Maechler
>> >  wrote:
>> >>> Henrik Bengtsson  on
>> >>> Fri, 11 Dec 2015 08:20:55 -0800 writes:
>> >>
>> >> > On Fri, Dec 11, 2015 at 8:10 AM, David Winsemius
>> >>  wrote:
>> >> >>
>> >> >>> On Dec 11, 2015, at 5:38 AM, Dario Beraldi
>> >>  wrote:
>> >> >>>
>> >> >>> Hi All,
>> >> >>>
>> >> >>> I'd like to understand the reason why
>> >> stopifnot(logical(0) == x) doesn't >>> (never?) throw an
>> >> exception, at least in these cases:
>> >> >>
>> >> >> The usual way to test for a length-0 logical object is
>> >> to use length():
>> >> >>
>> >> >> x <- logical(0)
>> >> >>
>> >> >> stopifnot( !length(x) & mode(x)=="logical" )
>> >>
>> >> > I found
>> >>
>> >> > stopifnot(!length(x), mode(x) == "logical")
>> >>
>> >> > more helpful when troubleshooting, because it will tell
>> >> you whether > it's !length(x) or mode(x) == "logical"
>> >> that is FALSE.  It's as if you > wrote:
>> >>
>> >> > stopifnot(!length(x)) > stopifnot(mode(x) == "logical")
>> >>
>> >> > /Henrik
>> >>
>> >> Yes, indeed, thank you Henrik --- and Jeff Newmiller
>> >> who's nice humorous reply added other relevant points.
>> >>
>> >> As author stopifnot(), I do agree with Dario's "gut
>> >> feeling" that stopifnot() "somehow ought to do the right
>> >> thing" in cases such as
>> >>
>> >> stopifnot(dim(x) == c(3,4))
>> >>
>> >> which is really subtle version of his cases {But the gut
>> >> feeling is wrong, as I argue from now on}.
>>
>> > Personally, I think the problem there is that people
>> > forget that == is vectorised, and for a non-vectorised
>> > equality check you really should use identical:
>>
>> > stopifnot(identical(dim(x), c(3,4)))
>>
>> You are right "in theory"  but practice is less easy:
>> identical() tends to be  too subtle for many users ... even
>> yourself (;-), not really of course!),  Hadley, in the above case:
>>
>> Your stopifnot() would *always* stop, i.e., signal an error
>> because typically all dim() methods return integer, and c(3,4)
>> is double.
>> So, if even Hadley gets it wrong so easily, I wonder if its good
>> to advertize to always use  identical() in such cases.
>> I indeed would quite often use identical() in such tests, and
>> you'd too and would quickly find and fix the "trap" of course..
>> So you are mostly right also in my opinion...
>
> Ooops, yes - but you would discover this pretty quickly if you weren't
> coding in a email client ;)
>
> I wonder if R is missing an equality operator for this case. Currently:
>
> * == is suboptimal because it's vectorised
> * all.equal is suboptimal because it returns TRUE or a text string
> * identical is suboptimal because it doesn't do common coercions
>
> Do we need another function (equals()?) that uses the same coercion
> rules as == but isn't vectorised? (Like == it would only work with
> vectors, so you'd still need identical() for (e.g.) comparing
> environments)
>
> Hadley
>
> --
> http://had.co.nz/
>
> __
>

Re: [R] stopifnot with logical(0)

2015-12-14 Thread Duncan Murdoch

On 14/12/2015 11:45 AM, Hadley Wickham wrote:

>> I wonder if R is missing an equality operator for this case. Currently:
>>
>> * == is suboptimal because it's vectorised
>> * all.equal is suboptimal because it returns TRUE or a text string
>> * identical is suboptimal because it doesn't do common coercions
>>
>> Do we need another function (equals()?) that uses the same coercion
>> rules as == but isn't vectorised? (Like == it would only work with
>> vectors, so you'd still need identical() for (e.g.) comparing
>> environments)
>
> I don't think so.  We already have all(), so all(x == y) would do what you
> want.

But that recycles, which is what we're trying to avoid here.

I think this is too special a case to need a function.  Usually all(x == 
y) is the test you want, because so many R functions will recycle.  It's 
not so when x is the dim of an array.  So I could see a weak argument 
for an equalDim() function, but I think it's better just to use

stopifnot(length(dim) == 2, all(dim == c(3,4)))

I know the all() in the second arg isn't needed, but I think it makes 
the intention clearer.  I think there would be less confusion if 
stopifnot() required its args to be single logical values, so I usually 
try to use it that way.

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] stopifnot with logical(0)

2015-12-14 Thread Hadley Wickham

>> I wonder if R is missing an equality operator for this case. Currently:
>>
>> * == is suboptimal because it's vectorised
>> * all.equal is suboptimal because it returns TRUE or a text string
>> * identical is suboptimal because it doesn't do common coercions
>>
>> Do we need another function (equals()?) that uses the same coercion
>> rules as == but isn't vectorised? (Like == it would only work with
>> vectors, so you'd still need identical() for (e.g.) comparing
>> environments)
>
> I don't think so.  We already have all(), so all(x == y) would do what you
> want.

But that recycles, which is what we're trying to avoid here.

Hadley

-- 
http://had.co.nz/

__
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] Extract survival probabilities from an Extended Cox model

2015-12-14 Thread David Winsemius

> On Dec 14, 2015, at 2:01 AM, Justine Nasejje  wrote:
> 
> Dear All,
>   How can I extract survival probabilities from a fitted extended
> Cox model in R? Thank you in advance.

Perhaps you should define the term "extended" in the context of Cox models?

Ideally you would offer one of a) reference to a specific help page in a 
specific package demonstating the constructions of such a model or b) use R 
code to construct an applicable example. 

>   [[alternative HTML version deleted]]

And when you do construct that example, do not post in HTML.

-- 

David Winsemius
Alameda, CA, USA

__
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] Random selection of a fixed number of values by interval

2015-12-14 Thread David Winsemius


> On Dec 14, 2015, at 12:46 PM, David L Carlson  wrote:
> 
> There are lots of ways to do this. For example,

Another method with mapply:

 mapply(function( n, vals) {sample(vals$id, n)} ,   # no replacement is the 
default for sample
vals= split(data, findInterval(data$value, 0:5) )[1:5] , # drops 
the values at 5 or above
n=c(10,7,5,5,3)   )

$`1`
 [1] 12 64 10 60 70 58 33 50 57 68

$`2`
[1] 43  8 15 26 19  3 20

$`3`
[1] 55  9 62 17 67

$`4`
[1] 61 21 31 24 48

$`5`
[1] 44 13 36

> 
>> groups <- cut(data$value, include.lowest = T, right = FALSE,
> +  breaks = 0:ceiling(max(data$value)))
>> grp <- c("[0,1)", "[1,2)", "[2,3)", "[3,4)", "[4,5)")
>> size <- c(10, 7, 5, 5, 3)
>> set.seed(42)
>> samples <- lapply(1:5, function(x) sample(data$id[groups==grp[x]],
> +  size[x]))
>> names(samples) <- grp
>> samples
> $`[0,1)`
> [1] 69 68 33 63 56 46 65 12 50 58
> 
> $`[1,2)`
> [1] 20 34 43  8 15 52 19
> 
> $`[2,3)`
> [1]  7 22 62 28  2
> 
> $`[3,4)`
> [1] 61 53  5 25 21
> 
> $`[4,5)`
> [1] 59 35 40
> 
>> 
>> groups <- cut(data$value, include.lowest = T, right = FALSE,
> +  breaks = 0:ceiling(max(data$value)))
>> grp <- c("[0,1)", "[1,2)", "[2,3)", "[3,4)", "[4,5)")
>> size <- c(10, 7, 5, 5, 3)
>> set.seed(42)
>> samples <- lapply(1:5, function(x) sample(data$id[groups==grp[x]],
> +  size[x]))
>> names(samples) <- grp
>> samples
> $`[0,1)`
> [1] 69 68 33 63 56 46 65 12 50 58
> 
> $`[1,2)`
> [1] 20 34 43  8 15 52 19
> 
> $`[2,3)`
> [1]  7 22 62 28  2
> 
> $`[3,4)`
> [1] 61 53  5 25 21
> 
> $`[4,5)`
> [1] 59 35 40
> 
> 
> -
> David L Carlson
> Department of Anthropology
> Texas A University
> College Station, TX 77840-4352
> 
> 
> 
> -Original Message-
> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Frank S.
> Sent: Monday, December 14, 2015 2:02 PM
> To: r-help@r-project.org
> Subject: [R] Random selection of a fixed number of values by interval
> 
> Dear R users,
> 
> I'm writing to this list because I must get a random sample (without 
> replacement) from a given vector, but the clue is that I need to extract a 
> fixed number of values by each prespecified 1-unit interval. As an example I 
> try to say, I have a data frame that looks like this (my real dataframe is 
> bigger):
> 
> data <- data.frame(id = 1:70, value=  c(0.68, 2.96, 1.93, 5.63, 3.08, 3.10, 
> 2.99, 1.79, 2.96, 0.85, 11.79, 0.06, 4.31, 0.64, 1.43, 0.88, 2.79, 4.67,
>  1.23, 1.43, 3.05, 2.44, 2.55, 3.82, 3.55, 1.56, 7.25, 2.75, 9.64, 5.14, 
> 3.54, 3.12, 0.17, 1.07, 4.08, 4.47, 5.58, 7.41, 0.85, 4.30, 7.58,
>  0.58, 1.40, 4.74, 5.04, 0.14, 1.14, 3.28, 7.84, 0.07, 3.97, 1.02, 3.47, 
> 0.66, 2.38, 0.06, 0.67, 0.48, 4.48, 0.12, 3.82, 2.27, 0.93, 0.30, 
>  0.73, 0.33, 2.91, 0.81, 0.18, 0.42))
> 
> And I would like to select, in a random manner:
> 
> 10 id's whose value belongs to [0,1) interval
> 7 id's whose value belongs to [1,2)
> 5 id's whose value belongs to [2,3)
> 5 id's whose value belongs to [3,4)
> 3 id's whose value belongs to [4,5)
> 
> # I have the following values by each 1-unit interval:
> table(cut(data$value, include.lowest = T, right = FALSE, breaks = 
> 0:ceiling(max(data$value
> 
> and the size vector:
> size <- c(10, 7, 5, 5, 3) 
> 
> But I'm not able to get it by using sample function. Does anyone have some 
> idea?
> 
> Thank you very much for any suggestions!!
> 
> Frank S.
> 
> 
> 
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

__
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] Make a box-whiskers plot in R with 5 variables, color coded.

2015-12-14 Thread Dmitri Leybman

I have a spreadsheet with five different columns standing for five
different variables:

Variable 1 Variable 2 Variable 3 Variable 4 Variable 5 0 0 7 1 0 0 0 7 0 0 1
1 8 2 0 5 5 8 0 0 1 4 8 1 0 4 5 8 0 0 0 1 7 2 1
I am trying to create five box and whiskers plots on a single graph with a
five x-label ticks named for each one of
the variables along with color coding. The names for the x-label would be
"Meeting"[ pertains to Variable1] "Meeting2"[pertains to Variable 2]
Meeting3[pertains to Variable 3], Meeting4[pertains to Variable4],
Meeting5[pertains to Variable5].

I have tried:

boxplot(data, las = 2, col =
c("red", "blue", "black", "aquamarine1", "darkorange3")
, at = c(1, 2, 3, 4, 5), par(mar = c(12, 5, 4, 2) + 0.1),
names = c("Meeting 1", "'Meeting2", "Meeting3", "Meeting4","Meeting5")


Error: unexpected string constant in:
"boxplot(data, las=2, col= c('red', 'blue', 'red', 'red', red')
boxplot(data, las=2, col= c('"

[[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] Make a box-whiskers plot in R with 5 variables, color coded.

2015-12-14 Thread David Winsemius


> On Dec 14, 2015, at 1:34 PM, Dmitri Leybman  wrote:
> 
> I have a spreadsheet with five different columns standing for five
> different variables:
> 
> Variable 1 Variable 2 Variable 3 Variable 4 Variable 5 0 0 7 1 0 0 0 7 0 0 1
> 1 8 2 0 5 5 8 0 0 1 4 8 1 0 4 5 8 0 0 0 1 7 2 1
> I am trying to create five box and whiskers plots on a single graph with a
> five x-label ticks named for each one of
> the variables along with color coding. The names for the x-label would be
> "Meeting"[ pertains to Variable1] "Meeting2"[pertains to Variable 2]
> Meeting3[pertains to Variable 3], Meeting4[pertains to Variable4],
> Meeting5[pertains to Variable5].
> 
> I have tried:
> 
> boxplot(data, las = 2, col =
> c("red", "blue", "black", "aquamarine1", "darkorange3")
> , at = c(1, 2, 3, 4, 5), par(mar = c(12, 5, 4, 2) + 0.1),
> names = c("Meeting 1", "'Meeting2",
 ^
Extra single quote about here


> "Meeting3", "Meeting4","Meeting5")
> 
> 
> Error: unexpected string constant in:
> "boxplot(data, las=2, col= c('red', 'blue', 'red', 'red', red')
> boxplot(data, las=2, col= c('"
> 
>   [[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.

David Winsemius
Alameda, CA, USA

__
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] Random selection of a fixed number of values by interval

2015-12-14 Thread Frank S.

Dear R users,
 
I'm writing to this list because I must get a random sample (without 
replacement) from a given vector, but the clue is that I need to extract a 
fixed number of values by each prespecified 1-unit interval. As an example I 
try to say, I have a data frame that looks like this (my real dataframe is 
bigger):
 
data <- data.frame(id = 1:70, value=  c(0.68, 2.96, 1.93, 5.63, 3.08, 3.10, 
2.99, 1.79, 2.96, 0.85, 11.79, 0.06, 4.31, 0.64, 1.43, 0.88, 2.79, 4.67,
  1.23, 1.43, 3.05, 2.44, 2.55, 3.82, 3.55, 1.56, 7.25, 2.75, 9.64, 5.14, 
3.54, 3.12, 0.17, 1.07, 4.08, 4.47, 5.58, 7.41, 0.85, 4.30, 7.58,
  0.58, 1.40, 4.74, 5.04, 0.14, 1.14, 3.28, 7.84, 0.07, 3.97, 1.02, 3.47, 
0.66, 2.38, 0.06, 0.67, 0.48, 4.48, 0.12, 3.82, 2.27, 0.93, 0.30, 
  0.73, 0.33, 2.91, 0.81, 0.18, 0.42))
 
And I would like to select, in a random manner:
 
10 id's whose value belongs to [0,1) interval
7 id's whose value belongs to [1,2)
5 id's whose value belongs to [2,3)
5 id's whose value belongs to [3,4)
3 id's whose value belongs to [4,5)
 
# I have the following values by each 1-unit interval:
table(cut(data$value, include.lowest = T, right = FALSE, breaks = 
0:ceiling(max(data$value
 
and the size vector:
size <- c(10, 7, 5, 5, 3) 
 
But I'm not able to get it by using sample function. Does anyone have some idea?
 
Thank you very much for any suggestions!!
 
Frank S.

 
 
 
  
[[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] Print minified XML in tree format

2015-12-14 Thread Giorgio Garziano

I may suggest this quick guide:

http://gastonsanchez.com/work/webdata/getting_web_data_r4_parsing_xml_html.pdf

and the following link:

http://www.r-datacollection.com/


I apologize for not being more specific.


--
GG



[[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] Print minified XML in tree format

2015-12-14 Thread Archit Soni

Hi All,

I am new to wonderful world of R, I am stuck at one problem.

I need to print the minified XML string in tree format, i have tried below
code:

ex  <- 
"ToveJaniReminderDon't
forget me this weekend!"> > XML::xmlParse(ex)
  Tove
  Jani
  Reminder
  Don't forget me this weekend!


and it is working fine for this short XML but it doesnt work for the XML
like:

Gambardella,
MatthewXML Developer's
GuideComputer44.952000-10-01An
in-depth look at creating applications
  with XML.Ralls,
KimMidnight
RainFantasy5.952000-12-16A
former architect battles corporate zombies,
  an evil sorceress, and her own childhood to become queen
  of the world.Corets, EvaMaeve
AscendantFantasy5.952000-11-17After
the collapse of a nanotechnology
  society in England, the young survivors lay the
  foundation for a new society.Corets, EvaOberon's
LegacyFantasy5.952001-03-10In
post-apocalypse England, the mysterious
  agent known only as Oberon helps to create a new life
  for the inhabitants of London. Sequel to Maeve
  Ascendant.Corets,
EvaThe Sundered
GrailFantasy5.952001-09-10The
two daughters of Maeve, half-sisters,
  battle one another for control of England. Sequel to
  Oberon's Legacy.Randall, CynthiaLover
BirdsRomance4.952000-09-02When
Carla meets Paul at an ornithology
  conference, tempers fly as feathers get
ruffled.Thurman,
PaulaSplish
SplashRomance4.952000-11-02A
deep sea diver finds true love twenty
  thousand leagues beneath the sea.Knorr, StefanCreepy
CrawliesHorror4.952000-12-06An
anthology of horror stories about roaches,
  centipedes, scorpions  and other
insects.Kress,
PeterParadox LostScience
Fiction6.952000-11-02After
an inadvertant trip through a Heisenberg
  Uncertainty Device, James Salway discovers the problems
  of being quantum.O'Brien, TimMicrosoft .NET: The
Programming 
BibleComputer36.952000-12-09Microsoft's
.NET initiative is explored in
  detail in this deep programmer's
reference.O'Brien,
TimMSXML3: A Comprehensive
GuideComputer36.952000-12-01The
Microsoft MSXML3 parser is covered in
  detail, with attention to XML DOM interfaces, XSLT processing,
  SAX and more.Galos, MikeVisual Studio 7: A
Comprehensive 
GuideComputer49.952001-04-16Microsoft
Visual Studio 7 is explored in depth,
  looking at how Visual Basic, Visual C++, C#, and ASP+ are
  integrated into a comprehensive development

environment.


Please advise.


-- 
Regards
Archit

[[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] repeated measure with quantitative independent variable

2015-12-14 Thread Fox, John

Dear Cristiano,

> -Original Message-
> From: Cristiano Alessandro [mailto:cri.alessan...@gmail.com]
> Sent: Monday, December 14, 2015 2:11 PM
> To: Fox, John
> Cc: r-help@r-project.org
> Subject: Re: [R] repeated measure with quantitative independent variable
> 
> Dear John,
> 
> thanks for your reply! The reason why I did not want to factorize the
> within-subjects variable was to avoid increasing the Df of the model
> from 1 (continuous variable) to k-1 (where k is the number of levels of
> the factors). I am now confused, because you have factorized the
> variable (indeed using "factor"), but the Df of myfactor_nc seems to be
> 1. Could you explain that?

I defined only one (linear) contrast for the factor -- but I made a mistake in 
how I did it, see below.

> 
> Comparing the results obtained with the two methods I seem to get
> completely different results:
> 
> * aov()*
> 
> dv <- c(1,3,4,2,2,3,2,5,6,3,4,4,3,5,6);
> subject <-
> factor(c("s1","s1","s1","s2","s2","s2","s3","s3","s3","s4","s4","s4","s5
> ","s5","s5"));
> myfactor_nc <- c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3)
> mydata_nc <- data.frame(dv, subject, myfactor_nc)
> 
> am1_nc <- aov(dv ~ myfactor_nc + Error(subject/myfactor_nc),
> data=mydata_nc)
> summary(am1_nc)
> 
> Error: subject
>Df Sum Sq Mean Sq F value Pr(>F)
> Residuals  4   12.4 3.1
> 
> Error: subject:myfactor_nc
>  Df Sum Sq Mean Sq F value Pr(>F)
> myfactor_nc  1   14.414.4  16 0.0161 *
> Residuals43.6 0.9
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Error: Within
>Df Sum Sq Mean Sq F value Pr(>F)
> Residuals  5  1.333  0.2667
> 
> 
> *Anova()*
> 
> dvm <- with(mydata_nc, cbind(dv[myfactor_nc==1],dv[myfactor_nc==2],
> dv[myfactor_nc==3]))
> 
> mlm1 <- lm(dvm ~ 1)
> myfactor_nc <- factor(1:3)
> contrasts(myfactor_nc) <- matrix(-1:1, ncol=1)
> idata <- data.frame(myfactor_nc)
> Anova(mlm1, idata=idata, idesign=~myfactor_nc)
> Note: model has only an intercept; equivalent type-III tests
> substituted.
> 
> Type III Repeated Measures MANOVA Tests: Pillai test statistic
>  Df test stat approx F num Df den Df   Pr(>F)
> (Intercept)  1   0.93790   60.409  1  4 0.001477 **
> myfactor_nc  1   0.834787.579  2  3 0.067156 .
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Why is that?

There are two differences: (1) The repeated-measures analysis reported by 
default by Anova() is a MANOVA; if you want the univariate analysis comparable 
to aov(), you have to ask for it. (2) As I mentioned, I made a mistake in 
defining the contrasts for the factor; what I did implicitly included the 
quadratic component, which wasn't my intention. 

Here's the correct version, which, as you can see, produces the same result as 
aov():

> contrasts(myfactor_nc, 1) <- matrix(-1:1, ncol=1)
> contrasts(myfactor_nc)
  [,1]
1   -1
20
31

> idata <- data.frame(myfactor_nc)

> summary(Anova(mlm1, idata=idata, idesign=~myfactor_nc), multivariate=FALSE)
Note: model has only an intercept; equivalent type-III tests substituted.

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

SS num Df Error SS den Df  F   Pr(>F)   
(Intercept) 187.27  1 12.4  4 60.409 0.001477 **
myfactor_nc  14.40  1  3.6  4 16.000 0.016130 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> summary(am1_nc) # your result using aov()

Error: subject
  Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4   12.4 3.1   

Error: subject:myfactor_nc
Df Sum Sq Mean Sq F value Pr(>F)  
myfactor_nc  1   14.414.4  16 0.0161 *
Residuals43.6 0.9 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
  Df Sum Sq Mean Sq F value Pr(>F)
Residuals  5  1.333  0.2667 

My apologies for the confusion about the contrasts.

John
  
> 
> Thanks a lot
> Cristiano
> 
> 
> On 12/14/2015 05:25 PM, Fox, John wrote:
> > Dear Cristiano,
> >
> > If I understand correctly what you want to do, you should be able to
> use Anova() in the car package (your second question) by treating your
> numeric repeated-measures predictor as a factor and defining a single
> linear contrast for it.
> >
> > Continuing with your toy example:
> >
> >> myfactor_nc <- factor(1:3)
> >> contrasts(myfactor_nc) <- matrix(-1:1, ncol=1)
> >> idata <- data.frame(myfactor_nc)
> >> Anova(mlm1, idata=idata, idesign=~myfactor_nc)
> > Note: model has only an intercept; equivalent type-III tests
> substituted.
> >
> > Type III Repeated Measures MANOVA Tests: Pillai test statistic
> >  Df test stat approx F num Df den Df   Pr(>F)
> > (Intercept)  1   0.93790   60.409  1  4 0.001477 **
> > myfactor_nc  1   0.834787.579  2  3 0.067156 .
> > ---
> > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >
> > With just 3

Re: [R] Make a box-whiskers plot in R with 5 variables, color coded.

2015-12-14 Thread peter dalgaard


> On 14 Dec 2015, at 22:54 , David Winsemius  wrote:
> 
>> 
>> On Dec 14, 2015, at 1:34 PM, Dmitri Leybman  wrote:
>> 
>> I have a spreadsheet with five different columns standing for five
>> different variables:
>> 
>> Variable 1 Variable 2 Variable 3 Variable 4 Variable 5 0 0 7 1 0 0 0 7 0 0 1
>> 1 8 2 0 5 5 8 0 0 1 4 8 1 0 4 5 8 0 0 0 1 7 2 1
>> I am trying to create five box and whiskers plots on a single graph with a
>> five x-label ticks named for each one of
>> the variables along with color coding. The names for the x-label would be
>> "Meeting"[ pertains to Variable1] "Meeting2"[pertains to Variable 2]
>> Meeting3[pertains to Variable 3], Meeting4[pertains to Variable4],
>> Meeting5[pertains to Variable5].
>> 
>> I have tried:
>> 
>> boxplot(data, las = 2, col =
>> c("red", "blue", "black", "aquamarine1", "darkorange3")
>> , at = c(1, 2, 3, 4, 5), par(mar = c(12, 5, 4, 2) + 0.1),
>> names = c("Meeting 1", "'Meeting2",
> ^
> Extra single quote about here
> 

Hmm, at face value, that should just become part of the label...

> 
>> "Meeting3", "Meeting4","Meeting5")
>> 
>> 
>> Error: unexpected string constant in:
>> "boxplot(data, las=2, col= c('red', 'blue', 'red', 'red', red')
>> boxplot(data, las=2, col= c('"


...but the error message doesn't match the input, which is lacking a right 
parenthesis. So is there a previous incomplete command maybe, or are we just 
being shown random snippets of things that didn't work??

I think we need to see a full transcript


>> 
>>  [[alternative HTML version deleted]]
>> 


...and PLEASE in plain text, not HTML.

>> __
>> 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.
> 
> David Winsemius
> Alameda, CA, USA
> 
> __
> 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.

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

Re: [R] repeated measure with quantitative independent variable

2015-12-14 Thread Cristiano Alessandro


Dear John,

thanks for your reply! The reason why I did not want to factorize the 
within-subjects variable was to avoid increasing the Df of the model 
from 1 (continuous variable) to k-1 (where k is the number of levels of 
the factors). I am now confused, because you have factorized the 
variable (indeed using "factor"), but the Df of myfactor_nc seems to be 
1. Could you explain that?


Comparing the results obtained with the two methods I seem to get 
completely different results:


* aov()*

dv <- c(1,3,4,2,2,3,2,5,6,3,4,4,3,5,6);
subject <- 
factor(c("s1","s1","s1","s2","s2","s2","s3","s3","s3","s4","s4","s4","s5","s5","s5"));

myfactor_nc <- c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3)
mydata_nc <- data.frame(dv, subject, myfactor_nc)

am1_nc <- aov(dv ~ myfactor_nc + Error(subject/myfactor_nc), data=mydata_nc)
summary(am1_nc)

Error: subject
  Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4   12.4 3.1

Error: subject:myfactor_nc
Df Sum Sq Mean Sq F value Pr(>F)
myfactor_nc  1   14.414.4  16 0.0161 *
Residuals43.6 0.9
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
  Df Sum Sq Mean Sq F value Pr(>F)
Residuals  5  1.333  0.2667


*Anova()*

dvm <- with(mydata_nc, cbind(dv[myfactor_nc==1],dv[myfactor_nc==2], 
dv[myfactor_nc==3]))


mlm1 <- lm(dvm ~ 1)
myfactor_nc <- factor(1:3)
contrasts(myfactor_nc) <- matrix(-1:1, ncol=1)
idata <- data.frame(myfactor_nc)
Anova(mlm1, idata=idata, idesign=~myfactor_nc)
Note: model has only an intercept; equivalent type-III tests substituted.

Type III Repeated Measures MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df   Pr(>F)
(Intercept)  1   0.93790   60.409  1  4 0.001477 **
myfactor_nc  1   0.834787.579  2  3 0.067156 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Why is that?

Thanks a lot
Cristiano


On 12/14/2015 05:25 PM, Fox, John wrote:

Dear Cristiano,

If I understand correctly what you want to do, you should be able to use 
Anova() in the car package (your second question) by treating your numeric 
repeated-measures predictor as a factor and defining a single linear contrast 
for it.

Continuing with your toy example:


myfactor_nc <- factor(1:3)
contrasts(myfactor_nc) <- matrix(-1:1, ncol=1)
idata <- data.frame(myfactor_nc)
Anova(mlm1, idata=idata, idesign=~myfactor_nc)

Note: model has only an intercept; equivalent type-III tests substituted.

Type III Repeated Measures MANOVA Tests: Pillai test statistic
 Df test stat approx F num Df den Df   Pr(>F)
(Intercept)  1   0.93790   60.409  1  4 0.001477 **
myfactor_nc  1   0.834787.579  2  3 0.067156 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

With just 3 distinct levels, however, you could just make myfactor_nc an 
ordered factor, not defining the contrasts explicitly, and then you'd get both 
linear and quadratic contrasts.

I hope this helps,
  John

---
John Fox, Professor
McMaster University
Hamilton, Ontario, Canada
http://socserv.socsci.mcmaster.ca/jfox/




-Original Message-
From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of
Cristiano Alessandro
Sent: Monday, December 14, 2015 8:43 AM
To: r-help@r-project.org
Subject: [R] repeated measure with quantitative independent variable

Hi all,

I am new to R, and I am trying to set up a repeated measure analysis
with a quantitative (as opposed to factorized/categorical)
within-subjects variable. For a variety of reasons I am not using
linear-mixed models, rather I am trying to fit a General Linear Model (I
am aware of assumptions and limitations) to assess whether the value of
the within-subjects variable affects statistically significantly the
response variable. I have two questions. To make myself clear I propose
the following exemplary dataset (where myfactor_nc is the quantitative
within-subjects variable; i.e. each subject performs the experiment
three times -- nc_factor=1,2,3 -- and produces the response in variable
dv).

dv <- c(1,3,4,2,2,3,2,5,6,3,4,4,3,5,6);
subject <-
factor(c("s1","s1","s1","s2","s2","s2","s3","s3","s3","s4","s4","s4","s5
","s5","s5"));
myfactor_nc <- c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3)
mydata_nc <- data.frame(dv, subject, myfactor_nc)

*Question 1 (using function aov)*

Easily done...

am1_nc <- aov(dv ~ myfactor_nc + Error(subject/myfactor_nc),
data=mydata_nc)
summary(am1_nc)

Unlike the case when myfactor_nc is categorical, this produces three
error strata: Error: subject, Error: subject:myfactor_nc, Error: Within.
I cannot understand the meaning of the latter. How is that computed?

*Question 2 (using function lm)*

Now I would like to do the same with the functions lm() and Anova()
(from the car package). What I have done so far (please correct me if I
am mistaking) is the following:

# Unstack the dataset
dvm <- with(mydata_nc,

Re: [R] choropleth packages (US)

2015-12-14 Thread Adrian Waddell

Alaska and Hawaii can be found in the 'world' or 'world2' databases of
the 'maps' package. The following is a bit a hack but it works



library(maps)
library(scales)


mergeMaps <- function(...) {
  maps <- list(...)
  if (length(maps) < 2)
stop("need at least two maps")
  map <- maps[[1]]
  for (i in 2:length(maps)) {
map$x <- c(map$x, NA, maps[[i]]$x)
map$y <- c(map$y, NA, maps[[i]]$y)
map$names <- c(map$names, maps[[i]]$names)
  }
  map$range <- c(range(map$x, na.rm = TRUE), range(map$y, na.rm = TRUE))
  map
}

shiftMap <- function(map, xmin=-180) {
  sel <- !is.na(map$x)
  map$x <- (map$x - xmin) %% 360 + xmin
  map$range <- c(range(map$x, na.rm = TRUE), range(map$y, na.rm = TRUE))
  map
}

m <- shiftMap(mergeMaps(map('state', fill=TRUE, plot=FALSE),
map('world', 'USA:Alaska', fill=TRUE, plot=FALSE),
map('world', 'Hawaii', fill=TRUE, plot=FALSE)),
  xmin=0)

s_data <- tolower(rownames(USArrests))
s_map <- tolower(m$names)

mapping <- lapply(s_data, function(state) {
  which(grepl(state, s_map))
})
## check if the mapping is good!

col_pal <- col_numeric("Greens", domain=NULL, na.color = 'lightyellow')

cols <- rep('lightyellow', length(s_data))

Map(function(indices, col) {
  cols[indices] <<- col
}, mapping, col_pal(USArrests$UrbanPop))

map(m, col=cols, fill=TRUE)
# map.axexs()

## or with no borders
map(m, col=cols, fill=TRUE, border=NA)



Greetings,

Adrian



On Fri, Dec 11, 2015 at 1:22 AM, Benjamin Tyner  wrote:
> Very nice Adrian. Is there a straightforward way to add Alaska and Hawaii at
> the lower left? (without resorting to choroplethr package)
>
>
> On 12/10/2015 06:09 AM, Adrian Waddell wrote:
>>
>> Hi,
>>
>> You can also use the 'maps' package for the map data and the 'scales'
>> package for the color mapping.
>>
>> E.g.
>>
>> library(maps)
>> library(scales)
>>
>> m <- map('state', fill=TRUE, plot=FALSE)
>>
>> s_data <- tolower(rownames(USArrests))
>> s_map <- tolower(m$names)
>>
>> mapping <- lapply(s_data, function(state) {
>>which(grepl(state, s_map))
>> })
>> ## check if the mapping is good!
>>
>> col_pal <- col_numeric("Greens", domain=NULL, na.color = 'lightyellow')
>>
>> cols <- rep('lightyellow', length(s_data))
>>
>> Map(function(indices, col) {
>>cols[indices] <<- col
>> }, mapping, col_pal(USArrests$UrbanPop))
>>
>> map(m, col=cols, fill=TRUE)
>>
>>
>> Adrian
>>
>>
>>
>> On Mon, Dec 7, 2015 at 9:34 AM, Erich Neuwirth
>>  wrote:
>>>
>>> ggplot2 also can do this with
>>> fortify
>>> geom_polygon
>>>
>>> Von meinem iPad gesendet
>>>
 Am 06.12.2015 um 21:03 schrieb Benjamin Tyner :

 Hi

 I wish to draw a basic choropleth (US, by state) and am wondering if
 anyone has any recommendations? I've tried the following thus far:

 1. choroplethr: this works, but required installation of 30+
 dependencies. I would prefer something with fewer dependencies.
 2. tmap: this also seems promising, but most of the examples I saw were
 specific to European maps. Can it be adapted for US?
 3. statebins: doesn't draw true choropleths, but I liked that it doesn't
 have many dependencies.

 Regards
 Ben

 __
 R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
>>>
>>> __
>>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>
>

__
R-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] paneling spplot and hist

2015-12-14 Thread Debasish Pai Mazumder

Hi all,

I am new in R. I am trying to panel spplot and hist.

spplot(hspdf, "CDP", col = "white", col.regions = blue2red(20), sp.layout =
list(l2), at = seq(1,10,1), colorkey = list(space = "bottom", labels =
list(labels = paste(seq(1,10,1)), cex = 1.5)), sub = list("CDP", cex = 1.5,
font = 2))

hist(cdp.obsc, col="grey", border="grey", main="CDP", probability=T)
lines(c.breaks, obs.cdp.d, col="blue")
lines(c.breaks, obs.cdp.e, col="red")


I have tried with par(mfrow=c(2,1))and layout(matrix, both don't work.

Any advice?

with regards
-Pai

[[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] Random selection of a fixed number of values by interval

2015-12-14 Thread David L Carlson

There are lots of ways to do this. For example,

> groups <- cut(data$value, include.lowest = T, right = FALSE,
+  breaks = 0:ceiling(max(data$value)))
> grp <- c("[0,1)", "[1,2)", "[2,3)", "[3,4)", "[4,5)")
> size <- c(10, 7, 5, 5, 3)
> set.seed(42)
> samples <- lapply(1:5, function(x) sample(data$id[groups==grp[x]],
+  size[x]))
> names(samples) <- grp
> samples
$`[0,1)`
 [1] 69 68 33 63 56 46 65 12 50 58

$`[1,2)`
[1] 20 34 43  8 15 52 19

$`[2,3)`
[1]  7 22 62 28  2

$`[3,4)`
[1] 61 53  5 25 21

$`[4,5)`
[1] 59 35 40

> 
> groups <- cut(data$value, include.lowest = T, right = FALSE,
+  breaks = 0:ceiling(max(data$value)))
> grp <- c("[0,1)", "[1,2)", "[2,3)", "[3,4)", "[4,5)")
> size <- c(10, 7, 5, 5, 3)
> set.seed(42)
> samples <- lapply(1:5, function(x) sample(data$id[groups==grp[x]],
+  size[x]))
> names(samples) <- grp
> samples
$`[0,1)`
 [1] 69 68 33 63 56 46 65 12 50 58

$`[1,2)`
[1] 20 34 43  8 15 52 19

$`[2,3)`
[1]  7 22 62 28  2

$`[3,4)`
[1] 61 53  5 25 21

$`[4,5)`
[1] 59 35 40


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



-Original Message-
From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Frank S.
Sent: Monday, December 14, 2015 2:02 PM
To: r-help@r-project.org
Subject: [R] Random selection of a fixed number of values by interval

Dear R users,

I'm writing to this list because I must get a random sample (without 
replacement) from a given vector, but the clue is that I need to extract a 
fixed number of values by each prespecified 1-unit interval. As an example I 
try to say, I have a data frame that looks like this (my real dataframe is 
bigger):

data <- data.frame(id = 1:70, value=  c(0.68, 2.96, 1.93, 5.63, 3.08, 3.10, 
2.99, 1.79, 2.96, 0.85, 11.79, 0.06, 4.31, 0.64, 1.43, 0.88, 2.79, 4.67,
  1.23, 1.43, 3.05, 2.44, 2.55, 3.82, 3.55, 1.56, 7.25, 2.75, 9.64, 5.14, 
3.54, 3.12, 0.17, 1.07, 4.08, 4.47, 5.58, 7.41, 0.85, 4.30, 7.58,
  0.58, 1.40, 4.74, 5.04, 0.14, 1.14, 3.28, 7.84, 0.07, 3.97, 1.02, 3.47, 
0.66, 2.38, 0.06, 0.67, 0.48, 4.48, 0.12, 3.82, 2.27, 0.93, 0.30, 
  0.73, 0.33, 2.91, 0.81, 0.18, 0.42))

And I would like to select, in a random manner:

10 id's whose value belongs to [0,1) interval
7 id's whose value belongs to [1,2)
5 id's whose value belongs to [2,3)
5 id's whose value belongs to [3,4)
3 id's whose value belongs to [4,5)

# I have the following values by each 1-unit interval:
table(cut(data$value, include.lowest = T, right = FALSE, breaks = 
0:ceiling(max(data$value

and the size vector:
size <- c(10, 7, 5, 5, 3) 

But I'm not able to get it by using sample function. Does anyone have some idea?

Thank you very much for any suggestions!!

Frank S.




  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Random selection of a fixed number of values by interval

2015-12-14 Thread Bert Gunter

Yes.

May I suggest:

grp <- c("[0,1)", "[1,2)", "[2,3)", "[3,4)", "[4,5)")

can be obtained more simply as
grp <- levels(groups)[1:5]

 and one slight aesthetic change in the indexing:

from:
samples <- lapply(1:5, function(x) sample(data$id[groups==grp[x]], size[x]))

to:
samples <- lapply(1:5, function(x) sample(data[groups==grp[x],"id"],  size[x]))

(rows and columns in a data frame can be simultaneously indexed)

Cheers,
Bert
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, Dec 14, 2015 at 12:46 PM, David L Carlson  wrote:
> There are lots of ways to do this. For example,
>
>> groups <- cut(data$value, include.lowest = T, right = FALSE,
> +  breaks = 0:ceiling(max(data$value)))
>> grp <- c("[0,1)", "[1,2)", "[2,3)", "[3,4)", "[4,5)")
>> size <- c(10, 7, 5, 5, 3)
>> set.seed(42)
>> samples <- lapply(1:5, function(x) sample(data$id[groups==grp[x]],
> +  size[x]))
>> names(samples) <- grp
>> samples
> $`[0,1)`
>  [1] 69 68 33 63 56 46 65 12 50 58
>
> $`[1,2)`
> [1] 20 34 43  8 15 52 19
>
> $`[2,3)`
> [1]  7 22 62 28  2
>
> $`[3,4)`
> [1] 61 53  5 25 21
>
> $`[4,5)`
> [1] 59 35 40
>
>>
>> groups <- cut(data$value, include.lowest = T, right = FALSE,
> +  breaks = 0:ceiling(max(data$value)))
>> grp <- c("[0,1)", "[1,2)", "[2,3)", "[3,4)", "[4,5)")
>> size <- c(10, 7, 5, 5, 3)
>> set.seed(42)
>> samples <- lapply(1:5, function(x) sample(data$id[groups==grp[x]],
> +  size[x]))
>> names(samples) <- grp
>> samples
> $`[0,1)`
>  [1] 69 68 33 63 56 46 65 12 50 58
>
> $`[1,2)`
> [1] 20 34 43  8 15 52 19
>
> $`[2,3)`
> [1]  7 22 62 28  2
>
> $`[3,4)`
> [1] 61 53  5 25 21
>
> $`[4,5)`
> [1] 59 35 40
>
>
> -
> David L Carlson
> Department of Anthropology
> Texas A University
> College Station, TX 77840-4352
>
>
>
> -Original Message-
> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Frank S.
> Sent: Monday, December 14, 2015 2:02 PM
> To: r-help@r-project.org
> Subject: [R] Random selection of a fixed number of values by interval
>
> Dear R users,
>
> I'm writing to this list because I must get a random sample (without 
> replacement) from a given vector, but the clue is that I need to extract a 
> fixed number of values by each prespecified 1-unit interval. As an example I 
> try to say, I have a data frame that looks like this (my real dataframe is 
> bigger):
>
> data <- data.frame(id = 1:70, value=  c(0.68, 2.96, 1.93, 5.63, 3.08, 3.10, 
> 2.99, 1.79, 2.96, 0.85, 11.79, 0.06, 4.31, 0.64, 1.43, 0.88, 2.79, 4.67,
>   1.23, 1.43, 3.05, 2.44, 2.55, 3.82, 3.55, 1.56, 7.25, 2.75, 9.64, 5.14, 
> 3.54, 3.12, 0.17, 1.07, 4.08, 4.47, 5.58, 7.41, 0.85, 4.30, 7.58,
>   0.58, 1.40, 4.74, 5.04, 0.14, 1.14, 3.28, 7.84, 0.07, 3.97, 1.02, 3.47, 
> 0.66, 2.38, 0.06, 0.67, 0.48, 4.48, 0.12, 3.82, 2.27, 0.93, 0.30,
>   0.73, 0.33, 2.91, 0.81, 0.18, 0.42))
>
> And I would like to select, in a random manner:
>
> 10 id's whose value belongs to [0,1) interval
> 7 id's whose value belongs to [1,2)
> 5 id's whose value belongs to [2,3)
> 5 id's whose value belongs to [3,4)
> 3 id's whose value belongs to [4,5)
>
> # I have the following values by each 1-unit interval:
> table(cut(data$value, include.lowest = T, right = FALSE, breaks = 
> 0:ceiling(max(data$value
>
> and the size vector:
> size <- c(10, 7, 5, 5, 3)
>
> But I'm not able to get it by using sample function. Does anyone have some 
> idea?
>
> Thank you very much for any suggestions!!
>
> Frank S.
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-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] Hexbin: Counting Bins That Meet Certain Criteria

Re: [R] stopifnot with logical(0)

Re: [R] matrix which results singular but at the same time positive definite

[R] R: forest plot metafor

[R] [R-pkgs] New Package "radiomics": First and Second Order Matrix Statistics

[R] Problems when using estimable() (gmodels) on lme()-objects

[R] Problems using estimable() (gmodels) on lme()-objects

Re: [R] stopifnot with logical(0)

[R] repeated measure with quantitative independent variable

[R] Extract survival probabilities from an Extended Cox model

[R] R-Help

Re: [R] Problems using estimable() (gmodels) on lme()-objects

[R] Draw a dendrogram with ROCK clustering (cba package)

[R] [R-pkgs] randomForestSRC 2.0.0

Re: [R] repeated measure with quantitative independent variable

Re: [R] stopifnot with logical(0)

Re: [R] stopifnot with logical(0)

Re: [R] stopifnot with logical(0)

Re: [R] stopifnot with logical(0)

Re: [R] Extract survival probabilities from an Extended Cox model

Re: [R] Random selection of a fixed number of values by interval

[R] Make a box-whiskers plot in R with 5 variables, color coded.

Re: [R] Make a box-whiskers plot in R with 5 variables, color coded.

[R] Random selection of a fixed number of values by interval

Re: [R] Print minified XML in tree format

[R] Print minified XML in tree format

Re: [R] repeated measure with quantitative independent variable

Re: [R] Make a box-whiskers plot in R with 5 variables, color coded.

Re: [R] repeated measure with quantitative independent variable

Re: [R] choropleth packages (US)

[R] paneling spplot and hist

Re: [R] Random selection of a fixed number of values by interval

Re: [R] Random selection of a fixed number of values by interval

33 matches

Site Navigation

Mail list logo

Footer information