Re: [R] Trellis.par.set/ family/ global change font?
A radical solution is to edit the file .../etc/Rdevga and to change the definition of the first set of fonts, e.g. # TT Arial : plain # TT Arial : bold # TT Arial : italic # TT Arial : bolditalic TT Times New Roman : plain TT Times New Roman : bold TT Times New Roman : italic TT Times New Roman : bolditalic Best, Renaud 2006/6/30, Deepayan Sarkar [EMAIL PROTECTED]: On 6/29/06, McClatchie, Sam (PIRSA-SARDI) [EMAIL PROTECTED] wrote: Background: OS: Linux Ubuntu Dapper 6.06 release: R 2.3.1 editor: GNU Emacs 21.4.1 front-end: ESS 5.2.3 - Colleagues I have a rather complicated trellis plot that a journal editor has requested I edit and change all the fonts to times. I'd like to change all fonts globally for the plot, as in par(family=serif) for non-trellis plots. Various experiments with trellis.par.set after reading the help page have not solved the problem for me. Even doing the local change trellis.par.set(par.xlab.text=list(cex=1.5, family=serif)) does not change the font to times for xlab (I mean my syntax is wrong, not that there is a bug). So I'm obviously misreading the help page or just missing the meaning. Any suggestions? Lattice uses 'fontfamily' rather than 'family' (borrowed from grid, I suppose). I don't think there's a way to set the family globally. You might try trellis.par.set(grid.pars = list(fontfamily = serif)) but I'm not sure if that will work. -Deepayan __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Renaud LANCELOT Département Elevage et Médecine Vétérinaire (EMVT) du CIRAD Directeur adjoint chargé des affaires scientifiques CIRAD, Animal Production and Veterinary Medicine Department Deputy director for scientific affairs Campus international de Baillarguet TA 30 / B (Bât. B, Bur. 214) 34398 Montpellier Cedex 5 - France Tél +33 (0)4 67 59 37 17 Secr. +33 (0)4 67 59 39 04 Fax +33 (0)4 67 59 37 95 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Trellis.par.set/ family/ global change font?
On Fri, 30 Jun 2006, Renaud Lancelot wrote: A radical solution is to edit the file .../etc/Rdevga and to change the definition of the first set of fonts, e.g. # TT Arial : plain # TT Arial : bold # TT Arial : italic # TT Arial : bolditalic TT Times New Roman : plain TT Times New Roman : bold TT Times New Roman : italic TT Times New Roman : bolditalic He is on Linux, so that will not work. We have not been told the graphics device, but both postscript() and pdf() allow the device to be set when the device is opened, and that works for me: pdf(family=serif) example(xyplot) dev.off() Best, Renaud 2006/6/30, Deepayan Sarkar [EMAIL PROTECTED]: On 6/29/06, McClatchie, Sam (PIRSA-SARDI) [EMAIL PROTECTED] wrote: Background: OS: Linux Ubuntu Dapper 6.06 release: R 2.3.1 editor: GNU Emacs 21.4.1 front-end: ESS 5.2.3 - Colleagues I have a rather complicated trellis plot that a journal editor has requested I edit and change all the fonts to times. I'd like to change all fonts globally for the plot, as in par(family=serif) for non-trellis plots. Various experiments with trellis.par.set after reading the help page have not solved the problem for me. Even doing the local change trellis.par.set(par.xlab.text=list(cex=1.5, family=serif)) does not change the font to times for xlab (I mean my syntax is wrong, not that there is a bug). So I'm obviously misreading the help page or just missing the meaning. Any suggestions? Lattice uses 'fontfamily' rather than 'family' (borrowed from grid, I suppose). I don't think there's a way to set the family globally. You might try trellis.par.set(grid.pars = list(fontfamily = serif)) but I'm not sure if that will work. -Deepayan __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] cat and positioning of the output
Thanks at all for your carefully help. I decided to used both version of your proposals and to combine with a request of the operating system. All the best - J?rn Schulz. -- View this message in context: http://www.nabble.com/cat-and-positioning-of-the-output-tf1869521.html#a5115737 Sent from the R help forum at Nabble.com. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Empirical CDF
[EMAIL PROTECTED] writes: Good day everyone, I want to assess the error when fitting a Gram-Charlier CDF to some data 'ws', that is, I want to calculate: Err = |ecdf(ws) - GCh_ser(ws)| The problem is, I cannot get the F(x) values from the ecdf. 'Summary(ecdf())' returns some of the x-axis values, but how do you get the F(x) values? By applying the function returned by ecdf? Fn - ecdf(some.data) Fn(some.other.data) should work. Thank you for any help you can provide. Regards, Augusto Augusto Sanabria. MSc, PhD. Mathematical Modeller Risk Research Group Geospatial Earth Monitoring Division Geoscience Australia (www.ga.gov.au) Cnr. Jerrabomberra Av. Hindmarsh Dr. Symonston ACT 2609 Ph. (02) 6249-9155 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Cointegration Test in R
Hello Dennis, have you considered the function bh6lrtest() in the package urca? To my knowledge, there is no other package available that offers VECM functionalities. Best, Bernhard ps: As you migth be interested in VAR and SVAR too: I am currently working on such a package which should be submitted to CRAN after summer. Dr. Bernhard Pfaff Global Structured Products Group (Europe) Invesco Asset Management Deutschland GmbH Bleichstrasse 60-62 D-60313 Frankfurt am Main Tel: +49(0)69 29807 230 Fax: +49(0)69 29807 178 Email: [EMAIL PROTECTED] -Ursprüngliche Nachricht- Von: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Im Auftrag von Dennis Heidemann Gesendet: Donnerstag, 29. Juni 2006 21:20 An: R-help@stat.math.ethz.ch Betreff: [R] Cointegration Test in R Hello! I'm using the blrtest() function in the urca package to test cointegration relationships. Unfortunately, the hypothesis (restrictions on beta) specifies the same restriction on all cointegration vectors. Is there any possibility to specify different restrictions on the cointegration vectors? Are there any other packages in R using cointegration tests? Thanks and best regards. Dennis Heidemann [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html * Confidentiality Note: The information contained in this mess...{{dropped}} __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R server
There is the Rserve package which you can look at here: http://stats.math.uni-augsburg.de/Rserve/down.shtml JS --- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] smoothing splines and degrees of freedom
Steven, I cannot vouch for the behaviour of the function smooth.spline(), but the theoretical answer to your question is yes. If g = Sy is the transformation from data vector y to spline vector g, the equivalent degrees of freedom are usually defined as EDF = trace(S), where S is the n x n smoothing matrix: EDF = sum_i(1/(1+theta*lambda_i)), where lambda_1 to lambda_n are the eigenvalues of S. Two of these are zero, so EDF = 2 + sum(1/(1+theta*lambda_i)) the sum now over i=3 to n. Here theta is the smoothing parameter. Setting theta = 0 (no smoothing) gives EDF=n and produces the interpolating spline. Setting theta = infty gives EDF=2 and a straight line fit. See either Green and Silverman, Nonparametric regression and generalized linear models, (p37), or Hastie and Tibshirani, Generalized additive models, p52. On Sat, Jun 24, 2006 at 11:35:16AM -0400, Steven Shechter wrote: Hi, If I set df=2 in my smooth.spline function, is that equivalent to running a linear regression through my data? It appears that df=# of data points gives the interpolating spline and that df = 2 gives the linear regression, but I just want to confirm this. Thank you, Steven __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- *I.White * *University of Edinburgh * *Ashworth Laboratories, West Mains Road* *Edinburgh EH9 3JT * *Fax: 0131 650 6564 Tel: 0131 650 5490 * *E-mail: [EMAIL PROTECTED] * __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R server
Mike, may also want to check out Rpad. Its mainly used to delivery R programs to other users via a web browser, so the other users don't need R installed on their machine. If you make a nice HTML gui, the other users don't even have to know R. Each user gets his own environment, so one person's variables won't interfer with anothers. HTH, Roger On 6/30/06, john seers (IFR) [EMAIL PROTECTED] wrote: There is the Rserve package which you can look at here: http://stats.math.uni-augsburg.de/Rserve/down.shtml JS --- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] data extraction
Dear mailing list I have a data that have 20,000 rows and 20 columns. Io wonted to extract the 10th row only. Example the 10th, 20th, 30th 40th ..2 th. can you please help me how do I do that.Than kyou. Example is below. Inpute: AG GG GG AG CC CC CC CC CT CC CT CT GG GG GG GG CC CC CC CC GG GG GG GG CC CC CC CC GG CG CG GG GG GG GG GG *CC CC CC CC* AA AG AG AA AA AA AA AA GG AG AG GG GG AG AG GG GG GG GG GG TT TT TT TT AA AG AG AA CT CC CT TT AG AA AG GG *AA AA AG AG* AA AA AA AA CC CC CC CG GG GG GG AG CT TT CT CT AT AA AT AT GG GG GG AG CG CC CG GG GG GG GG AG CC CC CC CT GG GG GG GG *GG GG AG AG *CC CC CC CT TT TT TT CT AG GG AG AG GG GG GG GG AG AG AA AA AG GG AG AA CT TT CT CT CT CT CC CT *AC CC AC AC* output *CC CC CC CC AA AA AG AG GG GG GG GG AC CC AC AC* thankyou a lot inadvance! [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] data extraction
This should work: data[seq(1,nrow(data),10),] Andris On Piektdiena, 30. Jūnijs 2006 14:45, yohannes alazar wrote: Dear mailing list I have a data that have 20,000 rows and 20 columns. Io wonted to extract the 10th row only. Example the 10th, 20th, 30th 40th�..2 th. can you please help me how do I do that.Than kyou. Example is below. Inpute: AG GG GG AG CC CC CC CC CT CC CT CT GG GG GG GG CC CC CC CC GG GG GG GG CC CC CC CC GG CG CG GG GG GG GG GG *CC CC CC CC* AA AG AG AA AA AA AA AA GG AG AG GG GG AG AG GG GG GG GG GG TT TT TT TT AA AG AG AA CT CC CT TT AG AA AG GG *AA AA AG AG* AA AA AA AA CC CC CC CG GG GG GG AG CT TT CT CT AT AA AT AT GG GG GG AG CG CC CG GG GG GG GG AG CC CC CC CT GG GG GG GG *GG GG AG AG *CC CC CC CT TT TT TT CT AG GG AG AG GG GG GG GG AG AG AA AA AG GG AG AA CT TT CT CT CT CT CC CT *AC CC AC AC* output *CC CC CC CC AA AA AG AG GG GG GG GG AC CC AC AC* thankyou a lot inadvance! [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Extracting R plots from MS Word
Just a note to say what I did. I think that the results were OK but I have yet to hear from the journal. 1. I saved the Word document under another name. 2. I deleted all the contents of the document except the target graphic. 3. I printed to file yielding a .prn file. 4. I changed the extension to .ps. 5. I used Gsview PS to EPS. Murray Jorgensen __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] data extraction
mydf[seq(10,2,10),] yohannes alazar wrote: Dear mailing list I have a data that have 20,000 rows and 20 columns. Io wonted to extract the 10th row only. Example the 10th, 20th, 30th 40th…..2 th. can you please help me how do I do that.Than kyou. Example is below. Inpute: AG GG GG AG CC CC CC CC CT CC CT CT GG GG GG GG CC CC CC CC GG GG GG GG CC CC CC CC GG CG CG GG GG GG GG GG *CC CC CC CC* AA AG AG AA AA AA AA AA GG AG AG GG GG AG AG GG GG GG GG GG TT TT TT TT AA AG AG AA CT CC CT TT AG AA AG GG *AA AA AG AG* AA AA AA AA CC CC CC CG GG GG GG AG CT TT CT CT AT AA AT AT GG GG GG AG CG CC CG GG GG GG GG AG CC CC CC CT GG GG GG GG *GG GG AG AG *CC CC CC CT TT TT TT CT AG GG AG AG GG GG GG GG AG AG AA AA AG GG AG AA CT TT CT CT CT CT CC CT *AC CC AC AC* output *CC CC CC CC AA AA AG AG GG GG GG GG AC CC AC AC* thankyou a lot inadvance! [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] data extraction
you need something like: new.data - data[seq(10, nrow(data), 10), ] I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm - Original Message - From: yohannes alazar [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Friday, June 30, 2006 1:45 PM Subject: [R] data extraction Dear mailing list I have a data that have 20,000 rows and 20 columns. Io wonted to extract the 10th row only. Example the 10th, 20th, 30th 40th...2 th. can you please help me how do I do that.Than kyou. Example is below. Inpute: AG GG GG AG CC CC CC CC CT CC CT CT GG GG GG GG CC CC CC CC GG GG GG GG CC CC CC CC GG CG CG GG GG GG GG GG *CC CC CC CC* AA AG AG AA AA AA AA AA GG AG AG GG GG AG AG GG GG GG GG GG TT TT TT TT AA AG AG AA CT CC CT TT AG AA AG GG *AA AA AG AG* AA AA AA AA CC CC CC CG GG GG GG AG CT TT CT CT AT AA AT AT GG GG GG AG CG CC CG GG GG GG GG AG CC CC CC CT GG GG GG GG *GG GG AG AG *CC CC CC CT TT TT TT CT AG GG AG AG GG GG GG GG AG AG AA AA AG GG AG AA CT TT CT CT CT CT CC CT *AC CC AC AC* output *CC CC CC CC AA AA AG AG GG GG GG GG AC CC AC AC* thankyou a lot inadvance! [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Biobass, SAGx, and Jonckheere-Terpstra test
Horace, Last year I found another Jonckheere-Terpstra function at http://www.biostat.wustl.edu/archives/html/s-news/2000-10/msg00126.html Using that and the version in SAGx and wanting a few other options, I created the version below. Hope it is useful, Chris -- Christopher Andrews, PhD SUNY Buffalo, Department of Biostatistics # # # j.test - function(x, y, alternative = c(two.sided, decreasing, increasing), asymp=TRUE, correct=FALSE, perm=0, na.action=c(omit, fail), permgraph=FALSE, permreps=FALSE, ...) { #Jonckheere-Terpstra test #x = response #y = group membership #alternative hypothesis #asymptotic = asymptotic formula for variance or don't bother #correct = continuity correction #perm = number of repetitions for a permutation test #na.action = what to do with missing data #permgraph = draw a histogram of the permutations? #permreps = return the permutations? #... for graph alternative - match.arg(alternative) complete - complete.cases(x,y) nn - sum(complete) cat(nn, complete cases.\n) if (!all(complete) (na.action==fail)) { cat(option na.action is 'fail' and, sum(!complete), cases are not complete.\n) return(NULL) } response - x[complete] groupvec - y[complete, drop=TRUE] if(!is.ordered(groupvec)) { groupvec - as.ordered(groupvec) cat(y was not an ordered factor. Redefined to be one.\n) } lev - levels(groupvec) nlev - length(lev) if(nlev = 2) { if (nlev 2) { cat(Not enough groups (k =, nlev, ) for Jonckheere.\n) return(NULL) } else { cat(Two groups. Could use wilcox.test.\n) } } computestat - function(response, groupvec, n.groups) { pairwise - function(x, y) { length(x)*(length(y) + (1+length(x))/2) - sum(rank(c(x,y))[seq(along=x)]) } H-0 for (i in seq(n.groups-1)) for (j in seq(i+1, n.groups)) H - H+pairwise(response[groupvec == lev[i]], response[groupvec == lev[j]]) return(H) } retval - list(statistic=computestat(response, groupvec, nlev), alternative = paste(alternative, paste(levels(groupvec),collapse=switch(alternative, two.sided = , , decreasing= , increasing = )), sep=: ), method = Jonckheere, data.name = paste(deparse(substitute(x)), by, deparse(substitute(y if (asymp) { ns - tabulate(as.numeric(groupvec)) retval$EH - (nn * nn - sum(ns * ns))/4 retval$VH - if (!any(duplicated(response))) { (nn * nn * (2 * nn + 3) - sum(ns * ns * (2 * ns + 3)))/72 } else { ds - as.vector(table(response)) ((nn*(nn-1)*(2*nn+5) - sum(ns*(ns-1)*(2*ns+5)) - sum(ds*(ds-1)*(2*ds+5)))/72 + sum(ns*(ns-1)*(ns-2))*sum(ds*(ds-1)*(ds-2))/(36*nn*(nn-1)*(nn-2)) + sum(ns*(ns-1))*sum(ds*(ds-1))/(8*nn*(nn - 1))) } pp - if (!correct) { pnorm(retval$statistic, retval$EH, sqrt(retval$VH)) } else if (retval$statistic = retval$EH + 1) { pnorm(retval$statistic - 1, retval$EH, sqrt(retval$VH)) } else if (retval$statistic = retval$EH - 1) { pnorm(retval$statistic + 1, retval$EH, sqrt(retval$VH)) } else { .5 } retval$p.value - switch(alternative, two.sided = 2*min(pp,1-pp), decreasing = pp, increasing = 1-pp) } if (perm0) { reps - numeric(perm) for (i in seq(perm)) { reps[i] - computestat(response, sample(groupvec), nlev) } retval$EH.perm - mean(reps) retval$VH.perm - var(reps) ppp - if (!correct) { pnorm(retval$statistic, retval$EH.perm, sqrt(retval$VH.perm)) } else if (retval$statistic = retval$EH.perm + 1) { pnorm(retval$statistic - 1, retval$EH.perm, sqrt(retval$VH.perm)) } else if (retval$statistic = retval$EH.perm - 1) { pnorm(retval$statistic + 1, retval$EH.perm, sqrt(retval$VH.perm)) } else { .5 } retval$p.value.perm - switch(alternative, two.sided = 2*min(ppp,1-ppp), decreasing = ppp, increasing = 1-ppp) rr - rank(c(retval$statistic, reps))[1]/(perm+2) retval$p.value.rank - switch(alternative, two.sided = 2*min(rr,1-rr), decreasing = rr, increasing = 1-rr) if (permreps) retval$reps - reps if (permgraph) { hist(reps, xlab=Jonckheere Statistic, prob=TRUE, main=Histogram of Permutations, sub=paste(perm, permutations)) abline(v=retval$statistic, col=red) points((-3:3)*sqrt(retval$VH.perm)+retval$EH.perm, rep(0,7), col=blue, pch=as.character(c(3:1,0:3))) } } class(retval) - htest names(retval$statistic) - J return(retval) #returns list (of class htest) with the following components #Always: #statistic = value of J #alternative = same as input #method = Jonckheere #data.name = source of data based on call #Most of the time (i.e., when asymp=TRUE): #EH = expected value based on sample size #VH = variance (adjusted for ties if necessary) #p.value = asymptotic p value #Sometimes
[R] Customizing breaks for histograms of unequal ranges
Hi, I would really appreciate any suggestions on this (rather trivial) problem.. Say I have two vectors: v1=seq(1:10) v2=seq(1:15) For a set of common breaks I need to divide the density of v2 over v1. This means that I want to avoid having 0 counts for any v1 breakpoint. But (unsuprisingly) if I define my common breaks as those returned by calling hist for v1: v1.hist=hist(v1, plot=F) v2.hist=hist(v2, plot=F, breaks=v1.hist$breaks) I get an error: Error in hist.default(v1, plot = F, breaks = v2.hist$breaks) : some 'x' not counted; maybe 'breaks' do not span range of 'x' If I had used the combined vector (c(v1,v2)) to set my breaks I would end up with 0 for some v1 counts. So my question is: is there a way to define breaks that cover the whole range of v1 and v2 while avoiding having 0 for the shortest vector? Many thanks Eleni Rapsomaniki Birkbeck College, UK __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lme convergence
I think the following part of lme.formula if (numIter controlvals$maxIter) { stop(Maximum number of iterations reached without convergence.) } should be something like if (numIter controlvals$maxIter) { if (controlvals$returnObject) { warning(Maximum number of iterations reached without convergence.) break } else { stop(Maximum number of iterations reached without convergence.) } } Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm - Original Message - From: Spencer Graves [EMAIL PROTECTED] To: Ravi Varadhan [EMAIL PROTECTED] Cc: 'Pryseley Assam' [EMAIL PROTECTED]; 'R-Users' R-help@stat.math.ethz.ch; Douglas Bates [EMAIL PROTECTED] Sent: Friday, June 30, 2006 4:08 AM Subject: Re: [R] lme convergence Does anyone know how to obtain the 'returnObject' from an 'lme' run that fails to converge? An argument of this name is described on the 'lmeControl' help page as, a logical value indicating whether the fitted object should be returned when the maximum number of iterations is reached without convergence of the algorithm. Default is 'FALSE'. Unfortunately, I've so far been unable to get it to work, as witnessed by the following modification of an example from the '?lme' help page: library(nlme) fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1)) Error in lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1)) : iteration limit reached without convergence (9) fm1 Error: object fm1 not found fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1, + returnObject=TRUE)) Error in lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1, : iteration limit reached without convergence (9) fm1 Error: object fm1 not found I might be able to fix the problem myself, working through the 'lme' code line by line, e.g., using 'debug'. However, I'm not ready to do that just now. Best Wishes, Spencer Graves Ravi Varadhan wrote: Use try to capture error messages without breaking the loop. ?try -- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -- -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of Pryseley Assam Sent: Wednesday, June 28, 2006 12:18 PM To: R-Users Subject: [R] lme convergence Dear R-Users, Is it possible to get the covariance matrix from an lme model that did not converge ? I am doing a simulation which entails fitting linear mixed models, using a for loop. Within each loop, i generate a new data set and analyze it using a mixed model. The loop stops When the lme function does not converge for a simulated dataset. I want to inquire if there is a method to suppress the error message from the lme function, or better still, a way of going about this issue of the loop ending once the lme function does not converge. Thanks in advance, Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] GARCH message 75
For those who can find no documention on the @ operator, ?@ says it is used to Extract tbe contents of a slot in a object with a formal class structure. ?slot provides more information, including some critical references. The primary reference is Chambers (1998) Programming with Data (Springer). This describes the S 4 standard for object oriented programming in R. Additional information can be found in a chapter of Venables and Ripley (2000) S Programming (Springer). When I read Chambers a few years ago, I found it very difficult, partly because of the inconsistencies between what Chambers described and what was actually implemented in versions of S-Plus and to which I had access at the time. That may have improved since then. Hope this helps, Spencer Graves Joe Byers wrote: All, You might check out the @ operator on Classes and objects. I can find no documentation on this but if you look at code in fseries or fbasic in the method fARMA.summary. You will find where the object fit in the returned object is access using the @ operator. object - [EMAIL PROTECTED] ans$coefmat - cbind(format(object$coef,digits=digits), format(object$se.coef,digits=digits), format(tval,digits=digits), prob=format.pval(prob,digits=digits)) where x is the from x-armaFit(). This I believe would be the same for the GARCH results. Note the format.pval on the prob. This is important because if you do not do this you will get 0 for small pvalues. I learned this by looking at the code of the pretty printing functions so I can save only results that I want from multiple models runs. You also can then to wald test and LL ratio tests on the models. Good luck Joe __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Running R
Does any one have solution for R not running? I am getting the error message. Please Advice. Thanks 1. gunzip R-patched.tar.gz 2. tar -xvf R-patched.tar 3. changed the directory to the newly created directory R-patched 4. Typed ./configure 5. Typed make 6. Make check 7. make check-all 8. Typed make install *9. Typed R # R Fatal error: unable to open the base package __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] [R-pkgs] twang - Toolkit for Weighting and Analysis of Nonequivalent Groups
The Toolkit for Weighting and Analysis of Nonequivalent Groups (twang 1.0) has been released to CRAN. The package collects functions useful for computing propensity score weights for treatment effect estimation, developing nonresponse weights, and diagnosing the quality of those weights. The package includes a vignette containing some basic theory and walks through two examples. It is available by typing vignette(twang) at the R prompt after loading the library. Background on the methodology with an application to evaluating drug treatment programs is available in McCaffrey, D., G. Ridgeway, A. Morral (2004). Propensity score estimation with boosted regression for evaluating adolescent substance abuse treatment, Psychological Methods 9(4):403-425. This email message is for the sole use of the intended recip...{{dropped}} ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] [R-pkgs] gmp: bigintegers with matrix computation
Dear all, gmp package is now available on cran with version 0.4 Aim of gmp is to provide a lot a function to manipulate big integers, big rationals and last but not least big integers in Z/nZ (modulo n) We add in last version support for matrix computation (standards operators, + * - / %*% ...) and inversion matrix (solve): in Z/nZ if matrix defined in Z/nZ with rationals (absolute precision) if not. Best Regards, Antoine Lucas. ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] incomplete final line found by readLines on ...
Dear R-users I need to read some text files produced by some other software. I used readLines (with n = -1 ) command and then tried to find some numbers I liked to extract or some line numbers I like to know. The problem is that there is no empty line at the end of the text files. R gave me this warning below; In addition: Warning message: incomplete final line found by readLines on 'output.txt' I know this doesn't affect anything and just warning messages. Is there any way to prevent this warning message. Thanks Taka __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Query : Chi Square goodness of fit test
I want to calculate chi square test of goodness of fit to test, Sample coming from Poisson distribution. please copy this script in R run the script The R script is as follows ## start # No_of_Frouds- c(4,1,6,9,9,10,2,4,8,2,3,0,1,2,3,1,3,4,5,4,4,4,9,5,4,3,11,8,12,3,10,0,7) N - length(No_of_Frouds) # Estimation of Parameter lambda- sum(No_of_Frouds)/N lambda pmf - dpois(i, lambda, log = FALSE) step_function - ppois(i, lambda, lower.tail = TRUE, log.p = FALSE) # Chi-Squared Goodness of Fit Test # Ho: The data follow a Poisson distribution Vs H1: Not Ho Frauds - c(1:13) counts- c(2,3,3,5,7,2,1,1,2,3,2,1,1,0) # Observed frequency Expected -c(0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.7817821 03,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.2344006 79,0.095299266,0.035764993) chisq.test(counts, Expected, simulate.p.value =FALSE, correct = FALSE) # end The result of R is as follows Pearson's Chi-squared test data: counts and poisson_fit X-squared = 70, df = 65, p-value = 0.3135 Warning message: Chi-squared approximation may be incorrect in: chisq.test(counts, poisson_fit, simulate.p.value = FALSE, correct = FALSE) But I have done calculations in Excel. I am getting different answer. Observed = 2,3,3,5,7,2,1,1,2,3,2,1,1,0 Expected=0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.78 1782103,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.23 4400679,0.095299266,0.035764993 Estimated Parameter =4.878788 Chi square stat = 0.000113 My excel answer tally with the book which I have refer for excel. Please tell me the correct calculations in R. ## Awaiting your positive reply. Regards. Priti. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Running R
Have you checked and see if there is any error in steps 5 through 8? If you can't start R, I doubt steps 6 and 7 ran fine. Andy From: Pramod Anugu Does any one have solution for R not running? I am getting the error message. Please Advice. Thanks 1. gunzip R-patched.tar.gz 2. tar -xvf R-patched.tar 3. changed the directory to the newly created directory R-patched 4. Typed ./configure 5. Typed make 6. Make check 7. make check-all 8. Typed make install *9. Typed R # R Fatal error: unable to open the base package __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Running R
On Fri, Jun 30, 2006 at 07:56:41AM -0500, Pramod Anugu wrote: Does any one have solution for R not running? I am getting the error message. Please Advice. Yes, you have emailing the group fairly regularly about your troubles. And you have been told just about each and every time to *check the output of the configure step*. Unless that step concludes with all relevant pieces found, *there is no point in resuming with steps 5 to 9*. As you are emailing from an educational institution, consider getting local help about configuring and compiling software. Building R is *extensively documented in a dedicated manual* but the manual may need someone with more familiarity with the process than you currently have. Good luck, and consider to stop sending the *same* email to the list every couple of days. Dirk Thanks ? ? 1. gunzip R-patched.tar.gz 2. tar -xvf R-patched.tar ?? 3. changed the directory to the newly created directory R-patched ?? 4. Typed ./configure ? ?5. Typed make ?? 6. Make check ?? 7. make check-all ?? 8. Typed make install ?? *9. Typed R # R Fatal error: unable to open the base package __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Hell, there are no rules here - we're trying to accomplish something. -- Thomas A. Edison __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Random numbers from noncentral t-distribution
Hi there: I'd thought these two versions of noncentral t-distribution are essentially the same: qqplot(rt(1000,df=20,ncp=3),qt(runif(1000),df=20,ncp=3)) But, the scales of the x-axis and the y-axis are quite different according to the QQ-plot. Did I make any mistakes somewhere? Thanks, Long - Mp3·è¿ñËÑ-иèÈȸè¸ßËÙÏ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] incomplete final line found by readLines on ...
[Taka Matzmoto] Is there any way to prevent [this] warning message. Hi, Taka. The easiest might be using the suppressWarnings wrapper. See ?suppressWarnings for more information. -- François Pinard http://pinard.progiciels-bpi.ca __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Query : Chi Square goodness of fit test
chisq.test(counts, p=Expected/sum(Expected), simulate.p.value =FALSE, correct = FALSE) Chi-squared test for given probabilities data: counts X-squared = 40.5207, df = 13, p-value = 0.0001139 Warning message: l'approximation du Chi-2 est peut-être incorrecte in: chisq.test(counts, p = Expected/sum(Expected), simulate.p.value = FALSE, but the use of Chi2 test is incorrect since some of Expected frequencies are lower than 5. --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- priti desai a écrit : I want to calculate chi square test of goodness of fit to test, Sample coming from Poisson distribution. please copy this script in R run the script The R script is as follows ## start # No_of_Frouds- c(4,1,6,9,9,10,2,4,8,2,3,0,1,2,3,1,3,4,5,4,4,4,9,5,4,3,11,8,12,3,10,0,7) N - length(No_of_Frouds) # Estimation of Parameter lambda- sum(No_of_Frouds)/N lambda pmf - dpois(i, lambda, log = FALSE) step_function - ppois(i, lambda, lower.tail = TRUE, log.p = FALSE) # Chi-Squared Goodness of Fit Test # Ho: The data follow a Poisson distribution Vs H1: Not Ho Frauds - c(1:13) counts- c(2,3,3,5,7,2,1,1,2,3,2,1,1,0) # Observed frequency Expected -c(0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.7817821 03,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.2344006 79,0.095299266,0.035764993) chisq.test(counts, Expected, simulate.p.value =FALSE, correct = FALSE) # end The result of R is as follows Pearson's Chi-squared test data: counts and poisson_fit X-squared = 70, df = 65, p-value = 0.3135 Warning message: Chi-squared approximation may be incorrect in: chisq.test(counts, poisson_fit, simulate.p.value = FALSE, correct = FALSE) But I have done calculations in Excel. I am getting different answer. Observed = 2,3,3,5,7,2,1,1,2,3,2,1,1,0 Expected=0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.78 1782103,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.23 4400679,0.095299266,0.035764993 Estimated Parameter =4.878788 Chi square stat = 0.000113 My excel answer tally with the book which I have refer for excel. Please tell me the correct calculations in R. ## Awaiting your positive reply. Regards. Priti. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Way to convert data frame to matrix
I have a text file that I have imported into R. It contains 3 columns and 316940 rows. The first column is vegetation plot ID, the second species names and the third is a cover value (numeric). I imported using the read.table function. My problem is this. I need to reformat the information as a matrix, with the first column becoming the row labels and the second the column labels and the cover values as the matrix cell data. However, since the read.tablefunction imported the data as an indexed data frame, I can't use the columns as vectors. Is there a way around this, to convert the data frame as 3 separate vectors? I have been looking all over for a function, and my programming skills are not great. Thanks in advance [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Way to convert data frame to matrix
?as.matrix On Friday 30 June 2006 09:11, Wade Wall wrote: I have a text file that I have imported into R. It contains 3 columns and 316940 rows. The first column is vegetation plot ID, the second species names and the third is a cover value (numeric). I imported using the read.table function. My problem is this. I need to reformat the information as a matrix, with the first column becoming the row labels and the second the column labels and the cover values as the matrix cell data. However, since the read.tablefunction imported the data as an indexed data frame, I can't use the columns as vectors. Is there a way around this, to convert the data frame as 3 separate vectors? I have been looking all over for a function, and my programming skills are not great. Thanks in advance [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Michael W. Sears, Ph.D. Associate Editor, Herpetologica Assistant Professor Center for Ecology Department of Zoology Southern Illinois University Carbondale, IL 62901 phone: 618-453-4137 web:http://equinox.unr.edu/homepage/msears http://www.ecology.siu.edu Natural selection is a mechanism for generating an exceedingly high degree of improbability Sir Ronald A. Fisher (1890-1962) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Way to convert data frame to matrix
On 6/30/2006 10:11 AM, Wade Wall wrote: I have a text file that I have imported into R. It contains 3 columns and 316940 rows. The first column is vegetation plot ID, the second species names and the third is a cover value (numeric). I imported using the read.table function. My problem is this. I need to reformat the information as a matrix, with the first column becoming the row labels and the second the column labels and the cover values as the matrix cell data. However, since the read.tablefunction imported the data as an indexed data frame, I can't use the columns as vectors. Is there a way around this, to convert the data frame as 3 separate vectors? I have been looking all over for a function, and my programming skills are not great. Internally, dataframes are just lists with a class=dataframe attribute. This means you can extract the columns as if they were just lists. So if your columns are named A, B, and C, and the dataframe is dataf, you get them as vectors using dataf$A, dataf$B, and dataf$C Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Way to convert data frame to matrix
On Fri, Jun 30, 2006 at 10:11:15AM -0400, Wade Wall wrote: I have a text file that I have imported into R. It contains 3 columns and 316940 rows. The first column is vegetation plot ID, the second species names and the third is a cover value (numeric). I imported using the read.table function. My problem is this. I need to reformat the information as a matrix, with the first column becoming the row labels and the second the column labels and the cover values as the matrix cell data. I'm not 100% sure but I think you are looking for reshape(). cu Philipp -- Dr. Philipp PagelTel. +49-8161-71 2131 Dept. of Genome Oriented Bioinformatics Fax. +49-8161-71 2186 Technical University of Munich Science Center Weihenstephan 85350 Freising, Germany and Institute for Bioinformatics / MIPS Tel. +49-89-3187 3675 GSF - National Research Center Fax. +49-89-3187 3585 for Environment and Health Ingolstädter Landstrasse 1 85764 Neuherberg, Germany http://mips.gsf.de/staff/pagel __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Way to convert data frame to matrix
Hi Wade, On 6/30/06, Wade Wall [EMAIL PROTECTED] wrote: I have a text file that I have imported into R. It contains 3 columns and 316940 rows. The first column is vegetation plot ID, the second species names and the third is a cover value (numeric). I imported using the read.table function. My problem is this. I need to reformat the information as a matrix, with the first column becoming the row labels and the second the column labels and the cover values as the matrix cell data. Try crosstab(mydata$plotID, mydata$species, mydata$cover) from the ecodist package. Sarah -- Sarah Goslee __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Passing arguments to glm()
Hi there I want to pass arguments (i.e. the response variable and the subset argument) in a self-made function to glm. Here is one way I can do this: f.myglm - function(y,subfact,subval) { glm(d.mydata[,y]~d.mydata[,'x1'],family=binomial,subset=d.mydata[,subfact]==subval) } str(d.mydata) `data.frame':15806 obs. of 3 variables: $ y : Factor w/ 2 levels no,yes: 1 1 1 1 1 1 1 NA 1 1 ... $ x1: Factor w/ 2 levels no,yes: 2 2 1 2 2 2 2 2 2 2 ... $ x2: Factor w/ 2 levels no,yes: 1 1 1 1 1 2 2 1 2 2 ... f.myglm('y','x2','yes') But is there a way I can pass the arguments and use the data argument of glm()? In a naive way of thinking I'd like to something like this: f.myglm - function(y,sub) { glm(y~x1,family=binomial,data=d.mydata,subset=sub) } f.myglm(y=y,sub=x2=='yes') I know that's not possible, because the objects y and x2 are not defined in the user workspace. So, something like passing the arguments as an expression and evaluate it in the glm function should work, but I didn't manage to do it. I'd appreciate your advice. Christian R.version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] apply a function to several lists' components
Dear R-user I have 100 lists. Each list has several components. For example, data1 $a [1] 1 2 $b [1] 3 4 $c [1] 5 There are data1, data2,, data100. All lists have the same number and the same name of components. Is there any function I can use for applying to only a specific component across 100 lists? (e.g., taking mean of $c acorss 100 lists) or do I need to write my own function for that? Thank you. Taka, __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Passing arguments to glm()
f.myglm - function(y=y, subset=x2 == 'yes', data=d.d.mydata) eval(parse(text=glm(, deparse(substitute(y)), ~ x1, family=binomial, data=, deparse(substitute(data)), , subset =, subset, ))) f.myglm() --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- Christian Bieli a écrit : Hi there I want to pass arguments (i.e. the response variable and the subset argument) in a self-made function to glm. Here is one way I can do this: f.myglm - function(y,subfact,subval) { glm(d.mydata[,y]~d.mydata[,'x1'],family=binomial,subset=d.mydata[,subfact]==subval) } str(d.mydata) `data.frame':15806 obs. of 3 variables: $ y : Factor w/ 2 levels no,yes: 1 1 1 1 1 1 1 NA 1 1 ... $ x1: Factor w/ 2 levels no,yes: 2 2 1 2 2 2 2 2 2 2 ... $ x2: Factor w/ 2 levels no,yes: 1 1 1 1 1 2 2 1 2 2 ... f.myglm('y','x2','yes') But is there a way I can pass the arguments and use the data argument of glm()? In a naive way of thinking I'd like to something like this: f.myglm - function(y,sub) { glm(y~x1,family=binomial,data=d.mydata,subset=sub) } f.myglm(y=y,sub=x2=='yes') I know that's not possible, because the objects y and x2 are not defined in the user workspace. So, something like passing the arguments as an expression and evaluate it in the glm function should work, but I didn't manage to do it. I'd appreciate your advice. Christian R.version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Passing arguments to glm()
i forgot a paste(): f.myglm - function(y=y, subset=x2 == 'yes', data=d.d.mydata) eval(parse(text=paste(glm(, deparse(substitute(y)), ~ x1, family=binomial, data=, deparse(substitute(data)), , subset =, subset, ), sep=))) --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- Jacques VESLOT a écrit : f.myglm - function(y=y, subset=x2 == 'yes', data=d.d.mydata) eval(parse(text=glm(, deparse(substitute(y)), ~ x1, family=binomial, data=, deparse(substitute(data)), , subset =, subset, ))) f.myglm() --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- Christian Bieli a écrit : Hi there I want to pass arguments (i.e. the response variable and the subset argument) in a self-made function to glm. Here is one way I can do this: f.myglm - function(y,subfact,subval) { glm(d.mydata[,y]~d.mydata[,'x1'],family=binomial,subset=d.mydata[,subfact]==subval) } str(d.mydata) `data.frame':15806 obs. of 3 variables: $ y : Factor w/ 2 levels no,yes: 1 1 1 1 1 1 1 NA 1 1 ... $ x1: Factor w/ 2 levels no,yes: 2 2 1 2 2 2 2 2 2 2 ... $ x2: Factor w/ 2 levels no,yes: 1 1 1 1 1 2 2 1 2 2 ... f.myglm('y','x2','yes') But is there a way I can pass the arguments and use the data argument of glm()? In a naive way of thinking I'd like to something like this: f.myglm - function(y,sub) { glm(y~x1,family=binomial,data=d.mydata,subset=sub) } f.myglm(y=y,sub=x2=='yes') I know that's not possible, because the objects y and x2 are not defined in the user workspace. So, something like passing the arguments as an expression and evaluate it in the glm function should work, but I didn't manage to do it. I'd appreciate your advice. Christian R.version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] median of gamma distribution
Doese anyone know a R function to find the median of a gamma distribution? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] median of gamma distribution
Philip He [EMAIL PROTECTED] writes: Doese anyone know a R function to find the median of a gamma distribution? qgamma(0.5, ) -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] median of gamma distribution
On Fri, 30 Jun 2006, Philip He wrote: Doese anyone know a R function to find the median of a gamma distribution? It's not clear what you mean. If you know the parameters of a gamma distribution then qgamma() will give you any quantile. If you have data and want to estimate the median then it's hard to beat median(), but you could use mle() to estimate the parameters and then qgamma(). -thomas __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] median of gamma distribution
From: Philip He [EMAIL PROTECTED] Date: Fri Jun 30 10:30:28 CDT 2006 To: R help list r-help@stat.math.ethz.ch Subject: [R] median of gamma distribution someone might know a more elegant way but one way is to use the distribution ( i think it's dnorm for the normal but i get confused because ther is also pnorm, rnorm and one other but you can look that up to figure out which one is the one you need ) functions in R and just put in .50 as the probabiliy input and itr will kick back the result which is the median. mark Doese anyone know a R function to find the median of a gamma distribution? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] median of gamma distribution
On 30-Jun-06 Philip He wrote: Doese anyone know a R function to find the median of a gamma distribution? qgamma will do it. Test: -log(0.5) [1] 0.6931472 qgamma(0.5,1) [1] 0.6931472 Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 30-Jun-06 Time: 16:53:16 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] getting the smoother matrix from smooth.spline
smooth.matrix = function(x, df){ n = length(x); A = matrix(0, n, n); for(i in 1:n){ y = rep(0, n); y[i]=1; yi = predict(smooth.spline(x, y, df=df),x)$y; A[,i]= yi; } (A+t(A))/2; } - Simon Wood, Mathematical Sciences, University of Bath, Bath BA2 7AY - +44 (0)1225 386603 www.maths.bath.ac.uk/~sw283/ On Sat, 24 Jun 2006, Gregory Gentlemen wrote: Can anyone tell me the trick for obtaining the smoother matrix from smooth.spline when there are non-unique values for x. I have the following code but, of course, it only works when all values of x are unique. ## get the smoother matrix (x having unique values smooth.matrix = function(x, df){ n = length(x); A = matrix(0, n, n); for(i in 1:n){ y = rep(0, n); y[i]=1; yi = smooth.spline(x, y, df=df)$y; A[,i]= yi; } (A+t(A))/2; } Thanks for any assistance, Gregory - - Get a sneak peak at messages with a handy reading pane. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Random numbers from noncentral t-distribution
On Fri, 30 Jun 2006, Long Qu wrote: Hi there: I'd thought these two versions of noncentral t-distribution are essentially the same: qqplot(rt(1000,df=20,ncp=3),qt(runif(1000),df=20,ncp=3)) But, the scales of the x-axis and the y-axis are quite different according to the QQ-plot. Did I make any mistakes somewhere? No, I think we did. We have rt function (n, df, ncp = 0) { if (ncp == 0) .Internal(rt(n, df)) else rnorm(n, ncp)/(rchisq(n, df)/sqrt(df)) } and the rchisq() in the denominator should be inside the sqrt(). -thomas __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] getting the smoother matrix from smooth.spline
Perhaps this could be developed into a spline smooth method for model.matrix and included in R. On 6/30/06, Simon Wood [EMAIL PROTECTED] wrote: smooth.matrix = function(x, df){ n = length(x); A = matrix(0, n, n); for(i in 1:n){ y = rep(0, n); y[i]=1; yi = predict(smooth.spline(x, y, df=df),x)$y; A[,i]= yi; } (A+t(A))/2; } - Simon Wood, Mathematical Sciences, University of Bath, Bath BA2 7AY - +44 (0)1225 386603 www.maths.bath.ac.uk/~sw283/ On Sat, 24 Jun 2006, Gregory Gentlemen wrote: Can anyone tell me the trick for obtaining the smoother matrix from smooth.spline when there are non-unique values for x. I have the following code but, of course, it only works when all values of x are unique. ## get the smoother matrix (x having unique values smooth.matrix = function(x, df){ n = length(x); A = matrix(0, n, n); for(i in 1:n){ y = rep(0, n); y[i]=1; yi = smooth.spline(x, y, df=df)$y; A[,i]= yi; } (A+t(A))/2; } Thanks for any assistance, Gregory - - Get a sneak peak at messages with a handy reading pane. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Passing arguments to glm()
This is the same as glm except that - formula may be of class character in which case its regarded as the name of the response variable and the formula defaults to resp ~ Species for that response - the data frame defaults to iris - modify as appropriate for your case myglm - function (formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart, offset, control = glm.control(...), model = TRUE, method = glm.fit, x = FALSE, y = TRUE, contrasts = NULL, ...) { cl - match.call() cl[[1]] - as.name(glm) if (is.character(formula)) { fo - . ~ Species ### default formula fo[[2]] - as.name(formula) cl$formula - fo } if (missing(data)) cl$data - as.name(iris) # default data frame eval(cl, parent.frame()) } # test myglm(Sepal.Length, subset = Petal.Length mean(Petal.Length)) myglm(Sepal.Length ~ Petal.Length) myglm(Sepal.Length, subset = 1:100) On 6/30/06, Christian Bieli [EMAIL PROTECTED] wrote: Hi there I want to pass arguments (i.e. the response variable and the subset argument) in a self-made function to glm. Here is one way I can do this: f.myglm - function(y,subfact,subval) { glm(d.mydata[,y]~d.mydata[,'x1'],family=binomial,subset=d.mydata[,subfact]==subval) } str(d.mydata) `data.frame':15806 obs. of 3 variables: $ y : Factor w/ 2 levels no,yes: 1 1 1 1 1 1 1 NA 1 1 ... $ x1: Factor w/ 2 levels no,yes: 2 2 1 2 2 2 2 2 2 2 ... $ x2: Factor w/ 2 levels no,yes: 1 1 1 1 1 2 2 1 2 2 ... f.myglm('y','x2','yes') But is there a way I can pass the arguments and use the data argument of glm()? In a naive way of thinking I'd like to something like this: f.myglm - function(y,sub) { glm(y~x1,family=binomial,data=d.mydata,subset=sub) } f.myglm(y=y,sub=x2=='yes') I know that's not possible, because the objects y and x2 are not defined in the user workspace. So, something like passing the arguments as an expression and evaluate it in the glm function should work, but I didn't manage to do it. I'd appreciate your advice. Christian R.version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lme convergence
It looks like in the call to lme fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1)) you did not specify any random effects. Why not try: fm1 - lme(distance ~ age, random= ~1| groupID, data = Orthodont, +control=lmeControl(msMaxIter=1)) where groupID is some factor that can be used to stratify the data. Also, the Othodont data set is used in Pinheiro Bates book, and you may want to consult that book to see the models they use in connection with that data set. For the Orthodont data set the groupID would most likely be the subject ID (Subject variable). So a possible model would be: fm1 - lme(distance ~ age, random= ~1|Subject, data=Orthodont) summary(fm1) Linear mixed-effects model fit by REML Data: Orthodont AIC BIClogLik 455.0025 465.6563 -223.5013 Random effects: Formula: ~1 | Subject (Intercept) Residual StdDev:2.114724 1.431592 Fixed effects: distance ~ age Value Std.Error DF t-value p-value (Intercept) 16.76 0.8023952 80 20.5 0 age 0.660185 0.0616059 80 10.71626 0 Correlation: (Intr) age -0.845 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.66453932 -0.53507984 -0.01289591 0.48742859 3.72178465 Number of Observations: 108 Number of Groups: 27 So this runs fine. As, I said this data set and its analysis is discussed extensively in Pinheiro and Bates book Michael Jerosch-Herold Spencer Graves [EMAIL PROTECTED] 06/29/06 7:08 PM Does anyone know how to obtain the 'returnObject' from an 'lme' run that fails to converge? An argument of this name is described on the 'lmeControl' help page as, a logical value indicating whether the fitted object should be returned when the maximum number of iterations is reached without convergence of the algorithm. Default is 'FALSE'. Unfortunately, I've so far been unable to get it to work, as witnessed by the following modification of an example from the '?lme' help page: library(nlme) fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1)) Error in lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1)) : iteration limit reached without convergence (9) fm1 Error: object fm1 not found fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1, + returnObject=TRUE)) Error in lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1, : iteration limit reached without convergence (9) fm1 Error: object fm1 not found I might be able to fix the problem myself, working through the 'lme' code line by line, e.g., using 'debug'. However, I'm not ready to do that just now. Best Wishes, Spencer Graves Ravi Varadhan wrote: Use try to capture error messages without breaking the loop. ?try -- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -- -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of Pryseley Assam Sent: Wednesday, June 28, 2006 12:18 PM To: R-Users Subject: [R] lme convergence Dear R-Users, Is it possible to get the covariance matrix from an lme model that did not converge ? I am doing a simulation which entails fitting linear mixed models, using a for loop. Within each loop, i generate a new data set and analyze it using a mixed model. The loop stops When the lme function does not converge for a simulated dataset. I want to inquire if there is a method to suppress the error message from the lme function, or better still, a way of going about this issue of the loop ending once the lme function does not converge. Thanks in advance, Pryseley - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Random numbers from noncentral t-distribution
Thank you very much for your kind reply. It solved the problem of rt( ). :D But it seems that the qt( ) also have problems: I modified the rt( ) function as you suggested, rt - function (n, df, ncp = 0) { if (ncp == 0) .Internal(rt(n, df)) else rnorm(n, ncp)/sqrt(rchisq(n, df)/df) } Then I increase the number of random variables to 1, and made a QQ-plot: qqplot(rt(1,df=20,ncp=3),qt(runif(1),df=20,ncp=3)) I've got some spurious points at lower-left corner. It seems that qt( ) results were truncated. I also tried this with another df and ncp: pt(-.75,df=2,ncp=1) [1] 0.05726429 sum(qt(1:1/10001,df=2,ncp=1) -.75) [1] 0 where I'd expected the last number should be 550 or so, not 0. Thanks again, the modified rt( ) is now OK for my work. Long Thomas Lumley [EMAIL PROTECTED] wrote£º On Fri, 30 Jun 2006, Long Qu wrote: Hi there: I'd thought these two versions of noncentral t-distribution are essentially the same: qqplot(rt(1000,df=20,ncp=3),qt(runif(1000),df=20,ncp=3)) But, the scales of the x-axis and the y-axis are quite different according to the QQ-plot. Did I make any mistakes somewhere? No, I think we did. We have rt function (n, df, ncp = 0) { if (ncp == 0) .Internal(rt(n, df)) else rnorm(n, ncp)/(rchisq(n, df)/sqrt(df)) } and the rchisq() in the denominator should be inside the sqrt(). -thomas - ÇÀ×¢ÑÅ»¢Ãâ·ÑÓÊÏä-3.5GÈÝÁ¿£¬20M¸½¼þ£¡ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lme convergence
In the old version of lme, one could construct a grouped data object and this would alleviate the need to specify the random portion of the model. So, Spencer's call is equivalent to fm1 - lme(distance ~ age, random= ~age| Subject, data = Orthodont) This condition does not hold under lmer, however. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Michael Jerosch-Herold Sent: Friday, June 30, 2006 12:37 PM To: [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Subject: Re: [R] lme convergence It looks like in the call to lme fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1)) you did not specify any random effects. Why not try: fm1 - lme(distance ~ age, random= ~1| groupID, data = Orthodont, +control=lmeControl(msMaxIter=1)) where groupID is some factor that can be used to stratify the data. Also, the Othodont data set is used in Pinheiro Bates book, and you may want to consult that book to see the models they use in connection with that data set. For the Orthodont data set the groupID would most likely be the subject ID (Subject variable). So a possible model would be: fm1 - lme(distance ~ age, random= ~1|Subject, data=Orthodont) summary(fm1) Linear mixed-effects model fit by REML Data: Orthodont AIC BIClogLik 455.0025 465.6563 -223.5013 Random effects: Formula: ~1 | Subject (Intercept) Residual StdDev:2.114724 1.431592 Fixed effects: distance ~ age Value Std.Error DF t-value p-value (Intercept) 16.76 0.8023952 80 20.5 0 age 0.660185 0.0616059 80 10.71626 0 Correlation: (Intr) age -0.845 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.66453932 -0.53507984 -0.01289591 0.48742859 3.72178465 Number of Observations: 108 Number of Groups: 27 So this runs fine. As, I said this data set and its analysis is discussed extensively in Pinheiro and Bates book Michael Jerosch-Herold Spencer Graves [EMAIL PROTECTED] 06/29/06 7:08 PM Does anyone know how to obtain the 'returnObject' from an 'lme' run that fails to converge? An argument of this name is described on the 'lmeControl' help page as, a logical value indicating whether the fitted object should be returned when the maximum number of iterations is reached without convergence of the algorithm. Default is 'FALSE'. Unfortunately, I've so far been unable to get it to work, as witnessed by the following modification of an example from the '?lme' help page: library(nlme) fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1)) Error in lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1)) : iteration limit reached without convergence (9) fm1 Error: object fm1 not found fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1, + returnObject=TRUE)) Error in lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1, : iteration limit reached without convergence (9) fm1 Error: object fm1 not found I might be able to fix the problem myself, working through the 'lme' code line by line, e.g., using 'debug'. However, I'm not ready to do that just now. Best Wishes, Spencer Graves Ravi Varadhan wrote: Use try to capture error messages without breaking the loop. ?try -- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -- -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of Pryseley Assam Sent: Wednesday, June 28, 2006 12:18 PM To: R-Users Subject: [R] lme convergence Dear R-Users, Is it possible to get the covariance matrix from an lme model that did not converge ? I am doing a simulation which entails fitting linear mixed models, using a for loop. Within each loop, i generate a new data set and analyze it using a mixed model. The loop stops When the lme function does not converge for a simulated dataset. I want to inquire if there is a method to suppress the error message from the lme function, or better still, a way of going about this issue of the loop ending once the lme function does not converge. Thanks in advance, Pryseley -
Re: [R] Random numbers from noncentral t-distribution
On Sat, 1 Jul 2006, Long Qu wrote: Thank you very much for your kind reply. It solved the problem of rt( ). :D But it seems that the qt( ) also have problems: Yes, there does seem to be a problem near zero. A clearer version is curve(qt(x,df=20,ncp=3),from=0,to=0.004) curve(qt(10^x,df=20,ncp=3),from=-10,to=-2,n=1000) The fix is less obvious here. I'll file it as a bug. -thomas Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Running R
Just for the record. (One of the) problem was that configure picked up ATLAS, but had problem with it at link time for whatever reason. (This is on some version of Redhat, x86_64. I should think there are people who have similar setup and got both ATLAS and readline to work.) Andy From: Pramod Anugu IT WORKS!!! Thanks you very much for your help. So the below line fixed my installation #./configure --with-blas=no --with-readline=no __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] tkbutton command - how to know which button was clicked?
In the below code fragment, print(arg) always prints the last element of rekeningen$rekening. Is this because of lazy evaluation? I.e. arg is evaluated at the time the button is pressed? And, if so, how can I avoid this? I tried function() {force(arg); print(arg)} but that didn't work either. Thanks, Jeebee. for(rek in seq(1,nrow(rekeningen))) { arg - rekeningen$rekening[rek] tkgrid(tkbutton(frame.1, text=paste(Saldo historie, arg), command=function() print(arg)), sticky=news) } __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lme convergence
Harold is correct. The help page for 'Orthodont' includes the following example: formula(Orthodont) distance ~ age | Subject If 'random' is not specified, 'lme' sets random = formula(data). If that's NULL, the 'lme' help page says it Defaults to a formula consisting of the right hand side of 'fixed'. This will generally return an error, indicated by the following example: tstDF - data.frame(gp=rep(1:2, 2), y=1:4) lme(y~1, tstDF) Error in getGroups.data.frame(dataMix, groups) : Invalid formula for groups lme(y~1, tstDF, random=~1) Error in getGroups.data.frame(dataMix, groups) : Invalid formula for groups Hope this helps. Spencer Doran, Harold wrote: In the old version of lme, one could construct a grouped data object and this would alleviate the need to specify the random portion of the model. So, Spencer's call is equivalent to fm1 - lme(distance ~ age, random= ~age| Subject, data = Orthodont) This condition does not hold under lmer, however. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Michael Jerosch-Herold Sent: Friday, June 30, 2006 12:37 PM To: [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Subject: Re: [R] lme convergence It looks like in the call to lme fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1)) you did not specify any random effects. Why not try: fm1 - lme(distance ~ age, random= ~1| groupID, data = Orthodont, +control=lmeControl(msMaxIter=1)) where groupID is some factor that can be used to stratify the data. Also, the Othodont data set is used in Pinheiro Bates book, and you may want to consult that book to see the models they use in connection with that data set. For the Orthodont data set the groupID would most likely be the subject ID (Subject variable). So a possible model would be: fm1 - lme(distance ~ age, random= ~1|Subject, data=Orthodont) summary(fm1) Linear mixed-effects model fit by REML Data: Orthodont AIC BIClogLik 455.0025 465.6563 -223.5013 Random effects: Formula: ~1 | Subject (Intercept) Residual StdDev:2.114724 1.431592 Fixed effects: distance ~ age Value Std.Error DF t-value p-value (Intercept) 16.76 0.8023952 80 20.5 0 age 0.660185 0.0616059 80 10.71626 0 Correlation: (Intr) age -0.845 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.66453932 -0.53507984 -0.01289591 0.48742859 3.72178465 Number of Observations: 108 Number of Groups: 27 So this runs fine. As, I said this data set and its analysis is discussed extensively in Pinheiro and Bates book Michael Jerosch-Herold Spencer Graves [EMAIL PROTECTED] 06/29/06 7:08 PM Does anyone know how to obtain the 'returnObject' from an 'lme' run that fails to converge? An argument of this name is described on the 'lmeControl' help page as, a logical value indicating whether the fitted object should be returned when the maximum number of iterations is reached without convergence of the algorithm. Default is 'FALSE'. Unfortunately, I've so far been unable to get it to work, as witnessed by the following modification of an example from the '?lme' help page: library(nlme) fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1)) Error in lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1)) : iteration limit reached without convergence (9) fm1 Error: object fm1 not found fm1 - lme(distance ~ age, data = Orthodont, +control=lmeControl(msMaxIter=1, + returnObject=TRUE)) Error in lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1, : iteration limit reached without convergence (9) fm1 Error: object fm1 not found I might be able to fix the problem myself, working through the 'lme' code line by line, e.g., using 'debug'. However, I'm not ready to do that just now. Best Wishes, Spencer Graves Ravi Varadhan wrote: Use try to capture error messages without breaking the loop. ?try -- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -- -Original Message- From: [EMAIL PROTECTED] [mailto:r-help- [EMAIL PROTECTED] On Behalf Of Pryseley Assam Sent: Wednesday, June 28, 2006 12:18 PM To: R-Users Subject: [R] lme
Re: [R] apply a function to several lists' components
Maybe this helps ( data1 = list(a=c(1,2), b=c(3,4), c=c(5,6,7)) ) ( data2 = list(a=c(10,11), b=c(30,40), c=c(70,80)) ) cc - NULL for(data in ls(pattern=^data[0-9]+$)) { cc - c(cc, with(get(data), c)) } mean(cc) JeeBee. On Fri, 30 Jun 2006 09:50:51 -0500, Taka Matzmoto wrote: Dear R-user I have 100 lists. Each list has several components. For example, data1 $a [1] 1 2 $b [1] 3 4 $c [1] 5 There are data1, data2,, data100. All lists have the same number and the same name of components. Is there any function I can use for applying to only a specific component across 100 lists? (e.g., taking mean of $c acorss 100 lists) or do I need to write my own function for that? Thank you. Taka, __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] send output to printer
It has to be a simple thing, but I could not figure it out: How do I send the text output from object x to the printer? As a shell user I would expect a pipe to the printer... |kprinter or |lpr -Pmyprinter somehow. And yes, I'm on Linux. Thanks! __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] SAS Proc Mixed and lme
I am trying to use lme to fit a mixed effects model to get the same results as when using the following SAS code: proc mixed; class refseqid probeid probeno end; model expression=end logpgc / ddfm=satterth; random probeno probeid / subject=refseqid type=cs; lsmeans end / diff cl; run; There are 3 genes (refseqid) which is the large grouping factor, with 2 probeids nested within each refseqid, and 16 probenos nested within each of the probeids. I have specified in the SAS Proc Mixed procedure that the variance-covariance structure is to be compound symmetric. Therefore, the variance-covariance matrix is a block diagonal matrix of the form, V_1 0 0 0 V_2 0 00 V3 where each V_i represents a RefSeqID. Moreover, for each V_i the structure within the block is v_{11} v{12} v_{21} v{22} where v_{11} and v_{22} are different probeids nested within the refseqid, and so are correlated. The structure of both v_{11} and v_{22} are compound symmetric, and v_{12} and v{21} contain a constant for all elements of the matrix. I have tried to reproduce this using lme, but it is unclear from the documentation (and Pinheiro Bates text) how the pdBlocked and compound symmetric structure can be combined. fit.lme-lme(expression~End+logpgc,random=list(RefSeqID=pdBlocked(list (~1,~ProbeID-1),pdClass=pdSymm)),data=dataset,correlation=corCompSym m(form=~1|RefSeqID/ProbeID/ProbeNo)) The point estimates are essentially the same comparing R and SAS for the fixed effects, but the 95% confidence intervals are much shorter using lme(). In order to find the difference in the algorithms used by SAS and R I tried to extract the variance-covariance matrix to look at its structure. I used the getVarCov() command, but it tells me that this function is not available for nested structures. Is there another way to extract the variance-covariance structure for nested models? Does anyone know how I could get the var-cov structure above using lme? Kellie J. Archer, Ph.D. Assistant Professor, Department of Biostatistics Fellow, Center for the Study of Biological Complexity Virginia Commonwealth University 1101 East Marshall St., B1-066 Richmond, VA 23298-0032 phone: (804) 827-2039 fax: (804) 828-8900 e-mail: [EMAIL PROTECTED] website: www.people.vcu.edu/~kjarcher __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] lme and SAS Proc mixed
I am trying to use lme to fit a mixed effects model to get the same results as when using the following SAS code: proc mixed; class refseqid probeid probeno end; model expression=end logpgc / ddfm=satterth; random probeno probeid / subject=refseqid type=cs; lsmeans end / diff cl; run; There are 3 genes (refseqid) which is the large grouping factor, with 2 probeids nested within each refseqid, and 16 probenos nested within each of the probeids. I have specified in the SAS Proc Mixed procedure that the variance-covariance structure is to be compound symmetric. Therefore, the variance-covariance matrix is a block diagonal matrix of the form, V_1 0 0 0 V_2 0 00 V3 where each V_i represents a RefSeqID. Moreover, for each V_i the structure within the block is v_{11} v{12} v_{21} v{22} where v_{11} and v_{22} are different probeids nested within the refseqid, and so are correlated. The structure of both v_{11} and v_{22} are compound symmetric, and v_{12} and v{21} contain a constant for all elements of the matrix. I have tried to reproduce this using lme, but it is unclear from the documentation (and Pinheiro Bates text) how the pdBlocked and compound symmetric structure can be combined. fit.lme-lme(expression~End+logpgc,random=list(RefSeqID=pdBlocked(list (~1,~ProbeID-1),pdClass=pdSymm)),data=dataset,correlation=corCompSym m(form=~1|RefSeqID/ProbeID/ProbeNo)) The point estimates are essentially the same comparing R and SAS for the fixed effects, but the 95% confidence intervals are much shorter using lme(). In order to find the difference in the algorithms used by SAS and R I tried to extract the variance-covariance matrix to look at its structure. I used the getVarCov() command, but it tells me that this function is not available for nested structures. Is there another way to extract the variance-covariance structure for nested models? Does anyone know how I could get the var-cov structure above using lme? Kellie J. Archer, Ph.D. Assistant Professor, Department of Biostatistics Fellow, Center for the Study of Biological Complexity Virginia Commonwealth University 1101 East Marshall St., B1-066 Richmond, VA 23298-0032 phone: (804) 827-2039 fax: (804) 828-8900 e-mail: [EMAIL PROTECTED] website: www.people.vcu.edu/~kjarcher __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Creating Vectors
type count 0 20 1 15 0 10 1 35 0 28 I would like to create two vectors from the data above. For example, type1=c(15, 35) and type0 = c(20, 10, 28). Can any one help Raphael __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] tkbutton command - how to know which button was clicked?
JeeBee [EMAIL PROTECTED] writes: In the below code fragment, print(arg) always prints the last element of rekeningen$rekening. Is this because of lazy evaluation? I.e. arg is evaluated at the time the button is pressed? No and yes. Lazy evaluation has nothing to do with it, but the function passed to command= is not evaluated until the button is pressed. If you want a different arg variable for each button-function, you have to make sure that they have different environments or have the arg actually part of the function. There are various possible variations. One is mk.button - function(arg) tkgrid(tkbutton(frame.1, text=paste(Saldo historie, arg), command=function() print(arg)), sticky=news) for { mk.button(arg) } another is to wrap a local({arg - arg;.}) around the tkgrid call in your code. A 3rd one is .command=eval(substitute(function()print(arg)))... And, if so, how can I avoid this? I tried function() {force(arg); print(arg)} but that didn't work either. Thanks, Jeebee. for(rek in seq(1,nrow(rekeningen))) { arg - rekeningen$rekening[rek] tkgrid(tkbutton(frame.1, text=paste(Saldo historie, arg), command=function() print(arg)), sticky=news) } -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] useR! 2006: presentation slides
Dear useRs, the useR! 2006 conference took place in Vienna two weeks ago: it was an exciting and interesting meeting with about 400 useRs from all over the world and more than 150 presentations. Especially for those of you who could not make it to the conference, we have made 4up PDF versions of the presentations slides available. The slides for the keynote lectures are available at http://www.R-project.org/useR-2006/Keynotes/ and the user-contributed presentations at http://www.R-project.org/useR-2006/Presentations/ Finally, some materials for the panel discussion on Getting recognition for excellence in computational statistics are provided at http://www.R-project.org/useR-2006/PanelDisc/ including a summary, the short presentation slides of the panelists and the full discussion as a video. A big thank you to everyone who contributed to the conference and... best wishes from Vienna! The organizing team Achim, Torsten, David, Bettina, Fritz and Kurt. ___ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-announce __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Creating Vectors
Hi all. I have a factor variable distributed over time. I am looking for an elegant way to code duration of a state. Suppose, rainfall.shocks - factor(sample(c(1,2,3), size = 15, replace = TRUE, prob = unit.p), + label = c(Drought, Normal, High)) rainfall.shocks [1] Normal HighHighDrought Normal Normal HighNormal Drought [10] Normal Drought Normal Normal Normal Normal So capture the duration of say drought, I'd need a variable that is able to keep track of rainfall.shocks as well as its past values. I was wondering if there is any obvious way to do this. the Drought variable in this case would have values 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 many thanks for the suggestions you are likely to make. Alexander Nervedi __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Fwd: time series patterns
-- Forwarded message -- From: Alexander Nervedi [EMAIL PROTECTED] Date: Jun 30, 2006 4:56 PM Subject: time series patterns To: [EMAIL PROTECTED] Hi all. I have a factor variable distributed over time. I am looking for an elegant way to code duration of a state. Suppose, rainfall.shocks - factor(sample(c(1,2,3), size = 15, replace = TRUE, prob = unit.p), + label = c(Drought, Normal, High)) rainfall.shocks [1] Normal HighHighDrought Normal Normal HighNormal Drought [10] Normal Drought Normal Normal Normal Normal So capture the duration of say drought, I'd need a variable that is able to keep track of rainfall.shocks as well as its past values. I was wondering if there is any obvious way to do this. the Drought variable in this case would have values 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 many thanks for the suggestions you are likely to make. Alexander Nervedi _ Express yourself instantly with MSN Messenger! Download today - it's FREE! http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] {Spam?} {Virus?} Message could not be delivered
I am out of the office and will be back on Monday 31 July. Your mail will not be forwarded automatically. Your mail regarding {Spam?} {Virus?} Message could not be delivered will be read when I return In urgent cases please contact Mr Kholdoun Torki mailto:[EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] weird error message
Hi! In the example below why is x[[1]] == 0.2237724 false? Alexander Nervedi x - runif(10) x[[1]] [1] 0.2237724 x [1] 0.2237724 0.2678944 0.9375811 0.5963889 0.6180519 0.6449580 0.7308510 [8] 0.7347386 0.4837286 0.1416100 x[[1]] == 0.2237724 FALSE From: Alexander Nervedi [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Subject: Re: [R] Creating Vectors Date: Fri, 30 Jun 2006 21:57:35 + Hi all. I have a factor variable distributed over time. I am looking for an elegant way to code duration of a state. Suppose, rainfall.shocks - factor(sample(c(1,2,3), size = 15, replace = TRUE, prob = unit.p), + label = c(Drought, Normal, High)) rainfall.shocks [1] Normal HighHighDrought Normal Normal HighNormal Drought [10] Normal Drought Normal Normal Normal Normal So capture the duration of say drought, I'd need a variable that is able to keep track of rainfall.shocks as well as its past values. I was wondering if there is any obvious way to do this. the Drought variable in this case would have values 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 many thanks for the suggestions you are likely to make. Alexander Nervedi __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Creating Vectors
Does this do what you want? It creates a 'list' with the vectors: x - 'type count + 0 20 + 1 15 +0 10 + 1 35 + 0 28 + ' x - read.table(textConnection(x), header=TRUE) x type count 1020 2115 3010 4135 5028 type - split(x$count, x$type) type $`0` [1] 20 10 28 $`1` [1] 15 35 type[['0']] [1] 20 10 28 type[['1']] [1] 15 35 On 6/30/06, Raphael Fraser [EMAIL PROTECTED] wrote: type count 0 20 1 15 0 10 1 35 0 28 I would like to create two vectors from the data above. For example, type1=c(15, 35) and type0 = c(20, 10, 28). Can any one help Raphael __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Jim Holtman Cincinnati, OH +1 513 646 9390 (Cell) +1 513 247 0281 (Home) What is the problem you are trying to solve? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] weird error message
This is FAQ 7.31: http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f Also please do not piggy back on other threads since it makes the archives less useful. On 6/30/06, Alexander Nervedi [EMAIL PROTECTED] wrote: Hi! In the example below why is x[[1]] == 0.2237724 false? Alexander Nervedi x - runif(10) x[[1]] [1] 0.2237724 x [1] 0.2237724 0.2678944 0.9375811 0.5963889 0.6180519 0.6449580 0.7308510 [8] 0.7347386 0.4837286 0.1416100 x[[1]] == 0.2237724 FALSE From: Alexander Nervedi [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Subject: Re: [R] Creating Vectors Date: Fri, 30 Jun 2006 21:57:35 + Hi all. I have a factor variable distributed over time. I am looking for an elegant way to code duration of a state. Suppose, rainfall.shocks - factor(sample(c(1,2,3), size = 15, replace = TRUE, prob = unit.p), + label = c(Drought, Normal, High)) rainfall.shocks [1] Normal HighHighDrought Normal Normal HighNormal Drought [10] Normal Drought Normal Normal Normal Normal So capture the duration of say drought, I'd need a variable that is able to keep track of rainfall.shocks as well as its past values. I was wondering if there is any obvious way to do this. the Drought variable in this case would have values 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 many thanks for the suggestions you are likely to make. Alexander Nervedi __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] postscript file too large : maybe an R question
i created a postscipt file in R and then i downloaded a free version of ghostview to view it. unfortunately, i get the message fata error : dynamic memory exhausted when i try to view it. when i do a dir on windows xp, the file size is 149,034,475 and i know there about 17,000 graphs. is there a way of possibly viewing this size postscript file in R itself ? Thanks __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] polynomial expansion in R
Hi: I have two vectors of data, x and y and I want to get the polynomial expansion of (x+y)^p with any integer power p in R. Suppose p=2, then I want a matrix of five vectors, namely, x y x^2 y^2 x*y. The coefficient of the polynomial is not needed. I can write it manully if p is small. But I want it in the case of p=10 or even bigger, is there any function in R can do that automatically for me with certain choice of p? Thx a lot! liu - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] generate bi-variate normal data
Dear all, I would like to generate bi-variate normal data given that the first column of the data is known. for example: I first generate a set of data using the command, x - rmvnorm(10, c(0, 0), matrix(c(1, 0, 0, 1), 2)) then I would like to sum up the two columns of x: x.sum - apply(x, 1, sum) now with x.sum I would like to generate another column of data, say y, that makes cbind(x.sum, y) follow a bi-variate normal distribution with mean = c(0, 0) and sigma = matrix(c(1, 0, 0, 1),2) I will appreciate for all insights. David s. This email may contain confidential material.\ If you were n...{{dropped}} __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] generate bi-variate normal data
From: Shin, David [EMAIL PROTECTED] Date: Sat Jul 01 00:15:21 CDT 2006 To: 'r-help@stat.math.ethz.ch' r-help@stat.math.ethz.ch Subject: [R] generate bi-variate normal data it's an interesting question. someone else on this list can answer more explicitly but i think you have to use the result for the multivariate normal distribution ( bivariate case ) where , if the joint is normal , then the conditional is normal also with parameters a function of the 2 means and the elements of the covariance matrix. the result in any decent mathematical statistics such as casella berger. so, given the one column, generate the other column conditionally using the formula and then the joint dist will be bivariate normal. Dear all, I would like to generate bi-variate normal data given that the first column of the data is known. for example: I first generate a set of data using the command, x - rmvnorm(10, c(0, 0), matrix(c(1, 0, 0, 1), 2)) then I would like to sum up the two columns of x: x.sum - apply(x, 1, sum) now with x.sum I would like to generate another column of data, say y, that makes cbind(x.sum, y) follow a bi-variate normal distribution with mean = c(0, 0) and sigma = matrix(c(1, 0, 0, 1),2) I will appreciate for all insights. David s. This email may contain confidential material.\ If you were n...{{dropped}} __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html