from:"Helios de Rosario"

Re: [R] [FORGED] Re: Crete stats course

2017-03-20 Thread Helios de Rosario

The fact related to Rolf's vague recollection is here:
https://stat.ethz.ch/pipermail/r-help/2015-March/426799.html 

And the relevant post (Martin Maechler's answer):
https://stat.ethz.ch/pipermail/r-help/2015-March/426834.html 

Should that be a FAQ?

Helios

> On 16/03/17 03:57, Bert Gunter wrote:
>> Perhaps this has been asked and settled before, but while such
courses
>> certainly might be of interest to those who read this list, they
are
>> for profit, and therefore advertising them here does seem somewhat
>> inappropriate.
>>
[...]
>
>>> On 15/03/2017 20:13, Rolf Turner 
answered:
>
> I have a *vague* recollection that it *has* been asked before and
that 
> there was a consensus, or a pronouncement from R core (or a
combination 
> of the two; or something like that) that such announcements were OK
as 
> long as they were reasonably brief and not overly frequent.  Or 
> something like that.

DESCARGATE el ANUARIO del IBV 2015 desde 
http://www.ibv.org/anuario/Anuario2015/

INSTITUTO DE BIOMECÁNICA 
Universitat Politècnica de València - Edificio 9C 
Camino de Vera s/n - 46022 VALENCIA (ESPAÑA) 
Tel. +34 96170 +34 610567200 - Fax +34 96 387 91 69 
www.ibv.org 

Antes de imprimir este e-mail piense bien si es necesario hacerlo. En
cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección de
Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

__
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] Release of phia 0.2-0

2015-02-17 Thread Helios de Rosario

Dear R users,

I want to announce an update of the package phia, version 0.2-0, now
on CRAN:
http://cran.r-project.org/web/packages/phia/

phia (Post-Hoc Interaction Analysis) is aimed at the analysis of the
expected values and other terms of in linear, generalized, and mixed
linear models, on the basis of multiple comparisons of factor contrasts,
and is specially suited for the analysis of interaction effects. The
function testFactors provides a flexible user interface for defining
combinations of factor levels and covariates, to evaluate and test the
model, using the function linearHypothesis from package car.
testInteractions gives a simpler interface to test multiple
comparisons of simple effects, interaction residuals, interaction
contrasts, or user-defined contrasts. interactionMeans may be used to
explore the cell means of factorial designs, and plot main effects or
first-order interactions.

This update incorporates some minor bug fixes to previous versions, and
a new feature that had long been requested by some package users: the
report of the error associated to the adjusted values calculated by
functions like interactionMeans. Now the standard errors are printed
together with the adjusted values, and the plot method for interaction
tables adds error bars with the size of such standard errors, or
asymptotic confidence intervals based on such errors.

In addition, I would like to refer the package users to the development
repository in GitHub: https://github.com/heliosdrm/phia, where new
contributions, issues/bug reports, and other comments are welcome.

Best regards

Helios De Rosario

biomecanicamente.org
Conoce la actualidad on-line del IBV
__

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96170- +34 610567200 • Fax +34 96 387 91 69
www.ibv.org

Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

___
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] Making a trial based Package for 30-days.

2014-05-12 Thread Helios de Rosario

I guess you could write some instructions in the configure file
(configure.win for windows) at the root level of the package
subdirectory, to create a scheduled task to uninstall and delete the
package after a fixed period. That script would normally execute at the
moment of installation. See Package structure in Writing R
extensions.

Yes, easy to hack. But R is free software, so that's what you should
normally expect for any mechanism. In any event, I suppose that you
don't want a trial-based package for your own use alone, but you want
to distribute it. And if you want to distribute it by mainstream
repositories (e.g. CRAN), you should publish the source anyway, and
won't be allowed to include such malicious mechanisms.

Helios De Rosario


 El día 12/05/2014 a las 11:31, Ashis Deb ashisde...@gmail.com
escribió:
 Hi  all   ,
 
 
I have   a  GUI  package , which  I  want  to  make 
it
 work  for  a  certain period,say‑30days   ,after which  it should  be
self
 destructive  .
 
 Could it  be possible  ??
 
 
 ASHIS
 

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] about lm()

2014-03-31 Thread Helios de Rosario

 El día 30/03/2014 a las 15:23, Si Qi L. liusiqi.n...@gmail.com
 escribió:
 Hi
 
 I have a problem with linear regression. This is my codes:
 
 acc1‑ lm(data$acc ~ dummy + data$Reqloanamount + data$Code + 
 data$Code.1 +
 data$EmpNetMonthlyPay + data$Customer_Age + data$RTI)
 summary(acc1)
 
 These attributes are all numerical except the acc(factors), so how
can I
 fix the problem R showed me? can anyone please help me out? Many 
 thanks!
 

If you want to fit a model with a factor dependent variable, you should
not use lm(). Perhaps you are looking for ?glm (family binomial) if it
is binomial data, or ?polr (in package MASS), ?multinom (in nnet), or
?clm and others (in ordinal) if the variable has more than two levels.

Helios De Rosario

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] outliers for Likert scale data

2013-09-04 Thread Helios de Rosario

 El día 01/09/2013 a las 15:13, Helen Sawaya
helensaw...@hotmail.com
escribió:
 Dear R experts,
 I have data from a questionnaire that I would like to factor analyse.
It is 
 in a likert scale form (0-3). I would like to check first for
univariate and 
 multivariate outliers but the most common ways of doing so assume the
data is 
 continuous and normal- neither of which is the case here. I found an
article 
 discussing this (Outlier Detection in Test and Questionnaire Data by
Wobbe P. 
 Zijlstra, L. Andries van der Ark, and Klaas Sijtsma), but I was
wondering if 
 I could get the exact R code on how to implement the outlier
detection 
 analyses.

I have not found an exact implementation of that article, but one of
its authors (van den Ark) has published the mokken package with some
methods referred to in it:
https://sites.google.com/a/tilburguniversity.edu/avdrark/mokken

The ESD method for identifying outliers, also used in the paper to
handle outlier scores, is implemented (together with others) in the
package parody:
http://www.bioconductor.org/packages/release/bioc/html/parody.html

Hope it helps
Helios De Rosario


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] ask help!

2013-07-25 Thread Helios de Rosario



-- 
Helios de Rosario Martínez
 
 Researcher
 El día 25/07/2013 a las 11:44, mei_yuan
mei_y...@dragon.nchu.edu.tw
escribió:
 Hi,
 
 In the R console, I have the following:
 
 runif(10)
 Error in runif(10) : 
   '.Random.seed' is not an integer vector but of type 'list'

It seems you have overwritten the default value of the variable
.Random.seed. Delete it and try again:

rm(.Random.seed)
runif(10)




INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

[R] Change R options only inside a package

2013-07-15 Thread Helios de Rosario

Hi,

What should I do if I want to use non-default R-options for the
functions inside a package, but not affect the options of the rest of
the session? Specifically, I'm thinking of using sum contrasts instead
of the default, in functions like lm, etc. when they are called by
other functions in a package.

I can think on the following solutions:

1) Quick-and-dirty solution: use options() to modify the default
contrasts before calling my functions. If I don't want to affect the
rest of the session, I should reset the options before returning,
preferably with exception handling to avoid that the options are not
reset if my function fails.

2) Complex-but-tidier solution: Use function-specific arguments to tell
the contrasts of the factors that I'm using. E.g., I should create a
named list of contrasts and pass it as the argument contrasts every
time I call lm:

mod - lm(form, contrasts=list(factor1=contr.sum,
factor2=contr.sum))


3) Encapsulated solution: create package-specific (unexported) versions
of lm and other functions, that already implement #2, e.g.:

lm - function(...)
{
# Code that examines dots, to get the factors of the model, and create
the list of contrasts
# Then call the stats lm (be sure that dots did not already include a
contrasts element).
stats::lm(..., contrasts = list.of.factors.with.sum.contrasts)
}

I don't like #1, because changing the options back and forth may be
error-prone and lead to unexpected behaviour.

