Re: [R] (coxph, se) Obtaining standard errors of coefficients from coxph to store

2007-08-27 Thread joris . dewolf


David,

It would be helpful to give an example of what you would like to extract.

I guess you know how to extract elements from vectors and lists.
However, sometimes the objects returned by functions can be rather complex
(output of coxph() is...)
A general method to capture printed output is via capture.output(). Maybe
not fast, but if you have no other solution...

Joris

 a - rnorm(10,1,1)
 b - rnorm(10,1,1)
 mod - lm(a~b)
 smod - summary(mod)
 smod

Call:
lm(formula = a ~ b)

Residuals:
Min  1Q  Median  3Q Max
-1.7482 -0.5991  0.1211  0.8341  1.4975

Coefficients:
Estimate Std. Error t value Pr(|t|)
(Intercept)   1.6210 0.5332   3.040   0.0161 *
b-0.7667 0.5037  -1.522   0.1664
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.142 on 8 degrees of freedom
Multiple R-Squared: 0.2246, Adjusted R-squared: 0.1277
F-statistic: 2.317 on 1 and 8 DF,  p-value: 0.1664

 output - capture.output(print(smod))
 output
 [1] 
 [2] Call:
 [3] lm(formula = a ~ b)
 [4] 
 [5] Residuals:
 [6] Min  1Q  Median  3Q Max 
 [7] -1.7482 -0.5991  0.1211  0.8341  1.4975 
 [8] 
 [9] Coefficients:
[10] Estimate Std. Error t value Pr(|t|)  
[11] (Intercept)   1.6210 0.5332   3.040   0.0161 *
[12] b-0.7667 0.5037  -1.522   0.1664  
[13] ---
[14] Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
[15] 
[16] Residual standard error: 1.142 on 8 degrees of freedom
[17] Multiple R-Squared: 0.2246,\tAdjusted R-squared: 0.1277 
[18] F-statistic: 2.317 on 1 and 8 DF,  p-value: 0.1664 
[19] 





   
 David Lloyd 
 [EMAIL PROTECTED] 
 lloyd.com To 
 Sent by:  r-help@stat.math.ethz.ch  
 [EMAIL PROTECTED]  cc 
 at.math.ethz.ch   
   Subject 
   [R] (coxph,  se) Obtaining standard 
 16/08/2007 11:31  errors of coefficients from coxph   
   to store
   
   
   
   
   
   




Hi all,

I'm wanting to be able to find and store the z-score of coxph below: -

modz=coxph(Surv(TSURV,STATUS)~RAGE+DAGE+REG_WTIME_M+CLD_ISCH+POLY_VS,
data=kidneyT,method=breslow)


I know summary(modz) will give me this, but how do i extract the
standard error or z-score values in a similar way to obtaining the
coefficients by coef(modz) ? I think it must be something to do with
modz$var but I'm having a complete mental blank.

I need this info so I can write a function to use within a bootstrap so
I can record the number of times (proportion) each variable in the Cox
PH model is actually significant over all the bootstrap resamples.

Any assistance is greatly appreciated

DL


Click to find local singles for dating, romance and fun.
http://tagline.bidsystem.com/fc/Ioyw36XJJVs581mfqGSywy0Z69Mq8VM03oVytPu
8otqP84CBZmNX2G/



span id=m2wTlpfont face=Arial, Helvetica, sans-serif size=2
style=font-size:13.5
px___BRGet
the Free email that has everyone talking at a
href=http://www.mail2world.com target=newhttp://www.mail2world.com/abr
font color=#99Unlimited Email Storage #150; POP3 #150; Calendar
#150; SMS #150; Translator #150; Much More!/font/font/span
 [[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] Custom axis

2007-08-01 Thread joris . dewolf


x - 1:10
y - rnorm(10,10,1)
x2 - 3*x + 2
plot(y ~ x, xaxt = n)
axis(side=1,at = x, labels = x2)

Joris








   
 Florent Bresson   
 [EMAIL PROTECTED] 
 frTo 
 Sent by:  r-help@stat.math.ethz.ch
 [EMAIL PROTECTED]  cc 
 at.math.ethz.ch   
   Subject 
   [R] Custom axis 
 01/08/2007 11:49  
   
   
   
   
   




Dear R users,

I would like to draw a plot with a custom scale for the axis. More
precisely, instead of plotting  y on x, I want to plot y on a monotone
function of x (for instance a*x+b). Which command and/or package should I
use in order to get this result?

Thanks

Florent Bresson





_


__
R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] Re : Custom axis

