Re: [R] Package creation - require statement?

[Prof Brian Ripley]

On Fri, 18 Dec 2009, Jeff Breiwick wrote:


All,

I have used package.skeleton, edited the .Rd files and the DESCRIPTON file.
In one of my package functions I am using RODBC so I have a line of code
that reads: require(RODBC).

In my DESCRIPTION file I have the following relevant line:
Depends: R (>= 2.9.0), RODBC

However, when I do a check (Rcmd check my_package_name) I get this warning:
* checking for unstated dependencies in R code ... WARNING
'library' or 'require' calls not declared from:
 RODBC
See the information on DESCRIPTION files in the chapter 'Creating R
packages' of the 'Writing R Extensions' manual.

Does anyone know where I have gone astray? Thank you.


Using require(RODBC): the Depends declaration already loaded RODBC.


Jeff Breiwick


--
Brian D. Ripley,  rip...@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
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] model matrix with a spline

[David Hugh-Jones]

Hi all

I want to get the design matrix for a model, evaluated at a single value.
For example, if I pass in a data frame with a=2, b=2, y=3, and my
model is y ~ a+b+a:b, then I would like to get
the values 3, 2, 2, 4 out. I can do this with:

tmp <- model.matrix(fit, data=mydata)

or

tmp <- predict(fit, newdata=mydata, type="terms")

However, if the fit had a smoothing spline component, this fails. It
seems like the prediction function is trying to reevaluate the basis
for the spline, and as there is only one row in the new data, it can't
do that.

Is there a way I can get the value of the already-created spline? And
is there a simple way to do this programmatically so I don't need to
check each term of the formula individually?

Cheers,

David Hugh-Jones
Post-doctoral Researcher
Max Planck Institute of Economics, Jena
http://davidhughjones.googlepages.com

__
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] ggplot2 / lattice

[Gary Miller]

Thanks Felix, I got others except:

1. Is it possible to reduce the gaps between the boxplots, may be say "NO
GAP".

2. I tried col="red" option, but its not changing the border color of boxes.

R Code:

g <- rep.int(c("A", "B", "C", "D"), 125)
t <- rnorm(5000)
a <- sample(t, 500, replace=TRUE)
b <- sample(t, 500, replace=TRUE)
dta <- data.frame(val = sample(t,1000), g = gl(4, 250, labels=c("A", "B",
"C", "D")) , G2 = gl(2,1, labels=c("XX", "YY")))
library(lattice)
bwplot(val ~ g | G2, dta,
 panel = function(x,y,...)
 {
  panel.bwplot(x,y,..., pch="|", col="red", fill="red",
  box.width=.3)
  }
)

Thanks,
Gary

On Thu, Dec 17, 2009 at 11:01 PM, Felix Andrews  wrote:

> Look at ?panel.bwplot specifically 'fill', 'pch', and 'box.width'.
>
>
>
> 2009/12/18 Gary Miller :
>  > @ Gabor: Thanks, it worked.
> >
> > @ All: Which options I should explore to change the appearance of boxes.
> > Like,
> >
> > 1. Color change.
> > 2. Fill-in boxes with colors.
> > 3. Change dot to line (used to represent median)
> > 4. No gaps between boxplots in each category.
> >
> > # R Code:
> > g <- rep.int(c("A", "B", "C", "D"), 125)
> > t <- rnorm(5000)
> > a <- sample(t, 500, replace=TRUE)
> > b <- sample(t, 500, replace=TRUE)
> > dta <- data.frame(val = sample(t,1000), g = gl(4, 250, labels=c("A", "B",
> > "C", "D")) , G2 = gl(2,1, labels=c("XX", "YY")))
> > library(lattice)
> > bwplot(val ~ g | G2, dta)
> >
> > Thanks,
> > Gary
> >
> > On Thu, Dec 17, 2009 at 10:09 PM, Gabor Grothendieck <
> > ggrothendi...@gmail.com> wrote:
> >
> >> Try this:
> >>
> >> bwplot(val ~ g | G2, dta)
> >>
> >>
> >> On Thu, Dec 17, 2009 at 9:34 PM, Gary Miller  >
> >> wrote:
> >> > Hi R Users,
> >> >
> >> > Is there a equivalent function for the following R code in Lattice
> >> package.
> >> > I want to plot a grouped box plots (grouped by two factors).
> >> >
> >> > g <- rep.int(c("A", "B", "C", "D"), 125)
> >> > t <- rnorm(5000)
> >> > a <- sample(t, 500, replace=TRUE)
> >> > b <- sample(t, 500, replace=TRUE)
> >> > dta <- data.frame(val = sample(t,1000), g = gl(4, 250, labels=c("A",
> "B",
> >> > "C", "D")) , G2 = gl(2,1, labels=c("XX", "YY")))
> >> > # ggplot2 package required
> >> > p <- ggplot(dta, aes(factor(G2), val))
> >> > p + geom_boxplot(aes(fill = factor(g)))
> >> >
> >> > Any suggestion would be highly appreciated!
> >> > Gary
> >> >
> >> >[[alternative HTML version deleted]]
> >> >
> >> > __
> >> > 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.
> >> >
> >>
> >
> >[[alternative HTML version deleted]]
> >
> > __
> > 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.
> >
>
>
>
> --
> Felix Andrews / 安福立
> Postdoctoral Fellow
> Integrated Catchment Assessment and Management (iCAM) Centre
> Fenner School of Environment and Society [Bldg 48a]
> The Australian National University
> Canberra ACT 0200 Australia
> M: +61 410 400 963
> T: + 61 2 6125 4670
> E: felix.andr...@anu.edu.au
> CRICOS Provider No. 00120C
> --
> http://www.neurofractal.org/felix/
>

[[alternative HTML version deleted]]

__
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] Error when running Conditional Logit Model

[Charles C. Berry]

On Fri, 18 Dec 2009, Hien Nguyen wrote:


Thanks a lot for answering my questions.

I have tried to run the clogit for only 64 observations and 4 independent 
variables and the results are solved instantly. However, when I run the same 
command (with only 4 dependent variables) for the full data, it keeps running 
for 50 minutes now. :(


Thomas, what do you mean by "maximizing the unconditional likelihood is fine 
when the stratum sizes are large"? What I put in "strata (__)" is actually 
the possible choices (1-64). Each choices will be recored more than 4000 
times (which means I have more than 4000 values of 1, 4000 values of 2 and so 
on).

Does it sound right?


So you have 64 cases and more than 25 controls.

Large strata will really slow down clogit. But I think that that isn't 
your problem.


If the strata really matter - in the sense that the conditional 
distributions of covariates for controls vary a lot from stratum to 
stratum - then you really gain little by having more than a handful of 
controls for each case. If that is the situation you are in, sampling a 
couple of dozen controls from the stratum of each case will give you 
results that are very nearly as precise as those obtained from using all 
4000 of them:


plot( 1:100, (1 + 1/1:100), xlab='n of controls',
ylab='relative variance of coef' )


will give you rough idea of the impact of increasing the number of 
controls per case. The variance with 1 control per case is 2; at the 
asymptote it is 1.


So you can probably spend things up a lot by using fewer controls with 
little loss in accuracy.


With only 64 cases you cannot fit terribly complicated models. This holds 
whether you approach things conditionally using clogit or unconditionally 
using glm. Fourteen degrees of freedom for regression is probably pushing 
matters.  ridge() is helpful in taming overlarge regressor sets in clogit, 
but you'll need to use survival:::summary.coxph.penal() on the result (or 
tinker with the class attribute).


BTW, when you say 'strata(___)', I hope you mean that you use something 
like 'strata( stratvar )' where stravar is a factor that encodes the 
64 levels.


HTH,

Chuck



Thanks a lot

Hien

tlum...@u.washington.edu wrote:

 On Fri, 18 Dec 2009, Hien Nguyen wrote:

>  Dear Drs Winsemius and Berry,
> 
>  Thanks a lot for your comment and suggestions on running my model. I am 
>  not just new to R but new to CLM as well. :( With your suggestions, I 
>  figure out that I have huge misunderstandings on the model and data 
>  arrangement.
> 
>  After my finals, I have read again related materials on CLM and 
>  rearranged in an appropriate way before running the model in R. This 
>  time, I have a data of more than 250,000 observations (created from more 
>  than 4000 response) and a model of 15 predictors.
> 
>  My question is that how long should it takes for the clogit command to 
>  run because it has been running for more 10 hours on a quad-core 
>  computer and still doesn't show any sign of done or almost done. Is it 
>  OK or my command just does not work.


 If you have a lot of records with case=1 in a stratum, conditional
 logistic regression will be extremely slow.   And unnecessary: maximizing
 the unconditional likelihood is fine when the stratum sizes are large.

 Note that a quad-core computer won't help. Only one core will be used in
 the computations.

  -thomas




>  Thanks a lot for your response
> 
>  Hien
> 
> 
>  Charles C. Berry wrote:

> >  On Fri, 4 Dec 2009, David Winsemius wrote:
> > 
> > > 
> > >  On Dec 4, 2009, at 5:49 PM, Hien Nguyen wrote:
> > > 
> > > >  Dear Dr. Winsemius,
> > > > 
> > > >  Thank you very much for your reply.
> > > > 
> > > >  I have tried many possible combinations (even with the model of 
> > > >  only 2 predictors) but it produces the same message. With more 
> > > >  than 4000 observations, I think 14 predictors might not be too 
> > > >  many.
> > > 
> > >  It is what happens in the factor combinations that concern me. I am 
> > >  guessing that some of those predictors are factors. You really 
> > >  should not ask r-help questions without providing better 
> > >  descriptions of both the outcomes and the predictor variables.
> > > 
> > > > 
> > > >  Although my dependent variable (Pin) is not discrete  (it ranges 
> > > >  from 0 to 1), I do not think it will create problems to the 
> > > >  estimation but I'm not sure
> > > 
> > >  I would think it _would_ cause problems. As I understand it, 
> > >  conditional methods create contingency tables. Why are you using an 
> > >  outcome type that is not consistent with the fundamental regression 
> > >  assumptions of the clogit function?
> > > 
> > >  I do not get that particular error when I munge the infert dataset 
> > >  to have case be a random uniform value, but I do get an error.

> > > >   infert$case <- runif(nrow(infert))
> > > >   clogit(case~spontaneous+induced+strata(stratum),data=infert)

Re: [R] Error when running Conditional Logit Model

[David Winsemius]

On Dec 18, 2009, at 7:39 PM, Hien Nguyen wrote:


Thanks a lot for answering my questions.

I have tried to run the clogit for only 64 observations and 4  
independent variables and the results are solved instantly. However,  
when I run the same command (with only 4 dependent variables) for  
the full data, it keeps running for 50 minutes now. :(


Thomas, what do you mean by "maximizing the unconditional likelihood  
is fine when the stratum sizes are large"? What I put in "strata  
(__)" is actually the possible choices (1-64). Each choices will be  
recored more than 4000 times (which means I have more than 4000  
values of 1, 4000 values of 2 and so on).

Does it sound right?


I'm pretty sure he means glm( formula, family="binomial", ...)  and  
skip the strata specification.


--
David.


Thanks a lot

Hien

tlum...@u.washington.edu wrote:

On Fri, 18 Dec 2009, Hien Nguyen wrote:


Dear Drs Winsemius and Berry,

Thanks a lot for your comment and suggestions on running my model.  
I am not just new to R but new to CLM as well. :( With your  
suggestions, I figure out that I have huge misunderstandings on  
the model and data arrangement.


After my finals, I have read again related materials on CLM and  
rearranged in an appropriate way before running the model in R.  
This time, I have a data of more than 250,000 observations  
(created from more than 4000 response) and a model of 15 predictors.


My question is that how long should it takes for the clogit  
command to run because it has been running for more 10 hours on a  
quad-core computer and still doesn't show any sign of done or  
almost done. Is it OK or my command just does not work.


If you have a lot of records with case=1 in a stratum, conditional  
logistic regression will be extremely slow.   And unnecessary:  
maximizing the unconditional likelihood is fine when the stratum  
sizes are large.


Note that a quad-core computer won't help. Only one core will be  
used in the computations.


-thomas





Thanks a lot for your response

Hien


Charles C. Berry wrote:

On Fri, 4 Dec 2009, David Winsemius wrote:



On Dec 4, 2009, at 5:49 PM, Hien Nguyen wrote:


Dear Dr. Winsemius,

Thank you very much for your reply.

I have tried many possible combinations (even with the model of  
only 2 predictors) but it produces the same message. With more  
than 4000 observations, I think 14 predictors might not be too  
many.


It is what happens in the factor combinations that concern me. I  
am guessing that some of those predictors are factors. You  
really should not ask r-help questions without providing better  
descriptions of both the outcomes and the predictor variables.




Although my dependent variable (Pin) is not discrete  (it  
ranges from 0 to 1), I do not think it will create problems to  
the estimation but I'm not sure


I would think it _would_ cause problems. As I understand it,  
conditional methods create contingency tables. Why are you using  
an outcome type that is not consistent with the fundamental  
regression assumptions of the clogit function?


I do not get that particular error when I munge the infert  
dataset to have case be a random uniform value, but I do get an  
error.

infert$case <- runif(nrow(infert))
clogit(case~spontaneous+induced+strata(stratum),data=infert)

Error in Surv(rep(1, 248L), case) : Invalid status value



David, I think you were on the right track. I get this:

---
clogit(I(case*runif(length(case)))~spontaneous+induced 
+strata(ifelse(stratum>40,NA,stratum)),data=infert)


Error in fitter(X, Y, strats, offset, init, control, weights =  
weights,  :

 NA/NaN/Inf in foreign function call (arg 6)
In addition: Warning messages:
1: In Surv(rep(1, 248L), I(case * runif(length(case :
 Invalid status value, converted to NA
2: In fitter(X, Y, strats, offset, init, control, weights =  
weights,  :

 Ran out of iterations and did not converge





which looks pretty much the same as Hien's error msg

So Hien needs to create a logical status value.

Chuck

p.s.


sessionInfo()

R version 2.10.0 (2009-10-26)
i386-pc-mingw32

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

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

[8] base

other attached packages:
[1] survival_2.35-7

loaded via a namespace (and not attached):
[1] tools_2.10.0





So I certainly would not have proceeded to submit a full  
analysis to clogit if I could not get a test case to run under  
the situation you propose.


--
David



I have checked the collinearity among predictors and they are  
all < 0.5 (which I think is OK). Do you know what else could  
make this errors?


Thanks a lot

Hien Nguyen

David Winsemius wrote:
> > On Dec 4, 2009, at 9:22 AM, Hien Nguyen wrote:
> > > Dear R-helpers,
> > > > I am very new

Re: [R] Error when running Conditional Logit Model

[Hien Nguyen]

Thanks a lot for answering my questions.

I have tried to run the clogit for only 64 observations and 4 
independent variables and the results are solved instantly. However, 
when I run the same command (with only 4 dependent variables) for the 
full data, it keeps running for 50 minutes now. :(


Thomas, what do you mean by "maximizing the unconditional likelihood is 
fine when the stratum sizes are large"? What I put in "strata (__)" is 
actually the possible choices (1-64). Each choices will be recored more 
than 4000 times (which means I have more than 4000 values of 1, 4000 
values of 2 and so on).

Does it sound right?

Thanks a lot

Hien

tlum...@u.washington.edu wrote:

On Fri, 18 Dec 2009, Hien Nguyen wrote:


Dear Drs Winsemius and Berry,

Thanks a lot for your comment and suggestions on running my model. I 
am not just new to R but new to CLM as well. :( With your 
suggestions, I figure out that I have huge misunderstandings on the 
model and data arrangement.


After my finals, I have read again related materials on CLM and 
rearranged in an appropriate way before running the model in R. This 
time, I have a data of more than 250,000 observations (created from 
more than 4000 response) and a model of 15 predictors.


My question is that how long should it takes for the clogit command 
to run because it has been running for more 10 hours on a quad-core 
computer and still doesn't show any sign of done or almost done. Is 
it OK or my command just does not work.


If you have a lot of records with case=1 in a stratum, conditional 
logistic regression will be extremely slow.   And unnecessary: 
maximizing the unconditional likelihood is fine when the stratum sizes 
are large.


Note that a quad-core computer won't help. Only one core will be used 
in the computations.


 -thomas





Thanks a lot for your response

Hien


Charles C. Berry wrote:

On Fri, 4 Dec 2009, David Winsemius wrote:



On Dec 4, 2009, at 5:49 PM, Hien Nguyen wrote:


Dear Dr. Winsemius,

Thank you very much for your reply.

I have tried many possible combinations (even with the model of 
only 2 predictors) but it produces the same message. With more 
than 4000 observations, I think 14 predictors might not be too many.


It is what happens in the factor combinations that concern me. I am 
guessing that some of those predictors are factors. You really 
should not ask r-help questions without providing better 
descriptions of both the outcomes and the predictor variables.




Although my dependent variable (Pin) is not discrete  (it ranges 
from 0 to 1), I do not think it will create problems to the 
estimation but I'm not sure


I would think it _would_ cause problems. As I understand it, 
conditional methods create contingency tables. Why are you using an 
outcome type that is not consistent with the fundamental regression 
assumptions of the clogit function?


I do not get that particular error when I munge the infert dataset 
to have case be a random uniform value, but I do get an error.

 infert$case <- runif(nrow(infert))
 clogit(case~spontaneous+induced+strata(stratum),data=infert)

Error in Surv(rep(1, 248L), case) : Invalid status value



David, I think you were on the right track. I get this:

---
clogit(I(case*runif(length(case)))~spontaneous+induced+strata(ifelse(stratum>40,NA,stratum)),data=infert) 


Error in fitter(X, Y, strats, offset, init, control, weights = 
weights,  :

  NA/NaN/Inf in foreign function call (arg 6)
In addition: Warning messages:
1: In Surv(rep(1, 248L), I(case * runif(length(case :
  Invalid status value, converted to NA
2: In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge





which looks pretty much the same as Hien's error msg

So Hien needs to create a logical status value.

Chuck

p.s.


sessionInfo()

R version 2.10.0 (2009-10-26)
i386-pc-mingw32

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

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

other attached packages:
[1] survival_2.35-7

loaded via a namespace (and not attached):
[1] tools_2.10.0





So I certainly would not have proceeded to submit a full analysis 
to clogit if I could not get a test case to run under the situation 
you propose.


--
David



I have checked the collinearity among predictors and they are all 
< 0.5 (which I think is OK). Do you know what else could make this 
errors?


Thanks a lot

Hien Nguyen

David Winsemius wrote:
> > On Dec 4, 2009, at 9:22 AM, Hien Nguyen wrote:
> > > Dear R-helpers,
> > > > I am very new to R and trying to run the conditional logit 
model using

> > "clogit " command.
> > I have more than 4000 observations in my dataset and try to 
predict the
> > dependent variable from 14 independent varia

Re: [R] how to combine multiple indicator variables in a single factor

[Frank E Harrell Jr]

You might also look at the Hmisc package's score.binary function.
Frank

Daniel Nordlund wrote:

Thanks to Ehud Cohen and William Dunlap for their assistance.  Both solutions 
worked and were helpful.

Dan

Daniel Nordlund
Bothell, WA USA


-Original Message-
From: ehud cohen [mailto:ehudco.l...@gmail.com]
Sent: Friday, December 18, 2009 1:05 PM
To: Daniel Nordlund
Cc: r-help@r-project.org
Subject: Re: [R] how to combine multiple indicator variables in a single
factor

you can try:

df$f<-names(df)[apply(df,1,function(x) which(x==1))]

Ehud

On Fri, Dec 18, 2009 at 10:48 PM, Daniel Nordlund
 wrote:

Say I have a dataframe like this:

df <- data.frame(cbind(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0)))

names(df) <- c('a','b','c')

I would like to create a factor in a new column, where the factor values

are taken from the column names, like this:

df2

 a b c f
1 1 0 0 a
2 0 1 0 b
3 0 0 1 c
4 1 0 0 a

How would I do this?  Thanks,

Dan

Daniel Nordlund
Bothell, WA USA



--
Frank E Harrell Jr   Professor and Chair   School of Medicine
 Department of Biostatistics   Vanderbilt University

__
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] ?OT: Probabilistic Simulation

[Polwart Calum (County Durham and Darlington NHS Foundation Trust)]

Sorry this may well be defined as Off Topic.  I apologize in advance.

I am interested in performing what I think would be a probabilistic sensitivity 
simulation.  I've done some crude ones before in excel but I'm wondering if R 
can help me do it more effectively?

I have a set of theoretical variables for simplicity lets use (what I think) is 
an easier example:  I have a peg and a hole which I want to predict how often 
they would not fit together securely.

The peg should be 1.8mm in diameter.
But there are some variables which might affect if it is actually 1.8mm but I'd 
expect it to be normally distributed around 1.8mm with a known standard 
deviation.  There are also variables such as temperature which would influence 
if the peg was the correct size (there would be a known relationship to 
temperature, and a mean temperature with a standard deviation would be known)

The hole should be 1.8mm diameter as well.  Again there are variables that 
affect if it is, drill size, drilling method, substrate being drilled, 
temperature at time of drilling, temperature at time of fitting peg (would be 
same as above).  I'd be looking to model if the peg would fit in the hole, AND 
if it fitted how well it fitted.

I'd then want to run a simulation of say 500 scenarios, randomly picking 
temperature, hole characteristics etc for each simulation.From there I'd 
get the number of times in 500 samples the peg wouldn't fit.  I could then try 
adjusting a variable - say using a different drilling method and see if that 
'easier' but less reliable drilling method actually would affect the number of 
times the peg didn't fit the hole properly.

So from what I understand this is a MonteCarlo simulation of probability.  And 
this is where I go off topic!  Am I right?  I know its off topic - so what I'd 
like to know is can someone point me to where I can find out?

Then if it is a monte-carlo can someone point me to a good description of how 
to get R to model this?

Ta

Calum



This message may contain confidential information. If yo...{{dropped:21}}

__
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] Package creation - require statement?

[Jeff Breiwick]

All,

I have used package.skeleton, edited the .Rd files and the DESCRIPTON file. 
In one of my package functions I am using RODBC so I have a line of code 
that reads: require(RODBC).

In my DESCRIPTION file I have the following relevant line:
Depends: R (>= 2.9.0), RODBC

However, when I do a check (Rcmd check my_package_name) I get this warning:
* checking for unstated dependencies in R code ... WARNING
'library' or 'require' calls not declared from:
  RODBC
See the information on DESCRIPTION files in the chapter 'Creating R
packages' of the 'Writing R Extensions' manual.

Does anyone know where I have gone astray? Thank you.

Jeff Breiwick

__
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 combine multiple indicator variables in a single factor

[Daniel Nordlund]

Thanks to Ehud Cohen and William Dunlap for their assistance.  Both solutions 
worked and were helpful.

Dan

Daniel Nordlund
Bothell, WA USA

> -Original Message-
> From: ehud cohen [mailto:ehudco.l...@gmail.com]
> Sent: Friday, December 18, 2009 1:05 PM
> To: Daniel Nordlund
> Cc: r-help@r-project.org
> Subject: Re: [R] how to combine multiple indicator variables in a single
> factor
> 
> you can try:
> 
> df$f<-names(df)[apply(df,1,function(x) which(x==1))]
> 
> Ehud
> 
> On Fri, Dec 18, 2009 at 10:48 PM, Daniel Nordlund
>  wrote:
> > Say I have a dataframe like this:
> >
> > df <- data.frame(cbind(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0)))
> >
> > names(df) <- c('a','b','c')
> >
> > I would like to create a factor in a new column, where the factor values
> are taken from the column names, like this:
> >
> >> df2
> >  a b c f
> > 1 1 0 0 a
> > 2 0 1 0 b
> > 3 0 0 1 c
> > 4 1 0 0 a
> >
> > How would I do this?  Thanks,
> >
> > Dan
> >
> > Daniel Nordlund
> > Bothell, WA USA
> >
> > __
> > 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-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 combine multiple indicator variables in a single factor

[William Dunlap]

> -Original Message-
> From: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-project.org] On Behalf Of Daniel Nordlund
> Sent: Friday, December 18, 2009 12:49 PM
> To: r-help@r-project.org
> Subject: [R] how to combine multiple indicator variables in a 
> single factor
> 
> Say I have a dataframe like this:
> 
> df <- data.frame(cbind(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0)))
> 
> names(df) <- c('a','b','c')
> 
> I would like to create a factor in a new column, where the 
> factor values are taken from the column names, like this:
> 
> > df2
>   a b c f
> 1 1 0 0 a
> 2 0 1 0 b
> 3 0 0 1 c
> 4 1 0 0 a
> 
> How would I do this?  Thanks,

One possibility is
   > df$f <- factor(df$a + df$b*2 + df$c*3, labels=names(df)[1:3])
   > df
 a b c f
   1 1 0 0 a
   2 0 1 0 b
   3 0 0 1 c
   4 1 0 0 a


Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com  

> 
> Dan
> 
> Daniel Nordlund
> Bothell, WA USA
> 
> __
> 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-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] Legend for two plots

[Chris Campbell]

On Fri, Dec 18, 2009 at 12:03, Marc Giombetti  wrote:

> Dear R users,
>
> I am new to R and I couldn't figure out how to solve the following
> problem:
>
> I am trying to put a legend below two plots using the code below. The
> legend appears in the second plot,
> but I want the legend to appear below the two plots in the center of
> the total chart. At the moment the
> graphic looks like this: http://i48.tinypic.com/2h2fvhf.jpg
>
> [code]
> layout(matrix(1:2, nrow=1))
> #c(down,left,top,right)
> par(mar=c(4,4,3,1))
> plot(totalExp , totalDiffs,
>main="Years of experience",
>cex.main = 0.9,
>cex.lab=0.8,
>xlab="Years of experience",
>ylab="COCOMO II - expert estimate",
>pch=totalPch,
>col = totalColors
> )
> abline(0,0)
>
> plot(totalSoftwareEXP , totalDiffs,
>main="Years of prof. software experience",
>cex.main = 0.9,
>cex.lab=0.8,
>xlab="Years of experience",
>ylab="COCOMO II - expert estimate",
>pch=totalPch,
>col = totalColors
> )
> abline(0,0)
>
> legend("topleft",  c("Security","Usability") ,col=c('red','blue'),
> pch=c(8,9), cex=0.8,bg='#e0dcdc')
> [/code]
>
> Thanks a lot for your help! I really appreciate it!
>
> Marc
>
> __
> 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.
>

Marc,

There a few issues with your legend call.  Using "topleft" causes the legend
to be drawn in the top left corner _inside_ the current plot.  So you need
to specify x and y coordinates that are outside the plotting region.
 Second, you need to use the xpd=TRUE parameter to have the legend appear in
the margins.  This line will get you close:

legend(par()$usr[1],par()$usr[3], c("Security","Usability")
,col=c('red','blue'),pch=c(8,9), cex=0.8,bg="#e0dcdc",xpd=TRUE,xjust=1)

but there are still problems.  The legend box is not centered between the
plots, it could be clipped by the first plot depending on the size of your
plotting window, etc.  Using locator() may be the easiest way (if your
graphics device supports it):  run this line then click on the plot where
you want the legend to appear:

legend(locator(1), c("Security","Usability")
,col=c('red','blue'),pch=c(8,9), cex=0.8,bg="#e0dcdc",xpd=TRUE,xjust=0)

If your legend is getting clipped, try decreasing the right margin of the
first plot.

Hope that helps.

-Chris

[[alternative HTML version deleted]]

__
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 combine multiple indicator variables in a single factor

[ehud cohen]

you can try:

df$f<-names(df)[apply(df,1,function(x) which(x==1))]

Ehud

On Fri, Dec 18, 2009 at 10:48 PM, Daniel Nordlund
 wrote:
> Say I have a dataframe like this:
>
> df <- data.frame(cbind(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0)))
>
> names(df) <- c('a','b','c')
>
> I would like to create a factor in a new column, where the factor values are 
> taken from the column names, like this:
>
>> df2
>  a b c f
> 1 1 0 0 a
> 2 0 1 0 b
> 3 0 0 1 c
> 4 1 0 0 a
>
> How would I do this?  Thanks,
>
> Dan
>
> Daniel Nordlund
> Bothell, WA USA
>
> __
> 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-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 print to file?

[Charlie Sharpsteen]

On Fri, Dec 18, 2009 at 12:12 PM, Peng Yu  wrote:

> I don't find a function to print a string to file. Would somebody let
> me know what function I should use?
>

I generally use writeLines() to print a string or vector of strings.

-Charlie

[[alternative HTML version deleted]]

__
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 print to file?

[Don MacQueen]

Use the online help:

   help.search('file')

and then see:

  base::sink  Send R Output to a File

help.search('print') will get you there also, but not as directly.

The cat() function is good also, as Gray mentioned.

-Don

At 2:12 PM +1800 12/19/09, Peng Yu wrote:

I don't find a function to print a string to file. Would somebody let
me know what function I should use?

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



--
--
Don MacQueen
Environmental Protection Department
Lawrence Livermore National Laboratory
Livermore, CA, USA
925-423-1062

__
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] how to combine multiple indicator variables in a single factor

[Daniel Nordlund]

Say I have a dataframe like this:

df <- data.frame(cbind(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0)))

names(df) <- c('a','b','c')

I would like to create a factor in a new column, where the factor values are 
taken from the column names, like this:

> df2
  a b c f
1 1 0 0 a
2 0 1 0 b
3 0 0 1 c
4 1 0 0 a

How would I do this?  Thanks,

Dan

Daniel Nordlund
Bothell, WA USA

__
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] lattice, latticeExtra, panel.refline()

[Ivan Gregoretti]

Hello fellow listers,

Can anybody show (or point to) a simple EXAMPLE using

panel.refline()

to draw gridlines on a log-log xyplot?

The idea is to draw something like figure 8.5 in
http://lmdvr.r-forge.r-project.org/figures/figures.html

but with log-scaled grid lines.

Thank you,

Ivan



Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1592
Fax: 1-301-496-9878

__
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 print to file?

[Gray Calhoun]

cat

On Fri, Dec 18, 2009 at 2:12 PM, Peng Yu  wrote:
> I don't find a function to print a string to file. Would somebody let
> me know what function I should use?
>
> __
> 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.
>



-- 
Gray Calhoun

Assistant Professor of Economics
Iowa State University

__
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] forecasting

[Patrick Burns]

DispersionMap wrote:


What about validating forecasting results. I have 5 years of data and have
been forecasting the following three years i.e. 2010 to 2012.

How can i check my forecast.

Why types of tests are out there that can be implemented in R??


isbn 9780553213515

Sorry, couldn't resist,

Pat

__
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] How to print to file?

[Peng Yu]

I don't find a function to print a string to file. Would somebody let
me know what function I should use?

__
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] Principal component

[kayj]

Hi All,


I have a quick question regarding Principal component. I ran principal
components on two data sets cases and controls. I am interested in knowing
if there is a significant difference between the cases and the controls, I
was wondering if there is a test that test if there is a significant
difference between the Principal components of the cases and the controls?

I appreciate your help 

-- 
View this message in context: 
http://n4.nabble.com/Principal-component-tp974949p974949.html
Sent from the R help mailing list archive at Nabble.com.

__
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] Hello all, How can I get corss-validation MSE of SVM in e1071?

[Max Kuhn]

You can get this using the caret package. There are a few package
vignettes that come with the package and a JSS article

  http://www.jstatsoft.org/v28/i05/paper

about the package.

Max

On Fri, Dec 18, 2009 at 12:26 PM, bbslover  wrote:
>
> as known, svm need tune some parameters like  cost,gamma and epsilon to get
> better performance,but one question appear, how can i monitor the
> performance . generally speaking ,we chose the cross-validation MSE in the
> training set, but It seems svm can not return the cross-validation MSE
> value, we just get it from "summary model.svm", if I write a loop, have no
> idear call the cros-validate MSE, and no way to monitor this performance
> ,how can I do?
>
>
> --
> View this message in context: 
> http://n4.nabble.com/Hello-all-How-can-I-get-corss-validation-MSE-of-SVM-in-e1071-tp974942p974942.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> 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.
>



-- 

Max

__
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] Error when running Conditional Logit Model

[David Winsemius]

On Dec 18, 2009, at 2:46 PM, Hien Nguyen wrote:


Dear Drs Winsemius and Berry,

