[R] Extraction of rules from Support Vector Machines

2011-12-15 Thread Paul Smith
Dear All,

I am using the svm() function of package e1071 for creating Support
Vector Machines prediction models. As far as I understand, there is no
function in this package to extract rules of prediction. Is there some
other package with such a functionality?

Thanks in advance,

Paul

__
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] optim with simulated annealing SANN for combinatorial optimization

2011-12-15 Thread Dae-Jin Lee
Hi all

I am trying to solve a combinatorial optimization problem. Basically, I can
reduce my problem into the next problem:

1.- Given a NxN grid of points, with some values in each cell
2.- Find the combination of K points on the grid such that, the maximum
mean value is obtained


 I took the Travel SalesMan problem example in ?optim documentation. I am
not sure if I have understood correctly the SANN implementation in optim,
as I do not see how the acceptance probability comes out. And it looks like
I am only evaluating the criteria several times and keep the maximum at the
end.

Thanks in advance


Here is some example code in R

### toy example
N=5
K=6

new.points = expand.grid(1:N,1:N)  # grid points

set.seed(1234)

resp=rnorm(N^2)  # random values on each grid cell

###  function to generate the sequence of candidates
generate<-function(ind){

idx <- seq(1, nrow(new.points), by=1)   # index of 1 to N^2 grid cells
swap.points <- c(sample(ind,1),sample(idx[-c(ind)], size=1))  # swap
between points of the initial set and other candidate
ind[which(ind==swap.points[1])]<-idx[which(idx==swap.points[2])]

cat("set  =",ind,"\n")
cat("crit =",media(ind),"\n")
cat("swap =",swap.points,"\n")

plot(new.points[,1:2],col='black',xlim=c(range(new.points[,1])),ylim=c(range(new.points[,2])))
points(new.points[ind,1:2],col='blue',pch=19,xlim=c(range(new.points[,1])),ylim=c(range(new.points[,2])))
return(as.numeric(ind))
}

###  function to calculate the mean value of the points on the grid
media=function(sq){
med=mean(resp[sq])
return(as.numeric(med))
}

###  initial set of K candidates
init=sample(1:K,K)

###  run SANN
res <- optim(init, media ,generate, method="SANN",control =
list(maxit=2, temp=100,tmax=1000, trace=TRUE, REPORT=1, fnscale=-1))
new.points[sort(res$par),]
plot(new.points,cex=.1)
points(new.points[res$par,],col=3,lwd=3,cex=1.5)