2007-08-01 Thread joris . dewolf


What about the other way round?

x - 1:10
y - rnorm(10,10,1)
x2 - 3*x + 2
plot(y ~ x2, xaxt = n)
axis(side=1,at = x2, labels = x)






   
 Florent Bresson   
 [EMAIL PROTECTED] 
 frTo 
   [EMAIL PROTECTED] 
 01/08/2007 12:14   cc 
   r-help@stat.math.ethz.ch,   
   [EMAIL PROTECTED]
   Subject 
   Re : [R] Custom axis
   
   
   
   
   
   




Maybe I do not explain well what I would like to do. I do not want to
change the labels of the axis, but the scale. What I want is a general
procedure for changing the scale. Its like using a logarithmic scale on a
plot. Labels are the same, but the increases of x along the x-axis are
defined by a known monotone and continuous function.

Florent Bresson

- Message d'origine 
De : [EMAIL PROTECTED] [EMAIL PROTECTED]
À : Florent Bresson [EMAIL PROTECTED]
Cc : r-help@stat.math.ethz.ch; [EMAIL PROTECTED]
Envoyé le : Mercredi, 1 Août 2007, 12h01mn 58s
Objet : Re: [R] Custom axis



x - 1:10
y - rnorm(10,10,1)
x2 - 3*x + 2
plot(y ~ x, xaxt = n)
axis(side=1,at = x, labels = x2)

Joris









 Florent Bresson
 [EMAIL PROTECTED]
 frTo
 Sent by:  r-help@stat.math.ethz.ch
 [EMAIL PROTECTED]  cc
 at.math.ethz.ch
   Subject
   [R] Custom axis
 01/08/2007 11:49









Dear R users,

I would like to draw a plot with a custom scale for the axis. More
precisely, instead of plotting  y on x, I want to plot y on a monotone
function of x (for instance a*x+b). Which command and/or package should I
use in order to get this result?

Thanks

Florent Bresson





_



__
R-help@stat.math.ethz.ch 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.








_

Ne gardez plus qu'une seule adresse mail ! Copiez vos mails vers Yahoo!
Mail

__
R-help@stat.math.ethz.ch 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] lme or gls prediction intervals

2007-07-25 Thread joris . dewolf


Martin,


Have you checked ?intervals.gls


This intervals are approximate, but this would be the obvious starting
point to me.


Joris










   
 Martin Henry H.  
 Stevens  
 [EMAIL PROTECTED]  To 
 edu  R-Help r-help@stat.math.ethz.ch   
 Sent by:   cc 
 [EMAIL PROTECTED] 
 at.math.ethz.ch   Subject 
   [R] lme or gls prediction intervals 
   
 24/07/2007 19:22  
   
   
   
   




Hi folks,
I am trying to generate 95% confidence intervals for a gls model
using predict.nlme
with
R version 2.5.1 (2007-06-27)
. nlme: Linear  and Nonlinear Mixed Effects Models. R package version
   3.1-83.

I have looked in help, and I can do it for lm and glm models, and I
can generate simple predictions for lme models with various levels --
I am familiar with the basics.

Is there a way to get prediction intervals for gls models? My best
model uses varPower(), so I am reluctant to fall back on lm predictions.

Thank you,

Hank


Dr. Hank Stevens, Associate Professor
338 Pearson Hall
Botany Department
Miami University
Oxford, OH 45056

Office: (513) 529-4206
Lab: (513) 529-4262
FAX: (513) 529-4243
http://www.cas.muohio.edu/~stevenmh/
http://www.muohio.edu/ecology/
http://www.muohio.edu/botany/

E Pluribus Unum

__
R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] The 'REP' term in AMMI{agricolae}

2007-07-24 Thread joris . dewolf