Thanks a lot for your comment and suggestions on running my model. I  
am not just new to R but new to CLM as well. :( With your  
suggestions, I figure out that I have huge misunderstandings on the  
model and data arrangement.


After my finals, I have read again related materials on CLM and  
rearranged in an appropriate way before running the model in R. This  
time, I have a data of more than 250,000 observations (created from  
more than 4000 response) and a model of 15 predictors.


My question is that how long should it takes for the clogit command  
to run because it has been running for more 10 hours on a quad-core  
computer and still doesn't show any sign of done or almost done. Is  
it OK or my command just does not work.


Quad-core would not help speed unless you used a multi-core  
application, which base R is not. Memory and OS is also important to  
specify. Don't have experience with that function but I suspect your  
machine is hung. My ordinary logistic regressions and Cox models with  
4.5 million cases and 40,000 events and 12 df in the X side of the  
formula take minutes. I would terminate the session and then try the  
same code with a 1/100 sample; see what the system times are; and  
scale up to a full situation.


--
David


Thanks a lot for your response

Hien


Charles C. Berry wrote:

On Fri, 4 Dec 2009, David Winsemius wrote:



On Dec 4, 2009, at 5:49 PM, Hien Nguyen wrote:


Dear Dr. Winsemius,

Thank you very much for your reply.

I have tried many possible combinations (even with the model of  
only 2 predictors) but it produces the same message. With more  
than 4000 observations, I think 14 predictors might not be too  
many.


It is what happens in the factor combinations that concern me. I  
am guessing that some of those predictors are factors. You really  
should not ask r-help questions without providing better  
descriptions of both the outcomes and the predictor variables.




Although my dependent variable (Pin) is not discrete  (it ranges  
from 0 to 1), I do not think it will create problems to the  
estimation but I'm not sure


I would think it _would_ cause problems. As I understand it,  
conditional methods create contingency tables. Why are you using  
an outcome type that is not consistent with the fundamental  
regression assumptions of the clogit function?


I do not get that particular error when I munge the infert dataset  
to have case be a random uniform value, but I do get an error.

infert$case <- runif(nrow(infert))
clogit(case~spontaneous+induced+strata(stratum),data=infert)

Error in Surv(rep(1, 248L), case) : Invalid status value



David, I think you were on the right track. I get this:

---
clogit(I(case*runif(length(case)))~spontaneous+induced 
+strata(ifelse(stratum>40,NA,stratum)),data=infert)
Error in fitter(X, Y, strats, offset, init, control, weights =  
weights,  :

 NA/NaN/Inf in foreign function call (arg 6)
In addition: Warning messages:
1: In Surv(rep(1, 248L), I(case * runif(length(case :
 Invalid status value, converted to NA
2: In fitter(X, Y, strats, offset, init, control, weights =  
weights,  :

 Ran out of iterations and did not converge





which looks pretty much the same as Hien's error msg

So Hien needs to create a logical status value.

Chuck

p.s.


sessionInfo()

R version 2.10.0 (2009-10-26)
i386-pc-mingw32

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

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

[8] base

other attached packages:
[1] survival_2.35-7

loaded via a namespace (and not attached):
[1] tools_2.10.0





So I certainly would not have proceeded to submit a full analysis  
to clogit if I could not get a test case to run under the  
situation you propose.


--
David



I have checked the collinearity among predictors and they are all  
< 0.5 (which I think is OK). Do you know what else could make  
this errors?


Thanks a lot

Hien Nguyen

David Winsemius wrote:
> > On Dec 4, 2009, at 9:22 AM, Hien Nguyen wrote:
> > > Dear R-helpers,
> > > > I am very new to R and trying to run the conditional  
logit model using

> > "clogit " command.
> > I have more than 4000 observations in my dataset and try to  
predict the
> > dependent variable from 14 independent variables. My command  
is as > > follows

> > > > clmtest1 <-
> > clogit(Pin~Income+Bus+Pop+Urbpro+Health+Student+Grad+NE+NW+NCC 
+SCC+CH+SE+MRD+strata(IDD),data=clmdata)

> > > > > > However, it produces the following errors:
> > > > Error in fitter(X, Y, strats, offset, init, control,  
weights = weights, > > :

> > NA/NaN/Inf in foreign function call (arg 6)
> > In addition: Warning messages:
> > 1: In Surv(rep(1, 4096L), Pinmig) : Invalid status v

Re: [R] Error when running Conditional Logit Model

[tlumley]

On Fri, 18 Dec 2009, Hien Nguyen wrote:


Dear Drs Winsemius and Berry,

Thanks a lot for your comment and suggestions on running my model. I am not 
just new to R but new to CLM as well. :( With your suggestions, I figure out 
that I have huge misunderstandings on the model and data arrangement.


After my finals, I have read again related materials on CLM and rearranged in 
an appropriate way before running the model in R. This time, I have a data of 
more than 250,000 observations (created from more than 4000 response) and a 
model of 15 predictors.


My question is that how long should it takes for the clogit command to run 
because it has been running for more 10 hours on a quad-core computer and 
still doesn't show any sign of done or almost done. Is it OK or my command 
just does not work.


If you have a lot of records with case=1 in a stratum, conditional logistic 
regression will be extremely slow.   And unnecessary: maximizing the 
unconditional likelihood is fine when the stratum sizes are large.

Note that a quad-core computer won't help. Only one core will be used in the 
computations.

 -thomas





Thanks a lot for your response

Hien


Charles C. Berry wrote:

On Fri, 4 Dec 2009, David Winsemius wrote:



On Dec 4, 2009, at 5:49 PM, Hien Nguyen wrote:


Dear Dr. Winsemius,

Thank you very much for your reply.

I have tried many possible combinations (even with the model of only 2 
predictors) but it produces the same message. With more than 4000 
observations, I think 14 predictors might not be too many.


It is what happens in the factor combinations that concern me. I am 
guessing that some of those predictors are factors. You really should not 
ask r-help questions without providing better descriptions of both the 
outcomes and the predictor variables.




Although my dependent variable (Pin) is not discrete  (it ranges from 0 
to 1), I do not think it will create problems to the estimation but I'm 
not sure


I would think it _would_ cause problems. As I understand it, conditional 
methods create contingency tables. Why are you using an outcome type that 
is not consistent with the fundamental regression assumptions of the 
clogit function?


I do not get that particular error when I munge the infert dataset to have 
case be a random uniform value, but I do get an error.

 infert$case <- runif(nrow(infert))
 clogit(case~spontaneous+induced+strata(stratum),data=infert)

Error in Surv(rep(1, 248L), case) : Invalid status value



David, I think you were on the right track. I get this:

---
clogit(I(case*runif(length(case)))~spontaneous+induced+strata(ifelse(stratum>40,NA,stratum)),data=infert) 

Error in fitter(X, Y, strats, offset, init, control, weights = weights,  :
  NA/NaN/Inf in foreign function call (arg 6)
In addition: Warning messages:
1: In Surv(rep(1, 248L), I(case * runif(length(case :
  Invalid status value, converted to NA
2: In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge





which looks pretty much the same as Hien's error msg

So Hien needs to create a logical status value.

Chuck

p.s.


sessionInfo()

R version 2.10.0 (2009-10-26)
i386-pc-mingw32

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

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

other attached packages:
[1] survival_2.35-7

loaded via a namespace (and not attached):
[1] tools_2.10.0





So I certainly would not have proceeded to submit a full analysis to 
clogit if I could not get a test case to run under the situation you 
propose.


--
David



I have checked the collinearity among predictors and they are all < 0.5 
(which I think is OK). Do you know what else could make this errors?


Thanks a lot

Hien Nguyen

David Winsemius wrote:
> > On Dec 4, 2009, at 9:22 AM, Hien Nguyen wrote:
> > > Dear R-helpers,
> > > > I am very new to R and trying to run the conditional logit model 
using

> > "clogit " command.
> > I have more than 4000 observations in my dataset and try to predict 
the
> > dependent variable from 14 independent variables. My command is as > 
> follows

> > > > clmtest1 <-
> > 
clogit(Pin~Income+Bus+Pop+Urbpro+Health+Student+Grad+NE+NW+NCC+SCC+CH+SE+MRD+strata(IDD),data=clmdata) 
> > > > > > However, it produces the following errors:
> > > > Error in fitter(X, Y, strats, offset, init, control, weights = 
weights, > > :

> > NA/NaN/Inf in foreign function call (arg 6)
> > In addition: Warning messages:
> > 1: In Surv(rep(1, 4096L), Pinmig) : Invalid status value, converted 
to > > NA
> > 2: In fitter(X, Y, strats, offset, init, control, weights = weights, 
:

> > Ran out of iterations and did not converge
> > > > I search the error message from R forums but it does not say 
anything

> > for Conditio

Re: [R] Error when running Conditional Logit Model

[Hien Nguyen]

Dear Drs Winsemius and Berry,

Thanks a lot for your comment and suggestions on running my model. I am 
not just new to R but new to CLM as well. :( With your suggestions, I 
figure out that I have huge misunderstandings on the model and data 
arrangement.


After my finals, I have read again related materials on CLM and 
rearranged in an appropriate way before running the model in R. This 
time, I have a data of more than 250,000 observations (created from more 
than 4000 response) and a model of 15 predictors.


My question is that how long should it takes for the clogit command to 
run because it has been running for more 10 hours on a quad-core 
computer and still doesn't show any sign of done or almost done. Is it 
OK or my command just does not work.


Thanks a lot for your response

Hien


Charles C. Berry wrote:

On Fri, 4 Dec 2009, David Winsemius wrote:



On Dec 4, 2009, at 5:49 PM, Hien Nguyen wrote:


Dear Dr. Winsemius,

Thank you very much for your reply.

I have tried many possible combinations (even with the model of only 
2 predictors) but it produces the same message. With more than 4000 
observations, I think 14 predictors might not be too many.


It is what happens in the factor combinations that concern me. I am 
guessing that some of those predictors are factors. You really should 
not ask r-help questions without providing better descriptions of 
both the outcomes and the predictor variables.




Although my dependent variable (Pin) is not discrete  (it ranges 
from 0 to 1), I do not think it will create problems to the 
estimation but I'm not sure


I would think it _would_ cause problems. As I understand it, 
conditional methods create contingency tables. Why are you using an 
outcome type that is not consistent with the fundamental regression 
assumptions of the clogit function?


I do not get that particular error when I munge the infert dataset to 
have case be a random uniform value, but I do get an error.

 infert$case <- runif(nrow(infert))
 clogit(case~spontaneous+induced+strata(stratum),data=infert)

Error in Surv(rep(1, 248L), case) : Invalid status value



David, I think you were on the right track. I get this:

---
clogit(I(case*runif(length(case)))~spontaneous+induced+strata(ifelse(stratum>40,NA,stratum)),data=infert) 

Error in fitter(X, Y, strats, offset, init, control, weights = 
weights,  :

  NA/NaN/Inf in foreign function call (arg 6)
In addition: Warning messages:
1: In Surv(rep(1, 248L), I(case * runif(length(case :
  Invalid status value, converted to NA
2: In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge





which looks pretty much the same as Hien's error msg

So Hien needs to create a logical status value.

Chuck

p.s.


sessionInfo()

R version 2.10.0 (2009-10-26)
i386-pc-mingw32

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

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

other attached packages:
[1] survival_2.35-7

loaded via a namespace (and not attached):
[1] tools_2.10.0





So I certainly would not have proceeded to submit a full analysis to 
clogit if I could not get a test case to run under the situation you 
propose.


--
David



I have checked the collinearity among predictors and they are all < 
0.5 (which I think is OK). Do you know what else could make this 
errors?


Thanks a lot

Hien Nguyen

David Winsemius wrote:
> > On Dec 4, 2009, at 9:22 AM, Hien Nguyen wrote:
> > > Dear R-helpers,
> > > > I am very new to R and trying to run the conditional logit 
model using

> > "clogit " command.
> > I have more than 4000 observations in my dataset and try to 
predict the
> > dependent variable from 14 independent variables. My command is 
as > > follows

> > > > clmtest1 <-
> > 
clogit(Pin~Income+Bus+Pop+Urbpro+Health+Student+Grad+NE+NW+NCC+SCC+CH+SE+MRD+strata(IDD),data=clmdata) 


> > > > > > However, it produces the following errors:
> > > > Error in fitter(X, Y, strats, offset, init, control, weights 
= weights, > > :

> > NA/NaN/Inf in foreign function call (arg 6)
> > In addition: Warning messages:
> > 1: In Surv(rep(1, 4096L), Pinmig) : Invalid status value, 
converted to > > NA
> > 2: In fitter(X, Y, strats, offset, init, control, weights = 
weights,  :

> > Ran out of iterations and did not converge
> > > > I search the error message from R forums but it does not say 
anything

> > for Conditional Logit Model.
> > With that many predictors in a small dataset, you may have 
created matrix > singularities. Perhaps you created a stratum where 
all of the subjects > experience the event and others where none did 
so. The coefficients might > be driven to infinities. Try 
simplifying the model.
> > > > > > Please check for me what it says and what shoul

Re: [R] forecasting

[DispersionMap]

What about validating forecasting results. I have 5 years of data and have
been forecasting the following three years i.e. 2010 to 2012.

How can i check my forecast.

Why types of tests are out there that can be implemented in R??
-- 
View this message in context: 
http://n4.nabble.com/forecasting-tp974326p975010.html
Sent from the R help mailing list archive at Nabble.com.

__
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] Legend for two plots

[Bert Gunter]

Marc:

As you have not received a (public) reply yet, let me try. You may need to
take what I say with a spoonful of salt, though.

Basically, your request makes no sense. The nature of a legend is that it
provides information about a specific plot (e.g. the group labels associated
with point colors). So the legend code apparently does not contemplate
association and positioning with respect to several similar plots, which is
what you'd like to do, although maybe there is something I've overlooked and
this can be straightforwardly done. In which case, please post the solution
to the list.

So I think what you need to do is draw the legend manually using, e.g.,
?text or ?mtext, ?segments, etc.(You may particularly want to note the
"outer" argument in mtext and par).  One way to do this might be to use
layout to set up your regions, including the legend area, and then "draw"
the legend as a plot where you want it. 

This does seem rather laborious, so I would hope that there is a more
elegant solution. Alternatively, you may wish to examine the ggplot or
lattice plotting systems to see whether they can do this simply (lattice
certainly can using it's "key" argument; but I'm not sure whether your plots
fit easily into the lattice framework). Of course, there's an additional
learning curve there to contend with.

HTH.

Cheers,

Bert Gunter
Genentech Nonclinical Biostatistics
467-7374


-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Marc Giombetti
Sent: Friday, December 18, 2009 9:04 AM
To: r-help@r-project.org
Subject: [R] Legend for two plots

Dear R users,

I am new to R and I couldn't figure out how to solve the following
problem:

I am trying to put a legend below two plots using the code below. The
legend appears in the second plot,
but I want the legend to appear below the two plots in the center of
the total chart. At the moment the
graphic looks like this: http://i48.tinypic.com/2h2fvhf.jpg

[code]
layout(matrix(1:2, nrow=1))
#c(down,left,top,right)
par(mar=c(4,4,3,1))
plot(totalExp , totalDiffs,
main="Years of experience",
cex.main = 0.9,
cex.lab=0.8,
xlab="Years of experience",
ylab="COCOMO II - expert estimate",
pch=totalPch,
col = totalColors
)
abline(0,0)

plot(totalSoftwareEXP , totalDiffs,
main="Years of prof. software experience",
cex.main = 0.9,
cex.lab=0.8,
xlab="Years of experience",
ylab="COCOMO II - expert estimate",
pch=totalPch,
col = totalColors
)
abline(0,0)

legend("topleft",  c("Security","Usability") ,col=c('red','blue'),
pch=c(8,9), cex=0.8,bg='#e0dcdc')
[/code]

Thanks a lot for your help! I really appreciate it!

Marc

__
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-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] NLS-Weibull-ERROR

[Gabor Grothendieck]

This only uses 2 parameters instead of 3 and gives quite a good fit:

> DF <- structure(list(conc = c(0.077, 0.328, 0.882, 1.195, 1.884, 3.577,
+ 6.549, 13, 33.69, 52.22, 90.14, 166.05, 233.62, 346.89), vel =
1:14), .Names = c("conc",
+ "vel"), class = "data.frame", row.names = c(NA, -14L))
>
> DF.nls <- nls(conc ~ a * exp(b/vel), DF, start = c(a=1, b=1))
> DF.nls
Nonlinear regression model
  model:  conc ~ a * exp(b/vel)
   data:  DF
   ab
35773.97   -65.02
 residual sum-of-squares: 249.6

Number of iterations to convergence: 12
Achieved convergence tolerance: 1.410e-06
> plot(conc ~ vel, DF)
> lines(DF$vel, fitted(DF.nls))


On Fri, Dec 18, 2009 at 7:53 AM, ruchita gupta  wrote:
> Hello
>
> I was trying to estimate the weibull model using nls after putting OLS
> values as the initial inputs to NLS.
> I tried multiple times but still i m getting the same error of Error in
> nlsModel(formula, mf, start, wts) :
>  singular gradient matrix at initial parameter estimates.
>
> The Program is as below
>
>> vel <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
>> df <- data.frame(conc, vel)
>> df
>      conc vel
> 1    0.077   1
> 2    0.328   2
> 3    0.882   3
> 4    1.195   4
> 5    1.884   5
> 6    3.577   6
> 7    6.549   7
> 8   13.000   8
> 9   33.690   9
> 10  52.220  10
> 11  90.140  11
> 12 166.050  12
> 13 233.620  13
> 14 346.890  14
>> plot(df$vel, df$conc)
>> para0.st <- c(K=450,
> +       alpha=0.054,beta=3.398 )
>>  fit0 <- nls(
> +      conc~ K-(K*exp(-(vel/alpha)^beta)), df, start= para0.st,trace=T)
> Error in nlsModel(formula, mf, start, wts) :
>  singular gradient matrix at initial parameter estimates
>
>
> I will be highly thankful if some one can please let me know where is the
> mistake as i m unable to trace it.
>
> Thanks
> Ruchita
>
> On Wed, Dec 16, 2009 at 3:18 PM, ruchita gupta wrote:
>
>> Hello
>>
>> After performing NLS estimates for some sigmoid model(like logistic growth
>> model  and Gompertz growth models), how can we get the RMSE(root mean square
>> error) and MAPE(mean absolute percentage error) in  R statistical tool for
>> comparison between two models
>>
>> Thanks
>> Ruchita
>>
>
>        [[alternative HTML version deleted]]
>
> __
> 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-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] What is the fastest way to see what are in an RData file?

[Peng Yu]

On Fri, Dec 18, 2009 at 3:33 AM, Gustaf Rydevik
 wrote:
> On Thu, Dec 17, 2009 at 4:33 PM, Peng Yu  wrote:
>> On Thu, Dec 17, 2009 at 5:33 AM, Gustaf Rydevik
>>  wrote:
>>> On Wed, Dec 16, 2009 at 10:13 PM, Peng Yu  wrote:

 Currently, I load the RData file then ls() and str(). But loading the file
 takes too long if the file is big. Most of the time, I only interested what
 the variables are in the the file and the attributes of the variables (like
 if it is a data.frame, matrix, what are the colnames/rownames, etc.)

 I'm wondering if there is any facility in R to help me avoid loading the
 whole file.
>>>
>>>
>>> I thought this was interesting as well, so i did a bit of searching
>>> through the R-help list archives and found this answer by Simon
>>> Urbanek:
>>> https://stat.ethz.ch/pipermail/r-devel/2007-August/046724.html
>>> The link to a c-routine that does what you want still works, but for
>>> future reference I'm pasting the code below.
>>
>> It doesn't work for the RData file that I saved by "save(list='test',
>> file='test.RData')".
>>
>> $ rdcopy test.RData
>> Format version 3ec, R version = 23813.88.84, release = f9db1dba
>> Sorry, this tool supported RXDR version 2 format only
>>
>
>
> What happens if you remove the version check?
> I.e. this one:
>  if (ver != 2) {
>    XdrInTerm(d);
>    error(_("Sorry, this tool supported RXDR version 2 format only\n"));
>  }

I commented the above lines out. But it doesn't work.

$ rdcopy some.RData
Format version 314, R version = -25755.184.212, release = e6164617
ReadItem: unknown type 93, perhaps written by later version of Rp

> From what I can read on the hel page for ?save, there hasn't been a
> change in the file format since 1.4.0

__
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] Hello all, How can I get corss-validation MSE of SVM in e1071?

[bbslover]

as known, svm need tune some parameters like  cost,gamma and epsilon to get
better performance,but one question appear, how can i monitor the
performance . generally speaking ,we chose the cross-validation MSE in the
training set, but It seems svm can not return the cross-validation MSE
value, we just get it from "summary model.svm", if I write a loop, have no
idear call the cros-validate MSE, and no way to monitor this performance
,how can I do?  


-- 
View this message in context: 
http://n4.nabble.com/Hello-all-How-can-I-get-corss-validation-MSE-of-SVM-in-e1071-tp974942p974942.html
Sent from the R help mailing list archive at Nabble.com.

__
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] Which hist cell each value falls in?

[Doug Hill]

Thanks, guys! Doug


David Winsemius wrote:
> 
> 
> On Dec 18, 2009, at 2:57 AM, Jim Lemon wrote:
> 
>> On 12/18/2009 06:35 AM, Doug Hill wrote:
>>> Hi, all. I'm using hist() to obtain a vector of break values in an  
>>> interval.
>>> I then want to be able to identify which cell any value from  
>>> another vector
>>> falls in.
>>>
>>> E.g. applying
>>>
>>>
 breaks

>>>  [1] -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0
>>>
>>> to
>>>
>>>
 x

>>>  [1] -3.74519666 -0.38183630 -1.22884247 -0.20971824 -0.30533939  
>>> -0.36271207
>>> [7] -2.27513499 -2.23688653 -1.98827155 -1.48666274 -1.26155084  
>>> -0.15234555
>>> [13] -0.09497287  0.34488440
>>>
>>> would give
>>>
>>>
 xcells

>>>  [1] 1  8  6  8  8  8  4  4  5  6  6  8  8  9
>>>
>>> where:
>>>   x<= breaks[1] ->  cell 1
>>> breaks[1]<  x<= breaks[2] ->  cell 2
>>> breaks[2]<  x<= breaks[3] ->  cell 3, etc.
> 
>  > breaks <- scan()
> 1: -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0
> 13:
> Read 12 items
>  > xcells <- scan()
> 1: -3.74519666 -0.38183630 -1.22884247 -0.20971824 -0.30533939  
> -0.36271207
> 7: -2.27513499 -2.23688653 -1.98827155 -1.48666274 -1.26155084  
> -0.15234555
> 13:  -0.09497287  0.34488440
> 15:
> Read 14 items
> 
>  > breaks <- c(-Inf,breaks,Inf)
>  > findInterval(xcells, breaks)
>   [1] 1 8 6 8 8 8 4 4 5 6 6 8 8 9
> 
> -- 
> David
>>>
>>>
>> Hi Doug,
>> The function below does more or less what you want, except that the  
>> bins are numbered from zero (meaning that a value is below the range  
>> of bins as x[1] is) to length(breaks) (meaning that a value is above  
>> the range of bins).
>>
>> whichBin<-function(x,breaks,right=TRUE) {
>> lenx<-length(x)
>> # any leftovers must be out of range
>> wb<-rep(lenx,lenx)
>> if(right) wb[x<=breaks[1]]<-0
>> else wb[x> for(bin in 1:(length(breaks)-1)) {
>>  if(right) wb[x>breaks[bin] & x<=breaks[bin+1]]<-bin
>>  else wb[x>=breaks[bin] & x> }
>> return(wb)
>> }
>>
>> Jim
>>
>> __
>> 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.
> 
> David Winsemius, MD
> Heritage Laboratories
> West Hartford, CT
> 
> __
> 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.
> 
> 

-- 
View this message in context: 
http://n4.nabble.com/Which-hist-cell-each-value-falls-in-tp974316p974865.html
Sent from the R help mailing list archive at Nabble.com.

__
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] Testing equality of regression model on multiple groups

[Clara Yuan]

Hi Daniel,

Thanks for your thorough response. You are indeed correct that I was looking 
for (a), and the Chow test fits the bill exactly. However, I believe my method 
is equivalent to the Chow test. Rather than summing the errors across 
regressions on different datasets, I formulate a single specification with 
errors equivalent to that sum. (Except there would be one difference: the 
reduced and full models include dummies for the grouping variable, and the Chow 
test doesn't.) What interfered with the method was the fact that my factor was 
ordered and therefore was not treated as a dummy variables by R in my model 
specification. Thanks to David who helped clarify R's behaviour on this point.

Once I unordered my factor, I found that the coefficient estimates were the 
same.
> data.ex$t = factor(data.ex$t, ordered = F)
> coef(lm.together)
 (Intercept)   t2   t3   t4   x1   x2
 2.691272263 -0.975915716 -0.811480858 -0.763039039  0.107520721  0.054694784
   t2:x1t3:x1t4:x1t2:x2t3:x2t4:x2
-0.060271900 -0.063503572 -0.087683111 -0.015070147  0.005671774 -0.024606167
   x1:x2 t2:x1:x2 t3:x1:x2 t4:x1:x2
 0.002180749 -0.001815363 -0.002329575 -0.002921691

You can see that these are equivalent to the estimates from lm.separate. 
Demonstrating with the intercept:
> sapply(lm.separate, coef)[1,]
   1234
2.691272 1.715357 1.879791 1.928233
> coef(lm.together)[1] + c(0, coef(lm.together)[2:4])
   t2   t3   t4
2.691272 1.715357 1.879791 1.928233

For reference, the method I'm using is described here:
www.stat.wisc.edu/~mchung/teaching/stat324/324.24.pdf

Thanks very much for everyone's help,
Clara

From: Daniel Malter [dan...@umd.edu]
Sent: Friday, December 18, 2009 4:22 AM
To: Clara Yuan; 'r-h...@lists.r-project.org'
Subject: RE: [R] Testing equality of regression model on multiple groups

Hi, your question is unclear. It is not clear whether you want to compare
a.) whether two subsets of data have the same coefficients for an identical
model or b.) whether a couple of coefficients are different once you include
additional regressors. The reason for this confusion is that your
lm.separate and lm.together are not the same models (i.e., they do not
include the same regressors). The former is (y~x1*x2); the latter is
(y~t*x1*x2), which is obviously different.

If a.) is your goal, you want to run a Chow test.

For that, you run the regression (y~x1*x2) on the entire dataset and on the
two subsets of the data separately. If the model on the entire data produces
much larger errors (in sum) than the two regressions on the subsets of data
(in sum), then this is evidence that the data-generating processes for the
two subsets differ significantly. In other words, the coefficients for the
subsets will jointly be significantly different from each other. This is, as
you said, an F-Test. Look here for further details:

http://en.wikipedia.org/wiki/Chow_test

It just beckons me that you probably want to compare multiple groups. If you
want to run (y~x1*x2) interacted with the group indicator, i.e.
(y~group*x1*x2), you should code your group indicator as "treatment
contrasts."

The t.L, t.Q, and t.C indicate a.) that t is a factor variable and b.) that
this factor variable is coded as orthogonal polynomial contrasts. This tests
for linear, quadractic, cubic, etc. influence of t on the dependent
variable. This makes sense if t has constant differences between the factor
levels (like 10 cm, 15 cm, 20 cm, 25 cm). Otherwise, you probably want
dummy-variable coding, which is called "treatment contrasts" in R lingo.

HTH,
Daniel

-
cuncta stricte discussurus
-
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Clara Yuan
Sent: Thursday, December 17, 2009 5:15 PM
To: r-h...@lists.r-project.org
Subject: [R] Testing equality of regression model on multiple groups

Hello,

I'm trying to test for the joint equality of coefficients of the same model
across different subsets of the data (ie, is it necessary to estimate the
same model on these different populations, or can I just estimate the model
once on the whole dataset?).

My plan is to use the F-test on the reduced model and the full model. By
full model, I mean a specification that mimics my regressions on separate
subsets of data, but I have found that the full model's coefficient
estimates don't correspond to my original model's estimates. I was under the
impression that they would be identical.

Original model:
> lm.separate = by(data.ex, data.ex$t, function(x) lm(y ~ x1 * x2, data =
x))

Full model:
> lm.together = lm(y ~ t * x1 * x2, data = data.ex)

The data are grouped by t.

When I examine the coefficients, I find that they are roughly in the same
ballpark, but not nearly identical

Re: [R] Significant performance difference between split of adata.frame and split of vectors

[Matthew Dowle]

Thanks for suggesting data.table. It does have advantages in this example 
but it has to be used in a particular way.

What does Peng actually want to achieve?  I'll guess (but its only a guess) 
that he doesn't actually need to hold the entire table in memory in a split 
up format before doing something with each subset.  He only needs each 
subset at a time.  Most of the time I guess most users are similar in this 
regard.

data.table has a mechanism to do this, and do it in a compact syntax that is 
quicker to program. When you use the 'by' argument,  it evaluates the j 
expression within the first subset by the 'by' and then doesn't need that 
subset anymore. It then moves on to the next subset. So it operates in a 
smaller memory footprint than i) creating an entire copy of the data in a 
split up format, and then ii) lapply'ing through that new list doing 
something on the subsets.

By guessing at something realistic we might want to actually do on each 
subset,  here is the example again :

> m = 3
> n = 6
> k = 3000
> x = replicate(n,rnorm(m))
> f = sample(1:k,size=m,replace=TRUE)
> dt = data.table(x,f)
> dt[,sum(V1),by="f"]# this is the proper mechanism to split a 
> data.table and operate on each subset. column names can be used as 
> variables directly.
 f  V1
 11.82720825112107
 24.22721189592209
 3  -0.409096014477913
...[ snip ] ...
> system.time(dt[,sum(V1),by="f"])# same again but timing it this time
user  system elapsed
   0.130.000.12
>
>  system.time(split(as.data.frame(x),f))
   user  system elapsed
   1.550.001.55
>

So just splitting a data.frame is about 10 times slower than splitting the 
data.table (in the proper way), and even then the data.frame split still 
needs something (more program code) to loop through afterwards and do 
something useful on the split up copy.

Its not just speed,  but data.table may work when data.frame may fail. An 
example could be constructed where x takes 55% of available ram. I'd expect 
that splitting a data.frame should fail with an out of memory error (as it 
needs another 55% for the full copy).  A data.table 'by' in contrast should 
work fine since it only requires the memory for the largest subset at any 
one time (other than in the esoteric case of there being only one subset).

Having said the above I still consider 'by' in data.table to be slow,  but 
just relative to how fast it could be. See feature request #195 in the 
package's R-forge project. There are some other more complicated 
circumstances when the 'by' argument can result in slow performance 
(possibly relative to data.frame or matrix methods) but this example doesn't 
seem to be in that category,  at least with what we've seen so far.  The 
problem in this thread so far looks to be due to split() itself.

Please note there is no split.data.table method and the default method does 
not appear efficient. methods(split) shows there is no split method for 
data.table but data.frame has its own special split method. The following 
method has been added as a feature request in the data.table project :
split.data.table = function(...){stop("Use 'by' argument instead of 
split() on data.table. See ?'[.data.table'")}

Hope this helps,
Matthew


"David Winsemius"  wrote in message 
news:ce3e3aff-f2b5-4c80-9d00-00954130a...@comcast.net...
>
> On Dec 9, 2009, at 2:59 PM, Peng Yu wrote:
>
>> On Tue, Dec 8, 2009 at 10:37 PM, David Winsemius > > wrote:
>>>
>>> On Dec 8, 2009, at 11:28 PM, Peng Yu wrote:
>>>
 I have the following code, which tests the split on a data.frame and
 the split on each column (as vector) separately. The runtimes are of
 10 time difference. When m and k increase, the difference become  even
 bigger.

 I'm wondering why the performance on data.frame is so bad. Is it a  bug
 in R? Can it be improved?
>>>
>>> You might want to look at the data.table package. The author calinms
>>> significant speed improvements over dta.frames
>>
>> 'data.table' doesn't seem to help. You can try the other set of m,n,k.
>> In both case, using as.data.frame is faster than using as.data.table.
>>
>> Please let me know if I understand what you meant.
>
> I was only suggesting that you look at it because it appeared in other 
> situation to have efficiency advantages. As it turned out, that  structure 
> offered no advantage, when I tested it.
>
> --
> David.
>
>
>>
>>> m=10
>>> n=6
>>> k=3
>>>
>>> #m=30
>>> #n=6
>>> #k=3
>>>
>>> set.seed(0)
>>> x=replicate(n,rnorm(m))
>>> f=sample(1:k, size=m, replace=T)
>>>
>>> library(data.table)
>> Loading required package: ref
>> dim(refdata) and dimnames(refdata) no longer allow parameter ref=TRUE,
>> use dim(derefdata(refdata)), dimnames(derefdata(refdata)) instead
>>> system.time(split(as.data.frame(x),f))
>>   user  system elapsed
>>  0.000   0.000   0.003
>>> system.time(split(as.data.table(x),f))
>>   user  system elapsed
>>  0.010   0.000

Re: [R] error when using multcomp and lm

[Achim Zeileis]

On Fri, 18 Dec 2009, chrischizinski wrote:

Just in case anyone ever has similar issues, it appears that 'glht' does 
not work without an intercept.


That is not precise enough to be useful. I've had a look again and appears 
to be the case that mcp() (or rather the internal function mcp2matrix) 
does not work for models where there are multiple factors but no 
intercept.


As the posting guide asks us to do, I have (1) provided a _reproducible_ 
example below, and (2) included the maintainer of the package in the 
posting. Torsten, it would be great if you could have a look.


Best,
Z

## package
library("multcomp")

## various models with and without intercept
m1a <- lm(breaks ~ tension, data = warpbreaks)
m1b <- lm(breaks ~ 0 + tension, data = warpbreaks)
m2a <- lm(breaks ~ wool + tension, data = warpbreaks)
m2b <- lm(breaks ~ 0 + wool + tension, data = warpbreaks)

## these two are equivalent: one factor with/without intercept
summary(glht(m1a, linfct = mcp(tension = "Tukey")))
summary(glht(m1b, linfct = mcp(tension = "Tukey")))

## these two should be equivalent: two factors with/without intercept
## but the latter fails
summary(glht(m2a, linfct = mcp(tension = "Tukey")))
summary(glht(m2b, linfct = mcp(tension = "Tukey")))

__
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] Numerical Integration

[William Dunlap]

> -Original Message-
> From: Julio Rojas [mailto:jcredbe...@ymail.com] 
> Sent: Friday, December 18, 2009 9:06 AM
> To: William Dunlap; r-help@r-project.org
> Subject: RE: [R] Numerical Integration
> 
> Thanks a lot William. I'm sorry about the syntax problem. I 
> was working at the same time with R code and a document, so I 
> copied from the wrong document.
> 
> I saw that I had a problem with the conditions of the 
> integrand. A ">" instead of a ">=" on the second one. Again, thanks.
> 
> One last question: Is there a way to use "approx" as the integrand?

Try
   integrand2 <- approxfun(x=c(-1e20,fx,Inf),y=c(0,0,1,1,0,0))
instead of your
   Vector(integrand)

The -1e20 should logically be -Inf, but approxfun produces
bad code in that case so the resulting function produces NaN
instead of 0 for x 
> Best regards.
> 
> Julio
> 
> 
> --- El vie 18-dic-09, William Dunlap  escribió:
> 
> > De: William Dunlap 
> > Asunto: RE: [R] Numerical Integration
> > A: "Julio Rojas" 
> > Fecha: viernes, 18 diciembre, 2009, 4:03 pm
> > The immediate problem is that
> > integrand(0.5) is NULL
> > because none of the if clauses is TRUE if x==fx[2].
> > 
> > If you restructured the code such that it had no return
> > statements and instead assigned things to a variable
> > and returned that variable at the end it would be more
> > apparent what happened: you would get a 'variable not
> > found' error when returning.
> >    i <- function(x) {
> >       if (x >       else if (x>fx[1] &&
> > x<=fx[2]) retval <- 1
> >       else if (x>fx[2]) retval <- 0
> >       retval
> >    }
> >    > i(.4)
> >    [1] 1
> >    > i(.3)
> >    Error in i(0.3) : object 'retval' not
> > found
> > 
> > When you convert your code for production use, look
> > at the ifelse and approx functions for doing this
> > sort of thing.  They are much faster than
> > Vector(yourIntegrand).
> > 
> > Also, when writing for help from R-help, please use
> > R syntax when defining variables, like
> >    fx<-c(0.3,0.5,0.5,0.6)
> > and not some other language, like
> >    fx<-[0.3,0.5,0.5,0.6]
> > If we can copy and paste your example into R instead
> > of translating the syntax by hand we will
> > be much more inclined to try to fix things up.
> > 
> > Bill Dunlap
> > Spotfire, TIBCO Software
> > wdunlap tibco.com  
> > 
> > > -Original Message-
> > > From: r-help-boun...@r-project.org
> > 
> > > [mailto:r-help-boun...@r-project.org]
> > On Behalf Of Julio Rojas
> > > Sent: Friday, December 18, 2009 7:02 AM
> > > To: r-help@r-project.org
> > > Subject: [R] Numerical Integration
> > > 
> > > Dear @ll. I have to calculate numerical integrals for
> > 
> > > triangular and trapezoidal figures. I know you can
> > calculate 
> > > the exactly, but I want to do it this way to learn how
> > to 
> > > proceed with more complicated shapes. The code I'm
> > using is 
> > > the following:
> > > 
> > > integrand<-function(x) {
> > >    print(x)
> > >    if(x > >    if(x>=fx[1] && x > return((x-fx[1])/(fx[2]-fx[1]))
> > >    if(x>fx[2] && x<=fx[3])
> > return(1)
> > >    if(x>fx[3] && x<=fx[4])
> > return((x-fx[4])/(fx[3]-fx[4]))
> > >    if(x>fx[4]) return(0)
> > > 
> > >  }
> > > 
> > > fx<-data[i,j,]
> > > reltol<-1e-07
> > >
> > integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions
> > > =200)$value
> > > 
> > > It works for most cases, but then, I tried for the
> > triangle 
> > > fx<-[0.3,0.5,0.5,0.6] and the following error was
> > presented:
> > > 
> > > >
> > integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions=200)
> > > [1] 0.5
> > > [1] 0.01304674
> > > [1] 0.9869533
> > > [1] 0.06746832
> > > [1] 0.9325317
> > > [1] 0.1602952
> > > [1] 0.8397048
> > > [1] 0.2833023
> > > [1] 0.7166977
> > > [1] 0.4255628
> > > [1] 0.5744372
> > > [1] 0.002171418
> > > [1] 0.9978286
> > > [1] 0.03492125
> > > [1] 0.9650787
> > > [1] 0.1095911
> > > [1] 0.8904089
> > > [1] 0.2186214
> > > [1] 0.7813786
> > > [1] 0.3528036
> > > [1] 0.6471964
> > > Error in integrate(Vectorize(integrand), 0, 1, rel.tol
> > = reltol,
> > > subdivisions = 200) :
> > >  evaluation of function gave a result of wrong
> > type
> > > 
> > > Does anybody know what happened? Thanks in advance!!!
> > > 
> > > 
> > >       
> > >
> > __
> > > __
> > > ¡Obtén la mejor experiencia en la web!
> > > Descarga gratis el nuevo Internet Explorer 8. 
> > > http://downloads.yahoo.com/ieak8/?l=e1
> > > 
> > > __
> > > 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.
> > > 
> > 
> 
> 
>   
> __
> __
> ¡Obtén la mejor experiencia en la web!
> Descarga gratis el nuevo Internet Explo

[R] Is there a way not to export the functions in the global namespace?

[Peng Yu]

Suppose that I 'library()' a package, in the 'NAMESPACE' file of the
package, it export some functions to the global namespace. Is there a
way to 'library()' the package without export the functions in
'NAMESPACE' to the global namespace?

__
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] The RSQLite version of dbGetQuery drops colums

[Magnus Torfason]

Hi all,

I just noticed (the hard way of course) that when a query returns 0 
rows, the columns in the resulting data.frame get dropped as well. See 
the following example code (where conn is an active connection to an 
SQLite db):


> dbGetQuery(conn, "select 1 as hey, 2 as ho where 1")
  hey ho
1   1  2
> dbGetQuery(conn, "select 1 as hey, 2 as ho where 0")
data frame with 0 columns and 0 rows

I believe that the second query should return a 0x2 data.frame instead, 
that is, the same value as:


> dbGetQuery(conn, "select 1 as hey, 2 as ho where 1")[FALSE,]
[1] hey ho
<0 rows> (or 0-length row.names)

Any thoughts? Is this a bug, and are the developers of RSQLite reading this?

Best regards,
Magnus

__
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] Numerical Integration

[Julio Rojas]

Thanks a lot William. I'm sorry about the syntax problem. I was working at the 
same time with R code and a document, so I copied from the wrong document.

I saw that I had a problem with the conditions of the integrand. A ">" instead 
of a ">=" on the second one. Again, thanks.

One last question: Is there a way to use "approx" as the integrand?

Best regards.

Julio


--- El vie 18-dic-09, William Dunlap  escribió:

> De: William Dunlap 
> Asunto: RE: [R] Numerical Integration
> A: "Julio Rojas" 
> Fecha: viernes, 18 diciembre, 2009, 4:03 pm
> The immediate problem is that
> integrand(0.5) is NULL
> because none of the if clauses is TRUE if x==fx[2].
> 
> If you restructured the code such that it had no return
> statements and instead assigned things to a variable
> and returned that variable at the end it would be more
> apparent what happened: you would get a 'variable not
> found' error when returning.
>    i <- function(x) {
>       if (x       else if (x>fx[1] &&
> x<=fx[2]) retval <- 1
>       else if (x>fx[2]) retval <- 0
>       retval
>    }
>    > i(.4)
>    [1] 1
>    > i(.3)
>    Error in i(0.3) : object 'retval' not
> found
> 
> When you convert your code for production use, look
> at the ifelse and approx functions for doing this
> sort of thing.  They are much faster than
> Vector(yourIntegrand).
> 
> Also, when writing for help from R-help, please use
> R syntax when defining variables, like
>    fx<-c(0.3,0.5,0.5,0.6)
> and not some other language, like
>    fx<-[0.3,0.5,0.5,0.6]
> If we can copy and paste your example into R instead
> of translating the syntax by hand we will
> be much more inclined to try to fix things up.
> 
> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com  
> 
> > -Original Message-
> > From: r-help-boun...@r-project.org
> 
> > [mailto:r-help-boun...@r-project.org]
> On Behalf Of Julio Rojas
> > Sent: Friday, December 18, 2009 7:02 AM
> > To: r-help@r-project.org
> > Subject: [R] Numerical Integration
> > 
> > Dear @ll. I have to calculate numerical integrals for
> 
> > triangular and trapezoidal figures. I know you can
> calculate 
> > the exactly, but I want to do it this way to learn how
> to 
> > proceed with more complicated shapes. The code I'm
> using is 
> > the following:
> > 
> > integrand<-function(x) {
> >    print(x)
> >    if(x >    if(x>=fx[1] && x return((x-fx[1])/(fx[2]-fx[1]))
> >    if(x>fx[2] && x<=fx[3])
> return(1)
> >    if(x>fx[3] && x<=fx[4])
> return((x-fx[4])/(fx[3]-fx[4]))
> >    if(x>fx[4]) return(0)
> > 
> >  }
> > 
> > fx<-data[i,j,]
> > reltol<-1e-07
> >
> integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions
> > =200)$value
> > 
> > It works for most cases, but then, I tried for the
> triangle 
> > fx<-[0.3,0.5,0.5,0.6] and the following error was
> presented:
> > 
> > >
> integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions=200)
> > [1] 0.5
> > [1] 0.01304674
> > [1] 0.9869533
> > [1] 0.06746832
> > [1] 0.9325317
> > [1] 0.1602952
> > [1] 0.8397048
> > [1] 0.2833023
> > [1] 0.7166977
> > [1] 0.4255628
> > [1] 0.5744372
> > [1] 0.002171418
> > [1] 0.9978286
> > [1] 0.03492125
> > [1] 0.9650787
> > [1] 0.1095911
> > [1] 0.8904089
> > [1] 0.2186214
> > [1] 0.7813786
> > [1] 0.3528036
> > [1] 0.6471964
> > Error in integrate(Vectorize(integrand), 0, 1, rel.tol
> = reltol,
> > subdivisions = 200) :
> >  evaluation of function gave a result of wrong
> type
> > 
> > Does anybody know what happened? Thanks in advance!!!
> > 
> > 
> >       
> >
> __
> > __
> > ¡Obtén la mejor experiencia en la web!
> > Descarga gratis el nuevo Internet Explorer 8. 
> > http://downloads.yahoo.com/ieak8/?l=e1
> > 
> > __
> > 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.
> > 
> 


  

¡Obtén la mejor experiencia en la web!
Descarga gratis el nuevo Internet Explorer 8. 
http://downloads.yahoo.com/ieak8/?l=e1

__
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] Legend for two plots

[Marc Giombetti]

Dear R users,

I am new to R and I couldn't figure out how to solve the following
problem:

I am trying to put a legend below two plots using the code below. The
legend appears in the second plot,
but I want the legend to appear below the two plots in the center of
the total chart. At the moment the
graphic looks like this: http://i48.tinypic.com/2h2fvhf.jpg

[code]
layout(matrix(1:2, nrow=1))
#c(down,left,top,right)
par(mar=c(4,4,3,1))
plot(totalExp , totalDiffs,
main="Years of experience",
cex.main = 0.9,
cex.lab=0.8,
xlab="Years of experience",
ylab="COCOMO II - expert estimate",
pch=totalPch,
col = totalColors
)
abline(0,0)

plot(totalSoftwareEXP , totalDiffs,
main="Years of prof. software experience",
cex.main = 0.9,
cex.lab=0.8,
xlab="Years of experience",
ylab="COCOMO II - expert estimate",
pch=totalPch,
col = totalColors
)
abline(0,0)

legend("topleft",  c("Security","Usability") ,col=c('red','blue'),
pch=c(8,9), cex=0.8,bg='#e0dcdc')
[/code]

Thanks a lot for your help! I really appreciate it!

Marc

__
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] Covariate adjusted survival curves

[Terry Therneau]

> We are using frailty models to estimate risk of one year death. Is
> there a way to generate survival curves adjusted for covariates and
> also include frailty term?

  I don't know how to do this as yet -- so it's not in the survival
package.  (On my list of things to think hard about; I'd like to know
the answer too, both how to compute and what the result is worth.)

Terry Therneau

__
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] birats example with R2WinBugs

[M JH]

Hello,

I am trying, and failing, to do the birats example from the WinBugs manual
with R2Winbugs.

I can manage the rats example OK. Also, I can manage birats in WinBugs.

Here is the code I am running, the .bug program, and the error message. Any
help would be gratefully received.

Thank you,
Mike


library(R2WinBUGS)

data = list(x = c(8.0, 15.0, 22.0, 29.0, 36.0), N = 30, T = 5,
Omega = structure(.Data = c(200, 0, 0, 0.2), .Dim = c(2, 2)),   
mean = c(0,0),
prec = structure(.Data = c(1.0E-6, 0, 0, 1.0E-6), .Dim = c(2, 
2)),
Y = structure(
.Data =   c(151, 199, 246, 283, 320,
 145, 199, 249, 293, 
354,
 147, 214, 263, 312, 
328,
 155, 200, 237, 272, 
297,
 135, 188, 230, 280, 
323,
 159, 210, 252, 298, 
331,
 141, 189, 231, 275, 
305,
 159, 201, 248, 297, 
338,
 177, 236, 285, 350, 
376,
 134, 182, 220, 260, 
296,
 160, 208, 261, 313, 
352,
 143, 188, 220, 273, 
314,
 154, 200, 244, 289, 
325,
 171, 221, 270, 326, 
358,
 163, 216, 242, 281, 
312,
 160, 207, 248, 288, 
324,
 142, 187, 234, 280, 
316,
 156, 203, 243, 283, 
317,
 157, 212, 259, 307, 
336,
 152, 203, 246, 286, 
321,
 154, 205, 253, 298, 
334,
 139, 190, 225, 267, 
302,
 146, 191, 229, 272, 
302,
 157, 211, 250, 285, 
323,
 132, 185, 237, 286, 
331,
 160, 207, 257, 303, 
345,
 169, 216, 261, 295, 
333,
 157, 205, 248, 289, 
316,
 137, 180, 219, 258, 
291,
 153, 200, 244, 286, 
324),
.Dim = c(30,5)))
  

inits <- list(mu.beta = c(0,0), tauC = 1,
  beta = structure(
.Data = c(100,6,100,6,100,6,100,6,100,6,
   100,6,100,6,100,6,100,6,100,6,
   100,6,100,6,100,6,100,6,100,6,
   100,6,100,6,100,6,100,6,100,6,
   100,6,100,6,100,6,100,6,100,6,
   100,6,100,6,100,6,100,6,100,6),
.Dim = c(30, 2)),
  R = structure(.Data = c(1,0,0,1), .Dim = c(2, 2)))
  
  
  names(inits)
  inits$R
  is.matrix(inits$beta)
  is.matrix(inits$R)
  
mhbirats = bugs(data, inits, model.file = "birats.bug",   

  parameters.to.save = c("beta", "R", "tauC",  "mu.beta"),   

   n.chains = 1, n.iter=1, n.burnin=1000, n.thin=1,

   bugs.directory = "c:/Program Files/WinBUGS14/",  debug = TRUE)

### birats example from winbugs

model
{
for( i in 1 : N ) {
beta[i , 1:2] ~ dmnorm(mu.beta[], R[ , ])   
for( j in 1 : T ) {
Y[i , j] ~ dnorm(mu[i , j], tauC)   

mu[i , j] <- beta[i , 1] + beta[i , 2] * x[j]
}
}

mu.beta[1:2] ~ dmnorm(mean[],prec[ , ])
R[1:2 , 1:2] ~ dwish(Omega[ , ], 2) 
tauC ~ dgamma(0.001, 0.001)
sigma <- 1 / sqrt(tauC) 
}

### the error message


Error in bugs(data, inits, m

Re: [R] Problem reading binaries created with fortran

[Don MacQueen]

However, source code is available at
  http://water.usgs.gov/nrp/gwsoftware/modflow2005/modflow2005.html
so it would seem that the details are available.

-Don

At 4:35 PM -0500 12/17/09, kapo coulibaly wrote:

Duncan,
I couldn't find clear details about the way fortran writes the binaries. I
was hoping someone here has done it before. But the hexView package seems
like a great idea i'll give it a shot.

Thanks a bunch

On Thu, Dec 17, 2009 at 4:02 PM, Duncan Murdoch wrote:


 On 17/12/2009 3:48 PM, kapo coulibaly wrote:


 Is it possible to read fortran binaries with R? I tried unsucessfully and
 my
 understanding is that fortran write binaries with leading and trailing
 bytes. I get numbers but not the right ones.
 Thanks

 ps: the binary I'm interested in reading is a MODFLOW output with a mix of
 character, double and integers.



 One other thing:  the hexView package does a very nice job of displaying
 the file, so you can work out what the structure is if the documentation is
 unclear (or nonexistent).

 Duncan Murdoch



[[alternative HTML version deleted]]

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



--
--
Don MacQueen
Environmental Protection Department
Lawrence Livermore National Laboratory
Livermore, CA, USA
925-423-1062

__
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] write.csv and col.names=F

[Don MacQueen]

At 9:37 PM +0800 12/18/09, 
=?UTF-8?B?UmVleWFybl/mnY7mmbrmtItfMTA5MjgxMTM=?= wrote:

On Fri, Dec 18, 2009 at 7:52 AM, kayj  wrote:


 Hi All,

 I always have a problem with write.csv when I want the column names to be
 ignored, when I specify col.names=F, I get a header of V1 V2 V3 V4 etc.



I tried that and found the same problem, however, I found
  write.table(mydata, file="data.csv",col.names=F)
works.

write.csv calls write.table to save data, is there something wrong with it?


No.

If you read the help page you will find that it answers your question. It says:

 'write.csv' and 'write.csv2' provide convenience wrappers for
 writing CSV files.  They set 'sep', 'dec' and 'qmethod', and
 'col.names' to 'NA' if 'row.names = TRUE' and 'TRUE' otherwise.

Notice that it sets col.names for you.
If you want full control, use write.table().
If you want to write a table in a commonly used and pre-defined 
format (and without having to worry about various options) use 
write.csv().


That's why it is a called "convenience wrapper".

-Don


--
Best Regards,
Reeyarn T. Lee

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



--
--
Don MacQueen
Environmental Protection Department
Lawrence Livermore National Laboratory
Livermore, CA, USA
925-423-1062

__
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] error when using multcomp and lm

[Torsten Hothorn]

On Fri, 18 Dec 2009, Achim Zeileis wrote:


On Fri, 18 Dec 2009, chrischizinski wrote:

Just in case anyone ever has similar issues, it appears that 'glht' does 
not work without an intercept.


That is not precise enough to be useful. I've had a look again and appears to 
be the case that mcp() (or rather the internal function mcp2matrix) does not 
work for models where there are multiple factors but no intercept.


As the posting guide asks us to do, I have (1) provided a _reproducible_ 
example below, and (2) included the maintainer of the package in the posting. 
Torsten, it would be great if you could have a look.


Achim and Chris,

thanks for spotting this problem. The bug is now fixed in multcomp 
1.1-3 available from R-forge.


Torsten



Best,
Z

## package
library("multcomp")

## various models with and without intercept
m1a <- lm(breaks ~ tension, data = warpbreaks)
m1b <- lm(breaks ~ 0 + tension, data = warpbreaks)
m2a <- lm(breaks ~ wool + tension, data = warpbreaks)
m2b <- lm(breaks ~ 0 + wool + tension, data = warpbreaks)

## these two are equivalent: one factor with/without intercept
summary(glht(m1a, linfct = mcp(tension = "Tukey")))
summary(glht(m1b, linfct = mcp(tension = "Tukey")))

## these two should be equivalent: two factors with/without intercept
## but the latter fails
summary(glht(m2a, linfct = mcp(tension = "Tukey")))
summary(glht(m2b, linfct = mcp(tension = "Tukey")))




__
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] linear contrasts for trends in an anova

[ONKELINX, Thierry]

Dear Colm,

Use lm() instead of aov() and your problem is solved. The parameter of Group.L 
is the linear trend along your groups. Group.Q the quandratic trend.

HTH,

Thierry

summary(lm(Mean_1 ~ Group*eventType, data=doi)) 
Call:
lm(formula = Mean_1 ~ Group * eventType, data = doi)

Residuals:
 Min   1Q   Median   3Q  Max 
-2.88504 -0.74862 -0.01771  0.92154  2.62031 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)
(Intercept)  5.1360 0.3610  14.229 3.41e-13 ***
Group.L  2.7761 0.6252   4.440 0.000172 ***
Group.Q -0.6375 0.6252  -1.020 0.318091
eventTypestops  -2.0464 0.5105  -4.009 0.000515 ***
Group.L:eventTypestops  -2.3123 0.8842  -2.615 0.015175 *  
Group.Q:eventTypestops   1.9637 0.8842   2.221 0.036049 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.398 on 24 degrees of freedom
Multiple R-squared: 0.6357, Adjusted R-squared: 0.5598 
F-statistic: 8.375 on 5 and 24 DF,  p-value: 0.0001073 





ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium

Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium

tel. + 32 54/436 185
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 
Colm Connolly
Verzonden: vrijdag 18 december 2009 16:13
Aan: r-h...@stat.math.ethz.ch
Onderwerp: [R] linear contrasts for trends in an anova

Hi everybody,

I'm trying to construct contrasts for an ANOVA to determine if there is a 
significant trend in the means of my groups.

In the following example, based on the type of 2x3 ANOVA I'm trying to perform, 
does the linear polynomial contrast generated by contr.poly allow me to test 
for a linear trend across groups?

doi=data.frame(
  Group=c(
rep(1, 5),  rep(2, 5),  rep(3, 5), rep(1, 5),  rep(2, 5),  rep(3, 5)),
  Mean_1=c(
rnorm(5, mean=3, sd=1), rnorm(5, mean=5.7, sd=1.2), rnorm(5, mean=7, 
sd=1.3),
rnorm(5, mean=3.5, sd=1), rnorm(5, mean=3, sd=1.2), rnorm(5, mean=4.2, 
sd=1.7)),
  eventType=c( rep("errors", 15), rep("stops", 15))) 
doi$Group=factor(doi$Group, labels=c("Control", "Short", "Long"), ordered=T)

## construct contrasts
contrasts(doi$Group)=contr.poly(levels(doi$Group))

model1=aov(Mean_1 ~ Group*eventType, data=doi)
model1.summary=summary(model1)
print(model1.summary)
print(summary.lm(model1))

Thanks in advance,

Regards,
--
Dr Colm G. Connolly
Institute of Neuroscience
The Lloyd Building
University of Dublin
Trinity College, Dublin 2, Éire
Fax: +353-1-671-3183

Please note that electronic mail to, from or within the Trinity College Dublin, 
may be the subject of a request under the Freedom of Information Act.




Druk dit bericht a.u.b. niet onnodig af.
Please do not print this message unnecessarily.

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.

__
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] NLS-Weibull-ERROR

[Erwin Kalvelagen]

ruchita gupta  gmail.com> writes:

> 
> Hello
> 
> I was trying to estimate the weibull model using nls after putting OLS
> values as the initial inputs to NLS.
> I tried multiple times but still i m getting the same error of Error in
> nlsModel(formula, mf, start, wts) :
>   singular gradient matrix at initial parameter estimates.
> 
> The Program is as below
> 
> > vel <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
> > df <- data.frame(conc, vel)
> > df
>   conc vel
> 10.077   1
> 20.328   2
> 30.882   3
> 41.195   4
> 51.884   5
> 63.577   6
> 76.549   7
> 8   13.000   8
> 9   33.690   9
> 10  52.220  10
> 11  90.140  11
> 12 166.050  12
> 13 233.620  13
> 14 346.890  14
> > plot(df$vel, df$conc)
> > para0.st <- c(K=450,
> +   alpha=0.054,beta=3.398 )
> >  fit0 <- nls(
> +  conc~ K-(K*exp(-(vel/alpha)^beta)), df, start= para0.st,trace=T)
> Error in nlsModel(formula, mf, start, wts) :
>   singular gradient matrix at initial parameter estimates
> 
> I will be highly thankful if some one can please let me know where is the
> mistake as i m unable to trace it.
> 
> Thanks
> Ruchita
> 


This is a better fit:

> nls(conc~K+alpha*vel^beta,df,start<-c(K=-0.7,alpha=0.0003,beta=5.3),trace=T)
349.1687 :  -0.7000  0.0003  5.3000 
264.4474 :  -0.7507993362  0.0002808334  5.3172851705 
263.6331 :  -0.7559514801  0.0002816438  5.316930 
263.6331 :  -0.7558506044  0.0002816331  5.3169456595 
Nonlinear regression model
  model:  conc ~ K + alpha * vel^beta 
   data:  df 
 K  alpha   beta 
-0.7558506  0.0002816  5.3169457 
 residual sum-of-squares: 263.6

Number of iterations to convergence: 3 
Achieved convergence tolerance: 1.599e-06 
> 



Erwin Kalvelagen
Amsterdam Optimization Modeling Group
er...@amsterdamoptimization.com
http://amsterdamoptimization.com

__
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] TRIANGLE AW: lattice: shape of box around wireframe

[David Winsemius]

On Dec 18, 2009, at 10:02 AM, Thomas Roth wrote:


Sorry,

I mixed up some words...
I wrote rectangle but meant TRIANGLE, sorry for that.


Huh? You want a triangular bounding box? We are supposed to fill in  
the rest of the details, ... how?




As for the book, I own it, studied it as well as the help pages but  
didn't
come to a conclusion for my problem. Aspect can't help me with my  
problem. I
need to set up the axes and obviously make some transformations and  
that's

where I'm looking for some help...



User-level control of the display is at the level of (1) either the
panel function, which is somewhat non-trivial for the 3-D functions,
or (2) the 'panel.3d.wireframe' function (for wireframe), where you
can basically do anything after the axes and the transformations have
been set up.
---

Thomas





-Ursprüngliche Nachricht-
Von: David Winsemius [mailto:dwinsem...@comcast.net]
Gesendet: Freitag, 18. Dezember 2009 15:32
An: Thomas Roth
Cc: 'Deepayan Sarkar'; r-help@r-project.org
Betreff: Re: [R] lattice: shape of box around wireframe


On Dec 18, 2009, at 8:52 AM, Thomas Roth wrote:


Hi,

Thank you for the quick answer...

I made a small example:

require(lattice)
x = 1:100
x = numeric(0)
z = numeric(0)
y = numeric(0)
for(i in 1:100)
{
y = c(y,1:i)
x = c(x, rep(i,i))
z = x + y
}

#omitted bounding box
wireframe(z ~ x*y, par.box = c(col = "transparent"))

#needed: rectangular bounding box
wireframe(z ~ x*y)


now what i'm looking for is a rectangular bounding box but I don't
know where to start...


When I tried to use the Lattice book, I found that I could not locate
the code that produced the example I remembered (Figure  and could not
initially findo Figures 13.8 and 13.9 (because they were in the color
plates section in chapter 3) that appeared to be the result for the
section to which I was referred from chapter 6.5.  So I would suggest
going to the book's website:

http://lmdvr.r-forge.r-project.org/figures/figures.html

and looking at the code for 13.9 and 6.5, both of which use cloud but
I think the parameter you want is aspect, which is also in the
parameter list for wireframe. Further examination of the wireframe
help page examples suggests that you have not yet looked at that
resource, since the first one uses aspect to create a rectangular
bounding box.

--
David.

r


Thomas

-Ursprüngliche Nachricht-
Von: Deepayan Sarkar [mailto:deepayan.sar...@gmail.com]
Gesendet: Freitag, 18. Dezember 2009 14:29
An: Thomas Roth
Cc: r-help@r-project.org
Betreff: Re: [R] lattice: shape of box around wireframe

On Fri, Dec 18, 2009 at 5:19 AM, Thomas Roth
 wrote:

Hi,

I was wondering if there is a way to adjust the shape of the box
around the
wireframe. By default this box is always a cube. For instance, is it
possible to cut this cube into two halfs each half being a 3d
rectangle? Or
just plot a 3d rectangle with a wireframe inside and adjusted axes?

So far I've read "Lattice: Multivariate Data Visualization with R"
and
searched through the mailing list.


User-level control of the display is at the level of (1) either the
panel function, which is somewhat non-trivial for the 3-D functions,
or (2) the 'panel.3d.wireframe' function (for wireframe), where you
can basically do anything after the axes and the transformations have
been set up. There are examples both in the book and on the
r-help/r-devel lists. Let me know if you need help in figuring out  
any

specific details (such as how to omit the default bounding box by
setting its color to transparent).

-Deepayan

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


David Winsemius, MD
Heritage Laboratories
West Hartford, CT




David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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 define new operators

[Duncan Murdoch]

On 18/12/2009 10:22 AM, berga...@gmail.com wrote:

Dear R community

I try to create a new operator to build a special sum of two CashFlows. (my  
S4 Class)


I tried the following but this is actually not what I want.

setGeneric("++",function(e1,e2)standardGeneric("++"))

setMethod("+",signature=list("CashFlow","CashFlow"),function(e1,e2){
print("+")
})
setMethod("++",signature=list("CashFlow","CashFlow"),function(e1,e2){
print("++")
})

The problem here is that this work but it treats "++" as a general  
function.. the call is "++"(e1,e2).

I'm looking for somthing like e3 <- e1 ++ e2
The R parser limits the operators it recognizes to the predefined ones, 
plus user-defined ones of the form %op% (where op could be ++ if you 
liked).  You can't define ++ as an operator because the parser will see 
it as two plus signs, not as a single infix operator.


Duncan Murdoch

__
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] forecasting

[Ista Zahn]

Hi,
To save the output to a text file you can use

sink(file="results.txt")
summary(output_dat)
sink()

See ?sink for details.

To output to a table, you can use xtable or Hmisc::latex.

-Ista
On Fri, Dec 18, 2009 at 6:34 AM, DispersionMap  wrote:
>
>
> Hi thanks.
>
> Another question...
>
> When R prints the regression output onto the screen after i perform my
> regression:
>
> output_dat <- lm(
> results <- summary(output_dat)
> results
>
> How do i print the regression output to a file for inclusion into a
> document.
>
> I tried:  cat(summary(output_dat), file="results.txt"))
>
> but this doesnt work and it tells me type list (in the first argument)
> cannot be handled by cat.
>
>
> So how do i get the output of lm to the desktopas i need to put this in
> a document (in atble would be even better).
>
>
>
> Stephan Kolassa wrote:
>>
>> ?predict
>>
>> HTH
>> Stephan
>>
>> DispersionMap schrieb:
>>>
>>> I have some data that i ran a regression on and got the usual r output,
>>> with
>>> the a intercept and b coefficient for my independent variable.
>>>
>>> i want to forecast the number of future events using these parameters.
>>>
>>>
>>> What function / packages in R would you recommend to be used in good
>>> effect
>>> for these purposes??
>>
>> __
>> 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.
>>
>>
>
> --
> View this message in context: 
> http://n4.nabble.com/forecasting-tp974326p974704.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> 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.
>



-- 
Ista Zahn
Graduate student
University of Rochester
Department of Clinical and Social Psychology
http://yourpsyche.org

__
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] What is the fastest way to see what are in an RData file?

[Peng Yu]

On Sat, Dec 19, 2009 at 3:35 AM, Patrick Connolly
 wrote:
> On Thu, 17-Dec-2009 at 03:13PM +1800, Peng Yu wrote:
>
> |> Currently, I load the RData file then ls() and str(). But loading the file
> |> takes too long if the file is big. Most of the time, I only interested what
> |> the variables are in the the file and the attributes of the variables (like
> |> if it is a data.frame, matrix, what are the colnames/rownames, etc.)
> |>
> |> I'm wondering if there is any facility in R to help me avoid loading the
> |> whole file.
>
>
> I have a pretty nifty way of seeing what's in such a file, but I still
> have to load all of the binary file before I can do so.  If it's
> taking you such a long time, maybe you could keep a larger number of
> smaller RData files.

What is your 'nifty way'? How fast is it?

__
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] How to define new operators

[bergarog]

Dear R community

I try to create a new operator to build a special sum of two CashFlows. (my  
S4 Class)

I tried the following but this is actually not what I want.

setGeneric("++",function(e1,e2)standardGeneric("++"))

setMethod("+",signature=list("CashFlow","CashFlow"),function(e1,e2){
print("+")
})
setMethod("++",signature=list("CashFlow","CashFlow"),function(e1,e2){
print("++")
})

The problem here is that this work but it treats "++" as a general  
function.. the call is "++"(e1,e2).
I'm looking for somthing like e3 <- e1 ++ e2

Thanks a lot for your help.
Best regards,
Roger

[[alternative HTML version deleted]]

__
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] linear contrasts for trends in an anova

[Colm Connolly]

Hi everybody,

I'm trying to construct contrasts for an ANOVA to determine if there is a 
significant trend in the means of my groups.

In the following example, based on the type of 2x3 ANOVA I'm trying to perform, 
does the linear polynomial contrast generated by contr.poly allow me to test 
for a linear trend across groups?

doi=data.frame(
  Group=c(
rep(1, 5),  rep(2, 5),  rep(3, 5), rep(1, 5),  rep(2, 5),  rep(3, 5)),
  Mean_1=c(
rnorm(5, mean=3, sd=1), rnorm(5, mean=5.7, sd=1.2), rnorm(5, mean=7, 
sd=1.3),
rnorm(5, mean=3.5, sd=1), rnorm(5, mean=3, sd=1.2), rnorm(5, mean=4.2, 
sd=1.7)),
  eventType=c( rep("errors", 15), rep("stops", 15)))
doi$Group=factor(doi$Group, labels=c("Control", "Short", "Long"), ordered=T)

## construct contrasts
contrasts(doi$Group)=contr.poly(levels(doi$Group))

model1=aov(Mean_1 ~ Group*eventType, data=doi)
model1.summary=summary(model1)
print(model1.summary)
print(summary.lm(model1))

Thanks in advance,

Regards,
--
Dr Colm G. Connolly
Institute of Neuroscience
The Lloyd Building
University of Dublin
Trinity College, Dublin 2, Éire
Fax: +353-1-671-3183

Please note that electronic mail to, from or within the Trinity College Dublin, 
may be the
subject of a request under the Freedom of Information Act.



__
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] mutlidimensional in.convex.hull(wasmultidimensionalpoint.in.polygon??)

[Keith Jewell]

Syntax suggestions implemented.
Inhull's original author consulted.
Function submitted for potential inclusion in geometry package.

Seasons greetings to all.

Keith Jewell
-
"baptiste auguie"  wrote in message 
news:de4e29f50912180238m45c4b7c3g792113d5b6c4c...@mail.gmail.com...
Hi,

Excellent, thanks for doing this!

I had tried the 2D case myself but I was put off by the fact that
Octave's convhulln had a different ordering of the points to R's
geometry package (an improvement to Octave's convhulln was made after
it was ported to R). I'm not sure how you got around this but it's
good to see that the result was worth the effort!

Below are a few minor syntax suggestions,

# p <- dim(calpts)[2]   # columns in calpts
# and other similar lines could be replaced with
ncol(calpts)
nrow(testpts)
nrow(hull)

# length(degenflag[degenflag])
# can probably be written
sum(degenflag)

# center = apply(calpts, 2, mean)
# more efficient
colMeans(calpts)

Would you consider submitting this function to the maintainer of the
geometry package, after checking it's OK with inhull's original
author?

Best regards,

baptiste

__
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] Numerical Integration

[Julio Rojas]

Dear @ll. I have to calculate numerical integrals for triangular and 
trapezoidal figures. I know you can calculate the exactly, but I want to do it 
this way to learn how to proceed with more complicated shapes. The code I'm 
using is the following:

integrand<-function(x) {
   print(x)
   if(x=fx[1] && xfx[2] && x<=fx[3]) return(1)
   if(x>fx[3] && x<=fx[4]) return((x-fx[4])/(fx[3]-fx[4]))
   if(x>fx[4]) return(0)

 }

fx<-data[i,j,]
reltol<-1e-07
integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions=200)$value