#2 is better, but prone to programmer errors: I should be careful of
not forgetting to set the contrasts argument every time I call lm.

I prefer #3, because it avoids the disadvantages of the previous
solutions, but perhaps it is overkill, hence my question: Is there a
clean way of setting the R-session options in such a way that they
only affect the functions and data inside a package?

Thanks and best regards

Helios De Rosario

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] Change R options only inside a package

2013-07-15 Thread Helios de Rosario

 El día 15/07/2013 a las 13:30, Duncan Murdoch
murdoch.dun...@gmail.com escribió:
 On 13-07-15 5:27 AM, Helios de Rosario wrote:
 Hi,

 What should I do if I want to use non-default R-options for the
 functions inside a package, but not affect the options of the rest
of
 the session? Specifically, I'm thinking of using sum contrasts
instead
 of the default, in functions like lm, etc. when they are called
by
 other functions in a package.

[...]

 I'd recommend that in your own functions that need this, you put this
at 
 the beginning:
 
   saveOpts - options(contrasts = c(contr.sum, contr.sum))
   on.exit(options(saveOpts))
 
 This means that the options will be changed just for the duration of
the 
 call.  If there are situations where users might not want this new 
 default, make the new value a default of a function argument.

Thank you very much. I was not aware of the on.exit function.

Helios De Rosario

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] suppress startup messages from default packages

2013-07-15 Thread Helios de Rosario

 Hi all,
 
 several packages print messages during loading.  How do I avoid to
see
 them when the packages are in the defaultPackages?
 
 Here is an example.
 
 With this in ~/.Rprofile
 ,[ ~/.Rprofile ]
 | old - getOption(defaultPackages)
 | options(defaultPackages = c(old, filehash))
 | rm(old)
 `
 
 I get as last line when starting R:
 ,
 | filehash: Simple key-value database (2.2-1 2012-03-12)
 `
 
 Another package with (even more) prints during startup is
tikzDevice.
 
How can I avoid to get these messages?


There are several options in ?library to control the messages that are
displayed when loading packages. However, this does not seem be able to
supress all the messages. Some messages are defined by the package
authors, because they feel necessary that the user reads them.

Helios De Rosario

 Regards,
 Andreas


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] GLMM post- hoc comparisons

2013-01-09 Thread Helios de Rosario

 El día 08/01/2013 a las 12:40, Silvina Velez
sve...@mendoza-conicet.gob.ar
escribió:
 Hi All,
 I have data about seed predation (SP) in fruits of three differents
colors 
 (yellow, motted, dark) and in two fruiting seasons (2007, 2008). I
performed 
 a GLMM (lmer function, lme4 package) and the outcome showed that the

 interaction term (color:season) was significant, and some
combinations of 
 this interaction have significant Pr(|z|), but I don't think they
are the 
 right significant combinations, because when I look the bwplot, this

 combinations seems to be very different from the other ones. So, I
would like 
 to know if there is any test a posteriori to know the p-values for
each 
 combination of color:season, and thereby be able to know what
conbination/s 
 is/are really significant.
 
 m1=lmer(SP ~ color + season:color +(1|Site:tree), data=datosfl, 
 family=poisson)
 AIC   BIC logLik deviance
 178.3 196.6 -81.14162.3
 Random effects:
 Groups  NameVariance Std.Dev.
 obsBR   (Intercept) 0.064324 0.25362 
 Site:tree   (Intercept) 0.266490 0.51623 
 Number of obs: 73, groups: obsBR, 73; Site:tree, 37
 
 Estimate Std. Error z value Pr(|z|)
 (Intercept)2.5089 0.2750   9.125   2e-16 ***
 colorM-0.1140 0.3242  -0.352   0.7250
 colorD-0.6450 0.4178  -1.544   0.1227
 Season2008-0.7343 0.3104  -2.365   0.0180 *  
 colorM:Season2008  0.2505 0.4352   0.576   0.5648
 colorD:Season2008  1.1445 0.5747   1.992   0.0464 * 

Hi Silvina,

What do you exactly mean with what combination(s) is/are significant?
If you mean what combinations have significantly greater SP than the
baseline combination (yellow:2007), the table that you have copied may
be what you actually want. If you want to test other contrasts between
color:season combinations, perhaps you can use the function
testInteractions() from package phia. For instance:

testInteractions(m1)

will give you a test of all the pairwise contrasts between color and
season. You can also test simple main effects, or other specific
contrasts by adding further arguments (see the documentation and the
package vignette). Anyway, the calculation of p-values in mixed models
must always be taken with care.

Helios De Rosario-Martinez
Instituto de Biomecánica de Valencia



INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] ANOVA repeated measures and post-hoc

2012-12-30 Thread Helios de Rosario

A (late) update to this question:

On Fri Aug 17 07:33:29, Henrik Singmann wrote:
 Hi Diego,

 I am struggeling with this question also for some time and there does

 not seem to be an easy and general solution to this problem. At least
I 
 haven't found one.
 However, if you have just one repeated-measures factor, use the
solution 
 describe by me here: http://stats.stackexchange.com/a/15532/442 
 Furthermore, you might wanna check the package phia with accompanying

 vignette (but it uses multivariate tests).

There is a new version of phia that also works with mixed-effects
models - but notice the caveat about the calculation of p-values for
such models, as explained in the vignette.

For more issues about mixed-effects models in R, there is a specifice
SIG-list:
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

 For running the ANOVA you could also check my new package afex.

 Cheers,
 Henrik
 
 
 Diego Bucci schrieb:
 Hi,
  I performed an ANOVA repeated measures but I still can't find any
good news regarding the possibility to perform multiple comparisons.
  Can anyone help me?
  Thanks


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

[R] Two-way linear model with interaction but without one main effect

2012-06-12 Thread Helios de Rosario

Hi,

I know that the type of model described in the subject line violates
the principle of marginality and it is rare in practice, but there may
be some circumstances where it has sense. Let's take this imaginary
example (not homework, just a silly made-up case for illustrating the
rare situation):

I'm measuring the energy absorption of sports footwear in jumping. I
have three models (S1, S2, S3), that are known by their having the same
average value of this variable for different types of ground, but I want
to model the energy absorption for specific ground types (grass, sand,
and pavement).

To fit the model I take 90 independent measures (different shoes,
different users for each observation), with 10 samples per footwear
model and ground type.

# Example data:
shoe - gl(3,30,labels=c(S1,S2,S3))
ground - rep(gl(3,10,labels=c(grass,sand,pavement)),3)
Y - rnorm(90,120,20)

My model may include a main effect of the ground type, and the
interaction shoe:ground, but I think that in this peculiar case I could
neglect the main effect of shoe, since my initial hypothesis is that the
average energy absorption is the same for the three models.

My first thought was fitting the following model (with effect coding,
so that the interaction coeffs have zero mean.):

mod1 - lm(Y ~ ground + ground:shoe,
contrasts=list(shoe=contr.sum,ground=contr.sum))

But this model has the same number of coefficients as a full factorial,
and actually represents the same model subspace, isn't it? In fact, the
marginal means are not the same for the three types of shoes:

# Marginal means for my (random) example data
 tapply(predict(mod1),shoe,FUN=mean)
  S1   S2   S3 
116.3581 121.0858 118.3800

If I'm not mistaken, to create the model that I want I can start with
the full factorial model and remove the part associated to the main shoe
effect:

# Full model and its model matrix
mod1 - lm(Y~shoe*ground,
contrasts=list(shoe=contr.sum,ground=contr.sum))
X - model.matrix(mod1)
# Split X columns by terms
X1 - X[,1]
X.shoe - X[,2:3]
X.ground - X[,4:5]
X.interact - X[,6:9]
# New model without method main effect
mod2 - lm(Y~X.ground+X.interact)

For this model the marginal means do coincide:
 tapply(predict(mod2),shoe,FUN=mean)
 S1  S2  S3 
118.608 118.608 118.608 

My questions are:
Is this correct? And is there an easier way of doing this?

Thanks
Helios De Rosario

-- 
Helios de Rosario Martínez
 
 Researcher


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] Two-way linear model with interaction but without one main effect