The AMMI senso strictu part starts from the corrected data for genotype and
environment. In most cases where AMMI is applied (maybe also in the
agricolae package, I haven't checked), starts from the interaction effects
obtained through a general linear model anova.
It should be possible to replace the input by blups from your mixed model.
R news 7/1 (p14-19) gave some basic code for AMMI. There you should be able
to get ideas how to modify the code.


Joris








   
 CG Pettersson   
 [EMAIL PROTECTED] 
 e.slu.se  To 
 Sent by:  r-help@stat.math.ethz.ch
 [EMAIL PROTECTED]  cc 
 at.math.ethz.ch   
   Subject 
   [R] The 'REP' term in   
 23/07/2007 19:59  AMMI{agricolae} 
   
   
   
   
   
   




Dear all,

W2k, R 2.5.1

I am trying out the AMMI function in the agricolae package, to analyse the
dependence of environment for a certain cultivar. The function responds to
four basic variables:

ENV Environment
GEN Genotype
REP Replication
Y Response

My question concerns the 'REP' term. When I normally do an analysis of
variance, the replication would refer to repeated observations within the
each field trial. In the case I am analysing right now, I have five years
of observations for each 'ENV', in such a case it feels natural to use
year as the most important replication of the environment- especially as I
am interested in long term trends. This approach also seems to work
allright.

But the field trials are also structured in three main blocks, subdivided
into five 'lattice' blocks, a structure that is powerful in the analysis
of variance. (I use a random call in lme{nlme}).

Is it possible to use the block structure also in the AMMI analysis? If
that is possible, how should I code? I have tried to find out in the
documentation, but if it is stated there I do not understand it.

Thank you
/CG

--
CG Pettersson, PhD
Swedish University of Agricultural Sciences (SLU)
Dept. of Crop Production Ecology. Box 7043.
SE-750 07 Uppsala, Sweden
[EMAIL PROTECTED]

__
R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] extract from anova

2007-07-16 Thread joris . dewolf


Mustapha,

You should tell us what fm1 actually is (which class), before we can give
you a correct answer.
It is helpfull to use str() in this case.

str(fm1) will give you the structure of fm1, and with that information you
extract any value you want.

I suppose that this is what you need:

summary(fm1)[[1]][Residuals,Mean Sq]








   
 elyakhlifi
 mustapha  
 elyakhlifi_musta  To 
 [EMAIL PROTECTED] R-help@stat.math.ethz.ch
 Sent by:   cc 
 [EMAIL PROTECTED] 
 at.math.ethz.ch   Subject 
   [R] extract from anova  
   
 16/07/2007 15:01  
   
   
   
   




Hello,
I would like extract the Mean Sq of the Residuals how should I do  please?
thanks.


 summary(fm1)
Df  Sum Sq Mean Sq F valuePr(F)
groupe  20 300.987  15.049  22.853 3.369e-16 ***
Residuals   41  27.000   0.659
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