It works for most cases, but then, I tried for the triangle 
fx<-[0.3,0.5,0.5,0.6] and the following error was presented:

> integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions=200)
[1] 0.5
[1] 0.01304674
[1] 0.9869533
[1] 0.06746832
[1] 0.9325317
[1] 0.1602952
[1] 0.8397048
[1] 0.2833023
[1] 0.7166977
[1] 0.4255628
[1] 0.5744372
[1] 0.002171418
[1] 0.9978286
[1] 0.03492125
[1] 0.9650787
[1] 0.1095911
[1] 0.8904089
[1] 0.2186214
[1] 0.7813786
[1] 0.3528036
[1] 0.6471964
Error in integrate(Vectorize(integrand), 0, 1, rel.tol = reltol,
subdivisions = 200) :
 evaluation of function gave a result of wrong type

Does anybody know what happened? Thanks in advance!!!


  

¡Obtén la mejor experiencia en la web!
Descarga gratis el nuevo Internet Explorer 8. 
http://downloads.yahoo.com/ieak8/?l=e1

__
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] TRIANGLE AW: lattice: shape of box around wireframe

[Thomas Roth]

Sorry,

I mixed up some words... 
I wrote rectangle but meant TRIANGLE, sorry for that. 

As for the book, I own it, studied it as well as the help pages but didn't
come to a conclusion for my problem. Aspect can't help me with my problem. I
need to set up the axes and obviously make some transformations and that's
where I'm looking for some help...