2012-06-12 Thread Helios de Rosario

Thanks for the suggestion, Thierry.

Nevertheless, in this example I'm not considering shoe as a random,
nuisance factor with zero mean. I'm considering three specific shoe
models, and I'm interested in modelling how the output changes between
the different shoes for those grounds, given that the average output is
the same for all shoes. That's not the type of question addressed by a
mixed model, I'm afraid.

Helios

 El día 12/06/2012 a las 14:17, ONKELINX, Thierry
thierry.onkel...@inbo.be
escribió:
 Dear Helios,
 
 I think you rather want a mixed model with shoe as random effect.
 
 library(lme4)
 lmer(Y ~ Ground + (1|Shoe)) #the effect of shoe is independent of the
ground 
 effect
 or
 lmer(Y ~ Ground + (0 + Ground|Shoe)) #the effect of shoe is different
per 
 ground.
 
 Best regards,
 
 Thierry
 
 ir. Thierry Onkelinx
 Instituut voor natuur- en bosonderzoek / Research Institute for
Nature and 
 Forest
 team Biometrie  Kwaliteitszorg / team Biometrics  Quality
Assurance
 Kliniekstraat 25
 1070 Anderlecht
 Belgium
 + 32 2 525 02 51
 + 32 54 43 61 85
 thierry.onkel...@inbo.be 
 www.inbo.be 
 
 To call in the statistician after the experiment is done may be no
more than 
 asking him to perform a post-mortem examination: he may be able to
say what 
 the experiment died of.
 ~ Sir Ronald Aylmer Fisher
 
 The plural of anecdote is not data.
 ~ Roger Brinner
 
 The combination of some data and an aching desire for an answer does
not 
 ensure that a reasonable answer can be extracted from a given body of
data.
 ~ John Tukey
 
 
 -Oorspronkelijk bericht-
 Van: r-help-boun...@r-project.org
[mailto:r-help-boun...@r-project.org] 
 Namens Helios de Rosario
 Verzonden: dinsdag 12 juni 2012 13:35
 Aan: r-help@r-project.org 
 Onderwerp: [R] Two-way linear model with interaction but without one
main 
 effect
 
 Hi,
 
 I know that the type of model described in the subject line violates
the 
 principle of marginality and it is rare in practice, but there may be
some 
 circumstances where it has sense. Let's take this imaginary example
(not 
 homework, just a silly made-up case for illustrating the rare
situation):
 
 I'm measuring the energy absorption of sports footwear in jumping. I
have 
 three models (S1, S2, S3), that are known by their having the same
average 
 value of this variable for different types of ground, but I want to
model the 
 energy absorption for specific ground types (grass, sand, and
pavement).
 
 To fit the model I take 90 independent measures (different shoes,
different 
 users for each observation), with 10 samples per footwear model and
ground 
 type.
 
 # Example data:
 shoe - gl(3,30,labels=c(S1,S2,S3)) ground - 
 rep(gl(3,10,labels=c(grass,sand,pavement)),3)
 Y - rnorm(90,120,20)
 
 My model may include a main effect of the ground type, and the
interaction 
 shoe:ground, but I think that in this peculiar case I could neglect
the main 
 effect of shoe, since my initial hypothesis is that the average
energy 
 absorption is the same for the three models.
 
 My first thought was fitting the following model (with effect coding,
so 
 that the interaction coeffs have zero mean.):
 
 mod1 - lm(Y ~ ground + ground:shoe,
 contrasts=list(shoe=contr.sum,ground=contr.sum))
 
 But this model has the same number of coefficients as a full
factorial, and 
 actually represents the same model subspace, isn't it? In fact, the
marginal 
 means are not the same for the three types of shoes:
 
 # Marginal means for my (random) example data
 tapply(predict(mod1),shoe,FUN=mean)
   S1   S2   S3
 116.3581 121.0858 118.3800
 
 If I'm not mistaken, to create the model that I want I can start with
the 
 full factorial model and remove the part associated to the main shoe
 effect:
 
 # Full model and its model matrix
 mod1 - lm(Y~shoe*ground,
 contrasts=list(shoe=contr.sum,ground=contr.sum))
 X - model.matrix(mod1)
 # Split X columns by terms
 X1 - X[,1]
 X.shoe - X[,2:3]
 X.ground - X[,4:5]
 X.interact - X[,6:9]
 # New model without method main effect
 mod2 - lm(Y~X.ground+X.interact)
 
 For this model the marginal means do coincide:
 tapply(predict(mod2),shoe,FUN=mean)
  S1  S2  S3
 118.608 118.608 118.608
 
 My questions are:
 Is this correct? And is there an easier way of doing this?
 
 Thanks
 Helios De Rosario
 
 --
 Helios de Rosario Martínez
 
  Researcher
 
 
 INSTITUTO DE BIOMECÁNICA DE VALENCIA
 Universidad Politécnica de Valencia ● Edificio 9C Camino de Vera
s/n ● 46022 
 VALENCIA (ESPAÑA) Tel. +34 96 387 91 60 ● Fax +34 96 387 91 69
www.ibv.org 
 
   Antes de imprimir este e-mail piense bien si es necesario hacerlo.
 En cumplimiento de la Ley Orgánica 15/1999 reguladora de la
Protección de 
 Datos de Carácter Personal, le informamos de que el presente mensaje
contiene 
 información confidencial, siendo para uso exclusivo del destinatario
arriba 
 indicado. En caso de no ser usted el destinatario del mismo le
informamos que 
 su recepción no le autoriza a su divulgación

Re: [R] replacing NA for zero

2012-06-12 Thread Helios de Rosario

 Dear R users,

   I have a very basic query, but was unable to find a proper anwser.

   I have the following data.frame

   x y
   2   0.12
   3   0.25
   4   0.11
   6   0.16
   7   0.20

   and, due to further calculations, I need the data to be stored as

   x y
   1 0
   2   0.12
   3   0.25
   4   0.11
   5 0
   6   0.16
   7   0.20
   80

   How do I do the transformation?

Let's suppose your original data frame is DF1. Then:

 DF2 - data.frame(x=1:8,y=rep(0,8))
 DF2$y[DF1$x] - DF1$y

But I'd recommend you use NA for the missing values, unless you have a
good reason to code them as zero.

   Many many thjanks in advance.

You are welcome.
Helios De Rosario

-- 
Helios de Rosario Martínez
 
 Researcher


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] anovas ss typeI vs typeIII

2012-05-21 Thread Helios de Rosario

Hi, you should give more details of your problem (at least some output,
as Peter Daalgard says). But you are probably asking for something like
this:
http://www.r-bloggers.com/anova-%E2%80%93-type-ii-ss-explained/
or many other webpages that you may find if you Google or R-seek with
keywords like anova type i, type iii, etc.

If you have within-subjects factors, the same commands will not give
you a sensible result. If you want to use Anova() for that, you should
use the arguments idata, idesign, etc., as explained in the help
page of that function, or in the web appendix to the CAR book:
http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Multivariate-Linear-Models.pdf

Helios De Rosario

--
Helios de Rosario Martínez

Researcher

El día 20/05/2012 a las 0:58, jacaranda tree
myjacara...@yahoo.com escribió:
Hi all,
I have been struggling with ANOVAs on R. I am new to R, so I created
a
simple data frame, and I do some analyses on R just to learn R and
then check
them on SPSS to make sure that I am doing fine. Here is the problem
that I've
run into:

when we use the aov function, it uses SS Type I as default (on SPSS
it is
Type III). Then I used the Anova function under cars package using
the
command:

mod - lm(DV
~ IV1*IV2, data = mydata,
contrasts=list(IV1=contr.sum,
IV2=contr.sum))
Anova(mod, type=”3”)

Above, both of my IVs are between-SS variables. But still, results
from this
model do not match the results from SPSS (I have to say they are not
too
different either). But I was wondering if I am doing something wrong.
If what
I am doing is okay, then my next question is can I use the same set
of
commands (for Anova function) if one of my IVs was within-SS and the
other IV
was between-SS?