_


 [[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] A More efficient method?

2007-07-04 Thread joris . dewolf


or
Cat - c('a','a','a','b','b','b','a','a','b')
C1 - (Cat=='a')*1







   
 ONKELINX,
 Thierry  
 Thierry.ONKELINX  To 
 @inbo.be Keith Alan Chamberlain
 Sent by:  [EMAIL PROTECTED],   
 [EMAIL PROTECTED] r-help@stat.math.ethz.ch  
 at.math.ethz.chcc 
   
   Subject 
 04/07/2007 17:17  Re: [R] A More efficient method?
   
   
   
   
   
   




Cat - c('a','a','a','b','b','b','a','a','b')
C1 - ifelse(Cat == 'a', -1, 1)



ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
[EMAIL PROTECTED]
www.inbo.be

Do not put your faith in what statistics say until you have carefully
considered what they do not say.  ~William W. Watt
A statistical analysis, properly conducted, is a delicate dissection of
uncertainties, a surgery of suppositions. ~M.J.Moroney



 -Oorspronkelijk bericht-
 Van: [EMAIL PROTECTED]
 [mailto:[EMAIL PROTECTED] Namens Keith Alan
 Chamberlain
 Verzonden: woensdag 4 juli 2007 15:45
 Aan: r-help@stat.math.ethz.ch
 Onderwerp: [R] A More efficient method?

 Dear Rhelpers,

 Is there a faster way than below to set a vector based on
 values from another vector? I'd like to call a pre-existing
 function for this, but one which can also handle an
 arbitrarily large number of categories. Any ideas?

 Cat=c('a','a','a','b','b','b','a','a','b')   # Categorical
variable
 C1=vector(length=length(Cat))# New vector for numeric values

 # Cycle through each column and set C1 to corresponding value of Cat.
 for(i in 1:length(C1)){
if(Cat[i]=='a') C1[i]=-1 else C1[i]=1
 }

 C1
 [1] -1 -1 -1  1  1  1 -1 -1  1
 Cat
 [1] a a a b b b a a b

 Sincerely,
 KeithC.
 Psych Undergrad, CU Boulder (US)
 RE McNair Scholar

 __
 R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] unequal variance assumption for lme (mixed effect model)

2007-07-02 Thread joris . dewolf



gls() from the package nlme is similar to lm but is meant for models
without random effects.

Joris





   
 Spencer Graves
 [EMAIL PROTECTED] 
 df.comTo 
 Sent by:  shirley zhang   
 [EMAIL PROTECTED] [EMAIL PROTECTED] 
 at.math.ethz.chcc 
   R-help@stat.math.ethz.ch
   Subject 
 01/07/2007 23:30  Re: [R] unequal variance assumption 
   for lme (mixed effect model)
   
   
   
   
   
   




  The 'weights' argument on 'lm' is assumed to identify a vector of
the same length as the response, giving numbers that are inversely
proportional to the variance for each observation.

  However, 'lm' provides no capability to estimate weights.  If you
want to do that, the varFunc capabilities in the 'nlme' package is the
best tool I know for that purpose.

  If someone thinks there are better tools available for estimating
heterscedasticity, I hope s/he will enlighten us both.

  Hope this helps.
  Spencer Graves

shirley zhang wrote:
 Thanks for Spencer and Simon's help.  I've got very interesting
 results based on your suggestions.

 One more question,  how to handle unequal variance problme in lm()?
 Isn't the weights option also, which means weighted least squares,
 right?  Can you give me an example of setting this parameter in lm()
 to account for  different variance assumption in each group?

 Thanks again,
 Shirley


 On 6/29/07, Spencer Graves [EMAIL PROTECTED] wrote:
 comments in line

 shirley zhang wrote:
  Hi Simon,
 
  Thanks for your reply. Your reply reminds me that book. I've read it
  long time ago, but haven't  try the weights option in my projects
  yet:)
 
  Is the heteroscedastic test always less powerful because we have to
  estimate the within group variance from the given data?
 
 SG:  In general, I suspect we generally lose power when we estimate more
 parameters.

 SG:  You can check this using the 'simulate.lme' function, whose use is
 illustrated in the seminal work reported in sect. 2.4 of Pinheiro and
 Bates (2000) Mixed-Effects Models in S and S-Plus (Springer).
  Should we check whether each group has equal variance before using
  weights=varIdent()? If we should, what is the function for linear
  mixed model?
 
 SG:  The general advice I've seen is to avoid excessive
 overparameterization of heterscedasticity and correlations.  However,
 parsimonious correlation had heterscedasticity models would likely be
 wise.  Years ago, George Box expressed concern about people worrying too
 much about outliers, which are often fairly obvious and relatively easy
 to detect, while they worried too little, he thought, about dependence,
 especially serial dependence, which is generally more difficult to
 detect and creates bigger problems in inference than outliers.  He
 wrote, Why worry about mice when there are tigers about?

 SG:  Issues of this type can be fairly easily evaluated using
 'simulate.lme'.

  Hope this helps.
  Spencer Graves
  Thanks,
  Shirley
 
  On 6/27/07, Simon Blomberg [EMAIL PROTECTED] wrote:
 
  The default settings for lme do assume equal variances within groups.
  You can change that by using the various varClasses. see
 ?varClasses. A
  simple example would be to allow unequal variances across groups.
 So if
  your call to lme was:
 
  lme(...,random=~1|group,...)
 
  then to allow each group to have its own variance, use:
 
  lme(...,random=~1|group, weights=varIdent(form=~1|group),...)
 
  You really really should read Pinheiro  Bates (2000). It's all
 there.
 
  HTH,
 
  Simon.
 
  , On Wed, 2007-06-27 at 21:55 -0400, shirley zhang wrote:
 
  Dear Douglas and R-help,
 
  Does lme assume normal distribution AND equal variance among groups
  like anova() does? If it does, is there any method like unequal
  variance T-test (Welch T) in lme when each group has unequal
 variance
  in my data?
 
  Thanks,
  Shirley
 
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE 

