[R] Extraction of rules from Support Vector Machines
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
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
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
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
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
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
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?
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()
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)
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
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
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
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
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
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
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)
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
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
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"
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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?
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
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
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
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?
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
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
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
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()?
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()?
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()?
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
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
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
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
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()?
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()?
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)
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
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
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)
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
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)
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()?
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
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"
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"
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?
> -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.