User-level control of the display is at the level of (1) either the
 panel function, which is somewhat non-trivial for the 3-D functions,
 or (2) the 'panel.3d.wireframe' function (for wireframe), where you
 can basically do anything after the axes and the transformations have
 been set up.
---

Thomas



 

-Ursprüngliche Nachricht-
Von: David Winsemius [mailto:dwinsem...@comcast.net] 
Gesendet: Freitag, 18. Dezember 2009 15:32
An: Thomas Roth
Cc: 'Deepayan Sarkar'; r-help@r-project.org
Betreff: Re: [R] lattice: shape of box around wireframe


On Dec 18, 2009, at 8:52 AM, Thomas Roth wrote:

> Hi,
>
> Thank you for the quick answer...
>
> I made a small example:
>
> require(lattice)
> x = 1:100
> x = numeric(0)
> z = numeric(0)
> y = numeric(0)
> for(i in 1:100)
> {
>  y = c(y,1:i)
>  x = c(x, rep(i,i))
>  z = x + y
> }
>
> #omitted bounding box
> wireframe(z ~ x*y, par.box = c(col = "transparent"))
>
> #needed: rectangular bounding box
> wireframe(z ~ x*y)
>
>
> now what i'm looking for is a rectangular bounding box but I don't  
> know where to start...