Thank you very much!

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] MANOVA with random factor

2012-05-18 Thread Helios de Rosario

I suppose you mean this:
http://www.r-project.org/doc/Rnews/Rnews_2007-2.pdf

Another approach, less general but in my opinion easier to use (and
equivalent in many situations), is provided by Anova() in package car.
See:
http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Multivariate-Linear-Models.pdf

The package ez also has ezANOVA(), a wrapper to Anova() that may
simplify the task for full factorial designs.

Helios De Rosario



 El día 17/05/2012 a las 12:05, David Costantini
david.costant...@glasgow.ac.uk escribió:
 Dear All
 I would need to perform a MANOVA with both fixed (group, sex,
group*sex) and
 random (brood) effects. I wonder if this is at all possible and if R
does 
 that. 
 At the moment, I only know that I can run a classic MANOVA with R.
 
 Thank you
 David
 
 __ 
 
 David Costantini, PhD
 
 http://www.davidcostantini.it 
 NERC Postdoctoral research associate
 
 Institute of Biodiversity, Animal Health and Comparative Medicine
 School of Life Sciences
 College of Medical, Veterinary and Life Sciences
 University of Glasgow
 Graham Kerr Building, room 511
 Glasgow G12 8QQ, UK
 
 See also my association Ornis italica
 http://www.ornisitalica.com 
 http://www.birdcam.it 
 
 
 
 
 
   [[alternative HTML version deleted]]

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] MANOVA with random factor

2012-05-18 Thread Helios de Rosario

Hi, after re-reading I think that I misunderstood your question. You
don't provide many details, but I suppose that the brood effect is
nested within the fixed effects, so you don't mean a multivariate
approach for a split-plot or a repeated-measures design, but the
analysis of a multivariate mixed effects model. Do you?

In that case, perhaps the package MCMCglmm may help, although I have
never used it, so I can't tell anything for sure.

Helios

 El día 17/05/2012 a las 12:05, David Costantini
david.costant...@glasgow.ac.uk escribió:
 Dear All
 I would need to perform a MANOVA with both fixed (group, sex,
group*sex) and
 random (brood) effects. I wonder if this is at all possible and if R
does 
 that. 
 At the moment, I only know that I can run a classic MANOVA with R.
 
 Thank you
 David
 
 __ 
 
 David Costantini, PhD
 
 http://www.davidcostantini.it 
 NERC Postdoctoral research associate
 
 Institute of Biodiversity, Animal Health and Comparative Medicine
 School of Life Sciences
 College of Medical, Veterinary and Life Sciences
 University of Glasgow
 Graham Kerr Building, room 511
 Glasgow G12 8QQ, UK
 
 See also my association Ornis italica
 http://www.ornisitalica.com 
 http://www.birdcam.it 
 
 
 
 
 
   [[alternative HTML version deleted]]

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

[R] read.octave fails with data from Octave 3.2.X

2012-03-27 Thread Helios de Rosario

Hi,

I'm afraid that the function read.octave from package foreign has
some problems with the ASCII data format exported by new versions of
Octave (later than 3.2.X). It fails even for a simple case as:

[Octave code:]
octave:1 x=1;
octave:2 save -ascii testdata.mat x

[Now in R:]
 octavedata - read.octave('testdata.mat')
Mensajes de aviso perdidos
In read_octave_unknown(con, type) : cannot handle unknown type ''

In this simple case I guess that the problem is that new versions
Octave append two blank lines after each variable, and this confuses the
current implementation of read.octave()

The problem is worse if the saved variables include other types as
structs, or strings. The new syntax of the MAT files is not recognized
by read.octave().

Of course, it's always difficult to keep this kind of functions working
when the external program changes its specification for saving
variables, but if would be nice if the maintainers of foreign could at
least solve the issue of blank lines. That way, it would still be
possible to import simple data types as scalars and matrices.

Otherwise, I suppose that a workaround is saving the data in binary
(matlab) format, then load it with Octave 3.2.X, and save it in text
format from that version.

 sessionInfo()
R version 2.14.2 (2012-02-29)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252
[3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
[5] LC_TIME=Spanish_Spain.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] foreign_0.8-49




-- 
Helios de Rosario Martínez
 
 Researcher


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] read.octave fails with data from Octave 3.2.X

2012-03-27 Thread Helios de Rosario

I wrote in my previous message the following Octave code:

[Octave code:]
octave:1 x=1;
octave:2 save -ascii testdata.mat x

Forget the -ascii. It should be -text or nothing (-text is the
default).

By the way, read.octave() does not really fail (it does return a
value), but the result is somewhat corrupted: it contains the exported
x variable, plus other empty elements corresponding to the blank
lines, I think.

Helios
INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] Type II and III sum of squares (R and SPSS)

2012-03-22 Thread Helios de Rosario

 El día 21/03/2012 a las 11:27, Marco Tommasi
marco.tomm...@unich.it wrote
 To whom it may concern
 
 I made some analysis with R using the command Anova. However, I found

 some problmes with the output obtained by selecting type II o type
III 
 sum of squares.
 
 Briefly, I have to do a 2x3 mixed model anova, wherein the first
factor 
 is a between factor and the second factor is a within factor. I use
the 
 command Anova in the list below, because I want to obtain also the
sum 
 of squares of the linear and quadratic contrast between the levels of

 the within factor.
 
 Below I report the list of commands used in R (fattA is the
beteween 
 factor and fB is the within factor):
 
   a1b1-c(10,9,8,7)
   a1b2-c(7,6,4,5)
   a1b3-c(3,2,3,4)
   a2b1-c(9,9,8,7)
   a2b2-c(8,7,9,7)
   a2b3-c(7,8,8,6)
  
   M3-matrix(0,8,4)
   M3[,1]-cbind(a1b1,a2b1)
   M3[,2]-cbind(a1b2,a2b2)
   M3[,3]-cbind(a1b3,a2b3)
   M3[,4]-rep(c(1,2),each=4)
  
   colnames(M3)-c(b1,b2,b3,fattA)
  
   M3-as.data.frame(M3)
  
   require(car)
 Loading required package: car
 Loading required package: MASS
 Loading required package: nnet
   f5-lm(cbind(b1,b2,b3)~fattA,data=M3)

You are missing a couple of important points here. First, fattA is a
factor according to your description, but you have coded it as a numeric
variable. You must coerce it to factor:

M3$fattA - as.factor(M3$fattA)

Now, type III SS with Anova() only provides sensible results, if the
factor contrasts are orthogonal, e.g. contr.sum, contr.poly, or
contr.helmert. This is not the default in R, so you should make it
explicit:

f5-lm(cbind(b1,b2,b3)~fattA,data=M3,
 contrasts=list(fattA=contr.sum))

Now, Anova() gives the same results regardless of the SS type:

   fB-factor(c(1:3))
   d2-data.frame(fB)

Try:
summary(Anova(f5,idata=d2,idesign=~fB,type=2),multivariate=FALSE)
summary(Anova(f5,idata=d2,idesign=~fB,type=3),multivariate=FALSE)

I recommend you read Fox's and Weisberg's book that explain many
functions of car, including ANOVA. See specially pages 193 and
following. Type
 carWeb()
to go to its webpage.

Helios

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] Ggplot barchart drops factor levels: how to show them with zero counts?

2012-03-16 Thread Helios de Rosario

You can tell that levels must not be dropped with scale_x_discrete():

library(ggplot2)
mtcars$cyl-factor(mtcars$cyl)
plot1 - ggplot(mtcars[!mtcars$cyl==4,], aes(cyl))+geom_bar()
plot1 + scale_x_discrete(drop=FALSE)

Or explicitly set the values you want in the x axis with the argument
limits:

plot1 + scale_x_discrete(limits=levels(mtcars$cyl))


Regards,
Helios


 El día 15/03/2012 a las 16:47, Bart6114 bartsmeet...@gmail.com
escribió:
 Hello,
 
 When plotting a barchart with ggplot it drops the levels of the