Re: [R] Linear Mixed Models with nlme, more than one random effect

2007-06-20 Thread joris . dewolf


Guilio,

Have a look at Rnew volume 5/1 (http://cran.r-project.org/doc/Rnews/) where
Doug Bates explains this nicely. Condider using lme4 for your purpose.

Joris








   
 Giulio Di
 Giovanni 
 [EMAIL PROTECTED]  To 
 otmail.com   r-help@stat.math.ethz.ch
 Sent by:   cc 
 [EMAIL PROTECTED] 
 at.math.ethz.ch   Subject 
   [R] Linear Mixed Models with nlme,  
   more than one random effect 
 20/06/2007 14:09  
   
   
   
   
   





Hi, I' trying to learn how to use lme for Linear Mixed Models  and I have a

problem when I have to include more than one random effect in my model. I
know that this could be a stupid question to ask, but I'm not able to solve

it by myself... One example: if my model is

response = operator + block + day
with operator and block as fixed effects and day as random effect, I use

res.lme - lme(resp ~ oper + block , random=~1|day)

If I want to include also another random effect, as experiment, what I
should do ?
This effect doesn't have to be nested, at the and I would like to have the
COV matrix using (if I'm not wrong) getVarCov function.

Thanks in advance for any help or suggestions, I'm a beginner on this
field...

Davide

_
Cinema, Tv, Gossip e Orsoscopo…Tutto su MSN intrattenimento!

__
R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] Inverse BoxCox transformation

2007-06-18 Thread joris . dewolf


to backtransform 'estimate':




if (lambda == 0 ) {


log(estimate)


} else {



estimate^(1/lambda)

}











   
 Des Callaghan   
 [EMAIL PROTECTED] 
 er.co.uk  To 
 Sent by:  r-help@stat.math.ethz.ch  
 [EMAIL PROTECTED]  cc 
 at.math.ethz.ch   
   Subject 
   [R] Inverse BoxCox transformation   
 18/06/2007 09:32  
   
   
   
   
   




Hi,

I can't seem to find a function in R that will reverse a BoxCox
transformation. Can somebody help me locate one please? Thanks in advance.

Best wishes,
Des


 [[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] str() to extract components

2007-02-27 Thread joris . dewolf

mod - lmer(resp3~b$age+b$size+b$pcfat+(1|sex), data=b)
coef(mod)[1]$Subject[1,1]


[EMAIL PROTECTED] wrote on 27/02/2007 16:52:54:

 Hi,

 I have been dabbling with str() to extract values from outputs such as
 lmer etc and have found it very helpful sometimes.

 but only seem to manage to extract the values when the output is one
 simple table, any more complicated and I'm stumped :-(

 take this example of the extracted coeficients from a lmer analysis...

 using str(coef(lmer(resp3~b$age+b$size+b$pcfat+(1|sex), data=b))) yields

 Formal class 'lmer.coef' [package Matrix] with 3 slots
   ..@ .Data :List of 1
   .. ..$ :`data.frame': 2 obs. of  4 variables:
   .. .. ..$ (Intercept): num [1:2] 1.07 1.13
   .. .. ..$ b$age  : num [1:2] 0.00702 0.00702
   .. .. ..$ b$size : num [1:2] 0.0343 0.0343
   .. .. ..$ b$pcfat: num [1:2] 0.0451 0.0451
   ..@ varFac: list()
   ..@ stdErr: num(0)

 how do I get inside the first table to get the value 1.07 for instance?

 Any help much appreciated.


 Simon Pickett
 PhD student
 Centre For Ecology and Conservation
 Tremough Campus
 University of Exeter in Cornwall
 TR109EZ
 Tel 01326371852

 __
 R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] Hierarchical ANOVA

2007-02-13 Thread joris . dewolf

Romain,

Look for info on mixed models.
In R you do this either with the library nlme or lme4.

A good starting point is an article by Doug Bates in Rnews
http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf

Joris


[EMAIL PROTECTED] wrote on 13/02/2007 17:16:41:


 Hello,

 Does somebody could help me in the computation(formulation) of a
 hierarchical ANOVA using linear model in R?
 I'm working in a population biology study of an endangered species. My
aim
 is to see if I have effects of Density of individuals/m2 on several
 measured plants fitness traits.

 As independent variables I have:

 Humidity (continuous)
 Nutritive substance (continuous)
 LUX (continuous)
 Density of individuals/m2 (continuous)
 Population (categorical)
 Year (categorical)

 I want to test the 4 continous variables with Population as error term.
And
 the Population, the Year and interaction term, with the error of
Population
 x Year.

 Thank you for any help, best regards.

 Romain Mayor

 __
 R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] changes in the structure of mer objects?

2007-01-15 Thread joris . dewolf
Thanks Martin,

R --vanilla did the trick.
By the way, is there a way to check which version of a lme4 or Matrix I am
using in a certain instance of R?

We are running multiple versions of R on the same server until we are sure
all our operational code is behaving well under a new version of R or of a
package.
Now I start doubting if we ever have been testing what we were intendng to
test...

Joris


[EMAIL PROTECTED] wrote on 15/01/2007 14:00:12:

 Hi Joris,

 I suspect you somehow load an older version lme4 or Matrix than
 you think you are loading.
 Or then you have an lmer() function or a class definition
 {from a saved workspace }
 in your work space.

 example(lmer) *must* run correctly for the 'lme4' package to get
 onto CRAN at all,
 hence it must be unique to your setup.

 Maybe as a first step,
 run R as R --vanilla ?

 Regards,
 Martin Maechler, ETH Zurich

  joris == joris dewolf [EMAIL PROTECTED]
  on Sun, 14 Jan 2007 19:15:43 +0100 writes:

 joris Dear all,

 joris I try to run the example of lmer and get the following
 error message.

  library(lme4)
  example(lmer)
 lmer (fm1 - lmer(Reaction ~ Days + (Days | Subject), sleepstudy))
 joris [[1]]
 joris Error in get(x, envir, mode, inherits) : variable as.
 dpoMatrix was not
 joris found

 joris This error message is similar to what I get with
 joris other models. It looks like the mer class has a
 joris slightly different structure. Anybody an idea how to
 joris solve this?


 joris I am using R 2.4.1 under linux and the latest releases of
 lme4 and Matrix

 joris lme4_0.9975-10
 joris Matrix_0.9975-8

  version
 joris _
 joris platform   x86_64-unknown-linux-gnu
 joris arch   x86_64
 joris os linux-gnu
 joris system x86_64, linux-gnu
 joris status
 joris major  2
 joris minor  4.1
 joris year   2006
 joris month  12
 joris day18
 joris svn rev40228
 joris language   R
 joris version.string R version 2.4.1 (2006-12-18)

 __
 R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] changes in the structure of mer objects?

2007-01-14 Thread joris . dewolf

Dear all,

I try to run the example of lmer and get the following error message.

 library(lme4)
 example(lmer)
lmer (fm1 - lmer(Reaction ~ Days + (Days | Subject), sleepstudy))
[[1]]
Error in get(x, envir, mode, inherits) : variable as.dpoMatrix was not
found

This error message is similar to what I get with other models. It looks
like the mer class has a slightly different structure. Anybody an idea how
to solve this?


I am using R 2.4.1 under linux and the latest releases of lme4 and Matrix

lme4_0.9975-10
Matrix_0.9975-8

 version
   _
platform   x86_64-unknown-linux-gnu
arch   x86_64
os linux-gnu
system x86_64, linux-gnu
status
major  2
minor  4.1
year   2006
month  12
day18
svn rev40228
language   R
version.string R version 2.4.1 (2006-12-18)

Thanks

Joris De Wolf
Phone: +32 9 2429155, E-Mail: [EMAIL PROTECTED]
Postal Address: CropDesign N.V. Technologiepark 3, 9052 Gent Belgium

__
R-help@stat.math.ethz.ch 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 extract the variance componets from lme

2007-01-02 Thread joris . dewolf
I advice you strongly to use VarCorr(), but if you insist

tmp - as.matrix(m1$modelStruct$reStruct$MI)
c(sqrt(diag(tmp)), Residual = 1) * m1$sigma


Joris

[EMAIL PROTECTED] wrote on 02/01/2007 12:50:05:

 Here is a piece of code fitting a model to a (part) of a dataset, just
 for
 illustration. I can extract the random interaction and the residual
 variance
 in group meth==1 using VarCorr, but how do I get the other residual
 variance?

 Is there any way to get the other variances in numerical form directly -
 it
 seems a litte contraintuitive to use as.numeric when extracting
 estimates, it's
 a bit like good old days writing SAS-programs that reads the SAS output
 files...

 Bendix
 
 ---

 library( nlme )

 dfr -
 structure(list(meth = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1,
 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2), .Label = c(CO,
 pulse), class = factor), item = c(1, 1, 1, 2, 2, 2, 3, 3,
 3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4), repl = c(1,
 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1,
 2, 3), y = c(78, 76.4, 77.2, 68.7, 67.6, 68.3, 82.9, 80.1, 80.7,
 62.3, 65.8, 67.5, 71, 72, 73, 68, 67, 68, 82, 77, 77, 43, 69,
 77)), .Names = c(meth, item, repl, y), row.names = c(1,
 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 184,
 185, 186, 187, 188, 189, 190, 191, 192, 193,
 194, 195), class = data.frame)

 m1 -
 lme( y ~ factor( meth ) + factor( item ),
  random = ~1 | MI,
  weights = varIdent( form = ~1 | meth ),
  method =REML,
  data = cbind( dfr, MI=interaction( dfr$meth, dfr$item ) ) )

 m1

 # The MI std and the residual std for meth==1
 as.numeric(VarCorr(m1)[,2])

 __

 Bendix Carstensen
 Senior Statistician
 Steno Diabetes Center
 Niels Steensens Vej 2-4
 DK-2820 Gentofte
 Denmark
 +45 44 43 87 38 (direct)
 +45 30 75 87 38 (mobile)
 +45 44 43 73 13 (fax)
 [EMAIL PROTECTED]   http://www.biostat.ku.dk/~bxc

 __
 R-help@stat.math.ethz.ch 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-help@stat.math.ethz.ch 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] contrasts involving random terms in lme

2005-01-10 Thread Joris DeWolf
Hello,
Has somebody an idea of how to fit contrast involving random terms and 
obtain their standard errors with lme?

I am referring to
Welham, S.; Cullis, B.; Gogel, B.; Gilmour, A. and Thompson, R. (2004).
Prediction in linear mixed models.
Australian and New Zealand Journal of Statistics 46, 325-347
who develop the ideas of
McLean, R.A.; Sanders,W.L.; and Stroup,W.W. (1991).
A unified approach to linear mixed models.
American Statistician 45, 5464.
Any hint would be welcome.
Joris
--
==
Joris De Wolf
CropDesign N.V.
Plant Evaluation Group
Technologiepark 3
B-9052 Zwijnaarde
Belgium
Tel. : +32 9 242 91 55
Fax  : +32 9 241 91 73
==

confidentiality notice:
The information contained in this e-mail is confidential and...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Lattice: how to index in a custom panel function?

2004-08-06 Thread Joris DeWolf
Hi,
I have a lattice xyplot that contains panels according to FactorA, and 
curves for the 2 levels of Factor B within a panel.
I try to add text in the panels of a lattice graph. I suppose I have to 
write a custom function (panel.txt).
What I really would like is to adapt the text in the panel according to 
the levels of FactorA.
In the manuals, I find examples for the strips using which.given and 
which.panel. But this does not work for the panels...

Can Anybody give a hint? Thanks
Joris
I work with R 1.9.1 under Linux
panel.txt = function(sometxt,x,y,...){
   grid.text(sometxt,x,y)
}
xyplot(data = Pdata, P ~ DAS | FactorA,
  groups = FactorB,
  type =s,
  col = c(red,blue),
 panel = function(x,y,...){
   panel.abline(h = 100, lty = 5, lwd =0.5, col = darkgrey)
   panel.txt(mytext, 0.2, 0.8)
}


--
== 
Joris De Wolf
CropDesign N.V. 
Plant Evaluation Group
Technologiepark 3 
B-9052 Zwijnaarde 
Belgium 
Tel. : +32 9 242 91 55
Fax  : +32 9 241 91 73
== 


confidentiality notice:
The information contained in this e-mail is confidential and...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] difference in order() between Linux and Windows with mixtures of caps and normal letters

2004-06-25 Thread Joris DeWolf
Could anybody explain this difference in the function order() between R 
under windows and R under Linux ?
It seems that under Linux the order in of character dsitinguishes 
between caps and normal letters and sorts them starting with the 
capitals (see below).
How can I avoid this (I mean: get the same results as I get under Windows)?
Thanks
Joris

Under Windows R 1.6.2
 mstring - c(b,c,a,t,d)
 mstring[order(mstring)]
[1] a b c d t
 Mstring - c(b,c,a,T,D)
 Mstring[order(Mstring)]
[1] a b c D T
the same under Linux R.1.6.2
 mstring - c(b,c,a,t,d)
 mstring[order(mstring)]
[1] a b c d t
 Mstring - c(b,c,a,T,D)
 Mstring[order(Mstring)]
[1] D T a b c


--
== 
Joris De Wolf
CropDesign N.V. 
Plant Evaluation Group
Technologiepark 3 
B-9052 Zwijnaarde 
Belgium 
Tel. : +32 9 242 91 55
Fax  : +32 9 241 91 73
== 


confidentiality notice:
The information contained in this e-mail is confidential and...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Problems with postscript output

2004-03-24 Thread Joris DeWolf
This is not a PS problem. The barplot you made has its horizontl axis 
suppressed by default. See argument axis.lty in ?barplot
Joris

Frank Gerrit Zoellner wrote:

Hi all!

I have a little problem with saving plots to file.
I use the command postscript() followed by the plotting command and a dev.off().
When I then look at the resulting image saved to disk, some of the axis labels are missing (see attached image). Is there a way to fix this. 

Yours,
 



__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
--

== 
Joris De Wolf
CropDesign N.V. 
Plant Evaluation Group
Technologiepark 3 
B-9052 Zwijnaarde 
Belgium 
Tel. : +32 9 242 91 55
Fax  : +32 9 241 91 73
== 

confidentiality notice:
The information contained in this e-mail is confidential and...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Adding data.frames together

2004-03-10 Thread Joris DeWolf
check rbind
Joris
John Sweval wrote:

I have a series of data frames that are identical structurally, i.e. -
made with the same code, but I need to add them together so that they
become one, longer, data frame, i.e. - each of the slot vectors are
increased in length by the length of the added data frame vectors.
So if I have df1 with a slot A so that length(df1$A) = 100 and I have
df2 with a slot A so that length(df2$A)=200 then I need a method to
create df3 its slot A is the df1$A plus df2$A such that length(df3$A) =
300.
It does not appear that if you use data.frame to join two data frames it
just adds the slots of both sources to the destination data frame and
that is not what I want.
In my finally solution, I need to do this with multiple data.frames that
slot-wise are identical, but each slot length is different between data
frames.
Seems like there should be an easy solution, but I just have not
stumbled across it in the documentation.
Thanks,
John C. Sweval
 Database Architect
 Illumigen Biosciences, Inc.
   Email: [EMAIL PROTECTED]
   Phone: 206-378-0400
   Fax: 206-378-0408
	[[alternative HTML version deleted]]

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
 

--

== 
Joris De Wolf
CropDesign N.V. 
Plant Evaluation Group
Technologiepark 3 
B-9052 Zwijnaarde 
Belgium 
Tel. : +32 9 242 91 55
Fax  : +32 9 241 91 73
== 

confidentiality notice:
The information contained in this e-mail is confidential and...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html