When I tried to use the Lattice book, I found that I could not locate  
the code that produced the example I remembered (Figure  and could not  
initially findo Figures 13.8 and 13.9 (because they were in the color  
plates section in chapter 3) that appeared to be the result for the  
section to which I was referred from chapter 6.5.  So I would suggest  
going to the book's website:

http://lmdvr.r-forge.r-project.org/figures/figures.html

and looking at the code for 13.9 and 6.5, both of which use cloud but  
I think the parameter you want is aspect, which is also in the  
parameter list for wireframe. Further examination of the wireframe  
help page examples suggests that you have not yet looked at that  
resource, since the first one uses aspect to create a rectangular  
bounding box.

-- 
David.

r
>
> Thomas
>
> -Ursprüngliche Nachricht-
> Von: Deepayan Sarkar [mailto:deepayan.sar...@gmail.com]
> Gesendet: Freitag, 18. Dezember 2009 14:29
> An: Thomas Roth
> Cc: r-help@r-project.org
> Betreff: Re: [R] lattice: shape of box around wireframe
>
> On Fri, Dec 18, 2009 at 5:19 AM, Thomas Roth
>  wrote:
>> Hi,
>>
>> I was wondering if there is a way to adjust the shape of the box  
>> around the
>> wireframe. By default this box is always a cube. For instance, is it
>> possible to cut this cube into two halfs each half being a 3d  
>> rectangle? Or
>> just plot a 3d rectangle with a wireframe inside and adjusted axes?
>>
>> So far I've read "Lattice: Multivariate Data Visualization with R"  
>> and
>> searched through the mailing list.
>
> User-level control of the display is at the level of (1) either the
> panel function, which is somewhat non-trivial for the 3-D functions,
> or (2) the 'panel.3d.wireframe' function (for wireframe), where you
> can basically do anything after the axes and the transformations have
> been set up. There are examples both in the book and on the
> r-help/r-devel lists. Let me know if you need help in figuring out any
> specific details (such as how to omit the default bounding box by
> setting its color to transparent).
>
> -Deepayan
>
> __
> 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.

David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] error when using multcomp and lm

[chrischizinski]

Just in case anyone ever has similar issues, it appears that 'glht' does not
work without an intercept.


-- 
View this message in context: 
http://n4.nabble.com/error-when-using-multcomp-and-lm-tp964705p974823.html
Sent from the R help mailing list archive at Nabble.com.

__
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] lattice: shape of box around wireframe

[David Winsemius]

On Dec 18, 2009, at 8:52 AM, Thomas Roth wrote:


Hi,

Thank you for the quick answer...

I made a small example:

require(lattice)
x = 1:100
x = numeric(0)
z = numeric(0)
y = numeric(0)
for(i in 1:100)
{
 y = c(y,1:i)
 x = c(x, rep(i,i))
 z = x + y
}

#omitted bounding box
wireframe(z ~ x*y, par.box = c(col = "transparent"))

#needed: rectangular bounding box
wireframe(z ~ x*y)


now what i'm looking for is a rectangular bounding box but I don't  
know where to start...


When I tried to use the Lattice book, I found that I could not locate  
the code that produced the example I remembered (Figure  and could not  
initially findo Figures 13.8 and 13.9 (because they were in the color  
plates section in chapter 3) that appeared to be the result for the  
section to which I was referred from chapter 6.5.  So I would suggest  
going to the book's website:


http://lmdvr.r-forge.r-project.org/figures/figures.html

and looking at the code for 13.9 and 6.5, both of which use cloud but  
I think the parameter you want is aspect, which is also in the  
parameter list for wireframe. Further examination of the wireframe  
help page examples suggests that you have not yet looked at that  
resource, since the first one uses aspect to create a rectangular  
bounding box.


--
David.

r


Thomas

-Ursprüngliche Nachricht-
Von: Deepayan Sarkar [mailto:deepayan.sar...@gmail.com]
Gesendet: Freitag, 18. Dezember 2009 14:29
An: Thomas Roth
Cc: r-help@r-project.org
Betreff: Re: [R] lattice: shape of box around wireframe

On Fri, Dec 18, 2009 at 5:19 AM, Thomas Roth
 wrote:

Hi,

I was wondering if there is a way to adjust the shape of the box  
around the

wireframe. By default this box is always a cube. For instance, is it
possible to cut this cube into two halfs each half being a 3d  
rectangle? Or

just plot a 3d rectangle with a wireframe inside and adjusted axes?

So far I've read "Lattice: Multivariate Data Visualization with R"  
and

searched through the mailing list.


User-level control of the display is at the level of (1) either the
panel function, which is somewhat non-trivial for the 3-D functions,
or (2) the 'panel.3d.wireframe' function (for wireframe), where you
can basically do anything after the axes and the transformations have
been set up. There are examples both in the book and on the
r-help/r-devel lists. Let me know if you need help in figuring out any
specific details (such as how to omit the default bounding box by
setting its color to transparent).

-Deepayan

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


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] write.csv and col.names=F

[jim holtman]

In R.2.9.2 I get the following error message if setting col.names=FALSE:

> write.csv(x, '', col.names=FALSE)
"","a","b"
"1",1,1
"2",2,2
"3",3,3
"4",4,4
"5",5,5
"6",6,6
"7",7,7
"8",8,8
"9",9,9
"10",10,10
Warning message:
In write.csv(x, "", col.names = FALSE) : attempt to set 'col.names' ignored
You have to use write.table if you don't want the column names; it is on the
help page:

> write.table(x,sep=',', col.names=FALSE)
"1",1,1
"2",2,2
"3",3,3
"4",4,4
"5",5,5
"6",6,6
"7",7,7
"8",8,8
"9",9,9
"10",10,10
>


On Fri, Dec 18, 2009 at 8:37 AM, Reeyarn_李智洋_10928113 
wrote:

> On Fri, Dec 18, 2009 at 7:52 AM, kayj  wrote:
> >
> > Hi All,
> >
> > I always have a problem with write.csv when I want the column names to be
> > ignored, when I specify col.names=F, I get a header of V1 V2 V3 V4 etc.
> >
>
> I tried that and found the same problem, however, I found
>  write.table(mydata, file="data.csv",col.names=F)
> works.
>
> write.csv calls write.table to save data, is there something wrong with it?
>
> --
> Best Regards,
> Reeyarn T. Lee
>
> __
> 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.
>



-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem that you are trying to solve?

[[alternative HTML version deleted]]

__
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] mantel test and NAs

[Sarah Goslee]

You will most likely need to drop the samples that are causing the NA
values from both of your datasets. The function won't work with NA values.

Sarah

On Thu, Dec 17, 2009 at 12:33 PM, marciarocha  wrote:
>
> Dear R users,
>
> I am having a problem performing the mantel test (both with functions mantel
> {vegan} and mantel.test {ape}) due to I believe the presence of NAs on my
> distance matrices, which look like e.g.:
>
> NA
> 1    2
> 1    2   3
> NA  4   5  6
>
> and
>
> 1
> 1   2
> 2   3   4
> 5   6   7  8
>
> Would any of you have a solution for that?
>
> Thank you for much for your help!
>
> Marcia Rocha
>
>

-- 
Sarah Goslee
http://www.functionaldiversity.org

__
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] Doubts about ADE-4

[David Winsemius]