factor for
 which no counts are available.
 
 For example:
 
 library(ggplot)
 mtcars$cyl-factor(mtcars$cyl)
 ggplot(mtcars[!mtcars$cyl==4,], aes(cyl))+geom_bar()
 levels(mtcars[!mtcars$cyl==4,])
 
 This shows my problem. Because no counts are available for
factorlevel '4',
 the label 4 dissapears from the plot. However, I would still like it
to show
 up, but without a bar (zero observations).
 
 I would like to use this for the presentation of data with a
Likert-like
 scale.
 
 Thanks in advance!
 Bart
 
 --
 View this message in context: 

http://r.789695.n4.nabble.com/Ggplot-barchart-drops-factor-levels-how-to-show-the

 m-with-zero-counts-tp4475417p4475417.html
 Sent from the R help mailing list archive at Nabble.com.

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] Faceted bar plot shows wrong counts (ggplot2)

2012-03-13 Thread Helios de Rosario

Michael,

Thanks for the pointer to the discussion in the ggplot list. It seems
that the reason of this behaviour of facet_grid() is already known and
being discussed by the developers of ggplot2.

facet_grid() reduces the original data frame with unique() before
applying the stats.  If the data frame has any other column that
prevents duplicated rows, counts are correctly computed.

E.g.

diamonds25 - droplevels(diamonds[1:25,]) # Keep all columns

# Everything else as before:
base  -  ggplot(diamonds25,  aes(fill  =  cut))  +
 geom_bar(position  =  dodge)  +
 opts(legend.position  =  none)
 base  +  aes(x  =  cut)  +
 facet_grid(.  ~  color)


Helios

 El día 12/03/2012 a las 20:59, R. Michael Weylandt
michael.weyla...@gmail.com escribió:
 You get the good behavior with
 
 base + aes(x = cut) + facet_wrap(~ color, ncol = 5)
 
 so this seems buggy to me.
 
 If someone here doesn't step forward with more insight, I'd forward
it
 to the ggplot list to see if one of the developers there can give an
 explanation or possibly make the official call that it's a bug.
 
 There was another report of a possible bug in facet_grid() today
that
 could be related:

https://groups.google.com/group/ggplot2/browse_thread/thread/5213ac35da6b36d

 4
 
 Michael
 
 On Mon, Mar 12, 2012 at 7:16 AM, Helios de Rosario
 helios.derosa...@ibv.upv.es wrote:
 I have encountered a problem with faceted bar plots. I have tried
to
 create something like the example explained in the ggplot2 book (see
pp.
 126-128):

 library(ggplot2)
 mpg4  -  subset(mpg,  manufacturer  %in%
 c(audi,  volkswagen,  jeep))
 mpg4$manufacturer  -  as.character(mpg4$manufacturer)
 mpg4$model  -  as.character(mpg4$model)

 base  -  ggplot(mpg4,  aes(fill  =  model))  +
 geom_bar(position  =  dodge)  +
 opts(legend.position  =  none)
 base  +  aes(x  =  model)  +
 facet_grid(.  ~  manufacturer)

 That example works fine; the bar heights are just the same as the
 counts in the table:

 table(mpg4[,1:2])
  model
 manufacturer a4 a4 quattro a6 quattro grand cherokee 4wd gti jetta
new
 beetle
  audi7  8  3  0   0 0
0
  jeep0  0  0  8   0 0
0
  volkswagen  0  0  0  0   5 9
6
  model
 manufacturer passat
  audi0
  jeep0

 But in other cases this does not occur. For instance, take a small
 subset of data(diamonds):

 diamonds25 - droplevels(diamonds[1:25,2:3])
 table(diamonds25)
   color
 cut E F H I J
  Fair  1 0 0 0 0
  Good  1 0 0 1 4
  Very Good 1 0 3 1 4
  Premium   3 1 0 1 0
  Ideal 1 0 0 1 2

 And change the variables mapped in the previous plot:

 base  -  ggplot(diamonds25,  aes(fill  =  cut))  +
 geom_bar(position  =  dodge)  +
 opts(legend.position  =  none)
 base  +  aes(x  =  cut)  +
 facet_grid(.  ~  color)

 I see all bars with height = 1.
 I have ovserved this problem (wrong bar heights, but not always =
1),
 in other cases when all counts are very small or zero.
 What's wrong here?

 Regards,
 Helios

 sessionInfo()
 R version 2.14.2 (2012-02-29)
 Platform: i386-pc-mingw32/i386 (32-bit)

 locale:
 [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252
 [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
 [5] LC_TIME=Spanish_Spain.1252

 attached base packages:
 [1] stats graphics  grDevices utils datasets  methods  
base

 other attached packages:
 [1] ggplot2_0.9.0

 loaded via a namespace (and not attached):
  [1] colorspace_1.1-1   dichromat_1.2-4digest_0.5.1
 grid_2.14.2
  [5] MASS_7.3-17memoise_0.1munsell_0.3
 plyr_1.7.1
  [9] proto_0.3-9.2  RColorBrewer_1.0-5 reshape2_1.2.1
 scales_0.2.0
 [13] stringr_0.6



 INSTITUTO DE BIOMECÁNICA DE VALENCIA
 Universidad Politécnica de Valencia ● Edificio 9C
 Camino de Vera s/n ● 46022 VALENCIA (ESPAÑA)
 Tel. +34 96 387 91 60 ● Fax +34 96 387 91 69
 www.ibv.org 

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
 En cumplimiento de la Ley Orgánica 15/1999 reguladora de la
Protección
 de Datos de Carácter Personal, le informamos de que el presente
mensaje
 contiene información confidencial, siendo para uso exclusivo del
 destinatario arriba indicado. En caso de no ser usted el
destinatario
 del mismo le informamos que su recepción no le autoriza a su
divulgación
 o reproducción por cualquier medio, debiendo destruirlo de
inmediato,
 rogándole lo notifique al remitente.

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

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e

[R] Faceted bar plot shows wrong counts (ggplot2)

2012-03-12 Thread Helios de Rosario

I have encountered a problem with faceted bar plots. I have tried to
create something like the example explained in the ggplot2 book (see pp.
126-128):

library(ggplot2)
mpg4  -  subset(mpg,  manufacturer  %in%
c(audi,  volkswagen,  jeep))
mpg4$manufacturer  -  as.character(mpg4$manufacturer)
mpg4$model  -  as.character(mpg4$model)

base  -  ggplot(mpg4,  aes(fill  =  model))  +
geom_bar(position  =  dodge)  +
opts(legend.position  =  none)
base  +  aes(x  =  model)  +
facet_grid(.  ~  manufacturer)

That example works fine; the bar heights are just the same as the
counts in the table:

table(mpg4[,1:2])
  model
manufacturer a4 a4 quattro a6 quattro grand cherokee 4wd gti jetta new
beetle
  audi7  8  3  0   0 0 
0
  jeep0  0  0  8   0 0 
0
  volkswagen  0  0  0  0   5 9 
6
  model
manufacturer passat
  audi0
  jeep0

But in other cases this does not occur. For instance, take a small
subset of data(diamonds):

diamonds25 - droplevels(diamonds[1:25,2:3])
table(diamonds25)
   color
cut E F H I J
  Fair  1 0 0 0 0
  Good  1 0 0 1 4
  Very Good 1 0 3 1 4
  Premium   3 1 0 1 0
  Ideal 1 0 0 1 2

And change the variables mapped in the previous plot:

base  -  ggplot(diamonds25,  aes(fill  =  cut))  +
geom_bar(position  =  dodge)  +
opts(legend.position  =  none)
base  +  aes(x  =  cut)  +
facet_grid(.  ~  color)

I see all bars with height = 1.
I have ovserved this problem (wrong bar heights, but not always = 1),
in other cases when all counts are very small or zero.
What's wrong here?

Regards,
Helios

sessionInfo()
R version 2.14.2 (2012-02-29)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252
[3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
[5] LC_TIME=Spanish_Spain.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] ggplot2_0.9.0

loaded via a namespace (and not attached):
 [1] colorspace_1.1-1   dichromat_1.2-4digest_0.5.1  
grid_2.14.2
 [5] MASS_7.3-17memoise_0.1munsell_0.3   
plyr_1.7.1
 [9] proto_0.3-9.2  RColorBrewer_1.0-5 reshape2_1.2.1
scales_0.2.0
[13] stringr_0.6



INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] How to create interrupted boxplot

2012-03-12 Thread Helios de Rosario

Jianghong wrote:
 Hello,
 
 I have created two boxplots with following R code. There is one
outlier in
 B group.
 The outlier is 33. But the all other data are between 0 to 4.
 
 How can I skip y-axis around 5 to 25, and expand 0-4 for this case.
Also I
 want keep the outlier in my boxplot.
 I want my boxplot look like the second one, keep the outlier, and
make an
 interrupt of y-axis from 5 to 25.
 
 Thanks,
 Jianghong
 
 
 Example = c( 0.00,0.33,0.75,3.00,2.50,0.50,2.00,33.00)
 Grp_Example =c(A,A,A,B,B,B,B,B)
 
 Example_Data= cbind(Example,Grp_Example)
 attach(Example_Data)
 
 
 boxplot(Example ~ Grp_Example,main=paste(Boxplot of Example),
 pars = list(boxwex = 0.25, staplewex = 0.25, outwex = 0.25))
 
 boxplot(Example ~ Grp_Example,main=paste(Boxplot of
Example),ylim=c(0,4),
 pars = list(boxwex = 0.25, staplewex = 0.25, outwex = 0.25))

Try plotting in two regions of the same device. Quick and dirty
example:

par(mfrow=c(2,1))
boxplot(Example ~ Grp_Example,main=paste(Boxplot of Example),
pars = list(boxwex = 0.25, staplewex = 0.25, outwex = 0.25),
ylim=c(25,35),xaxt=n)
boxplot(Example ~ Grp_Example,
pars = list(boxwex = 0.25, staplewex = 0.25, outwex = 0.25),
ylim=c(0,4))

See  ?layout --instead of par(mfrow) -- and ?par for further
adjustments of region sizes and margins.

Helios

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] Loop