[[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] Odp: Problem with package compilation

2011-12-15 Thread Petr PIKAL
Hi

> 
> Hi,
> 
> I have a R package with some functions made all of then only with R 
> code. I use the command R CMD build to build a package that I can 
> install on linux, windows or mac, because all the code is only R code. 
> But I have some problems with R version. For each new R version I need 
> to rebuild the package to be install in this new version. It is possible 

> to make a package without this R Version dependence? I know that in a 
> little cases I need to adapt my package when R change some of this basic 

> functions, but it is very rare.
> 
> Any idea? Is only a tag to put on my package? Or this is the correct 
> way, maintain the R version dependence?

Not sure if this is "correct" way but on windows I usually only copy my 
package from one R version to another.

Regards
Petr

> 
> Thanks
> Ronaldo
> 
> -- 
> 6Ş lei - Sua produtividade varia como
>   (tempo produtivo gasto por dia)^1.000.
> 
> 7Ş lei - Sua produtividade também varia como
>   1/(seu atraso na análise dos dados obtidos)^1.000.
> 
>--Herman, I. P. 2007. Following the law. NATURE, Vol 445, p. 228.
> 
> >  Prof. Ronaldo Reis Júnior
> |  .''`. UNIMONTES/DBG/Lab. Ecologia Comportamental e Computacional
> | : :'  : Campus Universitário Prof. Darcy Ribeiro, Vila Mauricéia
> | `. `'` CP: 126, CEP: 39401-089, Montes Claros - MG - Brasil
> |   `- Fone: (38) 3229-8192 | ronaldo.r...@unimontes.br
> | http://www.ppgcb.unimontes.br/lecc | LinuxUser#: 205366
> 
> 
>[[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] modify the name of axis of an R function

2011-12-15 Thread Petr PIKAL
Hi

> 
> Thanks you very much!  The way plot(fit, main="main title", xlab="X-axis
> lable", ylab="y-axis label") seems to work quite well, I didn't notice 
that
> I could do this. 
> 
> I have in fact one more problem with it : the fact is that I have three
> plots that are called by the function. I can specify my three titles by
> doing "main=c("title1","title2","title3")" whithout any problem. But
> proceeding this way to specify the axes' titles doesn't work anymore. 
How
> can I specify three differents pairs of titles for these axes?

I am rather confused 

plot(1,1, main=c("title1", "title2"), xlab=c("x1","x2"), 
ylab=c("y1","y2"))

puts 2 titles and two axes labels into one plot but is it really what you 
want?

Regards
Petr 



> 
> Thank you again
> 
> --
> View this message in context: http://r.789695.n4.nabble.com/modify-the-
> name-of-axis-of-an-R-function-tp4198804p4200631.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] lattice key in blank panel

2011-12-15 Thread Deepayan Sarkar
On Fri, Dec 16, 2011 at 1:22 AM, Max Kuhn  wrote:
> Somewhere I've seen an example of an xyplot() where the key was placed
> in a location of a missing panel. For example, if there were 3
> conditioning levels, the panel grid would look like:
>
> 34
> 12
>
> In this (possibly imaginary) example, there were scatter plots in
> locations 1:3 and location 4 had no conditioning bar at the top, only
> the key.
>
> I can find examples of putting the legend outside of the panel
> locations (e.g to the right of locations 2 and 4 above), but that's
> not really what I'd like to do.

You are probably thinking of this example from ?splom:

splom(~iris[1:3]|Species, data = iris,
  layout=c(2,2), pscales = 0,
  varnames = c("Sepal\nLength", "Sepal\nWidth", "Petal\nLength"),
  page = function(...) {
  ltext(x = seq(.6, .8, length.out = 4),
y = seq(.9, .6, length.out = 4),
labels = c("Three", "Varieties", "of", "Iris"),
cex = 2)
  })

It's actually easier to do that with legends, as illustrated by this
(somewhat silly) modification:

splom(~iris[1:3]|Species, data = iris, groups = Species,
  layout=c(2,2), pscales = 0,
  varnames = c("Sepal\nLength", "Sepal\nWidth", "Petal\nLength"),
  auto.key = list(x = 0.75, y = 0.75, corner = c(0.5, 0.5)))

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


[R] xyplots with differentiated first data points

2011-12-15 Thread Russell Wyeth

Hi,

I'm trying to plot a series of lines without data point markers, except 
for the first data point in each line, which I want to also mark with an 
open circle.


The following code accomplishes this for a single line:

xyplot(yy ~ xx,
 panel=function(x, y){
panel.xyplot(x, y, col="black",type="l")
panel.xyplot(x[1], y[1], col="black",type="p")
})

However, I am trying to do this for every line for multiple lines in a 
panel, with several panels created by a 'by()' function (I do not want a 
grid of panels, I want separated panels).


expt1 is a dataframe with seven factor columns, and x,y data coordinates 
in columns 8 and 9.


I can appropriately create the grouped lines  (but without the single 
data markers) with following code:


by(expt1,expt1[,3:4],function(q) xyplot(q[,9] ~ q[,8], groups=q[,2], 
data=q, col="black", type="l"))


But I am stuck putting the two goals together because I can't figure out 
how to get the group information to be properly passed to the custom 
panel function above.


I suspect I may be going about this entirely the wrong way - I've been 
using R for all of 4 hours...


Help appreciated.

Russell


--

Russell Wyeth
Biology, St Francis Xavier University
P.O. Box 5000 Antigonish, NS B2G 2W5 Canada
Shipping: 1 West St. Antigonish, NS B2G 2W5 Canada
http://people.stfx.ca/rwyeth/
Ph: 9028673886 Fx: 9028672389
Cell: 9023180250

__
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] .nc files query

2011-12-15 Thread vaish
Hello Everyone, I have tried using open.ncdf(con,..) to open .nc files in R,
but somehow, its showing that R could not find function open.ncdf. I am new
to R, please help me out with this

Thanks

--
View this message in context: 
http://r.789695.n4.nabble.com/nc-files-query-tp4203048p4203048.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] fundamental guide to use of numerical optimizers?

2011-12-15 Thread cberry
Paul Johnson  writes:

> I was in a presentation of optimizations fitted with both MPlus and
> SAS yesterday.  In a batch of 1000 bootstrap samples, between 300 and
> 400 of the estimations did not converge.  The authors spoke as if this
> were the ordinary cost of doing business, and pointed to some
> publications in which the nonconvergence rate was as high or higher.
>
> I just don't believe that's right, and if some problem is posed so
> that the estimate is not obtained in such a large sample of
> applications, it either means the problem is badly asked or badly
> answered.  But I've got no traction unless I can actually do
> better

A few years back there was a brouhaha in which a too lax convergence
criterion in the Splus gam() function resulted in wrong results. 

See

http://www.ihapss.jhsph.edu/publications/Results/nmmaps_faq.htm

I think this was also reported in the lay press.

IIRC, at that time there was an assertion that gam() was buggy, but it
turned out that for the particular problem a more stringent tolerance
was needed than the default provided. The original report used results
that hadn't actually converged.

 The trouble is there are many instances of monkey-see, monkey-do
data analysis. It seems that some authors do not really want to dig into
their data if the story it tells is not simple and firmly supported. And
not understanding why many bootstrap samples do not converge seems like
an instance of sweeping data-dirt under the rug.

The questions you ask below full under the rubric of 'numerical
analysis'. You might look here to start:

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

Chuck

>
> Perhaps I can use this opportunity to learn about R functions like
> optim, or perhaps maxLik.
>
>>From reading r-help, it seems to me there are some basic tips for
> optimization, such as:
>
> 1. It is wise to scale the data so that all columns have the same
> range before running an optimizer.
>
> 2. With estimates of variance parameters, don't try to estimate sigma
> directly, instead estimate log(sigma) because that puts the domain of
> the solution onto the real number line.
>
> 3 With estimates of proportions, estimate instead the logit, for the
> same reason.
>
> Are these mistaken generalizations?  Are there other tips that
> everybody ought to know?
>
> I understand this is a vague question, perhaps the answers are just in
> the folklore. But if somebody has written them out, I would be glad to
> know.

-- 
Charles C. BerryDept of Family/Preventive Medicine
ccberry at ucsd edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

__
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] Incorrect Number of Dimensions in Zelig with setx()

2011-12-15 Thread Abraham Mathew
I'm running an ordered logit in R with the Zelig package and am trying to
calculate some predicted probabilities. However, I get the following error
message.

> x.low <- setx(mod, cars=1)Error in dta[complete.cases(mf), names(dta) %in% 
> vars, drop = FALSE] :
  incorrect number of dimensions


I googled this problem and couldn't find anything, minus a question by

me on this same problem from 1.5 years ago. Just don't remember what I

did to solve the problem.


Help!


-- 
*Abraham Mathew
Statistical Analyst
www.amathew.com
720-648-0108
@abmathewks*

[[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] (no subject)

2011-12-15 Thread Val
Hi all,

I am trying to extract the residuals from the following Weibull model..
The results I am getting are   a bit strange.

library(survival)
mod=survreg(Surv(time, delta)  ~  p1+p2+p3,  data=testd,  dist="weibull")
ores=resid(mod, type='response')

Are the commands correct? How do I get teh standardized residuals too?

Thanks
Val

[[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] Multicollinearty in logistic regression models

2011-12-15 Thread cberry
David Winsemius  writes:

> On Dec 15, 2011, at 11:34 AM, Mohamed Lajnef wrote:
>
>> Dear All,
>>
>> Is there a method to diagnostic multicollinearty in logistic
>> regression
>> models  like vif indicator in linear regression ( variance inflation
>> Factor ...) ?
>>
>
> Wouldn't matrix representation of the predictor "side" of the
> regression be the same? Couldn't you just use the same methods you
> employ for linear regression?

Trouble is that in logistic regression the Fisher Information for each
case has a factor of p[i]*(1-p[i]) (where 'p' is the vector of success
probabilites and 'i' indexes which case).

If the value of p[i] is very near one or zero, then the information
provided is scant. And this will happen if you have a really good
predictor in the mix.

Even with an orthogonal design, you can wind up with huge variances. And
you can have an ill-conditioned var-cov matrix for the coefficients
depending on how different cases get weighted. Thus, you could get the
equivalent of multicollinearity even with an orthogonal design.

And the diagnostics for linear regresson would not be all
that helpful if you have a good predictor.

OTOH, if the predictors were collectively pretty weak, the linear
regression diagnostics might be OK.

Mu advice: Google Scholar 'pregibon logistic regression', click where it
says 'cited by ...' and page through the results to find good leads on
this topic.

HTH,

Chuck

>
>> Thank you in advance
>> M
>>
>> --
>> 
>> Mohamed Lajnef,IE INSERM U955 eq 15#
>
>
> David Winsemius, MD
> West Hartford, CT
>

-- 
Charles C. BerryDept of Family/Preventive Medicine
cberry at ucsd  edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

__
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] MAximum location

2011-12-15 Thread Milan Bouchet-Valat
Le jeudi 15 décembre 2011 à 21:15 +0100, Trying To learn again a écrit :
> Hi all,
> 
> I have a matrix
> a<-c(2,3,4,Inf)
> 
> > b<-as.matrix(a)
>  [,1]
> [1,]2
> [2,]3
> [3,]4
> [4,]  Inf
> 
> > range(b, finite=TRUE)[2] (this is the maximum)
> [1] 4
> 
> There is a pre-def function to extract the location (in terms of rows) of
> the value in the matrix.
> 
> In my example would be
> 
> 3 (max is in the third row)
> 
> The maximum is in the position (row) 3.
Maybe using this:
> row(b)[b == range(b, finite=TRUE)[2]]
[1] 3
> col(b)[b == range(b, finite=TRUE)[2]]
[1] 1

Not very short, since in you case involving Inf you cannot use
which.max() directly, but it works.

Regards

__
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] Reordering a numeric variable

2011-12-15 Thread Bert Gunter
I believe this needs to be posted on general statistical help list
like stats.stackexchange.com
as you appear to have some confusion about how linear regression
works, especially with categorical variables.

?ordered

might help you with R's way of dealing with the issues that I believe
you're asking about, but it may not make much sense if I'm correct
about your underlying misunderstanding. If I'm wrong, that and perhaps
?C will probably do it for you.

-- Bert

On Thu, Dec 15, 2011 at 1:44 PM, Abraham Mathew  wrote:
> I'm running a linear model in R using the car package.
>
> I have a variable education, which i have recoded and regrouped to my
> wishes.
> However, R seems to place each element of that variable in alphabetical
> order.
>
> When I am running the model, don't I need the model order from lowest to
> highest to make an inference that
> a one unit change in one variable produced a one unit change in another.
>
> levels(educ)
> educ2 = NA
> educ2[educ %in% levels(educ)[c(4,7)]] <- "HS or Some College"
> educ2[educ %in% levels(educ)[1:2]] <- "College Degree"
> educ2[educ %in% levels(educ)[c(3,5)]] <- "Advanced Degree"
> educ2[educ %in% levels(educ)[c(6,8)]] <- "Other"
> educ2 = factor(educ2)
> levels(educ2)
>
> The above code is how I regrouped the variable. How can I regroup it so
> that it's levels
> are from lowest to highest. What if they're numeric values"
>
> --
> *Abraham Mathew
> Statistical Analyst
> This or That Media, Inc.
> abra...@thisorthat.com
> 720-648-0108
> @abmathewks
> www.amathew.com
> *
>
>        [[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.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

__
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] Reordering a numeric variable

2011-12-15 Thread Abraham Mathew
I'm running a linear model in R using the car package.

I have a variable education, which i have recoded and regrouped to my
wishes.
However, R seems to place each element of that variable in alphabetical
order.

When I am running the model, don't I need the model order from lowest to
highest to make an inference that
a one unit change in one variable produced a one unit change in another.

levels(educ)
educ2 = NA
educ2[educ %in% levels(educ)[c(4,7)]] <- "HS or Some College"
educ2[educ %in% levels(educ)[1:2]] <- "College Degree"
educ2[educ %in% levels(educ)[c(3,5)]] <- "Advanced Degree"
educ2[educ %in% levels(educ)[c(6,8)]] <- "Other"
educ2 = factor(educ2)
levels(educ2)

The above code is how I regrouped the variable. How can I regroup it so
that it's levels
are from lowest to highest. What if they're numeric values"

-- 
*Abraham Mathew
Statistical Analyst
This or That Media, Inc.
abra...@thisorthat.com
720-648-0108
@abmathewks
www.amathew.com
*

[[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] printing all htest class members

2011-12-15 Thread William Dunlap
I think you should subclass "htest" and write a print method
for that subclass.  E.g.,

  > print.htest_rb <- function(x, ...) {
  +NextMethod(x, ...)
  +cat("all.orders =", x$all.orders, "\n")
  +invisible(x)
  + }
  > ht2$all.orders <- "THIS IS THE 'all.orders' COMPONENT"
  > class(ht2) <- c("htest_rb", class(ht2))
  > ht2
  
  Test the 'print.htest' method
  
  data:  x 
  Q = 0.0129, df = 2, p-value = 0.9936
  alternative hypothesis: It doesn't print 'all.orders' 
  null values:
ord df  Q  p
  1   2  2 0.0129 0.9936
  2   3  3 0.0490 0.9972
  3   4  4 0.0684 0.9994
  
  all.orders = THIS IS THE 'all.orders' COMPONENT

If you want the 'all.order' information printed in the
middle of the htest printout then you will have to copy
the code in stats:::print.htest to your print.htest_rb
method and edit it to suit you.

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 Rui Barradas
> Sent: Thursday, December 15, 2011 12:12 PM
> To: r-help@r-project.org
> Subject: Re: [R] printing all htest class members
> 
> 
> Hello,
> 
> Once again, and as simple as possible,
> 
> res <- data.frame(ord=2:4, df=2:4, Q=c(0.0129, 0.049, 0.0684),
> p=c(0.9936, 0.9972, 0.9994))
> 
> ht2<-structure(
> list(statistic=c(Q=res$Q[1]),
> p.value=res$p[1],
> parameter=c(df=res$df[1]),
> alternative="It doesn't print 'all.orders'",
> method="Test the 'print.htest' method",
> data.name=deparse(substitute(x)),
> null.value=res  # this is printed
> #all.orders=res # this wouldn't be if uncommented
> ),
> .Names=c("statistic","p.value","parameter","alternative",
> "method","data.name","null.value"),
> class="htest"
> )
> ht2
> class(ht2)
> 
> Creation? Conversation in my head? Have you heard about the experimental
> method?
> It's a very usefull device created by one Galileo Galilei some time ago.
> 
> 
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/printing-all-htest-class-members-
> tp4200872p4201793.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] removing contractions for recode in car

2011-12-15 Thread Nicole Marie Ford
Dear John and David,

Thanks so much for the advice.  I have only come across this one other time in 
my research, so the code escaped me!  Much obliged!

~Nicole

Ph.D. Student 
University of Wisconsin Milwaukee 
c: 813.786.5715 
e: nmf...@uwm.edu 

- Original Message -
From: "John Fox" 
To: "Nicole Marie Ford" 
Cc: r-help@r-project.org
Sent: Thursday, December 15, 2011 3:16:58 PM
Subject: Re: [R] removing contractions for recode in car

Dear Nicole,

On Thu, 15 Dec 2011 14:49:04 -0600 (CST)
 Nicole Marie Ford  wrote:
> hello,
> 
> i need to recode a variable, however the contraction is causing problems.  i 
> had the code to change this written down somewhere and i just can't find it, 
> of course.
> 
> i am using the car library to recode.
> 
> it's only 5 levels when it should have 6... when i do levels(trust).
> 
> here is my recode: 
> 
> 
> trust <- recode(Poland$SN35B, " 'Strongly Agree' = 1; 'Agree' = 2; 'Neither 
> Agree nor Disagree' = 3; 'Disagree' = 4; 'Strongly Disagree' = 5; 'Can't 
> choose' = 6; else=NA")

Although it's awkward, you should be able to do what you want by "escaping" the 
quotation marks surrounding the level name; the following (untested) should 
work:

trust <- recode(Poland$SN35B, " 'Strongly Agree' = 1; 'Agree' = 2; 'Neither 
Agree nor Disagree' = 3; 'Disagree' = 4; 'Strongly Disagree' = 5; \"Can't 
choose\" = 6; else=NA")

I hope this helps,
 John


John Fox
Sen. William McMaster Prof. of Social Statistics
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

> 
> 
> thanks.
> 
> ~nicole
> - Original Message -
> From: "David Winsemius" 
> To: "Rui Barradas" 
> Cc: r-help@r-project.org
> Sent: Thursday, December 15, 2011 2:33:32 PM
> Subject: Re: [R] printing all htest class members
> 
> 
> On Dec 15, 2011, at 2:16 PM, Rui Barradas wrote:
> 
> > You're right, David,
> >
> > The first line is wrong, it should be
> >
> > ...  df=2:4  ...
> >
> > As for creating something, try
> >
> >> ht <- structure( ... etc ...
> >> ht
> >> class(ht)
> >
> > See what is printed and what function prints it.
> 
> Well, the function is stats:::print.htest.
> 
> (Do not expect any further replies to emails sent without context.)
> 
> #---#
> > 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
> 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-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] Monetary support to the R-project (Was: Re: Executable for Production Use)

2011-12-15 Thread Uwe Ligges

On 15.12.2011 04:09, Xu Wang wrote:

I am still interested in this. Is there no way to pay directly online? via
paypal or other?


No.

Uwe Ligges




Thanks,

Xu

--
View this message in context: 
http://r.789695.n4.nabble.com/Re-Monetary-support-to-the-R-project-Was-Re-Executable-for-Production-Use-tp1585186p4198369.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] removing contractions for recode in car

2011-12-15 Thread John Fox
Dear Nicole,

On Thu, 15 Dec 2011 14:49:04 -0600 (CST)
 Nicole Marie Ford  wrote:
> hello,
> 
> i need to recode a variable, however the contraction is causing problems.  i 
> had the code to change this written down somewhere and i just can't find it, 
> of course.
> 
> i am using the car library to recode.
> 
> it's only 5 levels when it should have 6... when i do levels(trust).
> 
> here is my recode: 
> 
> 
> trust <- recode(Poland$SN35B, " 'Strongly Agree' = 1; 'Agree' = 2; 'Neither 
> Agree nor Disagree' = 3; 'Disagree' = 4; 'Strongly Disagree' = 5; 'Can't 
> choose' = 6; else=NA")

Although it's awkward, you should be able to do what you want by "escaping" the 
quotation marks surrounding the level name; the following (untested) should 
work:

trust <- recode(Poland$SN35B, " 'Strongly Agree' = 1; 'Agree' = 2; 'Neither 
Agree nor Disagree' = 3; 'Disagree' = 4; 'Strongly Disagree' = 5; \"Can't 
choose\" = 6; else=NA")

I hope this helps,
 John


John Fox
Sen. William McMaster Prof. of Social Statistics
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

> 
> 
> thanks.
> 
> ~nicole
> - Original Message -
> From: "David Winsemius" 
> To: "Rui Barradas" 
> Cc: r-help@r-project.org
> Sent: Thursday, December 15, 2011 2:33:32 PM
> Subject: Re: [R] printing all htest class members
> 
> 
> On Dec 15, 2011, at 2:16 PM, Rui Barradas wrote:
> 
> > You're right, David,
> >
> > The first line is wrong, it should be
> >
> > ...  df=2:4  ...
> >
> > As for creating something, try
> >
> >> ht <- structure( ... etc ...
> >> ht
> >> class(ht)
> >
> > See what is printed and what function prints it.
> 
> Well, the function is stats:::print.htest.
> 
> (Do not expect any further replies to emails sent without context.)
> 
> #---#
> > 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
> 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-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] removing contractions for recode in car

2011-12-15 Thread David Winsemius
You can either use "\" to escape a character or you can mix two kinds  
of quotes. If you used double quotes around your text entries you can  
have singel quotes in the "interior of the text.


> 'test of bslash-single-quote \' continues'
[1] "test of bslash-single-quote ' continues"
> "test of isolated single-quote ' continues"
[1] "test of isolated single-quote ' continues"

On Dec 15, 2011, at 3:49 PM, Nicole Marie Ford wrote:


hello,

i need to recode a variable, however the contraction is causing  
problems.  i had the code to change this written down somewhere and  
i just can't find it, of course.


i am using the car library to recode.

it's only 5 levels when it should have 6... when i do levels(trust).

here is my recode:


trust <- recode(Poland$SN35B, " 'Strongly Agree' = 1; 'Agree' = 2;  
'Neither Agree nor Disagree' = 3; 'Disagree' = 4; 'Strongly  
Disagree' = 5; 'Can't choose' = 6; else=NA")



thanks.

~nicole
- Original Message -
From: "David Winsemius" 
To: "Rui Barradas" 
Cc: r-help@r-project.org
Sent: Thursday, December 15, 2011 2:33:32 PM
Subject: Re: [R] printing all htest class members


On Dec 15, 2011, at 2:16 PM, Rui Barradas wrote:


You're right, David,

The first line is wrong, it should be

...  df=2:4  ...

As for creating something, try


ht <- structure( ... etc ...
ht
class(ht)


See what is printed and what function prints it.


Well, the function is stats:::print.htest.

(Do not expect any further replies to emails sent without context.)

#---#

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
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-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
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] slight documentation error in "stats" package "arima"

2011-12-15 Thread Jan Theodore Galkowski
The documentation for the arima function in the package stats has
a slight error. It references:

Ripley, B. D. (2002) Time series in R 1.5.0. R News, 2/1,
2–7. [1]http://www.r-project.org/doc/Rnews/Rnews_2002-1.pdf

This should be:

Ripley, B. D. (2002) Time series in R 1.5.0. R News, 2/2,
2–7. [2]http://www.r-project.org/doc/Rnews/Rnews_2002-2.pdf

Anyone know who I should tell about this?

Thanks!

 - Jan

References

1. http://www.r-project.org/doc/Rnews/Rnews_2002-1.pdf
2. http://www.r-project.org/doc/Rnews/Rnews_2002-1.pdf
--

 Jan Theodore Galkowski
 Senior Systems Software Engineer
 Akamai Technologies
 Cambridge, MA 02142

 jgalk...@akamai.com
 bayesianlo...@acm.org

 607.239.1834 (m)
 607.239.1834 (h)
 617.444.4995 (w)




[[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] Random Forest Reading N/A's, I don't see them

2011-12-15 Thread Lost in R
Thanks Michael -  That was a help, i got rid of the "," in my numbers and the
"%" which were making many of the numeric variables FACTORS. It appears that
I made all of the those revisions, but still getting the same error.
Attached is the str() output if anyone could shed some light it would be
much appreciated.



Thanks,
Mike

http://r.789695.n4.nabble.com/file/n4201899/Str%28%29.docx Str%28%29.docx 

--
View this message in context: 
http://r.789695.n4.nabble.com/Random-Forest-Reading-N-A-s-I-don-t-see-them-tp4201546p4201899.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] MAximum location

2011-12-15 Thread Trying To learn again
Hi all,

I have a matrix
a<-c(2,3,4,Inf)

> b<-as.matrix(a)
 [,1]
[1,]2
[2,]3
[3,]4
[4,]  Inf

> range(b, finite=TRUE)[2] (this is the maximum)
[1] 4

There is a pre-def function to extract the location (in terms of rows) of
the value in the matrix.

In my example would be

3 (max is in the third row)

The maximum is in the position (row) 3.

[[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] printing all htest class members

2011-12-15 Thread Rui Barradas

Hello,

Once again, and as simple as possible,

res <- data.frame(ord=2:4, df=2:4, Q=c(0.0129, 0.049, 0.0684),
p=c(0.9936, 0.9972, 0.9994))

ht2<-structure(
list(statistic=c(Q=res$Q[1]),
p.value=res$p[1],
parameter=c(df=res$df[1]),
alternative="It doesn't print 'all.orders'",
method="Test the 'print.htest' method",
data.name=deparse(substitute(x)),
null.value=res  # this is printed
#all.orders=res # this wouldn't be if uncommented
),
.Names=c("statistic","p.value","parameter","alternative",
"method","data.name","null.value"),
class="htest"
)
ht2
class(ht2)

Creation? Conversation in my head? Have you heard about the experimental
method?
It's a very usefull device created by one Galileo Galilei some time ago.


--
View this message in context: 
http://r.789695.n4.nabble.com/printing-all-htest-class-members-tp4200872p4201793.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] removing contractions for recode in car

2011-12-15 Thread Nicole Marie Ford
hello,

i need to recode a variable, however the contraction is causing problems.  i 
had the code to change this written down somewhere and i just can't find it, of 
course.

i am using the car library to recode.

it's only 5 levels when it should have 6... when i do levels(trust).

here is my recode: 


trust <- recode(Poland$SN35B, " 'Strongly Agree' = 1; 'Agree' = 2; 'Neither 
Agree nor Disagree' = 3; 'Disagree' = 4; 'Strongly Disagree' = 5; 'Can't 
choose' = 6; else=NA")


thanks.

~nicole
- Original Message -
From: "David Winsemius" 
To: "Rui Barradas" 
Cc: r-help@r-project.org
Sent: Thursday, December 15, 2011 2:33:32 PM
Subject: Re: [R] printing all htest class members


On Dec 15, 2011, at 2:16 PM, Rui Barradas wrote:

> You're right, David,
>
> The first line is wrong, it should be
>
> ...  df=2:4  ...
>
> As for creating something, try
>
>> ht <- structure( ... etc ...
>> ht
>> class(ht)
>
> See what is printed and what function prints it.

Well, the function is stats:::print.htest.

(Do not expect any further replies to emails sent without context.)

#---#
> 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
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-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] printing all htest class members

2011-12-15 Thread David Winsemius


On Dec 15, 2011, at 2:16 PM, Rui Barradas wrote:


You're right, David,

The first line is wrong, it should be

...  df=2:4  ...

As for creating something, try


ht <- structure( ... etc ...
ht
class(ht)


See what is printed and what function prints it.


Well, the function is stats:::print.htest.

(Do not expect any further replies to emails sent without context.)

#---#

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
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] lme with nested factor and random effect

2011-12-15 Thread R. Michael Weylandt


On Dec 15, 2011, at 2:07 PM, Ben Bolker  wrote:

> Mari Pesek  gmail.com> writes:
> 
>> 
>> Hello all,
>> 
>> I'm having difficulty with setting up a mixed model using lme in the
>> nlme package. To summarize my study, I am testing for effects of
>> ornamentation on foraging behavior of wolf spiders. I tested spiders
>> at two different ages (penultimate vs. mature) and of two different
>> phenotypes (one species tested lacks ornamentation throughout life
>> [non-ornamented males] while the other acquires ornamentation upon
>> maturation [i.e. brush-legged males]). I tested a sample of
>> brush-legged and non-ornamented males (as both penultimates and
>> matures) in 2009, and an additional sample of brush-legged males in
>> 2010 (as both penultimates and matures again) because I had a very
>> small sample of brush-legged males in 2009.
>> 
>> I would like to set up my lme so the fixed effects are "age"
>> (penultimate vs mature), "phenotype" (non-ornamented vs brush-legged),
>> and "year" (2009 vs 2010) nested within "phenotype" to test for
>> differences between the two samples of brush-legged males.
>> Additionally I want to include "id" (a unique identification number
>> given to each spider tested) as a random factor to account for testing
>> each individual twice (once as a penultimate and once as a mature).
>> 
>> So far I have the following code: lme(behavior ~ age*phenotype,
>> random=~1|maturity/id, data)
>> 
>> But I don't know how to include the code to nest year within phenotype
>> while testing for all possible interactions. Any help would be greatly
>> appreciated.
> 
>   I have some thoughts on this.  I think your best bet is
> 
> lme(behavior~age*phenotype*year, random=~age|id, data)
> 
> or possibly
> 
> lme(behavior~age*phenotype + phenotype:year, random=~age|id, data)
> 
> ("crossing" for fixed effects is more or less equivalent to
> creating an interaction.  You should also make sure that you
> have converted 'year' to a factor rather than a numeric variable ...)
> 
>  but if you re-post this to the r-sig-mixed mod...@r-project.org list I will
> answer more fully ...

Note a hyphen got lost along the way (or at least it didn't make it to my 
machine): the email address is r-sig-mixed-mod...@r-project.org

M

> 
>  Ben Bolker
> 
> __
> 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] Random Forest Reading N/A's, I don't see them

2011-12-15 Thread R. Michael Weylandt
Use str() on your object and attach the result. For even faster help, use 
dput() on a *small* sample of your data to make the problem reproducible. 

My guess is that there are characters or, less likely, factors lurking about...

Michael

On Dec 15, 2011, at 2:39 PM, Lost in R  
wrote:

> After checking the original data in Excel for blanks and running Summary(cm3)
> to identify any null values in my data, I'm unable to identify an instances.
> Yet when I attempted to use the data in Random Forest, I get the following
> error. Is there something that Random Forest is reading as null which is not
> actually null? Is there a better way to check for this?
> 
>> library(randomForest)
>> system.time(
> + rf1 <- randomForest(as.matrix(cm3[,c(2:length(colnames(cm3)))]),
> + cm3[,1],data=cm3,ntree=50)
> + )
> *Error in randomForest.default(as.matrix(cm3[, c(2:length(colnames(cm3)))]), 
> : 
>  NA/NaN/Inf in foreign function call (arg 1)
> In addition: Warning message:
> In storage.mode(x) <- "double" : NAs introduced by coercion
> Timing stopped at: 1.33 0.01 1.35 *
> 
> 
> Thanks in advance,
> Mike
> 
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Random-Forest-Reading-N-A-s-I-don-t-see-them-tp4201546p4201546.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] printing all htest class members

2011-12-15 Thread Rui Barradas
You're right, David,

The first line is wrong, it should be 

...  df=2:4  ...

As for creating something, try

> ht <- structure( ... etc ...
> ht
> class(ht)

See what is printed and what function prints it.

Rui



--
View this message in context: 
http://r.789695.n4.nabble.com/printing-all-htest-class-members-tp4200872p4201395.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] Random Forest Reading N/A's, I don't see them

2011-12-15 Thread Lost in R
After checking the original data in Excel for blanks and running Summary(cm3)
to identify any null values in my data, I'm unable to identify an instances.
Yet when I attempted to use the data in Random Forest, I get the following
error. Is there something that Random Forest is reading as null which is not
actually null? Is there a better way to check for this?

> library(randomForest)
> system.time(
+ rf1 <- randomForest(as.matrix(cm3[,c(2:length(colnames(cm3)))]),
+ cm3[,1],data=cm3,ntree=50)
+ )
*Error in randomForest.default(as.matrix(cm3[, c(2:length(colnames(cm3)))]), 
: 
  NA/NaN/Inf in foreign function call (arg 1)
In addition: Warning message:
In storage.mode(x) <- "double" : NAs introduced by coercion
Timing stopped at: 1.33 0.01 1.35 *


Thanks in advance,
Mike

--
View this message in context: 
http://r.789695.n4.nabble.com/Random-Forest-Reading-N-A-s-I-don-t-see-them-tp4201546p4201546.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] Data Manipulation - make diagonal matrix of each element of a matrix

2011-12-15 Thread Rui Barradas
Hello,

I believe I can help, or at least, my code is simpler.
First, look at your first line:

idd <- length(diag(1,tt))   # length of intercept matrix
#
not needed: diag(tt) would do the job but it's not needed,
why call 2 functions, and one of them, 'diag', uses memory(*), if the
result is tt squared? It's much simpler!
(*)like you say, "larger and larger" amounts of it

My solution to your problem is as follows (as a function, and yours).

fun2 <- function(n, tt, numco){
M.Unit <- matrix(rep(diag(1,tt),n), ncol=tt, byrow=TRUE)
M <- NULL
for(i in 1:numco) M <- cbind(M, M.Unit*rep(x[,i], each=tt))
M
}

fun1 <- function(n, tt, numco){
idd <- length(diag(1,tt))# length of intercept matrix
X <- matrix(numeric(n*numco*idd),ncol=tt*numco)
for(i in 1:numco){
  X[,((i-1)*tt+1):(i*tt)] <- matrix(
c(matrix(rep(diag(1,tt),n),ncol=tt, byrow=TRUE))*
rep(rep(x[,i],each=tt),tt)
   , ncol=tt)
}
X
}

I' ve tested the two with larger values of 'n', 'tt' and 'numco'
using the following timing instructions


n  <- 1000
tt <- 50
numco <- 15
set.seed(1)
x <- matrix(round(rnorm(n*numco),2), ncol=numco)   # the actual covariates

Runs <- 10^1

t1 <- system.time(for(i in 1:Runs) a1 <- fun1(n, tt, numco))[c(1,3)]
t2 <- system.time(for(i in 1:Runs) a2 <- fun2(n, tt, numco))[c(1,3)]

rbind(t1, t2, t1/t2)

  user.self elapsed
t1 23.21   31.06
t2 14.97   22.54
 1.5504341.377995

As you can see, it's not a great speed improvement.
I hope it's at least usefull.

Rui Barradas


--
View this message in context: 
http://r.789695.n4.nabble.com/Data-Manipulation-make-diagonal-matrix-of-each-element-of-a-matrix-tp4200321p4201305.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] UseR! 2011 slides and videos - now online

2011-12-15 Thread larrydag
Tal,

Keep up the great work with r-bloggers.  I'm the organizer of the Dallas RUG
and if there is any content you wish to provide from the Dallas RUG meetup
site feel free to post it at www.r-bloggers.com/RUG.  We're still a fairly
young user group but we do have some content from our presentations.

http://www.meetup.com/Dallas-R-Users-Group/files/

Thanks for support the R user groups.

Larry

--
View this message in context: 
http://r.789695.n4.nabble.com/UseR-2011-slides-and-videos-now-online-tp4190222p4201280.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] lattice key in blank panel

2011-12-15 Thread Max Kuhn
Somewhere I've seen an example of an xyplot() where the key was placed
in a location of a missing panel. For example, if there were 3
conditioning levels, the panel grid would look like:

34
12

In this (possibly imaginary) example, there were scatter plots in
locations 1:3 and location 4 had no conditioning bar at the top, only
the key.

I can find examples of putting the legend outside of the panel
locations (e.g to the right of locations 2 and 4 above), but that's
not really what I'd like to do.

Thanks,

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] Bioconductor. MA plot for qPCR array

2011-12-15 Thread Juliet Hannah
You may find the following discussion helpful.

http://comments.gmane.org/gmane.science.biology.informatics.conductor/37388

On Sun, Dec 11, 2011 at 8:08 AM, ali_protocol
 wrote:
> Dear all,
>
> Is there anyway too generate MA plot for 2 qPCR assays (an array of 2x 400).
>
>
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Bioconductor-MA-plot-for-qPCR-array-tp4182805p4182805.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] Ratio of huge products

2011-12-15 Thread David Winsemius


On Dec 15, 2011, at 1:35 PM, Alberto Magni wrote:


Hello everybody,

I have to compute something in this form:

x = prod(a:b) / prod(c:d),where: a < c and b < d and obviously: a
< b and c < d

I cannot make assumptions on the relative position of c,b and a,d.

The problem is that a,b,c,d are large and the products are huge (R  
return Inf).


Well, R does have some limitations.


Their ratio is less than 1 but significantly higher than 0: it is a
non-tiny probability.

I need to find a way to simplify this ratio.


x <- exp( sum(log(a:b)) -sum(log(c:d)) )


The only way to solve this that I see is to decompose into prime
factors all the
numbers in the numerator and the denominator and to remove the ones  
in common


Ewww. That does sound painful.

--

David Winsemius, MD
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] lme with nested factor and random effect

2011-12-15 Thread Ben Bolker
Mari Pesek  gmail.com> writes:

> 
> Hello all,
> 
> I'm having difficulty with setting up a mixed model using lme in the
> nlme package. To summarize my study, I am testing for effects of
> ornamentation on foraging behavior of wolf spiders. I tested spiders
> at two different ages (penultimate vs. mature) and of two different
> phenotypes (one species tested lacks ornamentation throughout life
> [non-ornamented males] while the other acquires ornamentation upon
> maturation [i.e. brush-legged males]). I tested a sample of
> brush-legged and non-ornamented males (as both penultimates and
> matures) in 2009, and an additional sample of brush-legged males in
> 2010 (as both penultimates and matures again) because I had a very
> small sample of brush-legged males in 2009.
> 
> I would like to set up my lme so the fixed effects are "age"
> (penultimate vs mature), "phenotype" (non-ornamented vs brush-legged),
> and "year" (2009 vs 2010) nested within "phenotype" to test for
> differences between the two samples of brush-legged males.
> Additionally I want to include "id" (a unique identification number
> given to each spider tested) as a random factor to account for testing
> each individual twice (once as a penultimate and once as a mature).
> 
> So far I have the following code: lme(behavior ~ age*phenotype,
> random=~1|maturity/id, data)
> 
> But I don't know how to include the code to nest year within phenotype
> while testing for all possible interactions. Any help would be greatly
> appreciated.

   I have some thoughts on this.  I think your best bet is

lme(behavior~age*phenotype*year, random=~age|id, data)

or possibly

lme(behavior~age*phenotype + phenotype:year, random=~age|id, data)

("crossing" for fixed effects is more or less equivalent to
creating an interaction.  You should also make sure that you
have converted 'year' to a factor rather than a numeric variable ...)

  but if you re-post this to the r-sig-mixed mod...@r-project.org list I will
answer more fully ...

  Ben Bolker

__
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] printing all htest class members

2011-12-15 Thread David Winsemius


On Dec 15, 2011, at 12:26 PM, Rui Barradas wrote:


Hello,

I've posted a question about this subject yesterday, but since there  
was no

R code to comment,
no one did.

I'm trying to have the print method for class 'htest' print some extra
information common in some test, like the time series linearity  
related

tests. Many of them have an 'order' parameter, representing a lag or
embedding dimension, and it would be a nice feature to pass a vector  
of

orders and have all corresponding tests run.

It's easy to do this, but the print method doesn't print those extra  
values.

Here the code  goes:

#Some values from the McLeod-Li test, with x <- rnorm(100)
res <- data.frame(ord=2:4, df=ord, Q=c(0.0129, 0.049, 0.0684),
   p=c(0.9936, 0.9972, 0.9994))
attach(res)
nr <- nrow(res)



My problem reading this is that you are using the attach function  
which generally confuses discussions and then you are not actually  
creating any object with a class of htest. Print methods are  
dispatched by class value.



# print.htest prints everything but 'all.orders'
# but when it's named 'null.value' it works and it's easier to make
# it work, all what's needed is 'null.value=res', whithout the need
# for a second 'structure()'


You seem to be transcribing a conversation inside your head. Most of  
the four lines above has no connection to what you typed earlier.



structure(


This does nothing. No object is created.


list(statistic=c(Q=Q[1]),
p.value=p[1],
parameter=c(df=df[1]),
alternative="It doesn't print 'all.orders' and I find 
'null.value'
misleading",
method="Test the 'print.htest' method using McLeod-Li test 
values.",
data.name=deparse(substitute(x)),
all.orders=structure(
list(order=ord, df=df, Q=Q, p.value=p),
.Names=c("order", "df", "Q", "p.value"),
row.names=c(NA,-nr),
class="data.frame"
)
),
.Names=c("statistic","p.value","parameter","alternative",
"method","data.name","all.orders"),
class="htest"
)

Is there a way to have 'all.orders' printed by 'print.htest' ?
The problem is NOT the return values, at least the time wasn't wasted,
I can use them for whatever I'll do next.

If not, anyone has any suggestions?


Yes. Give the code you are actually using!





--

David Winsemius, MD
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] Ratio of huge products

2011-12-15 Thread Peter Langfelder
On Thu, Dec 15, 2011 at 10:35 AM, Alberto Magni
 wrote:
> Hello everybody,
>
> I have to compute something in this form:
>
> x = prod(a:b) / prod(c:d),    where: a < c and b < d and obviously: a
> < b and c < d
>
> I cannot make assumptions on the relative position of c,b and a,d.
>
> The problem is that a,b,c,d are large and the products are huge (R return 
> Inf).
> Their ratio is less than 1 but significantly higher than 0: it is a
> non-tiny probability.
>
> I need to find a way to simplify this ratio.
> The only way to solve this that I see is to decompose into prime
> factors all the
> numbers in the numerator and the denominator and to remove the ones in common
>
> Do you know a better way to do this ?

Yes, exp( sum(log(a:c))-sum(log(b:d)) ), which is mathematically
exactly equivalent, unless I made a typo.
Remeber that log of a product is a sum of the logs of the arguments.

Peter.

__
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] From Distance Matrix to 2D coordinates

2011-12-15 Thread Peter Langfelder
On Thu, Dec 15, 2011 at 10:47 AM, Lorenzo Isella
 wrote:
> Thanks a lot!
> Precisely what I had in mind.
> One last question (an extension of the previous one): can this be extended
> to points in 3D? Once again, given the distance matrix, can I reconstruct a
> set of coordinates (among many possible) for the points in three-dimensional
> space?
> Cheers
>
> Lorenzo

Yes, you need to specify argument 'k=3' to cmdscale which instructs it
to 'scale' the input into 3 dimensions (the default is 2).

nPoints = 10;
nDim = 3;

set.seed(10);
points = matrix(runif(nPoints * nDim), nPoints, nDim);

# Their distance:
dst = dist(points)

# Classical multidimensional scaling
mds = cmdscale(dst, k=nDim);

# Distance of the points calculated by mds
dst2 = dist(mds);

# The two distances are equal
all.equal(as.vector(dst), as.vector(dst2))

HTH

Peter

__
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] Ratio of huge products

2011-12-15 Thread R. Michael Weylandt
Use logs?

Michael

On Thu, Dec 15, 2011 at 1:35 PM, Alberto Magni
 wrote:
> Hello everybody,
>
> I have to compute something in this form:
>
> x = prod(a:b) / prod(c:d),    where: a < c and b < d and obviously: a
> < b and c < d
>
> I cannot make assumptions on the relative position of c,b and a,d.
>
> The problem is that a,b,c,d are large and the products are huge (R return 
> Inf).
> Their ratio is less than 1 but significantly higher than 0: it is a
> non-tiny probability.
>
> I need to find a way to simplify this ratio.
> The only way to solve this that I see is to decompose into prime
> factors all the
> numbers in the numerator and the denominator and to remove the ones in common
>
> Do you know a better way to do this ?
>
> Thank you,
> Alberto
>
> __
> 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.


[R] lme with nested factor and random effect

2011-12-15 Thread Mari Pesek
Hello all,

I'm having difficulty with setting up a mixed model using lme in the
nlme package. To summarize my study, I am testing for effects of
ornamentation on foraging behavior of wolf spiders. I tested spiders
at two different ages (penultimate vs. mature) and of two different
phenotypes (one species tested lacks ornamentation throughout life
[non-ornamented males] while the other acquires ornamentation upon
maturation [i.e. brush-legged males]). I tested a sample of
brush-legged and non-ornamented males (as both penultimates and
matures) in 2009, and an additional sample of brush-legged males in
2010 (as both penultimates and matures again) because I had a very
small sample of brush-legged males in 2009.

I would like to set up my lme so the fixed effects are "age"
(penultimate vs mature), "phenotype" (non-ornamented vs brush-legged),
and "year" (2009 vs 2010) nested within "phenotype" to test for
differences between the two samples of brush-legged males.
Additionally I want to include "id" (a unique identification number
given to each spider tested) as a random factor to account for testing
each individual twice (once as a penultimate and once as a mature).

So far I have the following code: lme(behavior ~ age*phenotype,
random=~1|maturity/id, data)

But I don't know how to include the code to nest year within phenotype
while testing for all possible interactions. Any help would be greatly
appreciated.

Thank you!
Mari Pesek

__
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] From Distance Matrix to 2D coordinates

2011-12-15 Thread Lorenzo Isella

Thanks a lot!
Precisely what I had in mind.
One last question (an extension of the previous one): can this be 
extended to points in 3D? Once again, given the distance matrix, can I 
reconstruct a set of coordinates (among many possible) for the points in 
three-dimensional space?

Cheers

Lorenzo

On 12/15/2011 07:22 PM, Peter Langfelder wrote:

On Thu, Dec 15, 2011 at 10:08 AM, Lorenzo Isella
  wrote:

Dear All,
I am struggling with the following problem: I am given a NxN symmetric
matrix P ( P[i,i]=0, i=1...N and P[i,j]>0 for i!=j) which stands for the
relative distances of N points.
I would like use it to get the coordinates of the N points in a 2D plane. Of
course, the solution is not unique (given one solution, I can translate or
rotate all the points by the same amount and generate another solution), but
any correct solution will do for me.
Any idea about how I can achieve that? Is there any clustering package that
can help me?
Many thanks.


If your matrix really corresponds to distances of points (in 2
dimensions), you can try multidimensional scaling, function
cmdscale().

This little code illustrates that cmdscale recovers the 2-dimensional
points used to generate a distance matrix, up to a shift and rotation:

# Generate 10 random points in 2 dimensions
nPoints = 10;
nDim = 2;

set.seed(10);
points = matrix(runif(nPoints * nDim), nPoints, nDim);

# Their distance:
dst = dist(points)

# Classical multidimensional scaling
mds = cmdscale(dst);

# Distance of the points calculated by mds
dst2 = dist(mds);

# The two distances are equal
all.equal(as.vector(dst), as.vector(dst2))

HTH,

Peter


__
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] Ratio of huge products

2011-12-15 Thread Alberto Magni
Hello everybody,

I have to compute something in this form:

x = prod(a:b) / prod(c:d),where: a < c and b < d and obviously: a
< b and c < d

I cannot make assumptions on the relative position of c,b and a,d.

The problem is that a,b,c,d are large and the products are huge (R return Inf).
Their ratio is less than 1 but significantly higher than 0: it is a
non-tiny probability.

I need to find a way to simplify this ratio.
The only way to solve this that I see is to decompose into prime
factors all the
numbers in the numerator and the denominator and to remove the ones in common

Do you know a better way to do this ?

Thank you,
Alberto

__
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] printing all htest class members

2011-12-15 Thread Rui Barradas
Hello, 

I've posted a question about this subject yesterday, but since there was no
R code to comment,
no one did.

I'm trying to have the print method for class 'htest' print some extra
information common in some test, like the time series linearity related
tests. Many of them have an 'order' parameter, representing a lag or
embedding dimension, and it would be a nice feature to pass a vector of
orders and have all corresponding tests run.

It's easy to do this, but the print method doesn't print those extra values.
Here the code  goes:

#Some values from the McLeod-Li test, with x <- rnorm(100)
res <- data.frame(ord=2:4, df=ord, Q=c(0.0129, 0.049, 0.0684),
p=c(0.9936, 0.9972, 0.9994))
attach(res)
nr <- nrow(res)

# print.htest prints everything but 'all.orders'
# but when it's named 'null.value' it works and it's easier to make
# it work, all what's needed is 'null.value=res', whithout the need
# for a second 'structure()'
structure(
list(statistic=c(Q=Q[1]),
p.value=p[1],
parameter=c(df=df[1]),
alternative="It doesn't print 'all.orders' and I find 
'null.value'
misleading",
method="Test the 'print.htest' method using McLeod-Li test 
values.",
data.name=deparse(substitute(x)),
all.orders=structure(
list(order=ord, df=df, Q=Q, p.value=p),
.Names=c("order", "df", "Q", "p.value"),
row.names=c(NA,-nr),
class="data.frame"
)
),
.Names=c("statistic","p.value","parameter","alternative",
"method","data.name","all.orders"),
class="htest"
)

Is there a way to have 'all.orders' printed by 'print.htest' ? 
The problem is NOT the return values, at least the time wasn't wasted,
I can use them for whatever I'll do next.

If not, anyone has any suggestions?

Thank you in advance,

Rui Barradas


--
View this message in context: 
http://r.789695.n4.nabble.com/printing-all-htest-class-members-tp4200872p4200872.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] From Distance Matrix to 2D coordinates

2011-12-15 Thread Peter Langfelder
On Thu, Dec 15, 2011 at 10:08 AM, Lorenzo Isella
 wrote:
> Dear All,
> I am struggling with the following problem: I am given a NxN symmetric
> matrix P ( P[i,i]=0, i=1...N and P[i,j]>0 for i!=j) which stands for the
> relative distances of N points.
> I would like use it to get the coordinates of the N points in a 2D plane. Of
> course, the solution is not unique (given one solution, I can translate or
> rotate all the points by the same amount and generate another solution), but
> any correct solution will do for me.
> Any idea about how I can achieve that? Is there any clustering package that
> can help me?
> Many thanks.

If your matrix really corresponds to distances of points (in 2
dimensions), you can try multidimensional scaling, function
cmdscale().

This little code illustrates that cmdscale recovers the 2-dimensional
points used to generate a distance matrix, up to a shift and rotation:

# Generate 10 random points in 2 dimensions
nPoints = 10;
nDim = 2;

set.seed(10);
points = matrix(runif(nPoints * nDim), nPoints, nDim);

# Their distance:
dst = dist(points)

# Classical multidimensional scaling
mds = cmdscale(dst);

# Distance of the points calculated by mds
dst2 = dist(mds);

# The two distances are equal
all.equal(as.vector(dst), as.vector(dst2))

HTH,

Peter

__
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] From Distance Matrix to 2D coordinates

2011-12-15 Thread Sarah Goslee
That's exactly what ordination is for (not clustering).

I'd try principal coordinates analysis, or non-metric multidimensional
scaling, depending on whether the dissimilarity you'v been given is
metric or nonmetric.

There are implementations of both in the ecodist package, and in
various other packages as well, so you have lots of choice.

Sarah

On Thu, Dec 15, 2011 at 1:08 PM, Lorenzo Isella
 wrote:
> Dear All,
> I am struggling with the following problem: I am given a NxN symmetric
> matrix P ( P[i,i]=0, i=1...N and P[i,j]>0 for i!=j) which stands for the
> relative distances of N points.
> I would like use it to get the coordinates of the N points in a 2D plane. Of
> course, the solution is not unique (given one solution, I can translate or
> rotate all the points by the same amount and generate another solution), but
> any correct solution will do for me.
> Any idea about how I can achieve that? Is there any clustering package that
> can help me?
> Many thanks.
>
> Lorenzo
>


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


[R] From Distance Matrix to 2D coordinates

2011-12-15 Thread Lorenzo Isella

Dear All,
I am struggling with the following problem: I am given a NxN symmetric 
matrix P ( P[i,i]=0, i=1...N and P[i,j]>0 for i!=j) which stands for the 
relative distances of N points.
I would like use it to get the coordinates of the N points in a 2D 
plane. Of course, the solution is not unique (given one solution, I can 
translate or rotate all the points by the same amount and generate 
another solution), but any correct solution will do for me.
Any idea about how I can achieve that? Is there any clustering package 
that can help me?

Many thanks.

Lorenzo

__
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] fundamental guide to use of numerical optimizers?

2011-12-15 Thread Greg Snow
This really depends on more than just the optimizer, a lot can depend on what 
the data looks like and what question is being asked.  In bootstrapping it is 
possible to get bootstrap samples for which there is no unique correct answer 
to converge to, for example if there is a category where there ends up being no 
data due to the bootstrap but you still want to estimate a parameter for that 
category then there are an infinite number of possible answers that are all 
equal in the likelihood so there will be a lack of convergence on that 
parameter.  A stratified bootstrap or semi-parametric bootstrap can be used to 
avoid this problem (but may change the assumptions being made as well), or you 
can just throw out all those samples that don't have a full answer (which could 
be what your presenter did).

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
> project.org] On Behalf Of Paul Johnson
> Sent: Thursday, December 15, 2011 9:38 AM
> To: R-help
> Subject: [R] fundamental guide to use of numerical optimizers?
> 
> I was in a presentation of optimizations fitted with both MPlus and
> SAS yesterday.  In a batch of 1000 bootstrap samples, between 300 and
> 400 of the estimations did not converge.  The authors spoke as if this
> were the ordinary cost of doing business, and pointed to some
> publications in which the nonconvergence rate was as high or higher.
> 
> I just don't believe that's right, and if some problem is posed so
> that the estimate is not obtained in such a large sample of
> applications, it either means the problem is badly asked or badly
> answered.  But I've got no traction unless I can actually do
> better
> 
> Perhaps I can use this opportunity to learn about R functions like
> optim, or perhaps maxLik.
> 
> >From reading r-help, it seems to me there are some basic tips for
> optimization, such as:
> 
> 1. It is wise to scale the data so that all columns have the same
> range before running an optimizer.
> 
> 2. With estimates of variance parameters, don't try to estimate sigma
> directly, instead estimate log(sigma) because that puts the domain of
> the solution onto the real number line.
> 
> 3 With estimates of proportions, estimate instead the logit, for the
> same reason.
> 
> Are these mistaken generalizations?  Are there other tips that
> everybody ought to know?
> 
> I understand this is a vague question, perhaps the answers are just in
> the folklore. But if somebody has written them out, I would be glad to
> know.
> 
> --
> Paul E. Johnson
> Professor, Political Science
> 1541 Lilac Lane, Room 504
> University of Kansas
> 
> __
> 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] modify the name of axis of an R function

2011-12-15 Thread David Winsemius


On Dec 15, 2011, at 11:30 AM, plocq wrote:

Thanks you very much!  The way plot(fit, main="main title", xlab="X- 
axis
lable", ylab="y-axis label") seems to work quite well, I didn't  
notice that

I could do this.


You are replying to a message on a mailing list (which most people do  
NOT see as a web page)  and the Postng Guide asks that you include  
context.


I have in fact one more problem with it : the fact is that I have  
three
plots that are called by the function. I can specify my three titles  
by

doing "main=c("title1","title2","title3")" whithout any problem. But
proceeding this way to specify the axes' titles doesn't work  
anymore. How

can I specify three differents pairs of titles for these axes?


The Posting Guide also asks that you provide a specific example of  
what you mean by "not working".




Thank you again

--
View this message in context: 
http://r.789695.n4.nabble.com/modify-the-name-of-axis-of-an-R-function-tp4198804p4200631.html


No, we don't really want to do that. Nabble user are guests and need  
to learn that Rhelp is not Nabble.



Sent from the R help mailing list archive at Nabble.com.


And Rhelp is not at Nabble, and it's not an archive, either.

--

David Winsemius, MD
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] Multicollinearty in logistic regression models

2011-12-15 Thread David Winsemius


On Dec 15, 2011, at 11:34 AM, Mohamed Lajnef wrote:


Dear All,

Is there a method to diagnostic multicollinearty in logistic  
regression

models  like vif indicator in linear regression ( variance inflation
Factor ...) ?



Wouldn't matrix representation of the predictor "side" of the  
regression be the same? Couldn't you just use the same methods you  
employ for linear regression?



Thank you in advance
M

--

Mohamed Lajnef,IE INSERM U955 eq 15#



David Winsemius, MD
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] modify the name of axis of an R function

2011-12-15 Thread plocq
Thanks you very much!  The way plot(fit, main="main title", xlab="X-axis
lable", ylab="y-axis label") seems to work quite well, I didn't notice that
I could do this. 

I have in fact one more problem with it : the fact is that I have three
plots that are called by the function. I can specify my three titles by
doing "main=c("title1","title2","title3")" whithout any problem. But
proceeding this way to specify the axes' titles doesn't work anymore. How
can I specify three differents pairs of titles for these axes?

Thank you again

--
View this message in context: 
http://r.789695.n4.nabble.com/modify-the-name-of-axis-of-an-R-function-tp4198804p4200631.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] fundamental guide to use of numerical optimizers?

2011-12-15 Thread Paul Johnson
I was in a presentation of optimizations fitted with both MPlus and
SAS yesterday.  In a batch of 1000 bootstrap samples, between 300 and
400 of the estimations did not converge.  The authors spoke as if this
were the ordinary cost of doing business, and pointed to some
publications in which the nonconvergence rate was as high or higher.

I just don't believe that's right, and if some problem is posed so
that the estimate is not obtained in such a large sample of
applications, it either means the problem is badly asked or badly
answered.  But I've got no traction unless I can actually do
better

Perhaps I can use this opportunity to learn about R functions like
optim, or perhaps maxLik.

>From reading r-help, it seems to me there are some basic tips for
optimization, such as:

1. It is wise to scale the data so that all columns have the same
range before running an optimizer.

2. With estimates of variance parameters, don't try to estimate sigma
directly, instead estimate log(sigma) because that puts the domain of
the solution onto the real number line.

3 With estimates of proportions, estimate instead the logit, for the
same reason.

Are these mistaken generalizations?  Are there other tips that
everybody ought to know?

I understand this is a vague question, perhaps the answers are just in
the folklore. But if somebody has written them out, I would be glad to
know.

-- 
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas

__
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] Multicollinearty in logistic regression models

2011-12-15 Thread Mohamed Lajnef
Dear All,

Is there a method to diagnostic multicollinearty in logistic regression 
models  like vif indicator in linear regression ( variance inflation 
Factor ...) ?

Thank you in advance
M

-- 

Mohamed Lajnef,IE INSERM U955 eq 15#
P?le de Psychiatrie#
H?pital CHENEVIER  #
40, rue Mesly  #
94010 CRETEIL Cedex FRANCE #
mohamed.laj...@inserm.fr   #
tel : 01 49 81 32 79   #
Sec : 01 49 81 32 90   #
fax : 01 49 81 30 99   #




[[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] Calculate AUC Using the Trapezoidal Method

2011-12-15 Thread Chris Campbell
The arguments time, id and dv take character strings.

AUC(Data, time = "Time", id = "Fraction", dv = "Variable")  
  Fraction  AUC
1   C1 4413.549

Regards, 

Chris Campbell
MANGO SOLUTIONS
Data Analysis that Delivers
+44 1249 767700

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of arivald
Sent: 15 December 2011 13:54
To: r-help@r-project.org
Subject: [R] Calculate AUC Using the Trapezoidal Method

Hello,
I want to use function AUC {MIfuns} but I have problem with arguments.
My data:
 Fraction Time Variable
C1 0  0.0
C1   15  20.95475
C130 28.55030
C160 36.33064
C190 48.80438
C1  120 60.18636
 
AUC(data, time = Time, id = Fraction, dv = Variable)  ##this is not working!



Thx!

--
View this message in context: 
http://r.789695.n4.nabble.com/Calculate-AUC-Using-the-Trapezoidal-Method-tp4200049p4200049.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.
LEGAL NOTICE
This message is intended for the use o...{{dropped:10}}

__
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] R error

2011-12-15 Thread verse123
Hi guys, 

I am new to R and I am bascially trying to load a library that I installed
and use external data that I have. When trying to use an R package called
cummeRbund (http://compbio.mit.edu/cummeRbund/), I am doing:

> library(cummeRbund)
Loading required package: RSQLite
Loading required package: DBI
Loading required package: reshape
Loading required package: plyr
> cuff <- readCufflinks(system.file("cds.diff", package="cummeRbund"))

I get this error message:

Creating database /cuffData.db
Error in sqliteNewConnection(drv, ...) :
  RS-DBI driver: (could not connect to dbname:
unable to open database file

Any ideas why?



--
View this message in context: 
http://r.789695.n4.nabble.com/R-error-tp4200447p4200447.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] Am I misunderstanding loop variable assignment or how to use print()?

2011-12-15 Thread Tony Stocker
On Thu, Dec 15, 2011 at 10:54, Sarah Goslee  wrote:
>
> print(get(x)[["Pr"]]) maybe. Do the get(), then do the subsetting.
>
>
>>>
>>> It's often neater and more efficient to store your anova objects in a
>>> list, though.
>>
>> anything since it's still a set of character strings.  Could you
>> elaborate a bit on what you mean by storing the anova objects as
>> lists?
>
> Yes: not the names, but the anova objects themselves. Rather than
> creating a bunch of individual objects, store them in a list when
> created:
>
> myanova <- list()
> myanova[["ag.m2529.az"]] <- anova(whatever)
> myanova[["ag.m2529.can"]] <- anova(whatever)
> ...
>
> Then you can quite elegantly use lapply() across all of the anovas at once,
> and don't have so many objects in your workspace.

Okay I'll give that a try and see how it works for me.  Thanks for the
suggestion.

__
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] Am I misunderstanding loop variable assignment or how to use print()?

2011-12-15 Thread Sarah Goslee
Hi,

On Thu, Dec 15, 2011 at 10:43 AM, Tony Stocker  wrote:
> On Thu, Dec 15, 2011 at 09:51, Sarah Goslee  wrote:
>> But "anova.ag.m2529.az"  is a character string that happens to be the
>> *name* of an anova object, but R has no way to know that unless you
>> specifically tell it that your character string is an object by using
>> get().
>>
>> Something like print(get(x)) would work.
>
> Sarah - Thanks very much!  That did indeed work great at printing the
> entire contents out.  I couldn't do print(get(x$Pr)), but I can live
> with that for now.

print(get(x)[["Pr"]]) maybe. Do the get(), then do the subsetting.


>>
>> It's often neater and more efficient to store your anova objects in a
>> list, though.
>
> So if I were to do:
>> is.list(an)
> [1] FALSE
>> alist<-list(an)
>> is.list(alist)
> [1] TRUE
>> alist
> [1] "anova.ag.m2529.az"   "anova.ag.m2529.can"   "anova.ag.m2529.fl"
>
> I would have created a list, but I'm assuming that you mean something
> different than that since I'm not sure how that functionally changed
> anything since it's still a set of character strings.  Could you
> elaborate a bit on what you mean by storing the anova objects as
> lists?

Yes: not the names, but the anova objects themselves. Rather than
creating a bunch of individual objects, store them in a list when
created:

myanova <- list()
myanova[["ag.m2529.az"]] <- anova(whatever)
myanova[["ag.m2529.can"]] <- anova(whatever)
...

Then you can quite elegantly use lapply() across all of the anovas at once,
and don't have so many objects in your workspace.

Sarah

-- 
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] Am I misunderstanding loop variable assignment or how to use print()?

2011-12-15 Thread Tony Stocker
On Thu, Dec 15, 2011 at 09:51, Sarah Goslee  wrote:
> But "anova.ag.m2529.az"  is a character string that happens to be the
> *name* of an anova object, but R has no way to know that unless you
> specifically tell it that your character string is an object by using
> get().
>
> Something like print(get(x)) would work.

Sarah - Thanks very much!  That did indeed work great at printing the
entire contents out.  I couldn't do print(get(x$Pr)), but I can live
with that for now.

>
> It's often neater and more efficient to store your anova objects in a
> list, though.

So if I were to do:
> is.list(an)
[1] FALSE
> alist<-list(an)
> is.list(alist)
[1] TRUE
> alist
[1] "anova.ag.m2529.az"   "anova.ag.m2529.can"   "anova.ag.m2529.fl"

I would have created a list, but I'm assuming that you mean something
different than that since I'm not sure how that functionally changed
anything since it's still a set of character strings.  Could you
elaborate a bit on what you mean by storing the anova objects as
lists?

Again, thanks for the help!

__
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] Chinese translation of: Beginner's Guide to R

2011-12-15 Thread Highland Statistics Ltd
Chinese members of this mailing list may be interested to know that the 
Chinese translation of 'A Beginner's Guide to R' is now available from 
amazon.cn. The full URL is:



http://www.amazon.cn/R%E8%AF%AD%E8%A8%80%E5%88%9D%E5%AD%A6%E8%80%85%E6%8C%87%E5%8D%97-%E9%98%BF%E5%85%B0-F-%E7%A5%96%E5%B0%94/dp/B005LTWESQ/ref=pd_sim_b_cnclic_2 



Kind regards,

Alain Zuur

__
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] Data Manipulation - make diagonal matrix of each element of a matrix

2011-12-15 Thread Clemontina Alexander
I'm sorry, the indices of my X matrix are wrong.
It should be:

 X = x11  0  0 x12  0  0
  0   x11  00   x12  0
  0  0   x110  0   x12
   x21  0  0 x22  0  0
  0   x21  00   x22  0
  0  0   x210  0   x22
   ...
   xn1  0  0 x52  0  0
  0   xn1  00   x52  0
  0  0   xn10  0   x52

or

 X = -0.630  0-0.82 0 0
  0   -0.630  0 0-0.82 0
  0  0   -0.630 0  0   -0.82
   0.180  0 0.49  0 0
  0 0.18  0 0  0.49 0
  0 0  0.18 0  0 0.49
   ...
0.33  0  0-0.31  0  0
  00.33  0  0-0.31  0
  00  0.33  0 0 -0.31

Sorry for the confusion.
Tina








On Thu, Dec 15, 2011 at 10:02 AM, Clemontina Alexander
 wrote:
> Dear R list,
> I have the following data:
>
> set.seed(1)
> n  <- 5     # number of subjects
> tt <- 3     # number of repeated observation per subject
> numco <- 2  # number of covariates
> x <- matrix(round(rnorm(n*numco),2), ncol=numco)   # the actual covariates
> x
>> x
>      [,1]  [,2]
> [1,] -0.63 -0.82
> [2,]  0.18  0.49
> [3,] -0.84  0.74
> [4,]  1.60  0.58
> [5,]  0.33 -0.31
>
> I need to form a matrix X such that
> X =      x11      0      0     x21      0      0
>              0   x11      0        0   x21      0
>              0      0   x11        0      0   x21
>           x12      0      0     x22      0      0
>              0   x12      0        0   x22      0
>              0      0   x12        0      0   x22
>                       ...
>           x15      0      0     x25      0      0
>              0   x15      0        0   x25      0
>              0      0   x15        0      0   x25
> where both tt and numco can change. (So if tt=5 and numco=4, then X
> needs to have 20 columns and n*tt rows. Each diagonal matrix should be
> 5x5 and there will be 4 of them for the 4 covariates.) I wrote this
> funky for loop:
>
> idd <- length(diag(1,tt))    # length of intercept matrix
> X <- matrix(numeric(n*numco*idd),ncol=tt*numco)
> for(i in 1:numco){
>      X[,((i-1)*tt+1):(i*tt)] <- matrix(
>        c(matrix(rep(diag(1,tt),n),ncol=tt, byrow=TRUE))   *
> rep(rep(x[,i],each=tt),tt)
>       , ncol=tt)
> }
> X
>
> It works fine, but is there an easier way when n, tt, and numco get
> larger and larger?
> Thanks,
> Tina
>
>
> --
> Clemontina Alexander
> Ph.D Student
> Department of Statistics
> NC State University
> Email: ckale...@ncsu.com



-- 
Clemontina Alexander
Ph.D Student
Department of Statistics
NC State University
Raleigh, NC 27695
Phone: (850) 322-6878
Email: ckale...@ncsu.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] Data Manipulation - make diagonal matrix of each element of a matrix

2011-12-15 Thread Clemontina Alexander
Dear R list,
I have the following data:

set.seed(1)
n  <- 5 # number of subjects
tt <- 3 # number of repeated observation per subject
numco <- 2  # number of covariates
x <- matrix(round(rnorm(n*numco),2), ncol=numco)   # the actual covariates
x
> x
  [,1]  [,2]
[1,] -0.63 -0.82
[2,]  0.18  0.49
[3,] -0.84  0.74
[4,]  1.60  0.58
[5,]  0.33 -0.31

I need to form a matrix X such that
X =  x11  0  0 x21  0  0
  0   x11  00   x21  0
  0  0   x110  0   x21
   x12  0  0 x22  0  0
  0   x12  00   x22  0
  0  0   x120  0   x22
   ...
   x15  0  0 x25  0  0
  0   x15  00   x25  0
  0  0   x150  0   x25
where both tt and numco can change. (So if tt=5 and numco=4, then X
needs to have 20 columns and n*tt rows. Each diagonal matrix should be
5x5 and there will be 4 of them for the 4 covariates.) I wrote this
funky for loop:

idd <- length(diag(1,tt))# length of intercept matrix
X <- matrix(numeric(n*numco*idd),ncol=tt*numco)
for(i in 1:numco){
  X[,((i-1)*tt+1):(i*tt)] <- matrix(
c(matrix(rep(diag(1,tt),n),ncol=tt, byrow=TRUE))   *
rep(rep(x[,i],each=tt),tt)
   , ncol=tt)
}
X

It works fine, but is there an easier way when n, tt, and numco get
larger and larger?
Thanks,
Tina


--
Clemontina Alexander
Ph.D Student
Department of Statistics
NC State University
Email: ckale...@ncsu.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] modify the name of axis of an R function

2011-12-15 Thread David Winsemius


On Dec 15, 2011, at 1:53 AM, plocq wrote:


Hi,

I use the function fpot of packages evd. If I call the fit that I  
obtain
"fit", I want to modify the name of the axis and the main title that  
is

produced by plot(fit1).


Usually this would be accomplished with

plot(fit, main="main title", xlab="X-axis lable", ylab="y-axis label")


To do this, I want to create a new function where I
would have the names modified.


Seems unlikely that this would be needed.


The problem is that I can't find the function
that produce these plots... I tried to see in plot.uvevd but it  
doesn't

seems to be the good one. Can anybody help me?
Many thanks!


The way to address this, if you are committed to this path, is to  
first determine the class of the fit-object and then to look for a  
plot method with methods(plot). If you see an S3 method you can call  
up the code with:


evd:::plot.fit-class  # when "fit-class" is the value you got with  
class(fit)


If it's an S4 method, then it's much more convoluted, and over time  
I've learned not to try.


--

David Winsemius, MD
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] Am I misunderstanding loop variable assignment or how to use print()?

2011-12-15 Thread Sarah Goslee
Hi,

An anova object is stored in anova.ag.m2529.az.

But "anova.ag.m2529.az"  is a character string that happens to be the
*name* of an anova object, but R has no way to know that unless you
specifically tell it that your character string is an object by using
get().

Something like print(get(x)) would work.

It's often neater and more efficient to store your anova objects in a
list, though.

Sarah

On Thu, Dec 15, 2011 at 9:30 AM, Tony Stocker  wrote:
> Given this interactive session:
>
>> an<-ls(pat="anova.ag.m2529")
>> an
>  [1] "anova.ag.m2529.az"   "anova.ag.m2529.can"   "anova.ag.m2529.fl"
>> print(anova.ag.m2529.az)
> Analysis of Variance Table
>
> Response: year
>                 Df Sum Sq Mean Sq F value   Pr(>F)
> time             1 14.823  14.8235    10.598 0.004392 **
> Residuals   18  25.177 1.3987
> ---
> Signif. codes:   0 '***'' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> How come the following just prints the names and not the contents of
> the existing variable/data frame with the name assigned to x?  Does it
> need to be dereferenced in some way to indicate that it's not a string
> but an actual variable living inside the workspace?:
>
>> for(x in an) {print(x)}
> [1] "anova.ag.m2529.az"
> [1] "anova.ag.m2529.can"
> [1] "anova.ag.m2529.fl"
>
> I'm trying to get this working interactively first, but ultimately I
> need to put it into an Rscript so that I can use dynamic listings and
> for loops to print everything.  I'd also be happy if I could pull just
> the Pr out, as I could interactively like so:
>
>> print(anova.ag.m2529.az$Pr)
> [1] 0.004392059               NA
>
> So what am I misunderstanding about how the language works?  I've been
> all over the web and through the usingR.pdf file but can't find an
> example that shows something like what I'm trying to do.  What exactly
> (data frame, table, function, character string, etc.) is stored in
> 'anova.ag.m2529.az' if the commands that created it were:
>
>> aov.m2529.az=aov(year~time,data=ag.m2529.az)
>> anova.ag.m2529.az=anova(aov.m2529.az)
>
> Much thanks in advance.
>


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


[R] Am I misunderstanding loop variable assignment or how to use print()?

2011-12-15 Thread Tony Stocker
Given this interactive session:

> an<-ls(pat="anova.ag.m2529")
> an
 [1] "anova.ag.m2529.az"   "anova.ag.m2529.can"   "anova.ag.m2529.fl"
> print(anova.ag.m2529.az)
Analysis of Variance Table

Response: year
 Df Sum Sq Mean Sq F value   Pr(>F)
time 1 14.823  14.823510.598 0.004392 **
Residuals   18  25.177   1.3987
---
Signif. codes:   0 '***'' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

How come the following just prints the names and not the contents of
the existing variable/data frame with the name assigned to x?  Does it
need to be dereferenced in some way to indicate that it's not a string
but an actual variable living inside the workspace?:

> for(x in an) {print(x)}
[1] "anova.ag.m2529.az"
[1] "anova.ag.m2529.can"
[1] "anova.ag.m2529.fl"

I'm trying to get this working interactively first, but ultimately I
need to put it into an Rscript so that I can use dynamic listings and
for loops to print everything.  I'd also be happy if I could pull just
the Pr out, as I could interactively like so:

> print(anova.ag.m2529.az$Pr)
[1] 0.004392059   NA

So what am I misunderstanding about how the language works?  I've been
all over the web and through the usingR.pdf file but can't find an
example that shows something like what I'm trying to do.  What exactly
(data frame, table, function, character string, etc.) is stored in
'anova.ag.m2529.az' if the commands that created it were:

> aov.m2529.az=aov(year~time,data=ag.m2529.az)
> anova.ag.m2529.az=anova(aov.m2529.az)

Much thanks in advance.

__
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] lm and R-squared (newbie)

2011-12-15 Thread David Winsemius


On Dec 15, 2011, at 8:35 AM, PtitBleu wrote:


Hello,

I've two data.frames (data1 and data4), dec="." and sep=";".
http://r.789695.n4.nabble.com/file/n4199964/data1.txt data1.txt
http://r.789695.n4.nabble.com/file/n4199964/data4.txt data4.txt

When I do
plot(data1$nx,data1$ny, col="red")
points(data4$nx,data4$ny, col="blue")
,  results seem very similar (at least to me) but the R-squared of
summary(lm(data1$ny ~ data1$nx))
and
summary(lm(data4$ny ~ data4$nx))
are very different (0.48 against 0.89).

Could someone explain me the reason?


Because you failed to do an adequate assessment of your data. Try this  
ploting exercsie and I think you will see the reason for the  
differences:


plot(data1$nx,data1$ny, col="red", xlim=range(c(data1$nx,data4$nx)),  
ylim=range(c(data1$ny,data4$ny)) )


--
David.



To be complete, I am looking for an simple indicator telling me if  
it is
worthwhile to keep the values provided by lm. I thought that R- 
squared could
do the job. For me, if R-squared is far from 1, the data are not  
good enough

to perform a linear fit.
It seems that I'm wrong.

Thanks for your explainations.
Ptit Bleu.





--
View this message in context: 
http://r.789695.n4.nabble.com/lm-and-R-squared-newbie-tp4199964p4199964.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
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] Calculate AUC Using the Trapezoidal Method

2011-12-15 Thread arivald
Hello,
I want to use function AUC {MIfuns} but I have problem with arguments.
My data:
 Fraction Time Variable
C1 0  0.0
C1   15  20.95475
C130 28.55030
C160 36.33064
C190 48.80438
C1  120 60.18636
 
AUC(data, time = Time, id = Fraction, dv = Variable)  ##this is not working!



Thx!

--
View this message in context: 
http://r.789695.n4.nabble.com/Calculate-AUC-Using-the-Trapezoidal-Method-tp4200049p4200049.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] Load Libraries from list

2011-12-15 Thread Barry Rowlingson
On Thu, Dec 15, 2011 at 1:19 PM, anaraster  wrote:
> Hi
>
> How can I load libraries from a list containing the library names?
>
> Something like this
>
> ListOfLibraries=c("foreign","survival")
>
> for (in in 1:length(ListOfLibraries)){library(as.name(ListOfLibraries[i]))}
>
> I already tried this  with no results.

Read the help for library and you'll see there's an option to help
you. Also, use lapply to iterate over lists. Also also, they're
PACKAGES!

lapply(ListOfLibraries,function(x){library(x,character.only=TRUE)})

Barry

__
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] lm and R-squared (newbie)

2011-12-15 Thread Gabor Grothendieck
On Thu, Dec 15, 2011 at 8:35 AM, PtitBleu  wrote:
> Hello,
>
> I've two data.frames (data1 and data4), dec="." and sep=";".
> http://r.789695.n4.nabble.com/file/n4199964/data1.txt data1.txt
> http://r.789695.n4.nabble.com/file/n4199964/data4.txt data4.txt
>
> When I do
> plot(data1$nx,data1$ny, col="red")
> points(data4$nx,data4$ny, col="blue")
> ,  results seem very similar (at least to me) but the R-squared of
> summary(lm(data1$ny ~ data1$nx))
> and
> summary(lm(data4$ny ~ data4$nx))
> are very different (0.48 against 0.89).
>
> Could someone explain me the reason?
>
> To be complete, I am looking for an simple indicator telling me if it is
> worthwhile to keep the values provided by lm. I thought that R-squared could
> do the job. For me, if R-squared is far from 1, the data are not good enough
> to perform a linear fit.
> It seems that I'm wrong.

The problem is the outliers. Try using a robust measure instead.  If
we replace Pearson correlations with Spearman (rank) correlations they
are much closer:

> # R^2 based on Pearson correlations
> cor(fitted(lm(ny ~ nx, data4)), data4$ny)^2
[1] 0.8916924
> cor(fitted(lm(ny ~ nx, data1)), data1$ny)^2
[1] 0.4868575
>
> # R^2 based on Spearman (rank) correlations
> cor(fitted(lm(ny ~ nx, data4)), data4$ny, method = "spearman")^2
[1] 0.8104026
> cor(fitted(lm(ny ~ nx, data1)), data1$ny, method = "spearman")^2
[1] 0.7266705

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.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] htest class members and print method

2011-12-15 Thread David L Lorenz
Rui,
  The answer to your last question is easy--you cannot add a new component 
to an object of class "htest" and have it printed by print.htest. But that 
does not mean that you cannot add a component to the output for your own 
use.
  You will need to decide what you want for output, both visually and in 
terms of ease of any post processing. A data.frame would be nice in that 
it would already for the table for any report, but it would require that 
the test you have designed have the same output for all orders. An 
alternative would be a list that contains the output from each test. 
Printing the list would create a nicely formatted report for each test.
Dave

> Hello,
>
> I am writing a test function, and would like to have it return an
> appropriate 'htest'
> object. This is very easy, but the test can be run for several orders, a
> frequent situation
> in time series tests.
> 
> It would be nice to return a data.frame with order, params, test 
statistic,
> p.value.
> 
> After seeing the 'htest' members in the help pages, I discovered that 
the
> class has a 'null.value' member that fits the job (it's what I'm using) 
but
> it's name may be misleading.
> 
> Is it possible to create a new member with a new, or at least different,
> name and have print.htest correctly print it?
> 
> 
> Rui Barradas
> 

[[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] lm and R-squared (newbie)

2011-12-15 Thread PtitBleu
Hello,

I've two data.frames (data1 and data4), dec="." and sep=";".
http://r.789695.n4.nabble.com/file/n4199964/data1.txt data1.txt 
http://r.789695.n4.nabble.com/file/n4199964/data4.txt data4.txt 

When I do
plot(data1$nx,data1$ny, col="red")
points(data4$nx,data4$ny, col="blue")
,  results seem very similar (at least to me) but the R-squared of
summary(lm(data1$ny ~ data1$nx))
and
summary(lm(data4$ny ~ data4$nx))
are very different (0.48 against 0.89).

Could someone explain me the reason?

To be complete, I am looking for an simple indicator telling me if it is
worthwhile to keep the values provided by lm. I thought that R-squared could
do the job. For me, if R-squared is far from 1, the data are not good enough
to perform a linear fit.
It seems that I'm wrong.

Thanks for your explainations.
Ptit Bleu.


 


--
View this message in context: 
http://r.789695.n4.nabble.com/lm-and-R-squared-newbie-tp4199964p4199964.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] corCompSymm in gamm()?

2011-12-15 Thread Anna Zakrisson
 Hi,

I have confirmed temporal correlation problems in my data. Is there a 
possibility to use corCompSymm for a gamm()?
I am an R-beginner.

I have very short time series. There are three years and within each year, 
there are 10 weeks. he 10 weeks are the same every year and have not unique 
values, I seem not to be able to use AR-1 (I assume that I have too little 
data for autoregression models of higher orders (ARMA)).  If I rename weeks 
as week 1-week 30 to get unique identifiers, I loose the seasonal and 
year-specific effect in the final model.

  M2c.gamm <- gamm(het ~ s(LN.DIN, k=3) + s(LN.totn, k=3) + s(Ncell, 
k=3) + s(LN.biom) + s(temp, k=3) + s(week, k=3)
 + fstation + fyear, method = "ML",
 weights = varIdent(form=~1 | fstation),
 data = data1,
 correlation = corCompSymm(form = ~ week|year))#seems not to work in 
a gamm()

Thank you for your time!
Anna Zakrisson Braeunlich
PhD Student

Department of Systems Ecology
Stockholm University
Svante Arrheniusv. 21A
SE-106 91 Stockholm

E-mail: a...@ecology.su.se
Tel work: +46 (0)8 161103
Mobile: +46-(0)700-525015
Web site: http://www.ecology.su.se/staff/personal.asp?id=163

><º>`•. . • `•. .• `•. . ><º>`•. . • `•. .• 
`•. .><º>

[[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] Load Libraries from list

2011-12-15 Thread anaraster
Hi

How can I load libraries from a list containing the library names?

Something like this

ListOfLibraries=c("foreign","survival")

for (in in 1:length(ListOfLibraries)){library(as.name(ListOfLibraries[i]))}

I already tried this  with no results.


Any help appreciated. 

Thanks

--
View this message in context: 
http://r.789695.n4.nabble.com/Load-Libraries-from-list-tp4199903p4199903.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] How to open files that contain "0"

2011-12-15 Thread Sarah Goslee
It doesn't have anything to do with 0 values.

read.table() is treating everything after the # (the default comment
character) as a comment.

Using comment.char = "" in read.table() will change that behavior.

Special characters in files cause all sorts of issues.

Sarah

On Thu, Dec 15, 2011 at 7:03 AM, Neotropical bat risk assessments
 wrote:
> Hi all,
>
> How can I set open files that contain values of Zero =0?
>
> These are valid values for the parameters I need to evaluate.
> I have tried CSV and tab formats.
> Trying XL Connect and/or XLConnectJars dies not seem to work to open
> Excel files so I am at a loss on how to get the data into a DF.
>
>
> Sample of data with 0 values:
>
> Filename        Dur     TBC     Fmax    Fmin    Fmean   Fc      S1      Sc    
>   Pmc
> g8221843.13#    5.06    0       38.93   36.2    37.96   36.45   -34.08  
> 192.69  6.8
> g8221843.13#    0.41    5.29    38.83   36.04   38.83   38.83   -261.93       
>   -513.05         0
> g8221843.13#    0.66    0.68    35.71   33.4    36.42   35.63   -238.04       
>   -392.06         0.2
> g8221843.13#    0.58    54.84   42.78   40.3    41.1    40.3    410     0     
>   6.2
>
>
> This message is displayed in the console:
>
>  > BASCTWS <- read.table("C:/=Bat data
> working/R/BASCTWS.txt",header=T,sep="\t",quote="")
> Error in scan(file, what, nmax, sep, dec, quote, skip, nlines,
> na.strings,  :
>   line 1 did not have 10 elements
>  >
> Obviously I am missing something basic, seems I was able to do this in
> the past.
>
> Bruce
>
-- 
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.


[R] How to open files that contain "0"

2011-12-15 Thread Neotropical bat risk assessments
Hi all,

How can I set open files that contain values of Zero =0?

These are valid values for the parameters I need to evaluate.
I have tried CSV and tab formats.
Trying XL Connect and/or XLConnectJars dies not seem to work to open 
Excel files so I am at a loss on how to get the data into a DF.


Sample of data with 0 values:

FilenameDur TBC FmaxFminFmean   Fc  S1  Sc  
Pmc
g8221843.13#5.060   38.93   36.237.96   36.45   -34.08  192.69  
6.8
g8221843.13#0.415.2938.83   36.04   38.83   38.83   -261.93 
-513.05 0
g8221843.13#0.660.6835.71   33.436.42   35.63   -238.04 
-392.06 0.2
g8221843.13#0.5854.84   42.78   40.341.140.3410 0   
6.2


This message is displayed in the console:

 > BASCTWS <- read.table("C:/=Bat data 
working/R/BASCTWS.txt",header=T,sep="\t",quote="")
Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, 
na.strings,  :
   line 1 did not have 10 elements
 >
Obviously I am missing something basic, seems I was able to do this in 
the past.

Bruce

-- 
Bruce W. Miller, Ph.D.
Conservation Ecologist
Neotropical Bat Project


office details
Gallon Jug, Belize
Mailing address
P.O. Box 37, Belize City
Belize, Central America
Phone +501-220-9002



-- 
Bruce W. Miller, PhD.
Neotropical bat risk assessments

If we lose the bats, we may lose much of the tropical vegetation and the lungs 
of the planet

Using acoustic sampling to map species distributions for>15 years.

Providing Interactive identification keys to the vocal signatures of New World 
Bats

For various project details see:

https://sites.google.com/site/batsoundservices/


[[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 draw random numbers from many categorical distributions quickly?

2011-12-15 Thread Nordlund, Dan (DSHS/RDA)
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
> project.org] On Behalf Of Sean Zhang
> Sent: Wednesday, December 14, 2011 10:07 PM
> To: r-help@r-project.org
> Subject: [R] how to draw random numbers from many categorical
> distributions quickly?
> 
> Dear R helpers,
> 
> I have a question about drawing random numbers from many categorical
> distributions.
> 
> Consider n individuals, each follows a categorical distribution defined
> over k categories.
> Consider a simple case in which n=4, k=3 as below
> 
> catDisMat <-
> rbind(c(0.1,0.2,0.7),c(0.2,0.2,0.6),c(0.1,0.2,0.7),c(0.1,0.2,0.7))
> 
> outVec <- rep(NA,nrow(catDisMat))
> for (i in 1:nrow(catDisMat)){
> outVec[i] <- sample(1:3,1, prob=catDisMat[i,], replace = TRUE)
> }
> 
> I can think of one way to potentially speed it up (in reality, my n is
> very
> large, so speed matters). The approach above only samples 1 value each
> time. I could have sampled two values for c(0.1,0.2,0.7) because it
> appears
> three times. so by doing some manipulation, I think I can have the
> idea,
> "sample(1:3, 3, prob=c(0.1,0.2,0.7), replace = TRUE)",  implemented to
> improve speed a bit. But, I wonder whether there is a better approach
> for
> speed?
> 
> Thanks in advance.
> 
> -Sean
> 

Sean,

How about something like this:

outVec <- apply(catDisMat,1, function(x)sample(1:3, 1, prob = x, replace = 
TRUE))


I created a catDisMat matrix with a million rows and apply crunched through it 
in approximately 8-9 seconds on my 2.67 GHz 64-bit Windows 7 box with 12 GB of 
ram.  Your code above was substantially slower. 

Hope this is helpful,

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204


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