On Dec 18, 2009, at 6:49 AM, Diogo B. Provete wrote:


Hi,
I'm having some troubles in using the ADE-4 package. I want to  
analyze data
for the Outlying Mean Index in the function 'niche'. I've got two  
matrixes:
one for environmental data and the other for species abundance data.  
But to
implement the analysis, as I understand it, I need to put tocheger  
these two
matrixes in one, as in the example provided with the package. Un  
fortunately

I don't know how to do this. Could anyone help me?


Not familiar with the application you are rather vaguely describing,  
but the two functions for "putting two matrices together" would be  
rbind and cbind depending on whether you want them bound "side by  
side" or "on top of each other".


(I am reasonably sure that there is no ADE-4 package. This is R ...  
capitalization and punctuation matter ... please spell package names  
precisely as they are named.)


For the life of me I cannot see how one uses cbind or rbind for  
arguments to "niche". It takes as arguments the results of functions  
that create duality diagrams. The example data set on the niche help  
page is a 3 element list of dataframes each with 30 rows. So you  
should should be asking how to create lists of dataframes (not  
matrices) in a manner that parallels the output of :


library(ade4)
data(doubs)
str(doubs)

You should also be providing a much more complete description of your  
R objects, at a minimum the results of str on them.




Thank you in advance.

--
Respectfully,
Diogo Borges Provete




David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] some help regarding combining columns from different files

[jim holtman]

In your function, you have

 temp <- read.table(fnames,header=T,sep="\t",stringsAsFactors=F,quote="\"")

I think you mean:

 temp <- read.table(i,header=T,sep="\t",stringsAsFactors=F,quote="\"")

Also 'files' is a parameter, but you are using 'fnames' in the 'for' loop;
shouldn't that be 'files'?

On Thu, Dec 17, 2009 at 3:51 PM, Harikrishnadhar wrote:

> Dear all,
>
> Here is my code which am using to combine 5th column from different data
> sets.
>
> Here is the function  to do my job
>
>
> genesymbol.append.file <-NULL
> gene.column <- NULL
> readGeneSymbol <- function(files,genesymbol.column=5){
> for(i in fnames){
>  temp <- read.table(fnames,header=T,sep="\t",stringsAsFactors=F,quote="\"")
>  gene.column<-cbind(gene.column,temp[,genesymbol.column])
>  genesymbol.append.file$genecolumns <- gene.column
>  genesymbol.append.file
>  }
> }
>
>
>
>
> test <- readGeneSymbol(fnames,genesymbol.column=5)
>
> Here is the warning message  am getting only the 5th column from the first
> column is taken
>
>
> Warning messages:
> 1: In file(file, "r") : only first element of 'description' argument used
> 2: In file(file, "r") : only first element of 'description' argument used
> >
>
> Please help me to solve this
>
>
>
>
>
>
>
> --
> Thanks
> Hari
> 215-385-4122
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> "If there is anyone out there who still doubts that America is a place
> where
> all things are possible"
>
>[[alternative HTML version deleted]]
>
> __
> 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.
>



-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem that you are trying to solve?

[[alternative HTML version deleted]]

__
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] lattice: shape of box around wireframe

[Thomas Roth]

Hi,

Thank you for the quick answer...

I made a small example:

require(lattice)
x = 1:100
x = numeric(0)
z = numeric(0)
y = numeric(0)
for(i in 1:100)
{
  y = c(y,1:i)
  x = c(x, rep(i,i))
  z = x + y
} 

#omitted bounding box
wireframe(z ~ x*y, par.box = c(col = "transparent"))

#needed: rectangular bounding box
wireframe(z ~ x*y)


now what i'm looking for is a rectangular bounding box but I don't know where 
to start...

Thomas

-Ursprüngliche Nachricht-
Von: Deepayan Sarkar [mailto:deepayan.sar...@gmail.com] 
Gesendet: Freitag, 18. Dezember 2009 14:29
An: Thomas Roth
Cc: r-help@r-project.org
Betreff: Re: [R] lattice: shape of box around wireframe

On Fri, Dec 18, 2009 at 5:19 AM, Thomas Roth
 wrote:
> Hi,
>
> I was wondering if there is a way to adjust the shape of the box around the
> wireframe. By default this box is always a cube. For instance, is it
> possible to cut this cube into two halfs each half being a 3d rectangle? Or
> just plot a 3d rectangle with a wireframe inside and adjusted axes?
>
> So far I've read "Lattice: Multivariate Data Visualization with R" and
> searched through the mailing list.

User-level control of the display is at the level of (1) either the
panel function, which is somewhat non-trivial for the 3-D functions,
or (2) the 'panel.3d.wireframe' function (for wireframe), where you
can basically do anything after the axes and the transformations have
been set up. There are examples both in the book and on the
r-help/r-devel lists. Let me know if you need help in figuring out any
specific details (such as how to omit the default bounding box by
setting its color to transparent).

-Deepayan

__
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] write.csv and col.names=F

[Reeyarn_李智洋_10928113]

On Fri, Dec 18, 2009 at 7:52 AM, kayj  wrote:
>
> Hi All,
>
> I always have a problem with write.csv when I want the column names to be
> ignored, when I specify col.names=F, I get a header of V1 V2 V3 V4 etc.
>

I tried that and found the same problem, however, I found
  write.table(mydata, file="data.csv",col.names=F)
works.

write.csv calls write.table to save data, is there something wrong with it?

--
Best Regards,
Reeyarn T. Lee

__
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] forecasting

[David Winsemius]

On Dec 18, 2009, at 6:34 AM, DispersionMap wrote:




Hi thanks.

Another question...

When R prints the regression output onto the screen after i perform my
regression:

output_dat <- lm(
results <- summary(output_dat)
results

How do i print the regression output to a file for inclusion into a
document.



?capture.output



I tried:  cat(summary(output_dat), file="results.txt"))

but this doesnt work and it tells me type list (in the first argument)
cannot be handled by cat.


So how do i get the output of lm to the desktopas i need to put  
this in

a document (in atble would be even better).



Stephan Kolassa wrote:


?predict

HTH
Stephan

DispersionMap schrieb:


I have some data that i ran a regression on and got the usual r  
output,

with
the a intercept and b coefficient for my independent variable.

i want to forecast the number of future events using these  
parameters.



What function / packages in R would you recommend to be used in good
effect
for these purposes??


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




--
View this message in context: 
http://n4.nabble.com/forecasting-tp974326p974704.html
Sent from the R help mailing list archive at Nabble.com.

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


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] Testing equality of regression model on multiple groups

[David Winsemius]

On Dec 17, 2009, at 5:14 PM, Clara Yuan wrote:


Hello,






(Also, why are the coefficients renamed to t.L, t.Q, etc instead of  
t.1, t.2?)


You have coded your t variable as an ordered factor.

?ordered

if you do not want the linear, quadratic, cubic coefficients for the  
polynomial contrasts, then remove the ordering with:


t <- factor(t, is.ordered=FALSE)

... and rerun the analysis.

--
David



What am I missing?

Thanks for the help,
Clara

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


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] lattice: shape of box around wireframe

[Deepayan Sarkar]

On Fri, Dec 18, 2009 at 5:19 AM, Thomas Roth
 wrote:
> Hi,
>
> I was wondering if there is a way to adjust the shape of the box around the
> wireframe. By default this box is always a cube. For instance, is it
> possible to cut this cube into two halfs each half being a 3d rectangle? Or
> just plot a 3d rectangle with a wireframe inside and adjusted axes?
>
> So far I've read "Lattice: Multivariate Data Visualization with R" and
> searched through the mailing list.

User-level control of the display is at the level of (1) either the
panel function, which is somewhat non-trivial for the 3-D functions,
or (2) the 'panel.3d.wireframe' function (for wireframe), where you
can basically do anything after the axes and the transformations have
been set up. There are examples both in the book and on the
r-help/r-devel lists. Let me know if you need help in figuring out any
specific details (such as how to omit the default bounding box by
setting its color to transparent).

-Deepayan

__
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] Which hist cell each value falls in?

[David Winsemius]

On Dec 18, 2009, at 2:57 AM, Jim Lemon wrote:


On 12/18/2009 06:35 AM, Doug Hill wrote:
Hi, all. I'm using hist() to obtain a vector of break values in an  
interval.
I then want to be able to identify which cell any value from  
another vector

falls in.

E.g. applying



breaks


 [1] -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0

to



x

 [1] -3.74519666 -0.38183630 -1.22884247 -0.20971824 -0.30533939  
-0.36271207
[7] -2.27513499 -2.23688653 -1.98827155 -1.48666274 -1.26155084  
-0.15234555

[13] -0.09497287  0.34488440

would give



xcells


 [1] 1  8  6  8  8  8  4  4  5  6  6  8  8  9

where:
  x<= breaks[1] ->  cell 1
breaks[1]<  x<= breaks[2] ->  cell 2
breaks[2]<  x<= breaks[3] ->  cell 3, etc.


> breaks <- scan()
1: -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0
13:
Read 12 items
> xcells <- scan()
1: -3.74519666 -0.38183630 -1.22884247 -0.20971824 -0.30533939  
-0.36271207
7: -2.27513499 -2.23688653 -1.98827155 -1.48666274 -1.26155084  
-0.15234555

13:  -0.09497287  0.34488440
15:
Read 14 items

> breaks <- c(-Inf,breaks,Inf)
> findInterval(xcells, breaks)
 [1] 1 8 6 8 8 8 4 4 5 6 6 8 8 9

--
David




Hi Doug,
The function below does more or less what you want, except that the  
bins are numbered from zero (meaning that a value is below the range  
of bins as x[1] is) to length(breaks) (meaning that a value is above  
the range of bins).


whichBin<-function(x,breaks,right=TRUE) {
lenx<-length(x)
# any leftovers must be out of range
wb<-rep(lenx,lenx)
if(right) wb[x<=breaks[1]]<-0
else wb[xbreaks[bin] & x<=breaks[bin+1]]<-bin
 else wb[x>=breaks[bin] & xhttps://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] lattice: shape of box around wireframe

[Thomas Roth]

Hi,

I was wondering if there is a way to adjust the shape of the box around the
wireframe. By default this box is always a cube. For instance, is it
possible to cut this cube into two halfs each half being a 3d rectangle? Or
just plot a 3d rectangle with a wireframe inside and adjusted axes?

So far I've read "Lattice: Multivariate Data Visualization with R" and
searched through the mailing list.

I'm grateful for any ideas or hints.

Thanks in advance,

Thomas Roth

__
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] NLS-Weibull-ERROR

[ruchita gupta]

Hello

I was trying to estimate the weibull model using nls after putting OLS
values as the initial inputs to NLS.
I tried multiple times but still i m getting the same error of Error in
nlsModel(formula, mf, start, wts) :
  singular gradient matrix at initial parameter estimates.

The Program is as below

> vel <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
> df <- data.frame(conc, vel)
> df
  conc vel
10.077   1
20.328   2
30.882   3
41.195   4
51.884   5
63.577   6
76.549   7
8   13.000   8
9   33.690   9
10  52.220  10
11  90.140  11
12 166.050  12
13 233.620  13
14 346.890  14
> plot(df$vel, df$conc)
> para0.st <- c(K=450,
+   alpha=0.054,beta=3.398 )
>  fit0 <- nls(
+  conc~ K-(K*exp(-(vel/alpha)^beta)), df, start= para0.st,trace=T)
Error in nlsModel(formula, mf, start, wts) :
  singular gradient matrix at initial parameter estimates


I will be highly thankful if some one can please let me know where is the
mistake as i m unable to trace it.

Thanks
Ruchita

On Wed, Dec 16, 2009 at 3:18 PM, ruchita gupta wrote:

> Hello
>
> After performing NLS estimates for some sigmoid model(like logistic growth
> model  and Gompertz growth models), how can we get the RMSE(root mean square
> error) and MAPE(mean absolute percentage error) in  R statistical tool for
> comparison between two models
>
> Thanks
> Ruchita
>

[[alternative HTML version deleted]]

__
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] arrow plots

[Deepayan Sarkar]

On Wed, Dec 9, 2009 at 10:08 AM, Cable, Samuel B Civ USAF AFMC
AFRL/RVBXI  wrote:
> Thanks, all, for the help.  Much obliged.  I realize now that I should
> have said that I am using lattice graphics.  The par() command has not
> been helpful in convincing lattice to plot outside of the default
> window.  Any other advice is appreciated.  Thanks again.

[Sorry for the late reply, I missed your mail the first time around.]

I think you need the 'legend' functionality. A naive example is:

arrowFun <- function(x = c(0.25, 0.75), ...)
{
require(grid)
linesGrob(x = x, y = 0.5, default.units = "npc",
  arrow = arrow(...))
}