2012-02-22 Thread Helios de Rosario



-- 
Helios de Rosario Martínez
 
 Researcher
 El día 22/02/2012 a las 11:32, Florian Weiler fweile...@jhubc.it
escribió:
 Dear all,
 
 I have a (probably very basic) question. I am imputing data with the
mice
 package, using 10 chains. I can then write out the 10 final values of
the
 chains simply by
 
 name1 - complete(imp, 1)
   :
   :
 name10 - complete(imp,10)
 
 Not a big deal, I just wanted to do that in a little loop as
follows:
 
 for (i in 1:10){
   set[i] - complete(imp,i)}
 
 Yet that doesn't work, I also tried other things like:
 for (i in 1:10){
   set[[i]] - complete(imp,i)}
 
 Again, no success. It only saves a couple of lines of code, but there
must
 be an easy solution, right?
 Thanks,
 Florian

Hi Florian,

Please give more context. There are various reasons why that could
fail. Is your list set defined before starting to assign values?

Try:
set - vector(list,10)
before the loop.



INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

[R] Get default values of a function argument

2012-01-20 Thread Helios de Rosario

Hi, a quick question:

Is there a way to retrieve the default value of a function argument  -
if it exists?
(I know I can see it if I type the function name, but I would like get
the value programaticaly.)

Thanks,


-- 
Helios de Rosario Martínez
 
 Researcher


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] Get default values of a function argument

2012-01-20 Thread Helios de Rosario

That's exactly what I wanted. Thanks!
Helios

 El día 21/01/2012 a las 0:33, David Winsemius
dwinsem...@comcast.net
escribió:

 On Jan 20, 2012, at 6:26 PM, Helios de Rosario wrote:
 
 Hi, a quick question:

 Is there a way to retrieve the default value of a function argument 
-
 if it exists?
 (I know I can see it if I type the function name, but I would like
get
 the value programaticaly.)
 
 ?formals


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

[R] Error in summary.mlm: formula not subsettable

2011-10-26 Thread Helios de Rosario

When I fit a multivariate linear model, and the formula is defined
outside the call to lm(), the method summary.mlm() fails.

This works well:
 y - matrix(rnorm(20),nrow=10)
 x - matrix(rnorm(10))
 mod1 - lm(y~x)
 summary(mod1)
...

But this does not:
 f - y~x
 mod2 - lm(f)
 summary(mod2)
Error en object$call$formula[[2L]] - object$terms[[2L]] -
as.name(ynames[i]) :
  objeto de tipo 'symbol' no es subconjunto

I would say that the problem is in the following difference:
 class(mod1$call$formula)
[1] call
 class(mod2$call$formula)
[1] name

As far as I understand, summary.mlm() creates a list of .lm objects
from the individual columns of the matrices in the .mlm object, and then
it tries to change the second element of object$call$formula, to present
the name of the corresponding column as the response variable. But if
the formula has been defined outside the call to lm(), that element
cannot be modifed that way.

A bug, perhaps?

 sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252