xyplot(1 ~ 1, legend = list(top = list(fun = arrowFun, args = list(x =
c(0.1, 0.5)
xyplot(1 ~ 1, legend = list(top = list(fun = arrowFun(type = "closed"


Your actual requirements would be more complicated because (1) your
grob will be more complicated (but see draw.key and draw.colorkey for
examples) and (2) your legend has to get resized in sync with the
plot. The second part is more complicated, and the simplest approach
would be to fix the absolute size of your panels (see ?print.trellis).
Of course the easiest solution would be to draw your arrow inside the
panel.

-Deepayan

__
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] write.csv and col.names=F

[Gabor Grothendieck]

Use write.table instead.

On Thu, Dec 17, 2009 at 6:52 PM, kayj  wrote:
>
> Hi All,
>
> I always have a problem with write.csv when I want the column names to be
> ignored, when I specify col.names=F, I get a header of V1 V2 V3 V4 etc.
>
> for example I tried
>
> write.csv(mydata, file="data.csv", quote=FALSE, row.names=F, col.names=F)
> Warning message:
> In write.csv(mydata, file = "data.csv", quote = FALSE,  :
>  attempt to set 'col.names' ignored
>>
>
> how can I get rid of this problem ? I do not want the header of V1 V2 V3...
> to appear in my output file,
>
> Thanks
> --
> View this message in context: 
> http://n4.nabble.com/write-csv-and-col-names-F-tp974477p974477.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> 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-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] Getting Rd pages right for redefined S3 generic

[Ken Knoblauch]

S Ellison  lgc.co.uk> writes:
> I wanted to define a cbind equivalent for an object that mostly behaves
> like a data frame. base::cbind dispatches to a data frame method if
> _any_ parameter is a data frame, so I defined a new S3 cbind and
> cbind.default to handle dispatch on first object only. Though I confess
> that redefining cbind leaves me a tad nervous, that all works OK so far.
 clipped ---
> Steve Ellison

Why don't you just define a method for your object class.  It should
dispatch on your method if all of the objects that you give it have
that class.  I did this in the MLDS package in a situation in which
I had a data frame that had additional attributes and I wanted
rbind (in my case) to concatenate the attributes as well as the
actual data frame.  I gave the data frame with these extra attributes
the class of c("mld.df", "data.frame") so that they would still
inherit other data frame methods.  My rbind.mlds.df works fine
with them, and I document it accordingly.

HTH.

Ken

-- 
Ken Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.html

__
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] Doubts about ADE-4

[Diogo B. Provete]

Hi,
I'm having some troubles in using the ADE-4 package. I want to analyze data
for the Outlying Mean Index in the function 'niche'. I've got two matrixes:
one for environmental data and the other for species abundance data. But to
implement the analysis, as I understand it, I need to put tocheger these two
matrixes in one, as in the example provided with the package. Un fortunately
I don't know how to do this. Could anyone help me?

Thank you in advance.

-- 
Respectfully,
Diogo Borges Provete

Biólogo
Mestrando - PPG Biologia Animal
Laboratório de Ecologia Animal
Departamento de Zoologia e Botânica
Instituto de Biociências, Letras e Ciências Exatas
Universidade Estadual Paulista - UNESP
São José do Rio Preto-SP
Brazil

Rua Cristóvão Colombo, 2265
Jardim Nazareth -  15054-000

Skype: diogoprovete
MSN: diogop...@yahoo.com.br

[[alternative HTML version deleted]]

__
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] Getting Rd pages right for redefined S3 generic

[Duncan Murdoch]

On 17/12/2009 11:08 PM, S Ellison wrote:

I'm writing a package, and would appreciate advice on controlling the
help documentation cross-references for a redefined generic.

I wanted to define a cbind equivalent for an object that mostly behaves
like a data frame. base::cbind dispatches to a data frame method if
_any_ parameter is a data frame, so I defined a new S3 cbind and
cbind.default to handle dispatch on first object only. Though I confess
that redefining cbind leaves me a tad nervous, that all works OK so far.

However, R cmd tells me (rightly) that I haven't documented the generic
and new default. But I don't want to; I want ?cbind to point to the base
package help pages, not to mine.


But then the documentation will be incorrect, since it mentions the 
non-standard handling of dataframes, which you've done away with.



Assuming that's not unwise (?), what do I have to do to tell R cmd that
it should not look for Rd docs for cbind and cbind.default in my package
whilst still having my generic handle dispatch correctly?

It looks like one answer _might_ be not to export my redefined cbind  -
but will a local generic still work properly if it's not exported?


I would say that what you've done is unwise.  I would rename the 
function to avoid future confusion, especially if you decide to export 
it.  Even if it's for internal use only, won't you want a reminder in 2 
years time that cbind(a,b) doesn't mean what the documentation says it 
means?


Duncan Murdoch

__
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] forecasting

[DispersionMap]

Hi thanks.

Another question...

When R prints the regression output onto the screen after i perform my
regression:

output_dat <- lm(
results <- summary(output_dat)
results 

How do i print the regression output to a file for inclusion into a
document.

I tried:  cat(summary(output_dat), file="results.txt"))

but this doesnt work and it tells me type list (in the first argument)
cannot be handled by cat.


So how do i get the output of lm to the desktopas i need to put this in
a document (in atble would be even better).
 


Stephan Kolassa wrote:
> 
> ?predict
> 
> HTH
> Stephan
> 
> DispersionMap schrieb:
>> 
>> I have some data that i ran a regression on and got the usual r output,
>> with
>> the a intercept and b coefficient for my independent variable.
>> 
>> i want to forecast the number of future events using these parameters.
>> 
>> 
>> What function / packages in R would you recommend to be used in good
>> effect
>> for these purposes??
> 
> __
> 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.
> 
> 

-- 
View this message in context: 
http://n4.nabble.com/forecasting-tp974326p974704.html
Sent from the R help mailing list archive at Nabble.com.

__
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] Ryacas: problem with "Sym"

[Gabor Grothendieck]

Read the troubleshooting section on the home page:
http://ryacas.googlecode.com .  Some common problems are that you did
not install yacas itself (issue yacasInstall() command without args
and it will do that for you) and that you are using the wrong version
of the XML package.

> library(Ryacas)
Loading required package: XML
> a <- Sym("a")
> a
[1] "Starting Yacas!"
expression(a)


On Thu, Dec 17, 2009 at 2:06 PM, Gavin Steininger  wrote:
> I have been trying to get the Ryacas package to work without any
> success. I do not seem to be able to the "Sym" function work, for
> example:
>
>> library(XML)
>> library(Ryacas)
>> a=Sym("a")
>> a
> Error in summary.connection(x) : invalid connection
>
> I am windows XP  with R 2.10
>
> Thanks for you help
>
> gavin amw steininger
>
> __
> 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-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] mutlidimensional in.convex.hull (wasmultidimensionalpoint.in.polygon??)

[baptiste auguie]

Hi,

Excellent, thanks for doing this!

I had tried the 2D case myself but I was put off by the fact that
Octave's convhulln had a different ordering of the points to R's
geometry package (an improvement to Octave's convhulln was made after
it was ported to R). I'm not sure how you got around this but it's
good to see that the result was worth the effort!

Below are a few minor syntax suggestions,

# p <- dim(calpts)[2]   # columns in calpts
# and other similar lines could be replaced with
ncol(calpts)
nrow(testpts)
nrow(hull)

# length(degenflag[degenflag])
# can probably be written
sum(degenflag)

# center = apply(calpts, 2, mean)
# more efficient
colMeans(calpts)

Would you consider submitting this function to the maintainer of the
geometry package, after checking it's OK with inhull's original
author?

Best regards,

baptiste




2009/12/18 Keith Jewell :
> Hi All,
>
> I couldn't resist doing this the "right" way!
> A colleague explained the vector algebra to me (thanks Martin!) and I've
> followed the structure of the Matlab code referenced below, but all errors
> are mine!
>
> I don't pretend to great R expertise, so this code may not be optimal (in
> time or the memory issue addressed by the Matlab code), but it is a lot
> faster than the other algorithm discussed in this thread. These timings are
> on the real example data described in my previous message in this thread.
>
>> system.time(inhull(xs,ps))  # the "right" way
>   user  system elapsed
>   1.34    0.07    1.41
>> system.time({phull <-  convhulln(ps) # the other algorithm
> + phull2 <- convhulln(rbind(ps,xs))
> + nrp <- nrow(ps)
> + nrx <- nrow(xs)
> + outside <- unique(phull2[phull2>nrp])-nrp
> + done <- FALSE
> + while(!done){
> +     phull3 <- convhulln(rbind(ps,xs[-(outside),]))
> +     also.outside <- (1:nrx)[-outside][unique(phull3[phull3>nrp])-nrp]
> +     outside <- c(outside,also.outside)
> +     done <- length(also.outside)==0
> + }})
>   user  system elapsed
>  15.09    0.09   15.22
>
> Although I really must move on now, if anyone has comments, criticisms,
> suggestions for improvement I would be interested.
>
> Enjoy!
>
> Keith Jewell
> 
> inhull <- function(testpts, calpts, hull=convhulln(calpts),
> tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps)) {
> #
> # R implementation of the Matlab code by John D'Errico 04 Mar 2006 (Updated
> 30 Oct 2006)
> # downloaded from
> http://www.mathworks.com/matlabcentral/fileexchange/10226-inhull
> # with some modifications and simplifications
> #
> # Efficient test for points inside a convex hull in n dimensions
> #
> #% arguments: (input)
> #%  testpts - nxp array to test, n data points, in p dimensions
> #%       If you have many points to test, it is most efficient to
> #%       call this function once with the entire set.
> #%
> #%  calpts - mxp array of vertices of the convex hull, as used by
> #%       convhulln.
> #%
> #%  hull - (OPTIONAL) tessellation (or triangulation) generated by convhulln
> #%       If hull is left empty or not supplied, then it will be
> #%       generated.
> #%
> #%  tol - (OPTIONAL) tolerance on the tests for inclusion in the
> #%       convex hull. You can think of tol as the distance a point
> #%       may possibly lie outside the hull, and still be perceived
> #%       as on the surface of the hull. Because of numerical slop
> #%       nothing can ever be done exactly here. I might guess a
> #%       semi-intelligent value of tol to be
> #%
> #%         tol = 1.e-13*mean(abs(calpts(:)))
> #%
> #%       In higher dimensions, the numerical issues of floating
> #%       point arithmetic will probably suggest a larger value
> #%       of tol.
> #%
> # In this R implementation default
> tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps)
> #       DEFAULT: tol = 1e-6
> #
> # VALUE: Matlab returns a vector of TRUE (inside/on) or FALSE (outside)
> #       This R implementation returns an integer vector of length n
> #       1 = inside hull
> #      -1 = inside hull
> #       0 = on hull (to precision indicated by tol)
> #
>   require(geometry, quietly=TRUE)  # for  convhulln
>   require(MASS, quietly=TRUE)      # for Null
> # ensure arguments are matrices (not data frames) and get sizes
>   calpts <- as.matrix(calpts)
>   testpts <- as.matrix(testpts)
>   p <- dim(calpts)[2]   # columns in calpts
>   cx <- dim(testpts)[1]  # rows in testpts
>   nt <- dim(hull)[1]    # number of simplexes in hull
> # find normal vectors to each simplex
>   nrmls <- matrix(NA, nt, p)         # predefine each nrml as NA,
> degenerate
>   degenflag <- matrix(TRUE, nt, 1)
>   for (i in  1:nt) {
>    nullsp <- t(Null(t(calpts[hull[i,-1],] -
> matrix(calpts[hull[i,1],],p-1,p, byrow=TRUE
>    if (dim(nullsp)[1] == 1) { nrmls[i,] <- nullsp
>       degenflag[i] <- FALSE}}
> # Warn of degenerate faces, and remove corresponding normals
>   if(length(degenflag[degenflag]) > 0)
> warning(len

Re: [R] Testing equality of regression model on multiple groups

[Daniel Malter]

Hi, your question is unclear. It is not clear whether you want to compare
a.) whether two subsets of data have the same coefficients for an identical
model or b.) whether a couple of coefficients are different once you include
additional regressors. The reason for this confusion is that your
lm.separate and lm.together are not the same models (i.e., they do not
include the same regressors). The former is (y~x1*x2); the latter is
(y~t*x1*x2), which is obviously different. 

If a.) is your goal, you want to run a Chow test.

For that, you run the regression (y~x1*x2) on the entire dataset and on the
two subsets of the data separately. If the model on the entire data produces
much larger errors (in sum) than the two regressions on the subsets of data
(in sum), then this is evidence that the data-generating processes for the
two subsets differ significantly. In other words, the coefficients for the
subsets will jointly be significantly different from each other. This is, as
you said, an F-Test. Look here for further details:

http://en.wikipedia.org/wiki/Chow_test

It just beckons me that you probably want to compare multiple groups. If you
want to run (y~x1*x2) interacted with the group indicator, i.e.
(y~group*x1*x2), you should code your group indicator as "treatment
contrasts."

The t.L, t.Q, and t.C indicate a.) that t is a factor variable and b.) that
this factor variable is coded as orthogonal polynomial contrasts. This tests
for linear, quadractic, cubic, etc. influence of t on the dependent
variable. This makes sense if t has constant differences between the factor
levels (like 10 cm, 15 cm, 20 cm, 25 cm). Otherwise, you probably want
dummy-variable coding, which is called "treatment contrasts" in R lingo.

HTH,
Daniel

-
cuncta stricte discussurus
-
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Clara Yuan
Sent: Thursday, December 17, 2009 5:15 PM
To: r-h...@lists.r-project.org
Subject: [R] Testing equality of regression model on multiple groups

Hello,

I'm trying to test for the joint equality of coefficients of the same model
across different subsets of the data (ie, is it necessary to estimate the
same model on these different populations, or can I just estimate the model
once on the whole dataset?).

My plan is to use the F-test on the reduced model and the full model. By
full model, I mean a specification that mimics my regressions on separate
subsets of data, but I have found that the full model's coefficient
estimates don't correspond to my original model's estimates. I was under the
impression that they would be identical.

Original model:
> lm.separate = by(data.ex, data.ex$t, function(x) lm(y ~ x1 * x2, data =
x))

Full model:
> lm.together = lm(y ~ t * x1 * x2, data = data.ex)

The data are grouped by t.

When I examine the coefficients, I find that they are roughly in the same
ballpark, but not nearly identical:
> sapply(lm.separate, coef)
  12 3 4
(Intercept) 2.691272263 1.7153565472  1.8797914048  1.9282332240
x1  0.107520721 0.0472488208  0.0440171489  0.0198376096
x2  0.054694784 0.0396246366  0.0603665574  0.0300886164
x1:x2   0.002180749 0.0003653858 -0.0001488267 -0.0007409421

> coef(lm.together)
  (Intercept)   t.L   t.Q   t.Cx1
 2.0536633597 -0.4750933962  0.5121787674 -0.2809269719  0.0546560750
   x2t.L:x1t.Q:x1t.C:x1t.L:x2
 0.0461936485 -0.0595422428  0.0180461803 -0.0174386682 -0.0118682844
   t.Q:x2t.C:x2 x1:x2 t.L:x1:x2 t.Q:x1:x2
-0.0076038969 -0.0194162097  0.0004140914 -0.0020749112  0.0006116237
t.C:x1:x2
-0.0003083657

(Also, why are the coefficients renamed to t.L, t.Q, etc instead of t.1,
t.2?)

What am I missing?

Thanks for the help,
Clara

__
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-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] mutlidimensional in.convex.hull (wasmultidimensionalpoint.in.polygon??)

[Keith Jewell]

Hi All,

I couldn't resist doing this the "right" way!
A colleague explained the vector algebra to me (thanks Martin!) and I've 
followed the structure of the Matlab code referenced below, but all errors 
are mine!

I don't pretend to great R expertise, so this code may not be optimal (in 
time or the memory issue addressed by the Matlab code), but it is a lot 
faster than the other algorithm discussed in this thread. These timings are 
on the real example data described in my previous message in this thread.

> system.time(inhull(xs,ps))  # the "right" way
   user  system elapsed
   1.340.071.41
> system.time({phull <-  convhulln(ps) # the other algorithm
+ phull2 <- convhulln(rbind(ps,xs))
+ nrp <- nrow(ps)
+ nrx <- nrow(xs)
+ outside <- unique(phull2[phull2>nrp])-nrp
+ done <- FALSE
+ while(!done){
+ phull3 <- convhulln(rbind(ps,xs[-(outside),]))
+ also.outside <- (1:nrx)[-outside][unique(phull3[phull3>nrp])-nrp]
+ outside <- c(outside,also.outside)
+ done <- length(also.outside)==0
+ }})
   user  system elapsed
  15.090.09   15.22

Although I really must move on now, if anyone has comments, criticisms, 
suggestions for improvement I would be interested.

Enjoy!

Keith Jewell

inhull <- function(testpts, calpts, hull=convhulln(calpts), 
tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps)) {
#
# R implementation of the Matlab code by John D'Errico 04 Mar 2006 (Updated 
30 Oct 2006)
# downloaded from 
http://www.mathworks.com/matlabcentral/fileexchange/10226-inhull
# with some modifications and simplifications
#
# Efficient test for points inside a convex hull in n dimensions
#
#% arguments: (input)
#%  testpts - nxp array to test, n data points, in p dimensions
#%   If you have many points to test, it is most efficient to
#%   call this function once with the entire set.
#%
#%  calpts - mxp array of vertices of the convex hull, as used by
#%   convhulln.
#%
#%  hull - (OPTIONAL) tessellation (or triangulation) generated by convhulln
#%   If hull is left empty or not supplied, then it will be
#%   generated.
#%
#%  tol - (OPTIONAL) tolerance on the tests for inclusion in the
#%   convex hull. You can think of tol as the distance a point
#%   may possibly lie outside the hull, and still be perceived
#%   as on the surface of the hull. Because of numerical slop
#%   nothing can ever be done exactly here. I might guess a
#%   semi-intelligent value of tol to be
#%
#% tol = 1.e-13*mean(abs(calpts(:)))
#%
#%   In higher dimensions, the numerical issues of floating
#%   point arithmetic will probably suggest a larger value
#%   of tol.
#%
# In this R implementation default 
tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps)
#   DEFAULT: tol = 1e-6
#
# VALUE: Matlab returns a vector of TRUE (inside/on) or FALSE (outside)
#   This R implementation returns an integer vector of length n
#   1 = inside hull
#  -1 = inside hull
#   0 = on hull (to precision indicated by tol)
#
   require(geometry, quietly=TRUE)  # for  convhulln
   require(MASS, quietly=TRUE)  # for Null
# ensure arguments are matrices (not data frames) and get sizes
   calpts <- as.matrix(calpts)
   testpts <- as.matrix(testpts)
   p <- dim(calpts)[2]   # columns in calpts
   cx <- dim(testpts)[1]  # rows in testpts
   nt <- dim(hull)[1]# number of simplexes in hull
# find normal vectors to each simplex
   nrmls <- matrix(NA, nt, p) # predefine each nrml as NA, 
degenerate
   degenflag <- matrix(TRUE, nt, 1)
   for (i in  1:nt) {
nullsp <- t(Null(t(calpts[hull[i,-1],] - 
matrix(calpts[hull[i,1],],p-1,p, byrow=TRUE
if (dim(nullsp)[1] == 1) { nrmls[i,] <- nullsp
   degenflag[i] <- FALSE}}
# Warn of degenerate faces, and remove corresponding normals
   if(length(degenflag[degenflag]) > 0) 
warning(length(degenflag[degenflag])," degenerate faces in convex hull")
   nrmls <- nrmls[!degenflag,]
   nt <- dim(nrmls)[1]
# find center point in hull, and any (1st) point in the plane of each 
simplex
   center = apply(calpts, 2, mean)
   a <- calpts[hull[!degenflag,1],]
# scale normal vectors to unit length and ensure pointing inwards
   nrmls <- nrmls/matrix(apply(nrmls, 1, function(x) sqrt(sum(x^2))), nt, p)
   dp <- sign(apply((matrix(center, nt, p, byrow=TRUE)-a) * nrmls, 1, sum))
   nrmls <- nrmls*matrix(dp, nt, p)
# if  min across all faces of dot((x - a),nrml) is
#  +ve then x is inside hull
#  0   then x is on hull
#  -ve then x is outside hull
# Instead of dot((x - a),nrml)  use dot(x,nrml) - dot(a, nrml)
   aN <- diag(a %*% t(nrmls))
   val <- apply(testpts %*% t(nrmls) - matrix(aN, cx, nt, byrow=TRUE), 
1,min)
# code  values inside 'tol' to zero, return sign as integer
   val[abs(val) < tol] <- 0
   as.integer(sign(val))
}

__
R-help@r-project.org mailing list
ht

Re: [R] What is the fastest way to see what are in an RData file?

[Patrick Connolly]

On Thu, 17-Dec-2009 at 03:13PM +1800, Peng Yu wrote:

|> Currently, I load the RData file then ls() and str(). But loading the file
|> takes too long if the file is big. Most of the time, I only interested what
|> the variables are in the the file and the attributes of the variables (like
|> if it is a data.frame, matrix, what are the colnames/rownames, etc.)
|> 
|> I'm wondering if there is any facility in R to help me avoid loading the
|> whole file.


I have a pretty nifty way of seeing what's in such a file, but I still
have to load all of the binary file before I can do so.  If it's
taking you such a long time, maybe you could keep a larger number of
smaller RData files.


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

__
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] mantel test and NAs

[Tal Galili]

Could you please supply with a reproducible self sufficient script ?


Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com/ (English)
--




On Thu, Dec 17, 2009 at 7:33 PM, marciarocha  wrote:

>
> Dear R users,
>
> I am having a problem performing the mantel test (both with functions
> mantel
> {vegan} and mantel.test {ape}) due to I believe the presence of NAs on my
> distance matrices, which look like e.g.:
>
> NA
> 12
> 12   3
> NA  4   5  6
>
> and
>
> 1
> 1   2
> 2   3   4
> 5   6   7  8
>
> Would any of you have a solution for that?
>
> Thank you for much for your help!
>
> Marcia Rocha
>
>
> --
> View this message in context:
> http://n4.nabble.com/mantel-test-and-NAs-tp966682p966682.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> 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.
>

[[alternative HTML version deleted]]

__
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] What is the fastest way to see what are in an RData file?

[Gustaf Rydevik]

On Thu, Dec 17, 2009 at 4:33 PM, Peng Yu  wrote:
> On Thu, Dec 17, 2009 at 5:33 AM, Gustaf Rydevik
>  wrote:
>> On Wed, Dec 16, 2009 at 10:13 PM, Peng Yu  wrote:
>>>
>>> Currently, I load the RData file then ls() and str(). But loading the file
>>> takes too long if the file is big. Most of the time, I only interested what
>>> the variables are in the the file and the attributes of the variables (like
>>> if it is a data.frame, matrix, what are the colnames/rownames, etc.)
>>>
>>> I'm wondering if there is any facility in R to help me avoid loading the
>>> whole file.
>>
>>
>> I thought this was interesting as well, so i did a bit of searching
>> through the R-help list archives and found this answer by Simon
>> Urbanek:
>> https://stat.ethz.ch/pipermail/r-devel/2007-August/046724.html
>> The link to a c-routine that does what you want still works, but for
>> future reference I'm pasting the code below.
>
> It doesn't work for the RData file that I saved by "save(list='test',
> file='test.RData')".
>
> $ rdcopy test.RData
> Format version 3ec, R version = 23813.88.84, release = f9db1dba
> Sorry, this tool supported RXDR version 2 format only
>


What happens if you remove the version check?
I.e. this one:
  if (ver != 2) {
XdrInTerm(d);
error(_("Sorry, this tool supported RXDR version 2 format only\n"));
  }


>From what I can read on the hel page for ?save, there hasn't been a
change in the file format since 1.4.0


/Gustaf
-- 
Gustaf Rydevik, M.Sci.
tel: +46(0)703 051 451
address:Essingetorget 40,112 66 Stockholm, SE
skype:gustaf_rydevik

__
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] Plot with a factor variable which has large values

[Jim Lemon]

On 12/18/2009 06:42 AM, changzhu wrote:

I wanted to do a simple plot: Death Rate vs Clinical Site, and then draw a
confidence interval of Death Rate at each clinical site, which is a factor
variable with 100 levels. I used the following code:

plot(Site, Rate, col="red", type="o")

However, the graph was a short line at each site instead of points, and the
color was not red at all. Seem to me for factor variable all options are not
working ! Could anybody help me on how  to draw points at each site, with
specific color? if possible, how to further draw confidence interval at each
site?
   

Hi Changzu,
Since I don't have Site or Rate, I'll have to be imaginative:

Site<-paste(sample(LETTERS,100,TRUE),
 sample(letters,100,TRUE),sep="")
SiteOrder<-order(Site)
Rate<-runif(100)
plot(as.numeric(factor(Site))[SiteOrder],
 Rate[SiteOrder],col="red",type="o")

Having performed that chicanery, I feel justified in telling you to use 
the dispersion function in the plotrix package for the CIs, although 
there are many other similar functions.


Jim

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