[3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
[5] LC_TIME=Spanish_Spain.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base




-- 
Helios de Rosario Martínez
  Researcher

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

Re: [R] How to Code Random Nested Variables within Two-way Fixed Model in lmer or lme

2011-09-30 Thread Helios de Rosario

Dave, your situation is clearer now. You wrote (see the full context at
the end of the message):

  From this, you will see that I have 4 control sites and 7 treatment

 sites that are measured each week.  All 13 locations have different 
 names, and Location is a random varaible.  Is Location nested within

 Habitat?  I thought it was, but maybe I am wrong.  Perhaps it is a 
 random variable that is not nested?

I think so, or at least that nesting is irrelevant. If all factors were
fixed, to specify a nesting of Location within Habitat, you should write
the term Habitat/Location, but that is equal to Habitat +
Habitat:Location, and since there cannot be two different values of
Habitat for the same Location, that's exactly the same as Habitat +
Location. So you can safely separate both terms.

Moreover, the different type of effect for Habitat and Location makes
the nesting notation unsuitable, because Habitat is a fixed effect,
whereas Habitat:Location (i.e. Location) is random. If you put
Habitat/Location in the fixed part of the formula, both terms would
be treated as fixed, and if you put it in the random part, both would
be treated as random, and you don't want that, I suppose.

 My main goal is to look for an effect of Habitat.  But if there is a

 significant Week x Habitat interaction, I would examine the effect of

 Habitat separately for each Week.
 
 Hopefully, the above helps to clarify my situation.  I should
re-state, 
 I would like to use an lmer or lme syntax to properly analyze these 
 data, especially given that they are counts, I would like to try
family 
 = poisson or quasipoisson.

If you need a generalized linear model, you can try lmer (in package
lme4), since it supports the family argument where you can specify the
type of error distribution. On the other hand, lme (in package nlme)
only considers normal errors.

Regarding the formula, I like building formulas term by term, asking
myself if each possible term has a potential effect, and including it in
the model if I can answer yes. From what you have said, I can infer
that you do think that there may be a Week:Habitat interaction, so that
term must be in your formula. Now:

1. Do you think that the habitat may influence your outcome (the
counts), regardless of the other factors? I guess you do, so let's
include Habitat as well.

2. Do you think that Week may influence the counts, regardless of the
other factors?

  - 2.1. If so, the fixed part of your formula would be Habitat*Week (=
Habitat + Week + Week:Habitat)
  
  - 2.2. Otherwise, it would be Habitat/Week (= Habitat +
Habitat:Week)

As already commented on, the random part would just be Location. In
theory, since each location is measured in various weeks, you might
consider that the (fixed) effect of Week could be influenced by the
random Location as well, and in that case you would have an additional
random term, the Location:Week interaction. (I.e., you could write the
random term as Location/Week.) However, in your data set there is only
one observation for each value of Location:Week, so it would be
impossible to distinguish that random term from the residual error, and
you may just omit it.

All in all, you can try:

m. - lmer(CO ~ Habitat*Week + (1|Location), family=poisson)

or if Week is only relevant for different types of habitat:

m. - lmer(CO ~ Habitat/Week + (1|Location), family=poisson)

I must admit that I'm not used to analyse generalized linear models, so
I don't know if that approach is correct, but I'd say that's the code to
do what you asked for.

Now, the bad news is that perhaps you are expecting to get p-values
from anova(m.), but you won't get it. Douglas Bates explained why here:
https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html 

On the other hand, you can get p-values from Anova (from package car),
instead of anova, but I don't entirely understand the calculations of
that approach.

Helios

 El día 29/09/2011 a las 20:49, Dave Robichaud drobich...@lgl.com
escribió:
 Hi again,
 
 Thank you very much for taking the time to respond to my question.  I
am 
 sorry that my explanation was confusing.  Please allow me to try to
clarify.
 
 First, please ignore my attempts to define a lmer model.  By putting

 forward my best first guess, which was clearly wrong, I have only
served 
 to confuse matters.  My goal here is to get advice on how to
formulate 
 the correct lmer model.  Hopefully someone can help with that.
 
 I should describe my data in more detail.  I have the following
columns:
 
 LocationHabitatWeek CO
 1   Control110
 2   Control112
 3   Control1 0
 4   Control1 5
 5   Treatment  110
 6 Treatment 1 7
 7 Treatment  1 8
 8 Treatment  1 6
 9 Treatment  1 0
 10 Treatment  1 5
 11 Treatment  1 3
 12 Treatment  1 12
 13 Treatment  1 0
 ...(9 weeks of data omitted

[R] F and Wald chi-square tests in mixed-effects models

2011-09-29 Thread Helios de Rosario

I have a doubt about the calculation of tests for fixed effects in
mixed-effects models.

I have read that, except in well-balanced designs, the F statistic that
is usually calculated for ANOVA tables may be far from being distributed
as an exact F distribution, and that's the reason why the anova method
on mer objects (calculated by lmer) do not calculate the denominator
df nor a p-value. --- See for instance Douglas Bates' long post on this
topic, in:
https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html

However, Anova does calculate p-values from Wald chi-square tests for
fixed terms from mer objects (as well as from lme objects, from
lme). I suppose that the key to understand the logic for this is in Fox
 Weisberg's commentary in An R Companion to Applied Regression (2nd
edition, p. 272), where they say: Likelihood ratio tests and F tests
require fitting more than one model to the data, while Wald tests do
not.

Unfortunately, that's too brief a commentary for me to understand why
and how the Wald test can overcome the deficiencies of F-tests in
mixed-effects models. The online appendix of An R Companion... about
mixed-effects models does not comment on hypothesis tests either.

I would appreciate if someone can give some clues or references to read
about this issue.

Thanks,
Helios

INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

Re: [R] How to Code Random Nested Variables within Two-way Fixed Model in lmer or lme

2011-09-29 Thread Helios de Rosario

Dear Dave, there are some inconsistencies in your explanation of the
problem. You said your variables are:

 CO is a continuous response variable,
 
 Week is a fixed categorical factor,
 
 Habitat is a fixed categorical factor, and
 
 Location is a random categorical factor nested within Habitat.

What does this last statement mean? Are the Locations identified with
the same names in both Habitats (e.g. Location=1,2,3,... for
Habitat=Control, and then Location=1,2,3... for Habitat=Treatment,
although the Locations of both Habitats have nothing to do with each
other? Or do all 13 Locations receive different names?

If Location is really nested within Habitat, you would be in the former
case, and then your random terms should include the interaction
Location:Habitat. In the latter case, the random term would just be
Location.

But then, your model with aov is:

 mCO = aov(CO ~ Week * Habitat + Error(Location/Week))

Since you don't include Habitat in the Error term, I would say that
Location is not really nested within Habitat. But then, why is Week
nested within Location? Do you mean that the effect of Week may be
affected by the random Location?

Anyway, your interpretation of the ANOVA table is misleading:

 Error: Location
 
Df Sum Sq Mean Sq F value   Pr(F)   
 
 Habitat 1 182566   182566   8.6519 0.01341 *
 
 Residuals 11 23211521101   
 
   
 
 Error: Location:Week
 
   Df Sum Sq Mean Sq F value Pr(F) 
 
 Week  10 59643159643 11.05347.5e-13 ***
 
 Week:Habitat 10 19634919635   3.6389 0.0003251 ***
 
 Residuals110 593551 5396   

This actually means:

For the F test of Habitat, the denominator MS is that for Location
 
For the F test of Week, the denominator MS is that for the Location x
Week
 
For the F test of Habitat x Week, the denominator MS is that for
Location x Week


And then, you wrote your attempt with lmer:

 m. = lmer(CO ~ Week * Habitat + (1|Habitat/Location))

The random term here (1|Habitat/Location) has nothing to do with the
Error term you used in aov.

If location is really nested within Habitat, perhaps you meant

m. = lmer(CO ~ Week * Habitat + (1|Habitat:Location)) 

(Habitat/Location means that Habitat has a random effect per se as
well, and I guess you don't mean that!)
Or if Location is not really nested,

m. = lmer(CO ~ Week * Habitat + (1|Location))

or if you really wanted the same model as with aov:

m. = lmer(CO ~ Week * Habitat + (1|Location/Week)) 

Please clarify your model. Otherwise it would be impossible to make any
comparison.

Helios





 El día 29/09/2011 a las 6:30, Dave Robichaud drobich...@lgl.com
escribió:
 Hi All,
 
 I am frustrated by mixed-effects model! I have searched the web for 
 hours, and found lots on the nested anova, but nothing useful on my 
 specific case, which is: a random factor (C) is nested within one of
the 
 fixed-factors (A), and a second fixed factor (B) is crossed with the

 first fixed factor:
 
 C/A
 
 A
 
 B
 
 A x B
 
 My question: I have a functioning model using the aov command (see 
 below), and I would now would like to recode it, using a more
flexible 
 command such as lme or lmer. Once I have the equivalent syntax down,
I 
 would ideally like to re-run my analysis using family = poisson, as
CO 
 is actually count data.
 
 I have a dataset including a response variable CO, measured once per

 Week (for 11 weeks) at 13 Locations. The 13 Locations are divided
into 2 
 habitat types (Control and Treatment).
 
 Thus:
 
 CO is a continuous response variable,
 
 Week is a fixed categorical factor,
 
 Habitat is a fixed categorical factor, and
 
 Location is a random categorical factor nested within Habitat.
 
 Here is my model in R:
 
 mCO = aov(CO ~ Week * Habitat + Error(Location/Week))
 
 summary(mCO)
 
 And the output:
 
 Error: Location
 
Df Sum Sq Mean Sq F value   Pr(F)   
 
 Habitat 1 182566   182566   8.6519 0.01341 *
 
 Residuals 11 23211521101   
 
   
 
 Error: Location:Week
 
   Df Sum Sq Mean Sq F value Pr(F) 
 
 Week  10 59643159643 11.05347.5e-13 ***
 
 Week:Habitat 10 19634919635   3.6389 0.0003251 ***
 
 Residuals110 593551 5396   
 
 Given that this is a mixed model, I believe the appropriate error
terms 
 are as follows:
 
 For the F test of Habitat, the denominator MS is that for
location/habitat;
 
 For the F test of Week, the denominator MS is the residual; and
 
 For the F test of Habitat x Week, the denominator MS is the
residual.
 
 My tinkering with lmer and lme have not produced results similar to
the 
 above
 
 For example,
 
 m. = lmer(CO ~ Week * Habitat + (1|Habitat/Location))
 
 anova(m.)
 
 produces:
 
 Analysis of Variance Table
 
Df Sum Sq Mean Sq F value
 
 Week   10 59643159643 11.0534
 
 Habitat1   2865228652   5.3100
 
 Week:Habitat 10 19634919635

Re: [R] How make a x,y dataset from a formula based entry

2011-09-23 Thread Helios de Rosario

To separate the parts of a formula, use as.character
(check the examples in ?character)

Helios

22 Sep 2011 16:14:05 -0400
From: Jean-Christophe BOU?TT? jcboue...@gmail.com
 Hello,
 You can check ?model.frame.
 I do not know however to extract only the right-hand of left-hand
part
 of a formula.
 
 JC
 
 2011/9/22 trekvana trekv...@aol.com:
 Hello all,

 So I am using the (formula entry) method for randomForests:

 randomForest(y~x1+x2+...+x39+x40,data=xxx,...) but the issue is that
some of
 the items in that package dont take a formula entry - you have to
explicitly
 state the y and x vector:

 randomForest(x=xxx[,c('x1','x2',...,'x40')],y=xxx[,'y'],...)

 Now my question is whether there is a function/way to tell R to take
a
 formula and make the two corresponding datasets [x,y] (that way I
dont have
 to create the x dataset manually with all 40 variables I have).

 There must be a more elegant way to do this than
 x=xxx[,c('x1','x2',...,'x40')]

 Thanks!
 George


INSTITUTO DE BIOMECÁNICA DE VALENCIA
Universidad Politécnica de Valencia • Edificio 9C
Camino de Vera s/n • 46022 VALENCIA (ESPAÑA)
Tel. +34 96 387 91 60 • Fax +34 96 387 91 69
www.ibv.org

  Antes de imprimir este e-mail piense bien si es necesario hacerlo.
En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección
de Datos de Carácter Personal, le informamos de que el presente mensaje
contiene información confidencial, siendo para uso exclusivo del
destinatario arriba indicado. En caso de no ser usted el destinatario
del mismo le informamos que su recepción no le autoriza a su divulgación
o reproducción por cualquier medio, debiendo destruirlo de inmediato,
rogándole lo notifique al remitente.

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

[R] Wrapper of linearHypothesis (car) for post-hoc of repeated measures ANOVA

2011-09-22 Thread Helios de Rosario

For some time I have been looking for a convenient way of performing
post-hoc analysis to Repeated Measures ANOVA, that would be acceptable
if sphericity is violated (i.e. leaving aside post-hoc to lme models).

The best solution I found was John Fox's proposal to similar requests
in R-help:
http://tolstoy.newcastle.edu.au/R/e2/help/07/09/26518.html 
http://tolstoy.newcastle.edu.au/R/e10/help/10/04/1663.html 

However, I think that using linearHypothesis() is not as
straightforward as I would desire for testing specific contrasts between
factor levels. The hypotheses must be defined as linear combinations of
the model coefficients (subject to response transformations according to
the intra-subjects design), which may need some involved calculations if
one is thinking on differences between this and that factor levels
(either between-subjects or intra-subjects), and the issue gets worse
for post-hoc tests on interaction effects.

For that reason, I have spent some time in writing a wrapper to
linearHypothesis() that might be helpful in those situations. I copy the
commented code at the end of this message, because although I have
successfully used it in some cases, I would like more knowledgeable
people  to put it to test (and eventually help me create a worthwile
contribution for other people that could find it useful).

This function (which I have called factorltest.mlm) needs the
multivariate linear model and the intrasubject-related arguments (idata,
idesign...) that would be passed on to Anova() (from car) for a repeated
measures analysis, or directly the Anova.mlm object returned by Anova()
instead of idata, idesign... (I have tried to explain it clearly in the
commentaries to the code.)

Moreover, it needs an argument levelcomb: a list that represents the
level combinations of factors to be tested. There are different ways of
representing those combinations (through names of factor levels, or
coefficient vectors/matrices), and depending on the elements of that
list the test is made for main effects, simple effects, interaction
contrasts, etc.

For instance, let me use an example with the OBrienKaiser data set (as
in the help documentation for Anova() and linearHypothesis()).

The calculation of the multivariate linear model and Anova is copied
from those help files:

 phase - factor(rep(c(pretest, posttest, followup), c(5, 5,
5)),
+ levels=c(pretest, posttest, followup))
 hour - ordered(rep(1:5, 3))
 idata - data.frame(phase, hour)
 mod.ok - lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, 
+post.1, post.2, post.3, post.4, post.5, 
+fup.1, fup.2, fup.3, fup.4, fup.5) ~ 
treatment*gender, 
+   data=OBrienKaiser)
 av.ok - Anova(mod.ok, idata=idata, idesign=~phase*hour)

Then, let's suppose that I want to test pairwise comparisons for the
significant main effect treatment (whose levels are
c(control,A,B)). For the specific contrast between treatment A
and the control group I can define levelcomb in the following
(equivalent) ways:

 levelcomb - list(treatment=c(A,control))
 levelcomb - list(treatment=c(A=1,control=-1))
 levelcomb - list(treatment=c(-1,1,0))

Now, let's suppose that I am interested in the (marginally) significant
interaction between treatment and phase. First I test the simple main
effect of phase for different levels of treament (e.g. for the control
group). To do this, levelcomb must have one variable for each
interacting factor (treatment and phase): levelcomb$treatment will
specify the treatment that I want to fix for the simple main effects
test (control), and levelcomb$phase will have a NA value to represent
that I want to test all orthogonal contrasts within that factor:

 levelcomb - list(treatment=control,phase=NA)

I could also use numeric vectors to define the levels of treatment
that I want to fix, as in the previous example, or if I want a more
complicated combination (e.g. if I want to test the phase effect for
pooled treatments A and B):

 levelcomb - list(treatment=c(A=1,B=1),phase=NA)

The NA value can be replaced by the specific contrast matrix that I
would like to use. For instance, the previous statement is equivalent to
the following one:

 levelcomb - list(treatment=c(0,1,1),phase=contr.sum(levels(phase)))

And then, let's see an example of interaction contrast: I want to see
if the difference between the pre-test phase and the following two,
itself differs between the control group and the two treatments. This
would be

 levelcomb - list(treatment=c(2,-1,-1),phase=c(2,-1,-1))
or
 levelcomb -
list(treatment=c(control=2,A=-1,B=-1),phase=c(pretest=2,posttest=-1,followup=-1))

etc.

Now, to perform the test I just use this levelcomb list object with
the model and ANOVA previously performed, in the following way:

 factorltest.mlm(mod.ok,levelcomb,av.ok)

Or if I want to use idata and idesign directly:

 factorltest.mlm(mod.ok,levelcomb,idata=idata,idesign=~phase*hour)

Of course, this function only performs one test at a

Re: [R] [FORGED] Re: Crete stats course

[R] [R-pkgs] Release of phia 0.2-0

Re: [R] Making a trial based Package for 30-days.

Re: [R] about lm()

Re: [R] outliers for Likert scale data

Re: [R] ask help!

[R] Change R options only inside a package

Re: [R] Change R options only inside a package

Re: [R] suppress startup messages from default packages

Re: [R] GLMM post- hoc comparisons

Re: [R] ANOVA repeated measures and post-hoc

[R] Two-way linear model with interaction but without one main effect

Re: [R] Two-way linear model with interaction but without one main effect

Re: [R] replacing NA for zero

Re: [R] anovas ss typeI vs typeIII

Re: [R] MANOVA with random factor

Re: [R] MANOVA with random factor

[R] read.octave fails with data from Octave 3.2.X

Re: [R] read.octave fails with data from Octave 3.2.X

Re: [R] Type II and III sum of squares (R and SPSS)

Re: [R] Ggplot barchart drops factor levels: how to show them with zero counts?

Re: [R] Faceted bar plot shows wrong counts (ggplot2)

[R] Faceted bar plot shows wrong counts (ggplot2)

Re: [R] How to create interrupted boxplot

Re: [R] Loop

[R] Get default values of a function argument

Re: [R] Get default values of a function argument

[R] Error in summary.mlm: formula not subsettable

Re: [R] How to Code Random Nested Variables within Two-way Fixed Model in lmer or lme

[R] F and Wald chi-square tests in mixed-effects models

Re: [R] How to Code Random Nested Variables within Two-way Fixed Model in lmer or lme

Re: [R] How make a x,y dataset from a formula based entry

[R] Wrapper of linearHypothesis (car) for post-hoc of repeated measures ANOVA

33 matches

Site Navigation

Mail list logo

Footer information