Re: [R] import csv file problem
Erik Iverson-3 wrote: On 09/22/2010 07:24 PM, sisxy wrote: R will search in its working directory for Q.csv. What is the working directory, use getwd() to find out. Thanks , by using the getwd() , i found the place that i should save in my laptop... then after i saved in the location that show by the getwd() and i type the code read.csv(file=xx.csv,header=TRUE) My db came out It's useful. thanks ^^ -- View this message in context: http://r.789695.n4.nabble.com/import-csv-file-problem-tp2551256p2551455.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] ergm
Dear colleagues, I have another question, which, I think cannot be answered easily by the manual. What is the effect of including both nodefactor(Gender) and nodematch(Gender,diff=TRUE)) for the same variable in the model? Judging from the output (please see below), you cant have estimates for both for boys and girls forthe nodematch command, but I thought that the nodematch measures for boys and girls were not related. Does the nodematch command show the relative frequency of matching boy-boy against girl-girl, so one of the two categories has to be a reference category? Thank you Formula: ng ~ edges + nodefactor(Gender) + nodematch(LEA, diff = TRUE) + nodematch(Gender, diff = TRUE) Newton-Raphson iterations: 5 Maximum Likelihood Results: Estimate Std. Error MCMC s.e. p-value edges -2.098842 0.033393NA 1e-04 *** nodefactor.Gender.Girl 0.612314 0.021956NA 1e-04 *** nodematch.LEA.1 0.018116 0.032228NA 0.5740 nodematch.LEA.3-0.002688 0.154005NA 0.9861 nodematch.LEA.4 0.103054 0.043836NA 0.0187 * nodematch.LEA.5 0.016844 0.033565NA 0.6158 nodematch.LEA.6-0.086837 0.100466NA 0.3874 nodematch.Gender.Boy0.593721 0.038632NA 1e-04 *** nodematch.Gender.Girl NA 0.00NA NA --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 For this model, the pseudolikelihood is the same as the likelihood. Null Deviance: 92092 on 66430 degrees of freedom Residual Deviance: 67948 on 66421 degrees of freedom Deviance: 24144 on 9 degrees of freedom AIC: 67966BIC: 68048 Dr. Iasonas Lamprianou Assistant Professor (Educational Research and Evaluation) Department of Education Sciences European University-Cyprus P.O. Box 22006 1516 Nicosia Cyprus Tel.: +357-22-713178 Fax: +357-22-590539 Honorary Research Fellow Department of Education The University of Manchester Oxford Road, Manchester M13 9PL, UK Tel. 0044 161 275 3485 iasonas.lampria...@manchester.ac.uk __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] accumulation curves
Hello Kyran, Some more details of your data would be helpful. For example, is it cumulative species count over time ? Michael On 23 September 2010 15:05, Kyran Staunton staunto...@gmail.com wrote: Hi, I am trying to fit a logarithmic trendline to a scatterplot of a species accumulation curve. I've tried abline, lines, curve and scatter.smooth but none of these work. Can anyone help please, Kyran __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] lattice centre a diverging colour scale
Dear list, I'm using lattice::levelplot to plot a coloured image of 3D data. The range of the z values goes from negative to positive, but is not exactly centred around 0. I would however like to map a diverging colour scale with white falling exactly at 0, and both extremes being symmetrical in the legend to better contrast the opposite change in colour saturation. The following dummy example illustrates my problem, library(lattice) d - transform(expand.grid(x=seq(0, 10, length=100), y=seq(0, 10, length=100)), z = sin(x/pi)*cos(0.5*y/pi) - 0.2) levelplot(z~x*y, data=d, panel=panel.levelplot.raster, cuts = 100, interpolate = TRUE) The colour scale goes from -0.3 to 0.9 with a middle (white) value of 0.3 approximately. I'd like it to be from -1 (most saturated blue) to 1 (most saturated pink), say, with 0 being white. I read the entry in ?levelplot and in ?level.colors but could not find a solution to this particular case. Sincerely, baptiste __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Help me pls
Dear All, I need to create eps file which is the required figure format of the journal that I want to submit a paper. I am able to create files in pdf or wmf format but not in eps format. Is there a way to convert pdf or wmf to eps? or alternatively, how can I create an eps file in R? Any help is deeply appreciated. Kind Regards Seyit Ali __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] accumulation curves
OK, assuming a trend in estimated spp richness vs log(effort) you could do this... plot(effort, richness) rich.lm - lm( richness ~ log(effort) ) smooth.effort - seq(1, 10, 0.1) # or whatever is appropriate lines( smooth.effort, predict(rich.lm, newdata=list(effort=smooth.effort)), col=red ) Is that the sort of thing that you're after ? Michael On 23 September 2010 17:07, Kyran Staunton staunto...@gmail.com wrote: Sorry Michael, Yes it surely is, Chao1 species richness increase per sampling effort run through fossil. cheers, Kyran On 23 September 2010 16:58, Michael Bedward michael.bedw...@gmail.com wrote: Hello Kyran, Some more details of your data would be helpful. For example, is it cumulative species count over time ? Michael On 23 September 2010 15:05, Kyran Staunton staunto...@gmail.com wrote: Hi, I am trying to fit a logarithmic trendline to a scatterplot of a species accumulation curve. I've tried abline, lines, curve and scatter.smooth but none of these work. Can anyone help please, Kyran __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Help me pls
The postscript function ? On 23 September 2010 17:12, Seyit Ali KAYIS ska...@selcuk.edu.tr wrote: Dear All, I need to create eps file which is the required figure format of the journal that I want to submit a paper. I am able to create files in pdf or wmf format but not in eps format. Is there a way to convert pdf or wmf to eps? or alternatively, how can I create an eps file in R? Any help is deeply appreciated. Kind Regards Seyit Ali __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] dnorm
Dear R-users Idea: Plot a dnorm line using specific mean/sd to complete a histogram (skewed). xs:range of y-values, ys: dnorm function Problem: I expected to multiply the ys function with the sample size (n=250-300). I was wondering about a factor between 12'000 and 30'000 to match the size of the dnorm line with the specific histogram. Thanks Sibylle hist(Biotree[Ld,]$Height2008, main=Larix decidua, ylim=c(0,50), xlab=Tree Height 2008 (cm),col=aquamarine, font.main=3, cex.axis=0.8) xs-0:650 ys-dnorm(xs, mean=397.8, sd=97.6) lines(xs,ys*12000) -- GMX DSL SOMMER-SPECIAL: Surf Phone Flat 16.000 für nur 19,99 Euro/mtl.!* __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Help me pls
From: Seyit Ali KAYIS ska...@selcuk.edu.tr
To: r-help@r-project.org
Sent: Thu, September 23, 2010 12:12:49 AM
Subject: [R] Help me pls
Dear All,
I need to create eps file which is the required figure format of the
journal that I want to submit a paper. I am able to create files in pdf or
wmf format but not in eps format. Is there a way to convert pdf or wmf to
eps? or alternatively, how can I create an eps file in R?
?postscript
In particular:
The postscript produced for a single R plot is EPS (_Encapsulated
PostScript_) compatible, and can be included into other documents,
e.g., into LaTeX, using ‘\includegraphics{filename}’. For use
in this way you will probably want to use ‘setEPS()’ to set the
defaults as ‘horizontal = FALSE, onefile = FALSE, paper =
special’. Note that the bounding box is for the _device_
region: if you find the white space around the plot region
excessive, reduce the margins of the figure region via
‘par(mar=)’.
-- Mike
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] dnorm
On 09/23/2010 05:42 PM, Sibylle Stöckli wrote: Dear R-users Idea: Plot a dnorm line using specific mean/sd to complete a histogram (skewed). xs:range of y-values, ys: dnorm function Problem: I expected to multiply the ys function with the sample size (n=250-300). I was wondering about a factor between 12'000 and 30'000 to match the size of the dnorm line with the specific histogram. Thanks Sibylle hist(Biotree[Ld,]$Height2008, main=Larix decidua, ylim=c(0,50), xlab=Tree Height 2008 (cm),col=aquamarine, font.main=3, cex.axis=0.8) xs-0:650 ys-dnorm(xs, mean=397.8, sd=97.6) lines(xs,ys*12000) Hi Sibylle, You can use one of the rescale functions (one is in the plotrix package) like this: Ld.hist-hist(Biotree[Ld,]$Height2008, main=Larix decidua, ylim=c(0,50), xlab=Tree Height 2008 (cm), col=aquamarine, font.main=3,cex.axis=0.8) xs-0:650 ys-dnorm(xs, mean=397.8, sd=97.6) lines(xs,rescale(ys,0:max(Ld.hist))) Jim __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] lattice centre a diverging colour scale
On Thu, Sep 23, 2010 at 12:33 PM, baptiste auguie baptiste.aug...@googlemail.com wrote: Dear list, I'm using lattice::levelplot to plot a coloured image of 3D data. The range of the z values goes from negative to positive, but is not exactly centred around 0. I would however like to map a diverging colour scale with white falling exactly at 0, and both extremes being symmetrical in the legend to better contrast the opposite change in colour saturation. The following dummy example illustrates my problem, library(lattice) d - transform(expand.grid(x=seq(0, 10, length=100), y=seq(0, 10, length=100)), z = sin(x/pi)*cos(0.5*y/pi) - 0.2) levelplot(z~x*y, data=d, panel=panel.levelplot.raster, cuts = 100, interpolate = TRUE) The colour scale goes from -0.3 to 0.9 with a middle (white) value of 0.3 approximately. I'd like it to be from -1 (most saturated blue) to 1 (most saturated pink), say, with 0 being white. I read the entry in ?levelplot and in ?level.colors but could not find a solution to this particular case. Specify the cut-points explicitly: levelplot(z~x*y, data=d, panel = panel.levelplot.raster, at = do.breaks(c(-1, 1), 100)) -Deepayan __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] lattice centre a diverging colour scale
Great, I don't know how I missed that, thanks! baptiste On Sep 23, 2010, at 10:39 AM, Deepayan Sarkar wrote: On Thu, Sep 23, 2010 at 12:33 PM, baptiste auguie baptiste.aug...@googlemail.com wrote: Dear list, I'm using lattice::levelplot to plot a coloured image of 3D data. The range of the z values goes from negative to positive, but is not exactly centred around 0. I would however like to map a diverging colour scale with white falling exactly at 0, and both extremes being symmetrical in the legend to better contrast the opposite change in colour saturation. The following dummy example illustrates my problem, library(lattice) d - transform(expand.grid(x=seq(0, 10, length=100), y=seq(0, 10, length=100)), z = sin(x/pi)*cos(0.5*y/pi) - 0.2) levelplot(z~x*y, data=d, panel=panel.levelplot.raster, cuts = 100, interpolate = TRUE) The colour scale goes from -0.3 to 0.9 with a middle (white) value of 0.3 approximately. I'd like it to be from -1 (most saturated blue) to 1 (most saturated pink), say, with 0 being white. I read the entry in ?levelplot and in ?level.colors but could not find a solution to this particular case. Specify the cut-points explicitly: levelplot(z~x*y, data=d, panel = panel.levelplot.raster, at = do.breaks(c(-1, 1), 100)) -Deepayan __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] accumulation curves
Hi Kyran, Please reply via the list - you'll get more answers that way :) My example had the data in vectors, you are using a data.frame. Try replacing the relevant bits of your code (below) with these lines and see if that helps... rich.lm ~ lm(Chao1 ~ log(N.obs), data=table1) smooth.effort - seq(0,54,0.1) lines(smooth.effort, predict(rich.lm, newdata=list(N.obs=smooth.effort)), col=red) Michael On 23 September 2010 18:06, Kyran Staunton staunto...@gmail.com wrote: Ok so using your additions my script is.. library(adehabitat) in.dir = in.dir=D:/R output/Richness/BK12SR.csv out.dir = D:/R output/Richness/ table1= read.csv(file =in.dir) y_range - range(1:15) x_range = range(0,60) png(filename=D:/R output/Richness/BK12SR.png, height=800, width=1500, bg=white) par(mfrow = c(1, 1), pty = m, ps=15,cex.main=0.9,mar=c(5,7,4,7)) plot(table1$Chao1~table1$N.obs, ann=FALSE, axes=TRUE, ylim=y_range, xlim=x_range, cex=1) rich.lm - lm(table1$Chao1~ log (table1$N.obs)) smooth.effort - seq(0,54,0.1) lines(smooth.effort, predict(rich.lm, newdata=list(table1$N.obs=smooth.effort)), col=red) title(xlab= Samples) title(ylab= Chao1 (mean)) box() dev.off() #I have 54 samples (N.obs) and about 14 species (Chao1) R won't run the line as it states.. Error: unexpected '=' in lines(smooth.effort, predict(rich.lm, newdata=list(table1$N.obs= Can you please decipher this error, Thank you very much for your time Kyran. On 23 September 2010 17:36, Michael Bedward michael.bedw...@gmail.com wrote: OK, assuming a trend in estimated spp richness vs log(effort) you could do this... plot(effort, richness) rich.lm - lm( richness ~ log(effort) ) smooth.effort - seq(1, 10, 0.1) # or whatever is appropriate lines( smooth.effort, predict(rich.lm, newdata=list(effort=smooth.effort)), col=red ) Is that the sort of thing that you're after ? Michael On 23 September 2010 17:07, Kyran Staunton staunto...@gmail.com wrote: Sorry Michael, Yes it surely is, Chao1 species richness increase per sampling effort run through fossil. cheers, Kyran On 23 September 2010 16:58, Michael Bedward michael.bedw...@gmail.com wrote: Hello Kyran, Some more details of your data would be helpful. For example, is it cumulative species count over time ? Michael On 23 September 2010 15:05, Kyran Staunton staunto...@gmail.com wrote: Hi, I am trying to fit a logarithmic trendline to a scatterplot of a species accumulation curve. I've tried abline, lines, curve and scatter.smooth but none of these work. Can anyone help please, Kyran __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] eps file
Dear All, I need to create eps file which is the required figure format of the journal that I want to submit a paper. I am able to create files in pdf or wmf format but not in eps format. Is there a way to convert pdf or wmf to eps? or alternatively, how can I create an eps file in R? Any help is deeply appreciated. Kind Regards Seyit Ali Dr. Seyit Ali KAYIS Selcuk University, Faculty of Agriculture Kampus/Konya, Turkey Tel: +90 332 223 2830 Mobile: +90 535 587 1139 Greetings from Konya, Turkey http://www.ziraat.selcuk.edu.tr/skayis/ [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] apply union function vectorially
Thank you very much, this solve the problem, but more generally is there a function that allow to apply a function to a list of object, applying recursively the function to each answer... mathematically for a 2 argument function f(u,v) you would like a function g doing g(u1,u2,u3,u4,u5) =f(f(f(f(u1,u2),u3),u4),u5) Thanks if you can help -- View this message in context: http://r.789695.n4.nabble.com/apply-union-function-vectorially-tp2550162p2551536.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] referencing last row in a column
I am trying to set u a limit for my plot window according to the range of a year column I have. The first year in the range is just simply the first row, referred to as data$year[1]. How can I find the last row in a similar way to set the last year in the range? What I want to accomplish is the following: plot(,xlim = c(data$year[1], data$data[lastrow]),...), is this possible? -- View this message in context: http://r.789695.n4.nabble.com/referencing-last-row-in-a-column-tp2551645p2551645.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] referencing last row in a column
You could do data$year[ nrow(data) ] On 23 September 2010 18:56, fugelpitch jo...@runtimerecords.net wrote: I am trying to set u a limit for my plot window according to the range of a year column I have. The first year in the range is just simply the first row, referred to as data$year[1]. How can I find the last row in a similar way to set the last year in the range? What I want to accomplish is the following: plot(,xlim = c(data$year[1], data$data[lastrow]),...), is this possible? -- View this message in context: http://r.789695.n4.nabble.com/referencing-last-row-in-a-column-tp2551645p2551645.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] referencing last row in a column
Thanks, but I have multiple entries for each year so the column reads: year 1877 1877 1877 1877 1878 1878 1878 1878 1879 1879 1879 1879 So I need to pick out the year from the last row rather than just counting rows. -- View this message in context: http://r.789695.n4.nabble.com/referencing-last-row-in-a-column-tp2551645p2551675.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] quantile curves
Dear List, I fitted a bivariate extreme value distribution to my matrix (data). Now I would like to add ‘quantile curves’ to the plot, similar to those calculated with qcbvnonpar() But I need the probability of the upper right corner instead of the lower left corner. (So that the curves in the plot below would look like they were turned by 180 degrees). Is there a package that could help me? Or is there an option within the ‘evd’ package that I didn’t see? My data: a-c(6.8044, 6.4422, 6.0900, 6.1778, 6.2778, 6.3978, 6.2156, 6.1512, 6.2712, 6.3112, 6.1590, 6.2368, 6.1968, 6.1746, 6.2202, 6.5836, 6.1970, 6.2870, 6.1126, 6.3460, 6.0616, 6.2016, 6.3294, 6.1494, 6.0350, 6.2006, 6.1106, 6.3506, 6.3706, 6.0484, 6.3140, 6.0618, 6.0418, 6.4618, 6.1996, 6.3196, 6.0996, 6.2174, 6.1574, 6.0774, 6.1874, 6.0230, 6.0608, 6.3386, 6.0686, 6.3042, 6.2042, 6.6242, 6.0720, 6.1276, 6.1410, 6.1010, 6.2688, 6.5788, 6.3566, 6.1766, 6.1444, 6.6200, 6.0178, 6.1856, 6.1634) b-c(20, 16, 12, 16, 10, 26, 11, 5, 12, 16, 6, 12, 9, 15, 13, 22, 10, 15, 10, 9, 5, 8, 30, 8, 5, 8, 7, 23, 18, 5, 8, 6, 5, 28, 6, 9, 7, 8, 6, 5, 6, 6, 6, 15, 6, 8, 10, 26, 7, 5, 8, 8, 6, 11, 8, 8, 11, 16, 6, 12, 12) data-cbind(a,b) bifit-fbvevd(data,model=log) plot(bifit,which=5) Thanks a lot, Aishe __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] eps file
On Thu, Sep 23, 2010 at 8:05 AM, KAYIS Seyit Ali s_a_ka...@yahoo.com wrote:
I need to create eps file which is the required figure format
of the journal that I want to submit a paper. I am able to
create files in pdf or wmf format but not in eps format. Is
there a way to convert pdf or wmf to eps? or alternatively, how
can I create an eps file in R?
see ?postscript
The postscript produced by R is EPS (_Encapsulated PostScript_)
compatible, and can be included into other documents, e.g., into
LaTeX, using '\includegraphics{filename}'. For use in this way you
will probably want to use setEPS() to set the defaults as
horizontal = FALSE, onefile = FALSE, paper = special.
dev.copy2eps: for copying from screen to EPS.
Xpdf (http://www.foolabs.com/xpdf/) includes a pdftops utility for
converting PDF files to (Encapsulated) PostScript.
HTH,
Rainer
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Doing operations by grouping variable
Another option for doing opertions by grouping variables is to the the
plyr package.
d - data.frame(x=1:10,
g1=LETTERS[rep(11:12,each=5)],
g2=letters[rep(21:23,c(3,3,4))]
)
library(plyr)
ddply(d, c(g1, g2), function(z){
z$x - z$x / max(z$x)
z
})
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium
Research Institute for Nature and Forest
team Biometrics Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
thierry.onkel...@inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey
-Oorspronkelijk bericht-
Van: r-help-boun...@r-project.org
[mailto:r-help-boun...@r-project.org] Namens William Dunlap
Verzonden: woensdag 22 september 2010 19:51
Aan: Seth W Bigelow; bill.venab...@csiro.au
CC: r-help@r-project.org
Onderwerp: Re: [R] Doing operations by grouping variable
-Original Message-
From: r-help-boun...@r-project.org
[mailto:r-help-boun...@r-project.org] On Behalf Of Seth W Bigelow
Sent: Tuesday, September 21, 2010 4:22 PM
To: bill.venab...@csiro.au
Cc: r-help@r-project.org
Subject: Re: [R] Doing operations by grouping variable
Aah, that is the sort of truly elegant solution I have been
seeking.
And it's wrapped up in a nice programming shortcut to boot
(i.e., the
within statement). I retract anything I may have said about tapply
being clunky.
Many thanks
--Seth
Dr. Seth W. Bigelow
Biologist, USDA-FS Pacific Southwest Research Station
1731 Research Park Drive, Davis California
bill.venab...@csiro.au
09/21/2010 03:15 PM
To sbige...@fs.fed.us
You left out the subscript. Why not just do
d - within(data.frame(group = rep(1:5, each = 5),
variable = rnorm(25)),
scaled - variable/tapply(variable, group, max)[group])
This approach can be tricky when there is more than one
grouping variable. E.g., suppose we have grouping variables
g1 and g2:
d - data.frame(x=1:10,
g1=LETTERS[rep(11:12,each=5)],
g2=letters[rep(21:23,c(3,3,4))])
d
x g1 g2
1 1 K u
2 2 K u
3 3 K u
4 4 K v
5 5 K v
6 6 L v
7 7 L w
8 8 L w
9 9 L w
10 10 L w
and we want to divide each x value by it max for each
g1*g2 group (6 possible groups, of which 4 are in the data).
You can extend Bill V.'s approach with
with(d, x/tapply(x, list(g1,g2), FUN=max)[cbind(g1,g2)])
[1] 0.333 0.667 1.000 0.800 1.000
[6] 1.000 0.700 0.800 0.900 1.000 That
would fail if g1 and g2 were not factors but were integer
vectors. Try it with
di - data.frame(x=1:10,
g1=rep(11:12,each=5),
g2=rep(21:23,c(3,3,4)))
with(di, x/tapply(x, list(g1,g2), FUN=max)[cbind(g1,g2)])
Error in tapply(x, list(g1, g2), FUN = max)[cbind(g1, g2)] :
subscript out of bounds
To avoid that problem you can call tapply with no FUN to get
the indices to subscript by
with(d, x/tapply(x, list(g1,g2), FUN=max)[tapply(x,
list(g1, g2))])
[1] 0.333 0.667 1.000 0.800 1.000
[6] 1.000 0.700 0.800 0.900 1.000
The misleadingly named ave() can avoid the need to do the
subscripting after tapply but has other problems
with(d, x/ave(x, g1, g2, FUN=max))
[1] 0.333 0.667 1.000 0.800 1.000
[6] 1.000 0.700 0.800 0.900 1.000
Warning messages:
1: In FUN(X[[6L]], ...) : no non-missing arguments to max;
returning -Inf
2: In FUN(X[[6L]], ...) : no non-missing arguments to max;
returning -Inf It gives the right answer but it is calling
FUN even for the empty interaction groups. For some FUN's this would
abort the call, not just give a warning. In any case it
is a waste of time.
In either case you can also use the interaction() function to
change the multiple grouping vectors into one:
d - within(d, interaction(g1, g2, drop=TRUE))
with(d, x/ave(x, g1g2, FUN=max))
[1] 0.333 0.667 1.000 0.800 1.000
[6] 1.000 0.700 0.800 0.900 1.000
with(d, x/tapply(x, g1g2, FUN=max)[g1g2])
K.u K.u K.u K.v K.v L.v
0.333 0.667 1.000 0.800 1.000 1.000
L.w L.w L.w L.w
Re: [R] referencing last row in a column
Sorry, but now I understand what you're doing, and it works as a stand-alone in the console! :) But why doesn't it work in the xlim? __ pheno.dt$year[1] [1] 1877 pheno.dt$year[nrow(pheno.dt)] [1] 1916 plot(x,y, xlim = c(pheno.dt$year[1]:pheno.dt$year[nrow(pheno.dt)]), ylim = 50:150, xlab = Year, ylab = Julian Day) Error in plot.window(...) : invalid 'xlim' value Calls: plot - plot.default - localWindow - plot.window __ -- View this message in context: http://r.789695.n4.nabble.com/referencing-last-row-in-a-column-tp2551645p2551712.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] referencing last row in a column
OK, got it to work now, just messed up the ,' and :. Thank you! -- View this message in context: http://r.789695.n4.nabble.com/referencing-last-row-in-a-column-tp2551645p2551721.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] referencing last row in a column
On Sep 23, 2010, at 11:44 , fugelpitch wrote: Sorry, but now I understand what you're doing, and it works as a stand-alone in the console! :) But why doesn't it work in the xlim? You have a colon where a comma would work. BTW: Wouldn't xlim=range(pheno.dt$year) be more expedient? __ pheno.dt$year[1] [1] 1877 pheno.dt$year[nrow(pheno.dt)] [1] 1916 plot(x,y, xlim = c(pheno.dt$year[1]:pheno.dt$year[nrow(pheno.dt)]), ylim = 50:150, xlab = Year, ylab = Julian Day) Error in plot.window(...) : invalid 'xlim' value Calls: plot - plot.default - localWindow - plot.window __ -- View this message in context: http://r.789695.n4.nabble.com/referencing-last-row-in-a-column-tp2551645p2551712.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Peter Dalgaard Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] referencing last row in a column
Hi, xlim and ylim should be given the extremes only: plot(x,y, xlim=c(pheno.dt$year[1],pheno.dt$year[nrow(pheno.dt)]), ylim=c(50,150), xlab=Year, ylab=Julian Day) ^^ ^^ Does it work now? Ivan Le 9/23/2010 11:44, fugelpitch a écrit : Sorry, but now I understand what you're doing, and it works as a stand-alone in the console! :) But why doesn't it work in the xlim? __ pheno.dt$year[1] [1] 1877 pheno.dt$year[nrow(pheno.dt)] [1] 1916 plot(x,y, xlim = c(pheno.dt$year[1]:pheno.dt$year[nrow(pheno.dt)]), ylim = 50:150, xlab = Year, ylab = Julian Day) Error in plot.window(...) : invalid 'xlim' value Calls: plot - plot.default - localWindow - plot.window __ -- Ivan CALANDRA PhD Student University of Hamburg Biozentrum Grindel und Zoologisches Museum Abt. Säugetiere Martin-Luther-King-Platz 3 D-20146 Hamburg, GERMANY +49(0)40 42838 6231 ivan.calan...@uni-hamburg.de ** http://www.for771.uni-bonn.de http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] referencing last row in a column
Yes, thanks! :) 2010-09-23 11:55, Ivan Calandra skrev: Hi, xlim and ylim should be given the extremes only: plot(x,y, xlim=c(pheno.dt$year[1],pheno.dt$year[nrow(pheno.dt)]), ylim=c(50,150), xlab=Year, ylab=Julian Day) ^^ ^^ Does it work now? Ivan Le 9/23/2010 11:44, fugelpitch a écrit : Sorry, but now I understand what you're doing, and it works as a stand-alone in the console! :) But why doesn't it work in the xlim? __ pheno.dt$year[1] [1] 1877 pheno.dt$year[nrow(pheno.dt)] [1] 1916 plot(x,y, xlim = c(pheno.dt$year[1]:pheno.dt$year[nrow(pheno.dt)]), ylim = 50:150, xlab = Year, ylab = Julian Day) Error in plot.window(...) : invalid 'xlim' value Calls: plot - plot.default - localWindow - plot.window __ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] apply union function vectorially
Hi: Look at the examples in ?Reduce HTH, Dennis On Thu, Sep 23, 2010 at 12:31 AM, statquant2 statqu...@gmail.com wrote: Thank you very much, this solve the problem, but more generally is there a function that allow to apply a function to a list of object, applying recursively the function to each answer... mathematically for a 2 argument function f(u,v) you would like a function g doing g(u1,u2,u3,u4,u5) =f(f(f(f(u1,u2),u3),u4),u5) Thanks if you can help -- View this message in context: http://r.789695.n4.nabble.com/apply-union-function-vectorially-tp2550162p2551536.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [R-pkgs] depmixS4 1.0-0 on CRAN vignette/paper on jstatsoft.org
depmixS4 has reached some form of maturity and therefore we have bumped its version number to 1.0-0 which is now on CRAN: http://cran.r-project.org/web/packages/depmixS4/index.html depmixS4 fits hidden (latent) Markov models of multivariate, mixed categorical and continuous data, otherwise known as dependent mixture models. Responses or observations can be modeled using GLMs, and additionally with multivariate normal or multinomial distributions. The package includes an example of how to add new response distributions. There is a paper introducing the package in the Journal of Statistical Software http://www.jstatsoft.org/v36/i07 which is also included as a vignette http://cran.r-project.org/web/packages/depmixS4/index.html The development version of depmixS4 can be found on the R-forge site: https://r-forge.r-project.org/projects/depmix/ best, Ingmar Visser Maarten Speekenbrink [[alternative HTML version deleted]] ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Web forum - should I make one?: an opinion
On Tue, 21-Sep-2010 at 12:55PM -0500, Vojtěch Zeisek wrote:
[...]
| Well, if You think the niche is filled, never mind, but I think R
| should have an official web forum. It can be stackoverflow (then
| I'd expect links from main R pages to it) or I (or someone else)
| can create it. Still, I think it is good idea to create one... More
| opinions?
Here's my opinion:
I detest web forums. I find them an incredibly slow way to find
anything useful. If a topic has more than 2 people who've made
comments, it's a research project itself to work out whose comment is
being commented on at what stage. I'm not interested in becoming an
investigative journalist. For years I've been puzzled why forums are
so popular when IMHO, a mailing list is infinitely simpler, faster and
more informative.
Then I realized that a huge difference is the fact that at least 90%
of internet users have probably never seen a good simple text based
email client that is capable of displaying posts in threads. Therefor
they've never experienced what I can do so easily. Were I constrained
to using something like Outlook to read mail, an infliction which many
people have, I wouldn't find it that useful.
One of the many advantages of what we use now is the fact that it's
easy for me to remove clutter by deleting threads of posts that I'm
not interested in following and saving those I wish to keep in
specific folders for future reference and searching. The archives can
be searched if there's something I haven't saved. Though they are not
as handy as having it on my own computer, it's still far better than
the forum interface.
|
| The R-* e-mail lists are the official venues and you can read/post via
| e-mail. There are also other means of interacting with the e-mail lists
| using Gmane and Nabble.
|
| Mailman is very good tool, but not for everything. And IMHO many people are
| more used to search, ask and work on forums than on mailinglists.
I've never had the experience of finding help that way which was
remotely as prompt or helpful as what the R-help list or other mailing
lists provide. (For example, try finding anything on VSN's Forum page
to do with ASReml.) But, sad to say, it probably *is* true that many
people are more used to doing things the difficult poke-and-hope way.
There are times when I feel as though I'm using electricity when many
people have never seen electricity and are using candles and water
wheels.
YMMV evidently.
--
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
___Patrick Connolly
{~._.~} Great minds discuss ideas
_( Y )_ Average minds discuss events
(:_~*~_:) Small minds discuss people
(_)-(_) . Eleanor Roosevelt
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Ordinal mixed model
Dear Patrick, You can fit such a model with the MCMCglmm package. Have a look at the vignette of that package. install.packages(MCMCglmm) vignette(CourseNotes, package = MCMCglmm) But I'm affraid that this will require some rockclimbing upon the learning curve if you are a R novice. HTH, Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek team Biometrie Kwaliteitszorg Gaverstraat 4 9500 Geraardsbergen Belgium Research Institute for Nature and Forest team Biometrics Quality Assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 thierry.onkel...@inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey -Oorspronkelijk bericht- Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Namens Patrick Walsh Verzonden: woensdag 22 september 2010 18:29 Aan: r-help@r-project.org Onderwerp: [R] Ordinal mixed model Hello, I am trying to build a generalised linear mixed model. My dependent variable is ordinal. I have a random factor (7 individuals), and a repeated measure where the dependent variable was measured three times for each of four attempts (so the repeats are nested). I also have a few covariates. I am a complete novice in R, being used to using SPSS. SPSS lets me build an ordinal model with repeated measures, but can't include a random factor. So that is really the hurdle, is this possible in R. And is there a way to do this that could be explained to someone who has no experience in R? Any help would be greatly appreciated. Kind Regards, ptwal __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Druk dit bericht a.u.b. niet onnodig af. Please do not print this message unnecessarily. Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document. The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] non-linear integer optimization?
darckeen darckeen at hotmail.com writes: Are there any packages that do non-linear integer otimization? Looked at lpSolve but i'm pretty sure it only works with linear programming and not non-linear, tried L-BFGS-B optim after flooring all my params, works somewhat but seems really inefficient. Anything else I should look at? The most general answer is: There is no R package for MINLP problems ! But: Would have been nice if you had done the following things: - Have a look into the 'optimization' task view - Provide an example of a function you would like to optimize (can it be reduced to a quadratic or convex problem?). - Tell whether you need to do it repeatedly or only a few times (e.g., utilizing an interface to the NEOS solvers). - How many integer variables, binary or really integer, etc. (could you do your own branch-and-bound, e.g.)? Providing as little information as above makes communicating difficult. Hans Werner P. S.: The R-Forge project Integer and Nonlinear Optimization in R (RINO) will (hopefully) make some of those COIN-OR solvers available. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Web forum - should I make one?
Siterer Yihui Xie x...@yihui.name: lists are easier to post messages, but I really believe they have too many disadvantages, e.g. (relatively) difficult to search, dull interface, HTML not welcome (I don't like HTML in emails, though), no lively images, no code highlighting, attachments often get chopped off, no RSS to subscribe, admins are still fighting with spams somehow manually (?), and imagine how painful it is for us to read the threaded web archives... Yeah, the lively images are a great help. And [using Firefox] I have such a big search problem. Seriously, remembering the good ol' days of UNIX and VMS where mailing lists were the only option, I shy away from web forums cos they tend to be flashy (Powerpoint on acid) and MORE difficult to search than my mailbox. I'm also a pretty lazy person, so I enjoy getting mails into my mailbox without having to think too hard, or open web browsers unnecessarily. Not that I have a lot to contribute to the list, but I have learned a great deal and see this list as one of the most important resources that exists (thanks to all list members who teach me new stuff every single day!!!) Just my $0.02. Siri. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Creating table from data frame
That's golden Henrique, thanks a lot! Worked like a charm even with large datasets. On Tue, Sep 21, 2010 at 2:56 PM, Henrique Dallazuanna www...@gmail.comwrote: Try this: d - data.frame(A = letters[1:10], B = sample(letters[11:20]), C = sample(10)) xtabs(C ~ A + B, d) On Tue, Sep 21, 2010 at 8:39 AM, ZeMajik zema...@gmail.com wrote: Hey, I have a dataset where two columns are factors and another column consists of values. Each combination of factors can only have a single value assigned to it. I'd like to represent this as a matrix or table where the rows are the first column factors and the columns the second column factors. So that each cell a_ij in the matrix represents the associated value for the factor combination ij. When no such value exists for the combination the value should be 0. I've tried playing around with tables to get this to work, but I can't seem to get it right. I've also had little luck when trying to find a solution to this. Any help would be much appreciated! Mike [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Henrique Dallazuanna Curitiba-Paraná-Brasil 25° 25' 40 S 49° 16' 22 O [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] hdf-files
Dear All, I have data in HDF file format and would like to read it into R. I have tried the package hdf5 without success. Any ideas and suggestions?? Kind regards, Katrin -- Katrin Fleischer Vrije Universiteit Amsterdam Faculty of Earth and Life Sciences Subdepartment Hydrolgy and Geo-Environmental Sciences Room E-360 De Boelelaan 1085 1081 HV AMSTERDAM Tel: +31 20 59 87391 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Simple categorical scatter plot
New to R. I am trying to create a simple xy plot wherein the line
segment color is determined by a categorical column
The following does not change colors for me, probably because I don't
quite have a handle on either functions or value mapping syntax.
--
time - c(1, 2, 3, 7,10,11,14,16,20)
pressure - c(0,10,20,20,50,18,60,65,90)
status - c(0, 0, 1, 1, 1, 0, 3, 3, 3)
measures - c(time,pressure,status)
attach(measures)
statusColor - function (x) {
if (x==0) return (green)
if (x==1) return (orange)
if (x==2) return (pink)
if (x==3) return (red)
}
par(mfrow=c(3,2))
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l,col=statusColor(status))
plot(time,pressure,type=l)
--
Warning message:
In if (x == 0) return(green) :
the condition has length 1 and only the first element will be used
TIA,
Richard A. DeVenezia
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] referencing last row in a column
Hi the other option to get last item is tail(pheno.dt$year,1) Regards Petr r-help-boun...@r-project.org napsal dne 23.09.2010 12:02:35: Yes, thanks! :) 2010-09-23 11:55, Ivan Calandra skrev: Hi, xlim and ylim should be given the extremes only: plot(x,y, xlim=c(pheno.dt$year[1],pheno.dt$year[nrow(pheno.dt)]), ylim=c(50, 150), xlab=Year, ylab=Julian Day) ^^^^ Does it work now? Ivan Le 9/23/2010 11:44, fugelpitch a écrit : Sorry, but now I understand what you're doing, and it works as a stand-alone in the console! :) But why doesn't it work in the xlim? __ pheno.dt$year[1] [1] 1877 pheno.dt$year[nrow(pheno.dt)] [1] 1916 plot(x,y, xlim = c(pheno.dt$year[1]:pheno.dt$year[nrow(pheno.dt)]), ylim = 50:150, xlab = Year, ylab = Julian Day) Error in plot.window(...) : invalid 'xlim' value Calls: plot - plot.default - localWindow - plot.window __ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Merging two data sets
Greetings I need to merge two data sets. The first data set contains information on individual organisms roosting in clusters,e.g., cluster one has 5 individuals: 3 females, two males, and one juvenile and so on for hundreds of clusters. The second data set contains temperature data on each of the plots from the first data set. The problem: the second data set only has one data point per cluster and I need to match or duplicate those temperature data with each individual organism in each cluster from the first data set. I hope this explanation is clear. How do I go about doing this? Thanks in advance. Cheers Kurt __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Odp: Simple categorical scatter plot
Hi
I am not sure but you probably need to use segments.
?segments
AFAIK for line there is only one colour for whole line possible.
and in your function statusColor you shall not use if but ifelse or better
statusColor - function (x) {
c(green,orange,pink,red)[x+1]
}
based on fact that x is vector of integers from 0 to 3.
Regards
Petr
r-help-boun...@r-project.org napsal dne 23.09.2010 13:38:43:
New to R. I am trying to create a simple xy plot wherein the line
segment color is determined by a categorical column
The following does not change colors for me, probably because I don't
quite have a handle on either functions or value mapping syntax.
--
time - c(1, 2, 3, 7,10,11,14,16,20)
pressure - c(0,10,20,20,50,18,60,65,90)
status - c(0, 0, 1, 1, 1, 0, 3, 3, 3)
measures - c(time,pressure,status)
attach(measures)
statusColor - function (x) {
if (x==0) return (green)
if (x==1) return (orange)
if (x==2) return (pink)
if (x==3) return (red)
}
par(mfrow=c(3,2))
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l,col=statusColor(status))
plot(time,pressure,type=l)
--
Warning message:
In if (x == 0) return(green) :
the condition has length 1 and only the first element will be used
TIA,
Richard A. DeVenezia
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Odp: Merging two data sets
Hi ?merge see all.x, all.y parameters Regards Petr r-help-boun...@r-project.org napsal dne 23.09.2010 13:59:18: Greetings I need to merge two data sets. The first data set contains information on individual organisms roosting in clusters,e.g., cluster one has 5 individuals: 3 females, two males, and one juvenile and so on for hundreds of clusters. The second data set contains temperature data on each of the plots from the first data set. The problem: the second data set only has one data point per cluster and I need to match or duplicate those temperature data with each individual organism in each cluster from the first data set. I hope this explanation is clear. How do I go about doing this? Thanks in advance. Cheers Kurt __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Simple categorical scatter plot
Something like this ?
time - c(1, 2, 3, 7,10,11,14,16,20)
pressure - c(0,10,20,20,50,18,60,65,90)
status - c(0, 0, 1, 1, 1, 0, 3, 3, 3)
statusColor - c(green, orange, pink, red)
plot(time, pressure)
for (i in 2:length(time)) lines(time[(i-1):i], pressure[(i-1):i],
col=statusColor[status[i]+1])
Michael
On 23 September 2010 21:38, Richard DeVenezia rdevene...@gmail.com wrote:
New to R. I am trying to create a simple xy plot wherein the line
segment color is determined by a categorical column
The following does not change colors for me, probably because I don't
quite have a handle on either functions or value mapping syntax.
--
time - c(1, 2, 3, 7,10,11,14,16,20)
pressure - c(0,10,20,20,50,18,60,65,90)
status - c(0, 0, 1, 1, 1, 0, 3, 3, 3)
measures - c(time,pressure,status)
attach(measures)
statusColor - function (x) {
if (x==0) return (green)
if (x==1) return (orange)
if (x==2) return (pink)
if (x==3) return (red)
}
par(mfrow=c(3,2))
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l,col=statusColor(status))
plot(time,pressure,type=l)
--
Warning message:
In if (x == 0) return(green) :
the condition has length 1 and only the first element will be used
TIA,
Richard A. DeVenezia
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] local linear and local constant kernel regression with np
Hi there, I ran into a weird problem using the np-package doing some local linear kernel regression. Whenever I run the function npregbw(...) with the option regtype=ll (local linear modelling), my optimal bandwidth is supposed to be 1278946. This is kind of funny, because my regressor data (189 data points) only runs from about 3.4 to about 5.9. So, what I get as a result is a nice and straight line, the same one I would get if I would run a normal linear regression. Whenever I set the option regtype=lc (i.e. local constant modelling), the optimal bandwidth is calculated as 0.795 - which sounds right to me. Any ideas / similar experiences? Thanks! Philipp __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Simple categorical scatter plot
Hi:
time - c(1, 2, 3, 7,10,11,14,16,20)
pressure - c(0,10,20,20,50,18,60,65,90)
status - c(0, 0, 1, 1, 1, 0, 3, 3, 3)
# You want to combine the vectors into a data frame instead
# of concatenating them into one long vector.
df - data.frame(time, pressure, status)
df
time pressure status
110 0
22 10 0
33 20 1
47 20 1
5 10 50 1
6 11 18 0
7 14 60 3
8 16 65 3
9 20 90 3
library(lattice)
mycols - c('green', 'orange', 'red')
xyplot(pressure ~ time, data = df, group = status, type = 'b', col =
mycols)
For this type of graphic, the lattice package works well. The package has
its own book and extensive help pages.
On Thu, Sep 23, 2010 at 4:38 AM, Richard DeVenezia rdevene...@gmail.comwrote:
New to R. I am trying to create a simple xy plot wherein the line
segment color is determined by a categorical column
The following does not change colors for me, probably because I don't
quite have a handle on either functions or value mapping syntax.
--
time - c(1, 2, 3, 7,10,11,14,16,20)
pressure - c(0,10,20,20,50,18,60,65,90)
status - c(0, 0, 1, 1, 1, 0, 3, 3, 3)
measures - c(time,pressure,status)
# A bad idea. Don't use attach(); use the function with(dataframe,
expression) instead; see ?with
attach(measures)
statusColor - function (x) {
if (x==0) return (green)
if (x==1) return (orange)
if (x==2) return (pink)
if (x==3) return (red)
}
par(mfrow=c(3,2))
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l,col=statusColor(status))
plot(time,pressure,type=l)
--
Warning message:
In if (x == 0) return(green) :
the condition has length 1 and only the first element will be used
HTH,
Dennis
TIA,
Richard A. DeVenezia
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Simple categorical scatter plot
On 09/23/2010 09:38 PM, Richard DeVenezia wrote:
New to R. I am trying to create a simple xy plot wherein the line
segment color is determined by a categorical column
The following does not change colors for me, probably because I don't
quite have a handle on either functions or value mapping syntax.
--
time- c(1, 2, 3, 7,10,11,14,16,20)
pressure- c(0,10,20,20,50,18,60,65,90)
status- c(0, 0, 1, 1, 1, 0, 3, 3, 3)
measures- c(time,pressure,status)
attach(measures)
statusColor- function (x) {
if (x==0) return (green)
if (x==1) return (orange)
if (x==2) return (pink)
if (x==3) return (red)
}
par(mfrow=c(3,2))
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l,col=statusColor(status))
plot(time,pressure,type=l)
--
Warning message:
In if (x == 0) return(green) :
the condition has length 1 and only the first element will be used
Hi Richard,
Look at the color.scale.lines function in the plotrix package.
Jim
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] NetCDF file: adding a variable
I am trying to create a NetCDF file with bathymetry (a matrix) and then to add a grid type (an integer) as variables. This is the relevant part of the code: library(ncdf) wfile=my_file.nc msv=-10 grdtp=2 # Dimensions d1=dim.def.ncdf(lon,degrees,as.double(lon)) d2=dim.def.ncdf(lat,degrees,as.double(lat)) # Variables bathymetry=var.def.ncdf(bathymetry,m,list(d1,d2),msv,longname=Bathymetry) ncw=create.ncdf(writefile,list(bathymetry,loncp,latcp,dlonp,dlatp)) put.var.ncdf(ncw,bathymetry,bat_matrix) close.ncdf(ncw) # Here I am trying to add another variable which should have a value of grdtp ncw_old=open.ncdf(wfile,write=TRUE) d3=dim.def.ncdf(grid_type,'',1:1,create_dimvar=FALSE) grid_t=var.def.ncdf(grid_type,type,d3,1.e30,longname=Grid Type) ncw_new=var.add.ncdf(ncw_old,grid_t) Here I stop because R gives an error message: Error in var.add.ncdf(ncw_old, grid_t) : Error in create.ncdf, defining var! Any ideas of what am I doing wrong? Thank you! [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Fitting to a sum
I've run into a problem with a fitting procedure I would like R to solve for me. Basically I have to fit some data to an equation which includes a sum within the formula e.g. Y(x;a,b,c) = f_i(x;a,b,c_i) + m*x where a,b are unknown and f_i(x) is the sum of another function over a known interval in which c_i is changing for each element of the sum. So in the end i will be fitting to an equation that will look something like: Y(x) = f_1(x,a,b,c_1) + f_2(x,a,b,c_2) + f_3(x,a,b,c_3) + ... + f_n(x,a,b,c_n) Any tips on how to approach this with R? In a symbolic language such as mathematica I could express Y(x) using a loop and then solve for a and b, but struggling how to do this in R. Thanks, Paul -- View this message in context: http://r.789695.n4.nabble.com/Fitting-to-a-sum-tp2551960p2551960.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Newey West and Singular Matrix + library(sandwich)
thank you, achim. I will try chol2inv.
sandwich is a very nice package, but let me make some short
suggestions. I am not a good econometrician, so I do not know what
prewhitening is, and the vignette did not explain it. ?coeftest did
not work after I loaded the library. automatic bandwidth selection
can be a good thing, but is not always.
as to my own little function, I like the idea of specifying my choice
of overlapping periods. for me, the need often arises in a case of
regressions in which the variables are measured over overlapping
intervals, so I know the number of periods to worry about. besides,
my function has an attractive simplicity to it. it is short.
assert does indeed not exist in R, but it should be. YMMV.
assert - function( cond, ... ) { if (!cond) cat(...,
file=stderr()); stopifnot(cond) }
thanks again.
/iaw
On Wed, Sep 22, 2010 at 7:41 PM, Achim Zeileis achim.zeil...@uibk.ac.at wrote:
On Wed, 22 Sep 2010, ivo welch wrote:
dear R experts: I am writing my own little newey-west standard error
function, with heteroskedasticity and arbitrary x period
autocorrelation corrections. including my function in this post here
may help others searching for something similar. it is working quite
well, except on occasion, it complains that
Error in solve.default(crossprod(x.na.omitted, x.na.omitted)) :
system is computationally singular: reciprocal condition number =
3.63797e-23
I know that lm can do the inversion, so I presume that there is a more
stable way than qr.solve . I looked into lm, then into lm.fit, and it
seems to invoke dqrls . is this the recommended way, or is there a
higher-level more stable matrix inversion routine that I could use?
I typically leverage chol2inv(). In strucchange I've written a small
convenience function solveCrossprod() that provides two different approaches
(plus the naive method). The man page has a few comments. The function
defintions are all one-liners, though.
help is, as always, appreciated. (also, if you see something else
silly in my code, let me know, please.)
1. There is no assert() function in base R.
2. The error message of se.neweywest() refers to se.white.
3. A much more flexible and powerful solution of this is included
in package sandwich, see the vignette for details. The code
sqrt(diag(NeweyWest(lm_object, lag = 0, prewhite = 0)))
replicates se.neweywest(lm_object) but has the following advantages:
it also does automatic bandwidth selection, it does not require setting
x = TRUE, it incorporates other kernel weighting functions, supports
prewhitening etc.
Best,
Z
regards,
/iaw
se.neweywest - function( lmobject.withxtrue, ar.terms =0 ) {
assert( (class(lmobject.withxtrue)==lm),
[se.white] works only on 'lm' objects, not on ,
class(lmobject.withxtrue), objects \n )
x.na.omitted - lmobject.withxtrue$x
assert( class(x.na.omitted)==matrix, [se.white] internal
error---have no X matrix. did you invoke with 'x=TRUE'?\n)
r.na.omitted - residuals( lmobject.withxtrue )
diagband.matrix - function( m, ar.terms ) {
nomit - m-ar.terms-1
mm - matrix( TRUE, nrow=m, ncol=m )
mm[1:nomit,(ncol(mm)-nomit+1):ncol(mm)] - (lower.tri(
matrix(TRUE, nrow=nomit, ncol=nomit) ))
mm[(ncol(mm)-nomit+1):ncol(mm),1:nomit] - (upper.tri(
matrix(TRUE, nrow=nomit, ncol=nomit) ))
mm
}
## V(b) = inv(X'X) X' diag(e^2) X inv(X'X)
invx - qr.solve( crossprod( x.na.omitted, x.na.omitted ) )
if (!ar.terms) resid.matrix - diag( r.na.omitted^2 ) else {
full - r.na.omitted %*% t(r.na.omitted)
## the following is not particularly good. see, we could zero out also
## items which are multiplications with a missing residual for example,
## if we do an ar1 correction, and obs 5 is missing, then the AR term on
## 4 and 6 could be set to 0. right now, we just adjust for an add'l
## term.
maskmatrix - diagband.matrix( length(r.na.omitted), ar.terms )
resid.matrix - full * maskmatrix
}
invx.x - invx %*% t(x.na.omitted)
vmat - invx.x %*% resid.matrix %*% t(invx.x)
sqrt(diag(vmat))
}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Error: attempt to apply non-function
This code worked fine for me, then did some cleaning up of formatting using ESS
(Emacs) and now I get this error, no idea what is causing it, all the
brackets/parentheses seem to be balanced. What have I done wrong?
Thanks
Jim
p0.trial01 - 0.25
TruOR01 - 0.80
num.patients.01 - 50
num.trials.01 - 5
LOR01.het.in - 0.00
num.sims - 1
simLOR01 - vector(length=num.trials.01)
simLORSE01 - vector(length=num.trials.01)
simOR01 - vector(length=num.trials.01)
trialnum01 - vector(length=num.trials.01)
x0count - vector(length = num.trials.01)
x1count - vector(length = num.trials.01)
## Trial 1, comparison of treatment 1 with treatment 0
for (i in 1:num.trials.01) {
het01 - rnorm(1,0,LOR01.het.in)
log.odds.ratio.01 - log(TruOR01) + het01
odds.ratio.01 - exp (log.odds.ratio.01)
p1.trial01 - (p0.trial01 * odds.ratio.01) / (1 - p0.trial01 +
(p0.trial01 * odds.ratio.01))
x0.trial01[i] - rbinom(1,num.patients.01,p0.trial01)
x1.trial01[i] - rbinom(1,num.patients.01,p1.trial01)
trialnum01[i] - paste ( c (trial01-), i, sep=)
simLOR01[i] - log (x1.trial01 * (num.patients.01 - x0.trial01) /
((num.patients.01 - x1.trial01) * x0.trial01))
simLORSE01[i] - sqrt ( (1/x0.trial01) + (1/(num.patients.01 -
x0.trial01)) (1 / x1.trial01)
+ (1 / (num.patients.01 -
x1.trial01)))
simOR01[i] - (x1.trial01 * (num.patients.01 - x0.trial01) /
((num.patients.01 - x1.trial01) * x0.trial01))
}
## Output all results to a data.frame called results
results - data.frame (tn1 = trialnum01, SimOR01 = simOR01,
SimLOR01 = simLOR01, SimLORSE01 = simLORSE01,
p1trial = p1.trial01, x0count, x1count )
set.seed(9321685)
results
rm(list=ls())
===
Dr. Jim Maas
University of East Anglia
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] hdf files
Dear All, I want to upload data in HDF-file format into R. I have tried the package 'hdf5' without success. I presume that my files are in the hdf4 format and therefore are not readable with this package. I understand that its possible to convert hdf4 into hdf5 format using C, but I was wondering whether there is an alternative option for importing hdf4 files into R? kind regards, Katrin -- Katrin Fleischer Vrije Universiteit Amsterdam Faculty of Earth and Life Sciences Subdepartment Hydrolgy and Geo-Environmental Sciences Room E-360 De Boelelaan 1085 1081 HV AMSTERDAM Tel: +31 20 59 87391 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] problem opening pdf device on Windows 7
Peter Dalgaard pdalgd at gmail.com writes: Then I tried: pdf() Error in pdf() : unable to start device pdf In addition: Warning message: In pdf() : cannot open 'pdf' file argument 'Rplots.pdf' And your working directory is writable for you? Otherwise, Change dir from the File menu to somewhere that is. Thanks, Peter. I was using Tinn-R and default working directory set for the latest version of R: C:\Program Files\R\R-2.11.1\bin Changing the working directory worked. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] merging multiple data frames
hi guys i have multiple data frames which i want to merge. there are four of them..eg pdf SampleID UVDose_J RepairHours Day_0 Day_45 Day_90 1SDM001 1.0 3 485.612 465.142 490.873 2SDM001 1.0 3 503.658 457.863 487.783 3SDM001 1.0 2 533.193 451.044 456.973 4SDM001 1.0 2 538.334 452.887 474.915 5SDM001 1.0 1 526.034 481.123 477.801 6SDM001 1.0 1 546.543 472.322 481.546 7SDM001 1.0 0 NA NA NA 8SDM001 1.0 0 NA NA NA 9SDM001 0.5 3 432.134 457.245 497.975 10 SDM001 0.5 3 432.605 450.184 489.468 11 SDM001 0.5 2 450.335 496.520 488.784 12 SDM001 0.5 2 439.590 474.371 470.182 13 SDM001 0.5 1 510.480 489.561 525.029 14 SDM001 0.5 1 487.934 467.258 488.784 15 SDM001 0.5 0 NA NA NA 16 SDM001 0.5 0 NA NA NA 20 SDM002 1.0 3 465.549 528.715 501.374 21 SDM002 1.0 3 458.168 505.480 489.244 22 SDM002 1.0 2 447.317 464.009 478.058 23 SDM002 1.0 2 452.020 438.446 470.996 24 SDM002 1.0 1 441.718 458.760 499.221 25 SDM002 1.0 1 447.017 402.616 548.797 26 SDM002 1.0 0 NA NA NA 27 SDM002 1.0 0 NA NA NA 28 SDM002 0.5 3 421.409 448.870 476.392 29 SDM002 0.5 3 404.089 446.413 477.080 30 SDM002 0.5 2 399.775 432.678 465.015 31 SDM002 0.5 2 427.157 443.418 477.048 32 SDM002 0.5 1 389.674 449.353 482.264 33 SDM002 0.5 1 418.147 455.983 495.486 34 SDM002 0.5 0 NA NA NA 35 SDM002 0.5 0 NA NA NA 39 SDM005 1.0 3 579.836 441.040 476.382 40 SDM005 1.0 3 578.525 443.875 472.867 41 SDM005 1.0 2 564.266 432.116 469.416 42 SDM005 1.0 2 571.045 447.658 458.233 43 SDM005 1.0 1 564.664 427.673 524.122 44 SDM005 1.0 1 568.182 458.039 477.237 45 SDM005 1.0 0 NA NA NA 46 SDM005 1.0 0 NA NA NA 47 SDM005 0.5 3 556.534 424.786 501.658 48 SDM005 0.5 3 474.027 441.418 507.635 49 SDM005 0.5 2 481.355 430.346 468.021 50 SDM005 0.5 2 478.922 466.933 471.025 51 SDM005 0.5 1 505.539 937.759 460.985 52 SDM005 0.5 1 497.913 457.932 493.152 53 SDM005 0.5 0 NA NA NA 54 SDM005 0.5 0 NA NA NA 58 SDM006 1.0 3 589.164 459.578 509.565 59 SDM006 1.0 3 608.477 480.233 519.785 60 SDM006 1.0 2 598.354 449.266 487.058 61 SDM006 1.0 2 617.823 456.908 507.467 62 SDM006 1.0 1 566.477 500.189 526.744 63 SDM006 1.0 1 622.170 462.463 550.675 64 SDM006 1.0 0 NA NA NA 65 SDM006 1.0 0 NA NA NA 66 SDM006 0.5 3 546.472 457.880 468.129 67 SDM006 0.5 3 525.069 444.575 505.154 68 SDM006 0.5 2 569.068 446.196 473.739 69 SDM006 0.5 2 534.205 470.366 476.570 bdf SampleID UVDose_J RepairHoursDay_0 Day_45 Day_90 17SDM001 0.5 B 88.6145 388.3575 198.467 36SDM002 0.5 B 100.0760 384.9505 234.740 55SDM005 0.5 B 121.9595 300.3650 241.832 74SDM006 0.5 B 174.7005 378.3435 291.272 93SDM007 0.5 B 319.7750 335.4390 110.306 112 SDM008 0.5 B 292.8400 323.0370 172.615 131 SDM010 0.5 B 112.0225 281.0390 562.459 150 SDM011 0.5 B 125.0440 331.4650 230.026 169 SDM012 0.5 B 229.1310 264.5580 234.231 188 SDM013 0.5 B 212.9690 524.0240 420.211 207 SDM014 0.5 B 366.3350 225.0610 203.588 226 SDM016 0.5 B 305.6080 381.5770 155.052 245 SDM017 0.5 B 223.6260 281.7830 182.374 264 SDM018 0.5 B 200.7720 401.6350 193.573 283 SDM019 0.5 B 164.2360 156.9960 175.896 302 SDM023 0.5 B 198.8900 210.0600 215.629 321 SDM024 0.5 B 351.8460 100.0980 185.388 340 SDM026 0.5 B 190.4560 132.8970 160.213 359 SDM028 0.5 B 252.9760 158.2910 120.425 378 SDM029 0.5 B 411.0690 134.1060 138.528 tdf SampleID UVDose_J RepairHours Day_0 Day_45 Day_90 19SDM001 0.5 T 642.3255 579.6635 537.581 38SDM002 0.5 T 531.2000
Re: [R] Newey West and Singular Matrix + library(sandwich)
On Thu, 23 Sep 2010, ivo welch wrote:
thank you, achim. I will try chol2inv.
sandwich is a very nice package, but let me make some short
suggestions. I am not a good econometrician, so I do not know what
prewhitening is,
In general it means, that you do not look at a (typically multivariate)
series Y, but at the residuals of a VAR(p) model for Y. In case of
HAC covariances, people do not use the empirical estimating functions
directly, but the residuals of a VAR(1) model.
and the vignette did not explain it. ?coeftest did
not work after I loaded the library.
Presumably because you did load sandwich but did not load the lmtest
package in which coeftest() is provided.
automatic bandwidth selection can be a good thing, but is not always.
As my short command in the previous mail showed, it can be conveniently
set to any other value, just by specifying the lag argument.
as to my own little function, I like the idea of specifying my choice
of overlapping periods. for me, the need often arises in a case of
regressions in which the variables are measured over overlapping
intervals, so I know the number of periods to worry about. besides,
my function has an attractive simplicity to it. it is short.
assert does indeed not exist in R, but it should be. YMMV.
But the code would be even shorter without it ;-) Just kidding...
hth,
Z
assert - function( cond, ... ) { if (!cond) cat(...,
file=stderr()); stopifnot(cond) }
thanks again.
/iaw
On Wed, Sep 22, 2010 at 7:41 PM, Achim Zeileis achim.zeil...@uibk.ac.at wrote:
On Wed, 22 Sep 2010, ivo welch wrote:
dear R experts: I am writing my own little newey-west standard error
function, with heteroskedasticity and arbitrary x period
autocorrelation corrections. including my function in this post here
may help others searching for something similar. it is working quite
well, except on occasion, it complains that
Error in solve.default(crossprod(x.na.omitted, x.na.omitted)) :
system is computationally singular: reciprocal condition number =
3.63797e-23
I know that lm can do the inversion, so I presume that there is a more
stable way than qr.solve . I looked into lm, then into lm.fit, and it
seems to invoke dqrls . is this the recommended way, or is there a
higher-level more stable matrix inversion routine that I could use?
I typically leverage chol2inv(). In strucchange I've written a small
convenience function solveCrossprod() that provides two different approaches
(plus the naive method). The man page has a few comments. The function
defintions are all one-liners, though.
help is, as always, appreciated. (also, if you see something else
silly in my code, let me know, please.)
1. There is no assert() function in base R.
2. The error message of se.neweywest() refers to se.white.
3. A much more flexible and powerful solution of this is included
in package sandwich, see the vignette for details. The code
sqrt(diag(NeweyWest(lm_object, lag = 0, prewhite = 0)))
replicates se.neweywest(lm_object) but has the following advantages:
it also does automatic bandwidth selection, it does not require setting
x = TRUE, it incorporates other kernel weighting functions, supports
prewhitening etc.
Best,
Z
regards,
/iaw
se.neweywest - function( lmobject.withxtrue, ar.terms =0 ) {
assert( (class(lmobject.withxtrue)==lm),
[se.white] works only on 'lm' objects, not on ,
class(lmobject.withxtrue), objects \n )
x.na.omitted - lmobject.withxtrue$x
assert( class(x.na.omitted)==matrix, [se.white] internal
error---have no X matrix. did you invoke with 'x=TRUE'?\n)
r.na.omitted - residuals( lmobject.withxtrue )
diagband.matrix - function( m, ar.terms ) {
nomit - m-ar.terms-1
mm - matrix( TRUE, nrow=m, ncol=m )
mm[1:nomit,(ncol(mm)-nomit+1):ncol(mm)] - (lower.tri(
matrix(TRUE, nrow=nomit, ncol=nomit) ))
mm[(ncol(mm)-nomit+1):ncol(mm),1:nomit] - (upper.tri(
matrix(TRUE, nrow=nomit, ncol=nomit) ))
mm
}
## V(b) = inv(X'X) X' diag(e^2) X inv(X'X)
invx - qr.solve( crossprod( x.na.omitted, x.na.omitted ) )
if (!ar.terms) resid.matrix - diag( r.na.omitted^2 ) else {
full - r.na.omitted %*% t(r.na.omitted)
## the following is not particularly good. see, we could zero out also
## items which are multiplications with a missing residual for example,
## if we do an ar1 correction, and obs 5 is missing, then the AR term on
## 4 and 6 could be set to 0. right now, we just adjust for an add'l
## term.
maskmatrix - diagband.matrix( length(r.na.omitted), ar.terms )
resid.matrix - full * maskmatrix
}
invx.x - invx %*% t(x.na.omitted)
vmat - invx.x %*% resid.matrix %*% t(invx.x)
sqrt(diag(vmat))
}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Plotting densities
Hi group, I am currently plotting two densities using the following code: x1 - c(1,2,1,3,5,6,6,7,7,8) x2 - c(1,2,1,3,5,6,5,7) plot(density(x1, na.rm = TRUE)) polygon(density(x2, na.rm = TRUE), border=blue) However, I would like to avoid bordering the second density as it adds a nasty bottom line which I would like to avoid. I would also rather have a dashed or dotted line for the second (currently blue) density but without the bottom part. Any idea how to do that? Best, Ralf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Error: attempt to apply non-function
On Sep 23, 2010, at 14:54 , Maas James Dr (MED) wrote:
This code worked fine for me, then did some cleaning up of formatting using
ESS (Emacs) and now I get this error, no idea what is causing it, all the
brackets/parentheses seem to be balanced. What have I done wrong?
Look for ) ( -- missing a +, presumably. As a general matter, try using
tools like traceback() or options(recover=...) to catch such things in the act.
They are quite hard to spot by visual inspection.
-pd
Thanks
Jim
p0.trial01 - 0.25
TruOR01 - 0.80
num.patients.01 - 50
num.trials.01 - 5
LOR01.het.in - 0.00
num.sims - 1
simLOR01 - vector(length=num.trials.01)
simLORSE01 - vector(length=num.trials.01)
simOR01 - vector(length=num.trials.01)
trialnum01 - vector(length=num.trials.01)
x0count - vector(length = num.trials.01)
x1count - vector(length = num.trials.01)
## Trial 1, comparison of treatment 1 with treatment 0
for (i in 1:num.trials.01) {
het01 - rnorm(1,0,LOR01.het.in)
log.odds.ratio.01 - log(TruOR01) + het01
odds.ratio.01 - exp (log.odds.ratio.01)
p1.trial01 - (p0.trial01 * odds.ratio.01) / (1 - p0.trial01 +
(p0.trial01 * odds.ratio.01))
x0.trial01[i] - rbinom(1,num.patients.01,p0.trial01)
x1.trial01[i] - rbinom(1,num.patients.01,p1.trial01)
trialnum01[i] - paste ( c (trial01-), i, sep=)
simLOR01[i] - log (x1.trial01 * (num.patients.01 - x0.trial01) /
((num.patients.01 - x1.trial01) * x0.trial01))
simLORSE01[i] - sqrt ( (1/x0.trial01) + (1/(num.patients.01 -
x0.trial01)) (1 / x1.trial01)
+ (1 / (num.patients.01 -
x1.trial01)))
simOR01[i] - (x1.trial01 * (num.patients.01 - x0.trial01) /
((num.patients.01 - x1.trial01) * x0.trial01))
}
## Output all results to a data.frame called results
results - data.frame (tn1 = trialnum01, SimOR01 = simOR01,
SimLOR01 = simLOR01, SimLORSE01 = simLORSE01,
p1trial = p1.trial01, x0count, x1count )
set.seed(9321685)
results
rm(list=ls())
===
Dr. Jim Maas
University of East Anglia
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk Priv: pda...@gmail.com
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Error: attempt to apply non-function
On 23/09/2010 8:54 AM, Maas James Dr (MED) wrote:
This code worked fine for me, then did some cleaning up of formatting using ESS
(Emacs) and now I get this error, no idea what is causing it, all the
brackets/parentheses seem to be balanced. What have I done wrong?
It would help a lot if you posted the actual error message, and the
results of a call to traceback() just after it. We can't run your code
(we don't have x0.trial01 and don't know how you created it).
Duncan Murdoch
Thanks
Jim
p0.trial01- 0.25
TruOR01- 0.80
num.patients.01- 50
num.trials.01- 5
LOR01.het.in- 0.00
num.sims- 1
simLOR01- vector(length=num.trials.01)
simLORSE01- vector(length=num.trials.01)
simOR01- vector(length=num.trials.01)
trialnum01- vector(length=num.trials.01)
x0count- vector(length = num.trials.01)
x1count- vector(length = num.trials.01)
## Trial 1, comparison of treatment 1 with treatment 0
for (i in 1:num.trials.01) {
het01- rnorm(1,0,LOR01.het.in)
log.odds.ratio.01- log(TruOR01) + het01
odds.ratio.01- exp (log.odds.ratio.01)
p1.trial01- (p0.trial01 * odds.ratio.01) / (1 - p0.trial01 +
(p0.trial01 * odds.ratio.01))
x0.trial01[i]- rbinom(1,num.patients.01,p0.trial01)
x1.trial01[i]- rbinom(1,num.patients.01,p1.trial01)
trialnum01[i]- paste ( c (trial01-), i, sep=)
simLOR01[i]- log (x1.trial01 * (num.patients.01 - x0.trial01) /
((num.patients.01 - x1.trial01) * x0.trial01))
simLORSE01[i]- sqrt ( (1/x0.trial01) + (1/(num.patients.01 -
x0.trial01)) (1 / x1.trial01)
+ (1 / (num.patients.01 -
x1.trial01)))
simOR01[i]- (x1.trial01 * (num.patients.01 - x0.trial01) /
((num.patients.01 - x1.trial01) * x0.trial01))
}
## Output all results to a data.frame called results
results- data.frame (tn1 = trialnum01, SimOR01 = simOR01,
SimLOR01 = simLOR01, SimLORSE01 = simLORSE01,
p1trial = p1.trial01, x0count, x1count )
set.seed(9321685)
results
rm(list=ls())
===
Dr. Jim Maas
University of East Anglia
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Odp: Plotting densities
Hi r-help-boun...@r-project.org napsal dne 23.09.2010 15:10:04: Hi group, I am currently plotting two densities using the following code: x1 - c(1,2,1,3,5,6,6,7,7,8) x2 - c(1,2,1,3,5,6,5,7) plot(density(x1, na.rm = TRUE)) polygon(density(x2, na.rm = TRUE), border=blue) However, I would like to avoid bordering the second density as it adds a nasty bottom line which I would like to avoid. I would also rather have a dashed or dotted line for the second (currently blue) density but without the bottom part. Any idea how to do that? Something like lines(density(x2, na.rm = TRUE), col=blue) Regards Petr Best, Ralf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] extending survival curves past the last event using plot.survfit
On Sep 22, 2010, at 8:15 PM, Krambrink, Amy M wrote:
Hello,
I'm using plot.survfit to plot cumulative incidence of an event.
Essentially, my code boils down to:
cox -coxph(Surv(EVINF,STATUS) ~ strata(TREAT) + covariates, data=dat)
surv - survfit(cox)
plot(surv,mark.time=F,fun=event)
Follow-up time extends to 54 weeks, but the last event occurs at week
30, and no more people are censored in between. Is there a direct way
to extend the curves with a horizontal line to the end of follow-up
(54
weeks), rather than stopping at the time of the last event (30 weeks)?
The survfit.coxph function only records the survival curve at the death
times, so there is no data in the surv object about time 54. (This
will change in my next bundle of updates, which are currently under
test.) You will need to update the survival curves' output by hand.
For now this would be the easiest code for multiple strata, at least
that I can think of offhand.
plot(surv, mark.time=F, fun='event', xlim=c(0, 54))
for (i in 1:length(surv$strata)) { #number of curves
temp - surv[i]
lines(c(max(temp$time), 54), 1- rep(min(temp$surv),2))
}
Terry Therneau
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Simple categorical scatter plot
Hi:
Evidently I misunderstood your intention, so let's try again. I tried to
find a lattice solution but came up empty. I think this works instead using
base graphics. The idea is to build up vectors of start points (x0, y0),
endpoints (x1, y1) and the color of the segment between each pair of points.
The last line is to call segments, which (fortunately) takes vector
arguments. The last element is removed from the vectors to get the start
points, while the first element is removed to get the end points; if it
doesn't make sense, cbind them and it should be clearer.
x0 - pressure[-length(pressure)]
y0 - time[-length(time)]
x0 - time[-length(time)]
y0 - pressure[-length(pressure)]
x1 - time[-1]
y1 - pressure[-1]
cols - c(rep('green', 2), rep('orange', 3), 'green', rep('red', 2))
# Set up plot area, but don't plot (yet)
plot(time, pressure, type = 'n')
segments(x0, y0, x1, y1, col = cols)
Hope I got it right this time,
Dennis
PS: Something like this ought to work in lattice, too. A message to that
effect is cited here:
http://r.789695.n4.nabble.com/forum/PrintPost.jtp?post=2165353
On Thu, Sep 23, 2010 at 5:14 AM, Dennis Murphy djmu...@gmail.com wrote:
Hi:
time - c(1, 2, 3, 7,10,11,14,16,20)
pressure - c(0,10,20,20,50,18,60,65,90)
status - c(0, 0, 1, 1, 1, 0, 3, 3, 3)
# You want to combine the vectors into a data frame instead
# of concatenating them into one long vector.
df - data.frame(time, pressure, status)
df
time pressure status
110 0
22 10 0
33 20 1
47 20 1
5 10 50 1
6 11 18 0
7 14 60 3
8 16 65 3
9 20 90 3
library(lattice)
mycols - c('green', 'orange', 'red')
xyplot(pressure ~ time, data = df, group = status, type = 'b', col =
mycols)
For this type of graphic, the lattice package works well. The package has
its own book and extensive help pages.
On Thu, Sep 23, 2010 at 4:38 AM, Richard DeVenezia
rdevene...@gmail.comwrote:
New to R. I am trying to create a simple xy plot wherein the line
segment color is determined by a categorical column
The following does not change colors for me, probably because I don't
quite have a handle on either functions or value mapping syntax.
--
time - c(1, 2, 3, 7,10,11,14,16,20)
pressure - c(0,10,20,20,50,18,60,65,90)
status - c(0, 0, 1, 1, 1, 0, 3, 3, 3)
measures - c(time,pressure,status)
# A bad idea. Don't use attach(); use the function with(dataframe,
expression) instead; see ?with
attach(measures)
statusColor - function (x) {
if (x==0) return (green)
if (x==1) return (orange)
if (x==2) return (pink)
if (x==3) return (red)
}
par(mfrow=c(3,2))
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l)
plot(time,pressure,type=l,col=statusColor(status))
plot(time,pressure,type=l)
--
Warning message:
In if (x == 0) return(green) :
the condition has length 1 and only the first element will be used
HTH,
Dennis
TIA,
Richard A. DeVenezia
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] dnorm
You need to take into account the bin width as well (hence the extra multiple you asked about), but it is simpler to just include prob=TRUE in the hist call, then you do not need to do any adjustment on the y-values of the reference distribution. -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- project.org] On Behalf Of Sibylle Stöckli Sent: Thursday, September 23, 2010 1:43 AM To: r-help@r-project.org Subject: [R] dnorm Dear R-users Idea: Plot a dnorm line using specific mean/sd to complete a histogram (skewed). xs:range of y-values, ys: dnorm function Problem: I expected to multiply the ys function with the sample size (n=250- 300). I was wondering about a factor between 12'000 and 30'000 to match the size of the dnorm line with the specific histogram. Thanks Sibylle hist(Biotree[Ld,]$Height2008, main=Larix decidua, ylim=c(0,50), xlab=Tree Height 2008 (cm),col=aquamarine, font.main=3, cex.axis=0.8) xs-0:650 ys-dnorm(xs, mean=397.8, sd=97.6) lines(xs,ys*12000) -- GMX DSL SOMMER-SPECIAL: Surf Phone Flat 16.000 für nur 19,99 Euro/mtl.!* __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Plotting densities
In your call to polygon(), include lty=dashed or dotted or whatever you want your line type to look like. take a good look at all the options in ?par For everything you can customize in plots. Alternatively, Paul Murrel's R Graphics book is the best reference I know for this sort of stuff. Mike On Thu, Sep 23, 2010 at 8:10 AM, Ralf B ralf.bie...@gmail.com wrote: Hi group, I am currently plotting two densities using the following code: x1 - c(1,2,1,3,5,6,6,7,7,8) x2 - c(1,2,1,3,5,6,5,7) plot(density(x1, na.rm = TRUE)) polygon(density(x2, na.rm = TRUE), border=blue) However, I would like to avoid bordering the second density as it adds a nasty bottom line which I would like to avoid. I would also rather have a dashed or dotted line for the second (currently blue) density but without the bottom part. Any idea how to do that? Best, Ralf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.htmlhttp://www.r-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] non-linear integer optimization?
This is an example of the type of problem, and how i'm currently using
optim() to solve it.
mydata - runif(500,-1,1)
myfunc - function(p,d)
{
print(p - floor(p))
ws - function(i,n,x) sum(x[i-n+1]:x[i])
ws1 - c(rep(NA,p[1]-1),sapply(p[1]:NROW(d),ws,p[1],d))
ws2 - c(rep(NA,p[2]-1),sapply(p[2]:NROW(d),ws,p[2],d))
ws3 - c(rep(NA,p[3]-1),sapply(p[3]:NROW(d),ws,p[3],d))
var(ws1+ws2+ws3,na.rm=TRUE)
}
opt - optim(c(25,50,150),myfunc,method=L-BFGS-B,
control=list(fnscale=-1,parscale=c(1,1,1),factr=1,ndeps=c(5,5,5)),
lower=t(c(1,51,101)),upper=t(c(50,100,200)),d=mydata)
print(floor(opt$par))
print(myfunc(opt$par,mydata))
So the parameters to the function to be optimized are parameters to
functions that only accept integer values. All of the paramters to be
optimized are integers that are subject to upper lower bound constraints.
This was the solution I came up with after checking CRAN, searching nabble
etc.
It runs but not very efficiently, and does not seem to really explore the
sample space very well before focusing on a local minimum. I also looked at
the constrOptim function but I couldn't figure out how to implement two
sided constraints with it.
--
View this message in context:
http://r.789695.n4.nabble.com/non-linear-integer-optimization-tp2551409p2552080.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] merging multiple data frames
Hi, On Thu, Sep 23, 2010 at 9:04 AM, rasanpreet rasanpreet.k...@gmail.com wrote: hi guys i have multiple data frames which i want to merge. there are four of them..eg Can you provide a (correct) example of what you want your merged data.frame to look like? What column do you want to use in your data.frame to merge against? I'm guessing SampleID(?), but then again, these aren't unique in your `pdf` data.frame. For instance, what would the row for SDM001 look like in your merged data.frame? -steve pdf SampleID UVDose_J RepairHours Day_0 Day_45 Day_90 1 SDM001 1.0 3 485.612 465.142 490.873 2 SDM001 1.0 3 503.658 457.863 487.783 3 SDM001 1.0 2 533.193 451.044 456.973 4 SDM001 1.0 2 538.334 452.887 474.915 5 SDM001 1.0 1 526.034 481.123 477.801 6 SDM001 1.0 1 546.543 472.322 481.546 7 SDM001 1.0 0 NA NA NA 8 SDM001 1.0 0 NA NA NA 9 SDM001 0.5 3 432.134 457.245 497.975 10 SDM001 0.5 3 432.605 450.184 489.468 11 SDM001 0.5 2 450.335 496.520 488.784 12 SDM001 0.5 2 439.590 474.371 470.182 13 SDM001 0.5 1 510.480 489.561 525.029 14 SDM001 0.5 1 487.934 467.258 488.784 15 SDM001 0.5 0 NA NA NA 16 SDM001 0.5 0 NA NA NA 20 SDM002 1.0 3 465.549 528.715 501.374 21 SDM002 1.0 3 458.168 505.480 489.244 22 SDM002 1.0 2 447.317 464.009 478.058 23 SDM002 1.0 2 452.020 438.446 470.996 24 SDM002 1.0 1 441.718 458.760 499.221 25 SDM002 1.0 1 447.017 402.616 548.797 26 SDM002 1.0 0 NA NA NA 27 SDM002 1.0 0 NA NA NA 28 SDM002 0.5 3 421.409 448.870 476.392 29 SDM002 0.5 3 404.089 446.413 477.080 30 SDM002 0.5 2 399.775 432.678 465.015 31 SDM002 0.5 2 427.157 443.418 477.048 32 SDM002 0.5 1 389.674 449.353 482.264 33 SDM002 0.5 1 418.147 455.983 495.486 34 SDM002 0.5 0 NA NA NA 35 SDM002 0.5 0 NA NA NA 39 SDM005 1.0 3 579.836 441.040 476.382 40 SDM005 1.0 3 578.525 443.875 472.867 41 SDM005 1.0 2 564.266 432.116 469.416 42 SDM005 1.0 2 571.045 447.658 458.233 43 SDM005 1.0 1 564.664 427.673 524.122 44 SDM005 1.0 1 568.182 458.039 477.237 45 SDM005 1.0 0 NA NA NA 46 SDM005 1.0 0 NA NA NA 47 SDM005 0.5 3 556.534 424.786 501.658 48 SDM005 0.5 3 474.027 441.418 507.635 49 SDM005 0.5 2 481.355 430.346 468.021 50 SDM005 0.5 2 478.922 466.933 471.025 51 SDM005 0.5 1 505.539 937.759 460.985 52 SDM005 0.5 1 497.913 457.932 493.152 53 SDM005 0.5 0 NA NA NA 54 SDM005 0.5 0 NA NA NA 58 SDM006 1.0 3 589.164 459.578 509.565 59 SDM006 1.0 3 608.477 480.233 519.785 60 SDM006 1.0 2 598.354 449.266 487.058 61 SDM006 1.0 2 617.823 456.908 507.467 62 SDM006 1.0 1 566.477 500.189 526.744 63 SDM006 1.0 1 622.170 462.463 550.675 64 SDM006 1.0 0 NA NA NA 65 SDM006 1.0 0 NA NA NA 66 SDM006 0.5 3 546.472 457.880 468.129 67 SDM006 0.5 3 525.069 444.575 505.154 68 SDM006 0.5 2 569.068 446.196 473.739 69 SDM006 0.5 2 534.205 470.366 476.570 bdf SampleID UVDose_J RepairHours Day_0 Day_45 Day_90 17 SDM001 0.5 B 88.6145 388.3575 198.467 36 SDM002 0.5 B 100.0760 384.9505 234.740 55 SDM005 0.5 B 121.9595 300.3650 241.832 74 SDM006 0.5 B 174.7005 378.3435 291.272 93 SDM007 0.5 B 319.7750 335.4390 110.306 112 SDM008 0.5 B 292.8400 323.0370 172.615 131 SDM010 0.5 B 112.0225 281.0390 562.459 150 SDM011 0.5 B 125.0440 331.4650 230.026 169 SDM012 0.5 B 229.1310 264.5580 234.231 188 SDM013 0.5 B 212.9690 524.0240 420.211 207 SDM014 0.5 B 366.3350 225.0610 203.588 226 SDM016 0.5 B 305.6080 381.5770 155.052 245 SDM017 0.5 B 223.6260 281.7830 182.374 264 SDM018 0.5 B 200.7720 401.6350 193.573 283 SDM019 0.5
[R] How to pass a model formula as argument to with.mids
Hello I would like to pass a model formula as an argument to the with.mids function from the mice package. The with.mids functon fits models to multiply imputed data sets. Here's a simple example library(mice) #Create multiple imputations on the nhanes data contained in the mice package. imp - mice(nahnes) #Fitting a linear model with each imputed data set the regular way works fine with(imp, lm(bmi~hyp+chl)) #Creating a formula object and than passing it as argument does not work: form.obj - formula(bmi~hyp+chl) with(imp, lm(form.obj)) #The following doesn't work either expr -lm(bmi~hyp+chl) with(imp, expr) Looking at the definition of with.mids reveals that the second argument is first substituted and than evaluated within each data.frame of the multiply imputed data sets. Is there a way to pass lm(bmi~hyp+chl) without having to change the definition of the with.mids function? Thanks in advance, Erich __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] hdf-files
Hi, On Thu, Sep 23, 2010 at 7:13 AM, Katrin Fleischer katrin.fleisc...@falw.vu.nl wrote: Dear All, I have data in HDF file format and would like to read it into R. I have tried the package hdf5 without success. What type of errors are you getting? Any ideas and suggestions?? Sure, from: http://cran.r-project.org/doc/manuals/R-data.html Packages hdf5, RNetCDF and ncdf on CRAN provide interfaces to NASA's HDF5 ... I'm not familiar with the differences between hdf5 and netCDF, but perhaps those other packages can help (I see that they've been more recently updated, anyway ...). -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] For loop with ifelse help
Hi Pele, On Wed, Sep 22, 2010 at 12:40 PM, Pele drdi...@yahoo.com wrote: Hi David - thanks for your suggestion, but I am trying to avoid doing any merging and sorting for this step because the real file I will be working with has about 20 million records. If I can get this loop or something similar to work will be good enough. If that's the case, you might consider looking at the sqldf or data.table packages. They both implement data.frame-like objects, but can do subsetting (and merging) rather quickly since they implement indexes over keys (columns) of the respective data.frame(s). Subsetting normal data.frames in this scenario you describe involves a linear search for every query through the column(s) you are querying against, which can get slow as the size of your data.frames get large. -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] merging multiple data frames
First, you might want to start by generating a new column to identify your 'pdf and bdf or whatever once it's merged. For the merging, see ?merge But as someone's already pointed out, it's not clear what you are trying to merge by. Also, as your example calculations show, you don't need to merge it to to the calcuations you want to do... On Thu, Sep 23, 2010 at 8:53 AM, Steve Lianoglou mailinglist.honey...@gmail.com wrote: Hi, On Thu, Sep 23, 2010 at 9:04 AM, rasanpreet rasanpreet.k...@gmail.com wrote: hi guys i have multiple data frames which i want to merge. there are four of them..eg Can you provide a (correct) example of what you want your merged data.frame to look like? What column do you want to use in your data.frame to merge against? I'm guessing SampleID(?), but then again, these aren't unique in your `pdf` data.frame. For instance, what would the row for SDM001 look like in your merged data.frame? -steve pdf SampleID UVDose_J RepairHours Day_0 Day_45 Day_90 1SDM001 1.0 3 485.612 465.142 490.873 2SDM001 1.0 3 503.658 457.863 487.783 3SDM001 1.0 2 533.193 451.044 456.973 4SDM001 1.0 2 538.334 452.887 474.915 5SDM001 1.0 1 526.034 481.123 477.801 6SDM001 1.0 1 546.543 472.322 481.546 7SDM001 1.0 0 NA NA NA 8SDM001 1.0 0 NA NA NA 9SDM001 0.5 3 432.134 457.245 497.975 10 SDM001 0.5 3 432.605 450.184 489.468 11 SDM001 0.5 2 450.335 496.520 488.784 12 SDM001 0.5 2 439.590 474.371 470.182 13 SDM001 0.5 1 510.480 489.561 525.029 14 SDM001 0.5 1 487.934 467.258 488.784 15 SDM001 0.5 0 NA NA NA 16 SDM001 0.5 0 NA NA NA 20 SDM002 1.0 3 465.549 528.715 501.374 21 SDM002 1.0 3 458.168 505.480 489.244 22 SDM002 1.0 2 447.317 464.009 478.058 23 SDM002 1.0 2 452.020 438.446 470.996 24 SDM002 1.0 1 441.718 458.760 499.221 25 SDM002 1.0 1 447.017 402.616 548.797 26 SDM002 1.0 0 NA NA NA 27 SDM002 1.0 0 NA NA NA 28 SDM002 0.5 3 421.409 448.870 476.392 29 SDM002 0.5 3 404.089 446.413 477.080 30 SDM002 0.5 2 399.775 432.678 465.015 31 SDM002 0.5 2 427.157 443.418 477.048 32 SDM002 0.5 1 389.674 449.353 482.264 33 SDM002 0.5 1 418.147 455.983 495.486 34 SDM002 0.5 0 NA NA NA 35 SDM002 0.5 0 NA NA NA 39 SDM005 1.0 3 579.836 441.040 476.382 40 SDM005 1.0 3 578.525 443.875 472.867 41 SDM005 1.0 2 564.266 432.116 469.416 42 SDM005 1.0 2 571.045 447.658 458.233 43 SDM005 1.0 1 564.664 427.673 524.122 44 SDM005 1.0 1 568.182 458.039 477.237 45 SDM005 1.0 0 NA NA NA 46 SDM005 1.0 0 NA NA NA 47 SDM005 0.5 3 556.534 424.786 501.658 48 SDM005 0.5 3 474.027 441.418 507.635 49 SDM005 0.5 2 481.355 430.346 468.021 50 SDM005 0.5 2 478.922 466.933 471.025 51 SDM005 0.5 1 505.539 937.759 460.985 52 SDM005 0.5 1 497.913 457.932 493.152 53 SDM005 0.5 0 NA NA NA 54 SDM005 0.5 0 NA NA NA 58 SDM006 1.0 3 589.164 459.578 509.565 59 SDM006 1.0 3 608.477 480.233 519.785 60 SDM006 1.0 2 598.354 449.266 487.058 61 SDM006 1.0 2 617.823 456.908 507.467 62 SDM006 1.0 1 566.477 500.189 526.744 63 SDM006 1.0 1 622.170 462.463 550.675 64 SDM006 1.0 0 NA NA NA 65 SDM006 1.0 0 NA NA NA 66 SDM006 0.5 3 546.472 457.880 468.129 67 SDM006 0.5 3 525.069 444.575 505.154 68 SDM006 0.5 2 569.068 446.196 473.739 69 SDM006 0.5 2 534.205 470.366 476.570 bdf SampleID UVDose_J RepairHoursDay_0 Day_45 Day_90 17SDM001 0.5 B 88.6145 388.3575 198.467 36SDM002 0.5 B 100.0760 384.9505 234.740 55SDM005 0.5 B 121.9595 300.3650 241.832 74SDM006 0.5 B 174.7005 378.3435 291.272 93SDM007 0.5 B 319.7750 335.4390 110.306 112 SDM008 0.5 B 292.8400 323.0370
Re: [R] hdf-files
Hi Katrin, (don't forget to reply-all to R-help messages -- by default they are only sent to the recipient, and not back to the list). On Thu, Sep 23, 2010 at 10:03 AM, Katrin Fleischer katrin.fleisc...@falw.vu.nl wrote: The error message is Error in hdf5load(Jan2006.HDF) : unable to open HDF file: Jan2006.HDF Dumb question, but is that the correct path to the file? If you call `dir()` in your R session, is Jan2006.HDF one of the elements returned? Maybe try the absolute path to the file if not. I have also seen these other packages, but Im afraid that the data I have might be in hdf4 (thats would be an explanation to why hdf5load isnt working), so I was wondering whether there is a package for R dealing with hdf4 files.. Perhaps you could try converting it to hdf5 format, if that's the case. It looks like there is a utility you can use to do so here: http://www.hdfgroup.org/h4toh5/ -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Length of vector without NA's
Hi, this following code: x-c(1,2,NA) length(x) returns 3, correctly counting numbers as well as NA's. How can I exclude NA's from this count? Ralf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Prediction plot for logistic regression output
How do I construct a figure showing predicted value plots for the dependent variable as a function of each explanatory variable (separately) using the results of a logistic regression? It would also be helpful to know how to show uncertainty in the prediction (95% CI or SE). Thanks- This email has been processed by SmoothZap - www.smoothwall.net __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] For loop with ifelse help
Hi Pele, I think this should work file1$state.sum - rowSums(file2[file1$state,6:10],0) On Thu, Sep 23, 2010 at 7:46 PM, Steve Lianoglou mailinglist.honey...@gmail.com wrote: Hi Pele, On Wed, Sep 22, 2010 at 12:40 PM, Pele drdi...@yahoo.com wrote: Hi David - thanks for your suggestion, but I am trying to avoid doing any merging and sorting for this step because the real file I will be working with has about 20 million records. If I can get this loop or something similar to work will be good enough. If that's the case, you might consider looking at the sqldf or data.table packages. They both implement data.frame-like objects, but can do subsetting (and merging) rather quickly since they implement indexes over keys (columns) of the respective data.frame(s). Subsetting normal data.frames in this scenario you describe involves a linear search for every query through the column(s) you are querying against, which can get slow as the size of your data.frames get large. -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contacthttp://cbio.mskcc.org/%7Elianos/contact __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Length of vector without NA's
try this sum(!is.na(x)) I hope it helps. Best, Dimitris On 9/23/2010 5:08 PM, Ralf B wrote: Hi, this following code: x-c(1,2,NA) length(x) returns 3, correctly counting numbers as well as NA's. How can I exclude NA's from this count? Ralf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Dimitris Rizopoulos Assistant Professor Department of Biostatistics Erasmus University Medical Center Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands Tel: +31/(0)10/7043478 Fax: +31/(0)10/7043014 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Length of vector without NA's
this following code: x-c(1,2,NA) length(x) returns 3, correctly counting numbers as well as NA's. How can I exclude NA's from this count? sum(!is.na(x)) Peter __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Length of vector without NA's
Hi Ralf, The usual way (as others have shown you), takes advantage of the fact that the logical values TRUE and FALSE are counted as 1 and 0, respectively. is.na() returns TRUE if the value is NA, so to find how many are not NA, the result is reversed using ' ! '. Similar logic can be used to find how many meet any logical condition (e.g., sum(1:10 5) ). Cheers, Josh On Thu, Sep 23, 2010 at 8:08 AM, Ralf B ralf.bie...@gmail.com wrote: Hi, this following code: x-c(1,2,NA) length(x) returns 3, correctly counting numbers as well as NA's. How can I exclude NA's from this count? Ralf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.com/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Prediction plot for logistic regression output
Hi: Try ?termplot. HTH, Dennis On Thu, Sep 23, 2010 at 8:08 AM, Jay Vaughan j...@fishsciences.net wrote: How do I construct a figure showing predicted value plots for the dependent variable as a function of each explanatory variable (separately) using the results of a logistic regression? It would also be helpful to know how to show uncertainty in the prediction (95% CI or SE). Thanks- This email has been processed by SmoothZap - www.smoothwall.net __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Length of vector without NA's
-Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- project.org] On Behalf Of Joshua Wiley Sent: Thursday, September 23, 2010 8:23 AM To: Ralf B Cc: r-help Mailing List Subject: Re: [R] Length of vector without NA's Hi Ralf, The usual way (as others have shown you), takes advantage of the fact that the logical values TRUE and FALSE are counted as 1 and 0, respectively. is.na() returns TRUE if the value is NA, so to find how many are not NA, the result is reversed using ' ! '. Similar logic can be used to find how many meet any logical condition (e.g., sum(1:10 5) ). Cheers, Josh On Thu, Sep 23, 2010 at 8:08 AM, Ralf B ralf.bie...@gmail.com wrote: Hi, this following code: x-c(1,2,NA) length(x) returns 3, correctly counting numbers as well as NA's. How can I exclude NA's from this count? Ralf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. Ralf, While I might use the sum() function as others have posted, if you want the code to clearly show your intent (i.e. to get the length of a vector) then another option is length(x[!is.na(x)]) Hope this is helpful, Dan Daniel J. Nordlund Washington State Department of Social and Health Services Planning, Performance, and Accountability Research and Data Analysis Division Olympia, WA 98504-5204 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Passing a function as a parameter...
Are you aware that you can pass the function, itself, as an argument
(R is mostly a functional programming language where language objects
are first class objects).
e.g
g - function(x,fun)fun(x)
g(2,function(x)x^2)
[1] 4
On Wed, Sep 22, 2010 at 4:06 PM, Jonathan Greenberg
greenb...@ucdavis.edu wrote:
R-helpers:
If I want to pass a character name of a function TO a function, and then
have that function executed, how would I do this? I want
an arbitrary version of the following, where any function can be used (e.g.
I don't want the if-then statement here):
apply_some_function - function(data,function_name)
{
if(function_name==mean)
{
return(mean(data))
}
if(function_name==min)
{
return(min(data))
}
}
apply_some_function(1:10,mean)
apply_some_function(1:10,min)
Basically, I want the character name of the function used to actually
execute that function. Thanks!
--j
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Bert Gunter
Genentech Nonclinical Biostatistics
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Prediction plot for logistic regression output
On 23-Sep-10 15:08:06, Jay Vaughan wrote:
How do I construct a figure showing predicted value plots for the
dependent variable as a function of each explanatory variable
(separately) using the results of a logistic regression? It would also
be helpful to know how to show uncertainty in the prediction (95% CI or
SE).
Thanks-
The basic tool is the predict() function -- see '?predict.glm'
for details of its implementation for logistic regression etc.
Your desire to construct a figure showing predicted value plots
for the dependent variable as a function of each explanatory
variable (separately) will somehow have to cope with the fact
that the predicted value will be a function of *all* of the
independent variables, not just the one you are using for the
plot. How you decide to cope with this will depend on the purpose
for which you want the plot.
One approach to this could be to hold all the other independent
variables constant (say at their mean values), while including
a sequance of values for the variable if interest, in the 'newdata'
which you submit to the predict() function.
Another approach, which has a number of problematic aspects,
is, for each value of the independent variable in the plot,
to use the data to estimate a suitable representative value
of each other independent variable, and build your 'newdata'
accordingly. In particular, this introduces the complication
that there will be uncertainty in such representative values
which will have to be incorporated into the uncertainty of
prediction of the dependent variable.
Whichever approach you adopt, it is important with logistic
regression to note that, while you may predict the response
(i.e. plot the probability of Y=1), the use of the SE returned
by predict() when 'type=response' will in general give
nonsensical results (probabilities 0 or 1). The good
method is to use 'type=link' and transform using the
logistic response function exp(x)/(1 + exp(x)), which will
always give non-nonsensical results. The following example
illustrates this.
First, a logistic regression is fitted. Then Prob(Y=1) is
obtained from the '$fit' component of predict(), using
'type=response, and plotted. So far so good. Then 95%
confidence bounds around this are obtained as +/- 1.96*SE,
where SE is the $se component of predict(). It can be seen
that these give bad results near the ends of the range (red curves).
Then an SE is obtained jusing predict() with 'type=link,
and the +/-95% limits on the link fuinction are obtained.
These are then transformed using the logistic response function,
and plotted. These can be seen to be sensible (green curves).
lg - function(x){ exp(x)/(1 + exp(x)) }
set.seed(54321)
X - 0.3*(-15:15)
P - lg(X)
Y - 1*(runif(31) = P)
GLM - glm(Y ~ X, family=binomial)
prGLM - predict(GLM,type=response,se=TRUE)
plot(X,prGLM$fit,type=l,ylim=c(0,1),col=blue)
lines(X,prGLM$fit+1.96*prGLM$se,col=red)
lines(X,prGLM$fit-1.96*prGLM$se,col=red)
prGLM - predict(GLM,type=link,se=TRUE)
lines(X,lg(prGLM$fit+1.96*prGLM$se),col=green)
lines(X,lg(prGLM$fit-1.96*prGLM$se),col=green)
Hoping this helps,
Ted.
E-Mail: (Ted Harding) ted.hard...@manchester.ac.uk
Fax-to-email: +44 (0)870 094 0861
Date: 23-Sep-10 Time: 17:14:11
-- XFMail --
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to pass a model formula as argument to with.mids
??? Use lm's data argument to do this, not the with() construction: x - 1:10 y - 2*x+5 +rnorm(10) form - formula(y~x) class(form) [1] formula lm(form) Call: lm(formula = form) Coefficients: (Intercept)x 5.7471.921 dat - data.frame(x=x,y=y) rm(x,y) lm(form,data=dat) Call: lm(formula = form, data = dat) Coefficients: (Intercept)x 5.7471.921 On Thu, Sep 23, 2010 at 6:53 AM, Erich Studerus erich.stude...@bli.uzh.ch wrote: Hello I would like to pass a model formula as an argument to the with.mids function from the mice package. The with.mids functon fits models to multiply imputed data sets. Here's a simple example library(mice) #Create multiple imputations on the nhanes data contained in the mice package. imp - mice(nahnes) #Fitting a linear model with each imputed data set the regular way works fine with(imp, lm(bmi~hyp+chl)) #Creating a formula object and than passing it as argument does not work: form.obj - formula(bmi~hyp+chl) with(imp, lm(form.obj)) #The following doesn't work either expr -lm(bmi~hyp+chl) with(imp, expr) Looking at the definition of with.mids reveals that the second argument is first substituted and than evaluated within each data.frame of the multiply imputed data sets. Is there a way to pass lm(bmi~hyp+chl) without having to change the definition of the with.mids function? Thanks in advance, Erich __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Bert Gunter Genentech Nonclinical Biostatistics __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] empirical df, cdf and crossing points
Hello, I am using survey data with two stage sampling for two sub-populations. I am looking for a package (or packages) that can do the following for a measure of weight. The sub-populations are M (male) and F (female). (1) empirical df and cdf for weight, and compare that across two sub- populations. (i) get the x,y values where the plotted curves cross (more than one x in my data), (2) test for dominance, (3) give x value, and know which of the csf curves first or second order dominate at that value (using area under the curve). How can I do this? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] plotting multiple animal tracks against Date/Time
Dear list, I would like to create a time series plot in which the paths of several individuals are stacked above each other, with the x-axis being the total observation period of three years ( 1.1.2004 to 31.12.2007) and the y-axis being some defined range[min,max]. My data consist of Date/Time information and the paths of 45 individual as the distance from the location of release. An example data set for 2 individuals is given below.The observation period and frequency of observations varies between individuals. I believe stackplot() may be able to do this task, but I am not sure how to handle the variable time period and frequency of observations for different individuals. Could someone advise if stackplot() is suitable or if there is a better approach or package ? Thank you very much for your time and best wishes, Juliane Individual 1 DateDistance [m] 2006-08-18 22:05:15 1815.798 2006-08-18 22:06:35 1815.798 2006-08-18 22:08:33 1815.798 2006-08-18 22:09:49 1815.798 2006-08-18 22:12:50 1815.798 2006-08-18 22:16:26 1815.798 Individual 2 Date Distance [m] 2006-08-18 09:53:20 0.0 2006-08-18 09:59:07 0.0 2006-08-18 10:09:20 0.0 2006-08-18 10:21:14 0.0 2006-08-18 10:34:18 0.0 2006-08-18 10:36:44100.2 2 Date Distance 6 2006-08-18 09:53:20 0.0 7 2006-08-18 09:59:07 0.0 8 2006-08-18 10:09:20 0.0 9 2006-08-18 10:21:14 0.0 10 2006-08-18 10:34:18 0.0 11 2006-08-18 10:36:44100.2 006-03-1 22:05:15 1815.798 2006-03-18 22:06:35 1815.798 2006-03-18 22:08:33 1815.798 2006-03-18 22:09:49 1815.798 2006-03-18 22:12:50 1815.798 2006-03-18 22:16:26 1815.798 Dr. Juliane Struve Imperial College London Department of Life Sciences Silwood Park Campus Buckhurst Road Ascot, Berkshire, SL5 7PY, UK Tel: +44 (0)20 7594 2527 Fax: +44 (0)1344 874 957 http://www.aquaticresources.org __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Plotting multiple animal tracks against Date/Time
Sorry for posting this questions twice, but my previous question was accidentally sent unfinished. Dear list, I would like to create a time series plot in which the paths of several individuals are stacked above each other, with the x-axis being the total observation period of three years ( 1.1.2004 to 31.12.2007) and the y-axis being some defined range[min,max]. My data consist of Date/Time information and the paths of 45 individual as the distance from the location of release. An example data set for 2 individuals is given below.The observation period and frequency of observations varies between individuals. I believe stackplot() may be able to do this task, but I am not sure how to handle the variable time period and frequency of observations for different individuals. Could someone advise if stackplot() is suitable or if there is a better approach or package ? Thank you very much for your time, Juliane Individual 1 DateDistance [m] 2005-07-18 22:05:15 1815.798 2005-07-18 22:06:35 1815.798 2005-07-18 22:08:33 1815.798 2005-07-18 22:09:49 1815.798 2005-07-18 22:12:50 1815.798 2005-07-18 22:16:26 1815.798 Individual 2 Date Distance [m] 2006-08-18 09:53:20 0.0 2006-08-18 09:59:07 0.0 2006-08-18 10:09:20 0.0 2006-08-18 10:21:14 100.5 Dr. Juliane Struve Imperial College London Department of Life Sciences Ascot, Berkshire, SL5 7PY, UK __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] merging multiple data frames
hi guys..thx for help i have to perform a calculation P-B/T-B where P is the values in pdf and B is the values in bdf and T in tdf On Thu, Sep 23, 2010 at 7:49 PM, Mike Rennie mikerenni...@gmail.com wrote: First, you might want to start by generating a new column to identify your 'pdf and bdf or whatever once it's merged. For the merging, see ?merge But as someone's already pointed out, it's not clear what you are trying to merge by. Also, as your example calculations show, you don't need to merge it to to the calcuations you want to do... On Thu, Sep 23, 2010 at 8:53 AM, Steve Lianoglou mailinglist.honey...@gmail.com wrote: Hi, On Thu, Sep 23, 2010 at 9:04 AM, rasanpreet rasanpreet.k...@gmail.com wrote: hi guys i have multiple data frames which i want to merge. there are four of them..eg Can you provide a (correct) example of what you want your merged data.frame to look like? What column do you want to use in your data.frame to merge against? I'm guessing SampleID(?), but then again, these aren't unique in your `pdf` data.frame. For instance, what would the row for SDM001 look like in your merged data.frame? -steve pdf SampleID UVDose_J RepairHours Day_0 Day_45 Day_90 1SDM001 1.0 3 485.612 465.142 490.873 2SDM001 1.0 3 503.658 457.863 487.783 3SDM001 1.0 2 533.193 451.044 456.973 4SDM001 1.0 2 538.334 452.887 474.915 5SDM001 1.0 1 526.034 481.123 477.801 6SDM001 1.0 1 546.543 472.322 481.546 7SDM001 1.0 0 NA NA NA 8SDM001 1.0 0 NA NA NA 9SDM001 0.5 3 432.134 457.245 497.975 10 SDM001 0.5 3 432.605 450.184 489.468 11 SDM001 0.5 2 450.335 496.520 488.784 12 SDM001 0.5 2 439.590 474.371 470.182 13 SDM001 0.5 1 510.480 489.561 525.029 14 SDM001 0.5 1 487.934 467.258 488.784 15 SDM001 0.5 0 NA NA NA 16 SDM001 0.5 0 NA NA NA 20 SDM002 1.0 3 465.549 528.715 501.374 21 SDM002 1.0 3 458.168 505.480 489.244 22 SDM002 1.0 2 447.317 464.009 478.058 23 SDM002 1.0 2 452.020 438.446 470.996 24 SDM002 1.0 1 441.718 458.760 499.221 25 SDM002 1.0 1 447.017 402.616 548.797 26 SDM002 1.0 0 NA NA NA 27 SDM002 1.0 0 NA NA NA 28 SDM002 0.5 3 421.409 448.870 476.392 29 SDM002 0.5 3 404.089 446.413 477.080 30 SDM002 0.5 2 399.775 432.678 465.015 31 SDM002 0.5 2 427.157 443.418 477.048 32 SDM002 0.5 1 389.674 449.353 482.264 33 SDM002 0.5 1 418.147 455.983 495.486 34 SDM002 0.5 0 NA NA NA 35 SDM002 0.5 0 NA NA NA 39 SDM005 1.0 3 579.836 441.040 476.382 40 SDM005 1.0 3 578.525 443.875 472.867 41 SDM005 1.0 2 564.266 432.116 469.416 42 SDM005 1.0 2 571.045 447.658 458.233 43 SDM005 1.0 1 564.664 427.673 524.122 44 SDM005 1.0 1 568.182 458.039 477.237 45 SDM005 1.0 0 NA NA NA 46 SDM005 1.0 0 NA NA NA 47 SDM005 0.5 3 556.534 424.786 501.658 48 SDM005 0.5 3 474.027 441.418 507.635 49 SDM005 0.5 2 481.355 430.346 468.021 50 SDM005 0.5 2 478.922 466.933 471.025 51 SDM005 0.5 1 505.539 937.759 460.985 52 SDM005 0.5 1 497.913 457.932 493.152 53 SDM005 0.5 0 NA NA NA 54 SDM005 0.5 0 NA NA NA 58 SDM006 1.0 3 589.164 459.578 509.565 59 SDM006 1.0 3 608.477 480.233 519.785 60 SDM006 1.0 2 598.354 449.266 487.058 61 SDM006 1.0 2 617.823 456.908 507.467 62 SDM006 1.0 1 566.477 500.189 526.744 63 SDM006 1.0 1 622.170 462.463 550.675 64 SDM006 1.0 0 NA NA NA 65 SDM006 1.0 0 NA NA NA 66 SDM006 0.5 3 546.472 457.880 468.129 67 SDM006 0.5 3 525.069 444.575 505.154 68 SDM006 0.5 2 569.068 446.196 473.739 69 SDM006 0.5 2 534.205 470.366 476.570 bdf SampleID UVDose_J RepairHoursDay_0 Day_45 Day_90 17SDM001 0.5 B 88.6145 388.3575 198.467 36SDM002 0.5 B 100.0760 384.9505 234.740 55SDM005 0.5
Re: [R] Determine area between two density plots
I wonder what the best way is to access those values. I am using the following code: x1 - c(1,2,1,3,5,6,6,7,7,8) x2 - c(1,2,1,3,5,6,5,3,8,7) d1 - density(x1, na.rm = TRUE) d2 - density(x2, na.rm = TRUE) plot(d1, lwd=3, main=bla) lines(d2, lty=2, lwd=3) d1[1] d1[2] The last two lines allow me to access 1000 values, but I don't know if this is the right approach. I also don't know why they are in two columns. Does density have a saver way to get to those values? Ralf Ralf On Wed, Sep 22, 2010 at 5:25 PM, Peter Alspach peter.alsp...@plantandfood.co.nz wrote: Tena koe Ralf If you save the results of density() x1Den - density(x1) you get the x and y values of the line which is plotted. Similarly for x2 - you can then use these to shade the joint area and find the area. Tinkering with the arguments of density to make the x values for each the same will make this process easier. Let me know if you'd like more details. HTH Peter Alspach -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- project.org] On Behalf Of Ralf B Sent: Thursday, 23 September 2010 8:55 a.m. To: r-help Mailing List Subject: [R] Determine area between two density plots Hi group, I am creating two density plots as shown in the code below: x1 - c(1,4,5,3,2,3,4,5,6,5,4,3,2,1,1,1,2,3) x2 - c(1,4,5,3,5,7,4,5,6,1,1,1,2,1,1,1,2,3) plot(density(x1, na.rm = TRUE)) polygon(density(x2, na.rm = TRUE), border=blue) How can I determine the area that is covered between the two plots as a number and how can I grey (or highlight with a pattern) the area that lies between the two lines? Thanks, Ralf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. The contents of this e-mail are confidential and may be subject to legal privilege. If you are not the intended recipient you must not use, disseminate, distribute or reproduce all or any part of this e-mail or attachments. If you have received this e-mail in error, please notify the sender and delete all material pertaining to this e-mail. Any opinion or views expressed in this e-mail are those of the individual sender and may not represent those of The New Zealand Institute for Plant and Food Research Limited. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] help in density estimation
Hi, guys, I'm using kernel density estimation. But how can I return to a density estimation for a fixed point? For example, b-runif(1000,0,1) f-density(b) How can I get the value of density(b) at b=0.5? Your help is extremely appreciated. Thanks. Jay -- View this message in context: http://r.789695.n4.nabble.com/help-in-density-estimation-tp2552264p2552264.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] plotting multiple animal tracks against Date/Time
By range on the y-axis, do you mean distance? It would have to be if time is on your x? Or am I misreading this? You could just plot() with the data for your first individual, and then add additional individuals after that using lines(), specifying a different colour and/or line type for each individual, and could even plot points as well to identify your actual data collection points. Alternatively, you could do a bunch of multi-panel plots with the same axes, and individual data is given in each plot. It sounds like this is what you want to do by referencing stackplot() (based on it's description, I've never used it)? If so, see ?par for details, specifically mfrow, mfcol. By doing it with calls to par(), I think you'll have more control over the appearance of the plot than with stackplot(). Mike On Thu, Sep 23, 2010 at 8:50 AM, Struve, Juliane j.str...@imperial.ac.ukwrote: Dear list, I would like to create a time series plot in which the paths of several individuals are stacked above each other, with the x-axis being the total observation period of three years ( 1.1.2004 to 31.12.2007) and the y-axis being some defined range[min,max]. My data consist of Date/Time information and the paths of 45 individual as the distance from the location of release. An example data set for 2 individuals is given below.The observation period and frequency of observations varies between individuals. I believe stackplot() may be able to do this task, but I am not sure how to handle the variable time period and frequency of observations for different individuals. Could someone advise if stackplot() is suitable or if there is a better approach or package ? Thank you very much for your time and best wishes, Juliane Individual 1 DateDistance [m] 2006-08-18 22:05:15 1815.798 2006-08-18 22:06:35 1815.798 2006-08-18 22:08:33 1815.798 2006-08-18 22:09:49 1815.798 2006-08-18 22:12:50 1815.798 2006-08-18 22:16:26 1815.798 Individual 2 Date Distance [m] 2006-08-18 09:53:20 0.0 2006-08-18 09:59:07 0.0 2006-08-18 10:09:20 0.0 2006-08-18 10:21:14 0.0 2006-08-18 10:34:18 0.0 2006-08-18 10:36:44100.2 2 Date Distance 6 2006-08-18 09:53:20 0.0 7 2006-08-18 09:59:07 0.0 8 2006-08-18 10:09:20 0.0 9 2006-08-18 10:21:14 0.0 10 2006-08-18 10:34:18 0.0 11 2006-08-18 10:36:44100.2 006-03-1 22:05:15 1815.798 2006-03-18 22:06:35 1815.798 2006-03-18 22:08:33 1815.798 2006-03-18 22:09:49 1815.798 2006-03-18 22:12:50 1815.798 2006-03-18 22:16:26 1815.798 Dr. Juliane Struve Imperial College London Department of Life Sciences Silwood Park Campus Buckhurst Road Ascot, Berkshire, SL5 7PY, UK Tel: +44 (0)20 7594 2527 Fax: +44 (0)1344 874 957 http://www.aquaticresources.org __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.htmlhttp://www.r-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help in density estimation
On 23/09/2010 11:42 AM, wangguojie2006 wrote: b-runif(1000,0,1) f-density(b) f is a list of things, including x values where the density is computed, and y values for the density there. So you could do it by linear interpolation using approx or approxfun. For example b - runif(1000,0,1) flist - density(b) f - approxfun(flist$x, flist$y) f(0.2) [1] 0.9717893 f(-1) [1] NA If you don't like the NA for an out-of-range argument, then choose something different from the default for the rule argument to approxfun. Duncan Murdoch __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] looking for a faster way to compare two columns of a matrix
Please consider this matrix:
x - structure(c(5, 4, 3, 2, 1, 6, 3, 2, 1, 0, 3, 2, 1, 0, 0, 2, 1,
1, 0, 0, 2, 0, 0, 0, 0), .Dim = c(5L, 5L))
For each pair of columns, I want to calculate the proportion of entries
different than 0 in column j (i j) that have lower values than the entries
in the same row in column i:
x[, 1:2]
sum((x[,1] x[,2]) (x[,2] 0))/sum(x[,2] 0)
Thus, for this pair, 3 of the 4 entries in the second column are
lower than the entries in the same row in the first column.
When both columns of a given pair have the same number of cells different than
0, the value of the metric is 0.
x[, 3:4]
colSums(x[, 3:4] 0)
The same if column j has more valid ( 0) entries.
I've been doing this using this idea:
combinations - combn(1:ncol(x), 2)
values - numeric(ncol(combinations))
for (i in 1:ncol(combinations)) {
pair - combinations[,i]
first - x[, pair[1]]
second - x[, pair[2]]
if (sum(first 0) = sum(second 0)) next
values[i] - sum(first - second 0 second 0) / sum(second 0)
}
values
Anyway, I was wondering if there is a faster/better way. I've tried
putting the code from
the for loop into a function and passing it to combn but, as expected, it didn't
help much. Any pointers to functions that I should be looking into will be
greatly appreciated.
Thank you very much,
Gustavo.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] For loop with ifelse help
Hi Sayan, This is exactly what I was looking for - it worked perfectly. Many thanks!! Also, thanks to everyone else for their suggestions. Pele -- View this message in context: http://r.789695.n4.nabble.com/For-loop-with-ifelse-help-tp2550547p2552388.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Determine area between two density plots
Am 23.Sep.2010 um 18:27 schrieb Ralf B: I wonder what the best way is to access those values. I am using the following code: x1 - c(1,2,1,3,5,6,6,7,7,8) x2 - c(1,2,1,3,5,6,5,3,8,7) d1 - density(x1, na.rm = TRUE) d2 - density(x2, na.rm = TRUE) plot(d1, lwd=3, main=bla) lines(d2, lty=2, lwd=3) d1[1] d1[2] The last two lines allow me to access 1000 values, but I don't know if this is the right approach. I also don't know why they are in two columns. Does density have a saver way to get to those values? Do you mean d1$x, d1$y Benno __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] extracting random effects from model formula
I've found using the arm package to be very useful. require(arm) then use ranef(Full_model) and fixef(Full_model) On Wed, Sep 22, 2010 at 05:39:41PM +1000, Lorenzo Cattarino wrote: Hi again, Sorry, probably I was not clear enough. I was wondering if there is a way in R to identify (and extract) only the random effects, which, because I am using the lmer function, are the terms with the symbol | on the left side of the grouping variable (SITE in my example). Thanks Lorenzo From: Lorenzo Cattarino Sent: Wednesday, 22 September 2010 5:23 PM To: r-help@r-project.org Subject: extracting random effects from model formula Hi R-users I would like to extract the random effects (1|SITE, 1+SPECIES|SITE and BA|SITE) from this model formula: Full_model - formula (VAR ~ (1|SITE) + (1+SPECIES|SITE) + (BA|SITE) + HEIGHT + COND + NN_DIST) I tried: terms(Full_model) labels(terms(Full_model)) but I could not distinguish between random and fixed effects. thanks Lorenzo [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] NetCDF file: adding a variable
Hi Gennadi, I think this is fixed in the latest version of the ncdf package, which you can access by going to CRAN, then 'packages', then 'ncdf', and clicking on 'URL'. Or here is a direct link: http://cirrus.ucsd.edu/~pierce/ncdf Give that a try and let me know if you still have a problem, --Dave Gennadi Lessin wrote: I am trying to create a NetCDF file with bathymetry (a matrix) and then to add a grid type (an integer) as variables. This is the relevant part of the code: library(ncdf) wfile=my_file.nc msv=-10 grdtp=2 # Dimensions d1=dim.def.ncdf(lon,degrees,as.double(lon)) d2=dim.def.ncdf(lat,degrees,as.double(lat)) # Variables bathymetry=var.def.ncdf(bathymetry,m,list(d1,d2),msv,longname=Bathymetry) ncw=create.ncdf(writefile,list(bathymetry,loncp,latcp,dlonp,dlatp)) put.var.ncdf(ncw,bathymetry,bat_matrix) close.ncdf(ncw) # Here I am trying to add another variable which should have a value of grdtp ncw_old=open.ncdf(wfile,write=TRUE) d3=dim.def.ncdf(grid_type,'',1:1,create_dimvar=FALSE) grid_t=var.def.ncdf(grid_type,type,d3,1.e30,longname=Grid Type) ncw_new=var.add.ncdf(ncw_old,grid_t) Here I stop because R gives an error message: Error in var.add.ncdf(ncw_old, grid_t) : Error in create.ncdf, defining var! Any ideas of what am I doing wrong? Thank you! [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. --- David W. Pierce Division of Climate, Atmospheric Science, and Physical Oceanography Scripps Institution of Oceanography (858) 534-8276 (voice) / (858) 534-8561 (fax)dpie...@ucsd.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help in density estimation
On 23-Sep-10 16:52:09, Duncan Murdoch wrote: On 23/09/2010 11:42 AM, wangguojie2006 wrote: b-runif(1000,0,1) f-density(b) f is a list of things, including x values where the density is computed, and y values for the density there. So you could do it by linear interpolation using approx or approxfun. For example b - runif(1000,0,1) flist - density(b) f - approxfun(flist$x, flist$y) f(0.2) [1] 0.9717893 f(-1) [1] NA If you don't like the NA for an out-of-range argument, then choose something different from the default for the rule argument to approxfun. Duncan Murdoch Or, perhaps more transparently (and more explicitly modifiable): b-runif(1000,0,1) f - density(b, from=0, to=1, n=512) plot(f$x, f$y, type=l, col=blue, xlim=c(0,1), ylim=c(0,1.5)) ## Plot the density estimate x0 - 0.5## Target value of x i0 - max(which(f$x = x0)) i1 - min(which(f$x x0)) u0 - f$x[i0] ; v0 - f$y[i1] u1 - f$x[i1] ; v1 - f$y[i1] y0 - v0 + (v1-v0)*(x0-u0)/(u1-u0) ## Linear interpolation points(x0, y0, pch=+, col=red) ## Add interpolated point Ted. E-Mail: (Ted Harding) ted.hard...@manchester.ac.uk Fax-to-email: +44 (0)870 094 0861 Date: 23-Sep-10 Time: 18:05:41 -- XFMail -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] merging multiple data frames
On Thu, Sep 23, 2010 at 12:18 PM, rasanpreet kaur suri rasanpreet.k...@gmail.com wrote: hi guys..thx for help i have to perform a calculation P-B/T-B where P is the values in pdf and B is the values in bdf and T in tdf You have 69 rows in your pdf data.frame, and something like 20 rows in bdf and tdf, so my original question stands: What column do you want to use in your data.frame to merge against? I'm guessing SampleID(?), but then again, these aren't unique in your `pdf` data.frame. For instance, what would the row for SDM001 look like in your merged data.frame? You know what I mean? How do you intend to group the rows across your pdf,bdf,tdf to do this calculation? -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] plotting multiple animal tracks against Date/Time
On Thu, Sep 23, 2010 at 9:50 AM, Struve, Juliane j.str...@imperial.ac.uk wrote: Dear list, I would like to create a time series plot in which the paths of several individuals are stacked above each other, with the x-axis being the total observation period of three years ( 1.1.2004 to 31.12.2007) and the y-axis being some defined range[min,max]. My data consist of Date/Time information and the paths of 45 individual as the distance from the location of release. An example data set for 2 individuals is given below.The observation period and frequency of observations varies between individuals. I believe stackplot() may be able to do this task, but I am not sure how to handle the variable time period and frequency of observations for different individuals. Could someone advise if stackplot() is suitable or if there is a better approach or package ? Thank you very much for your time and best wishes, Juliane Try this: Lines1 - DateDistance [m] 2006-08-18 22:05:15 1815.798 2006-08-18 22:06:35 1815.798 2006-08-18 22:08:33 1815.798 2006-08-18 22:09:49 1815.798 2006-08-18 22:12:50 1815.798 2006-08-18 22:16:26 1815.798 Lines2 - Date Distance [m] 2006-08-18 09:53:20 0.0 2006-08-18 09:59:07 0.0 2006-08-18 10:09:20 0.0 2006-08-18 10:21:14 0.0 2006-08-18 10:34:18 0.0 2006-08-18 10:36:44100.2 library(chron) dt - function(date, time) as.chron(paste(date, time)) library(zoo) library(chron) # read in data dt - function(date, time) as.chron(paste(date, time)) z1 - read.zoo(textConnection(Lines1), skip = 1, index = list(1, 2), FUN = dt) z2 - read.zoo(textConnection(Lines2), skip = 1, index = list(1, 2), FUN = dt) # create single zoo object z - na.approx(cbind(z1, z2), na.rm = FALSE) # plot -- remove screen=1 if you want separate panels plot(z, screen = 1) -- Statistics Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help in density estimation
You could do: b - runif(1, 0, 1) tmp - density(b, from=0.5, to=0.5, n=1) tmp$y As one direct approach. You could also look at the logspline package for an alternative for density estimation that provides you with a density function (and also allows for finite domains). -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- project.org] On Behalf Of wangguojie2006 Sent: Thursday, September 23, 2010 9:43 AM To: r-help@r-project.org Subject: [R] help in density estimation Hi, guys, I'm using kernel density estimation. But how can I return to a density estimation for a fixed point? For example, b-runif(1000,0,1) f-density(b) How can I get the value of density(b) at b=0.5? Your help is extremely appreciated. Thanks. Jay -- View this message in context: http://r.789695.n4.nabble.com/help-in- density-estimation-tp2552264p2552264.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] scatterplot 3d equal axis sequence length limitation
I was wondering if anyone has a way out of the limitation that you must use equal length x,y, and z sequence lengths. For instance, x-seq(1,100) y-seq(1,100) z-rnorm(100) scatterplot3d(z,x,y) works fine. However, if I get some results that has a different y subset length, such as x-seq(1,100) y-seq(1,300) z-rnorm(100) scatterplot3d(z,x,y) I get the following error: Error in xyz.coords(x = x, y = y, z = z, xlab = xlabel, ylab = ylabel, : 'x', 'y' and 'z' lengths differ I have found one solution is to pad the values with n*0, where n is the number of excess variables of y over x and z. Unfortunately, the visual appearance is not that appealing. The situation is very practical as there are cases where you might limit the x axis variable length to some value, and have many more runs of y (monte carlo sweeps for instance). Ideally, rather than pad, If I cannot limit the length of one axis to the same length of another, I would just like the color to be transparent for those values (edges and vertices). thanks, rtist -- View this message in context: http://r.789695.n4.nabble.com/scatterplot-3d-equal-axis-sequence-length-limitation-tp2552476p2552476.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Contour Plot - water bodies
Hi All I have been trying to plot irregular XYZ data on world / country map, where X and Y are llongitude and latitude coordinates and Z is a variable. My problem is I get contour lines over water bodies. The method I follow is: library(akima) library(maps) library(mapdata) library(gpclib) library(maptools) ori - read.table (interp2.txt, header=T) ori.li - interp(ori$x, ori$y, ori$z) ind- readShapePoly(IND_adm0.shp) plot(ori$y ~ ori$x, main = ori data) plot(ind, add=TRUE)contour(ori.li, col= topo.colors (1), lwd=2, add=TRUE) Can someone please help me to draw the contour plot without the contour lines falling over water bodies like sea. I thank in advance for the time and help TrulyG. Arun Kumar Junior Research Fellow Dept of Immunology School of Biological Sciences Madurai Kamaraj University Madurai - 625 021 An eye for an eye will only make the entire world blind.-Mohandas Karamchand Gandhi [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help in density estimation
You guys are really good. Thanks. -- View this message in context: http://r.789695.n4.nabble.com/help-in-density-estimation-tp2552264p2552484.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] scatterplot 3d equal axis sequence length limitation
On 23/09/2010 2:10 PM, rtist wrote: I was wondering if anyone has a way out of the limitation that you must use equal length x,y, and z sequence lengths. For instance, x-seq(1,100) y-seq(1,100) z-rnorm(100) scatterplot3d(z,x,y) works fine. However, if I get some results that has a different y subset length, such as x-seq(1,100) y-seq(1,300) z-rnorm(100) scatterplot3d(z,x,y) I get the following error: Error in xyz.coords(x = x, y = y, z = z, xlab = xlabel, ylab = ylabel, : 'x', 'y' and 'z' lengths differ I have found one solution is to pad the values with n*0, where n is the number of excess variables of y over x and z. Unfortunately, the visual appearance is not that appealing. The situation is very practical as there are cases where you might limit the x axis variable length to some value, and have many more runs of y (monte carlo sweeps for instance). I don't understand what you would want to plot: scatterplot3d is for plotting triplets (x,y,z). If you have more y's than x's and z's, how do you form the triplets? Duncan Murdoch Ideally, rather than pad, If I cannot limit the length of one axis to the same length of another, I would just like the color to be transparent for those values (edges and vertices). thanks, rtist __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Ordinal mixed model
Patrick-- One other option in addition to Thierry's suggestion, within R you might also consider the ordinal package, which handles random-intercept models. That said, if you are used to SPSS, I suspect this will be a titanic pain trying to move to R (part. if just for this one analysis...). Having been taught SPSS in grad school (and still working amongst many SPSS users), the two programs are very different. SPSS will fit GEE models, which I would guess is a defensible approach for your data (given the very little that I know about it...). You might explore that before scaling the steep learning curve of R. Hope that helps. cheers, Dave Dear Patrick, You can fit such a model with the MCMCglmm package. Have a look at the vignette of that package. install.packages(MCMCglmm) vignette(CourseNotes, package = MCMCglmm) But I'm affraid that this will require some rockclimbing upon the learning curve if you are a R novice. HTH, Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek team Biometrie Kwaliteitszorg Gaverstraat 4 9500 Geraardsbergen Belgium Research Institute for Nature and Forest team Biometrics Quality Assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 Thierry.Onkelinx op inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey -Oorspronkelijk bericht- Van: r-help-bounces op r-project.org [mailto:r-help-bounces op r-project.org] Namens Patrick Walsh Verzonden: woensdag 22 september 2010 18:29 Aan: r-help op r-project.org Onderwerp: [R] Ordinal mixed model Hello, I am trying to build a generalised linear mixed model. My dependent variable is ordinal. I have a random factor (7 individuals), and a repeated measure where the dependent variable was measured three times for each of four attempts (so the repeats are nested). I also have a few covariates. I am a complete novice in R, being used to using SPSS. SPSS lets me build an ordinal model with repeated measures, but can't include a random factor. So that is really the hurdle, is this possible in R. And is there a way to do this that could be explained to someone who has no experience in R? Any help would be greatly appreciated. Kind Regards, ptwal __ R-help op r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Dave Atkins, PhD Research Associate Professor Department of Psychiatry and Behavioral Science University of Washington datk...@u.washington.edu Center for the Study of Health and Risk Behaviors (CSHRB) 1100 NE 45th Street, Suite 300 Seattle, WA 98105 206-616-3879 http://depts.washington.edu/cshrb/ (Mon-Wed) Center for Healthcare Improvement, for Addictions, Mental Illness, Medically Vulnerable Populations (CHAMMP) 325 9th Avenue, 2HH-15 Box 359911 Seattle, WA 98104 http://www.chammp.org (Thurs) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help in density estimation
There was a typo error in my code below. See the inserted correction. On 23-Sep-10 17:05:45, Ted Harding wrote: On 23-Sep-10 16:52:09, Duncan Murdoch wrote: On 23/09/2010 11:42 AM, wangguojie2006 wrote: b-runif(1000,0,1) f-density(b) f is a list of things, including x values where the density is computed, and y values for the density there. So you could do it by linear interpolation using approx or approxfun. For example b - runif(1000,0,1) flist - density(b) f - approxfun(flist$x, flist$y) f(0.2) [1] 0.9717893 f(-1) [1] NA If you don't like the NA for an out-of-range argument, then choose something different from the default for the rule argument to approxfun. Duncan Murdoch Or, perhaps more transparently (and more explicitly modifiable): b-runif(1000,0,1) f - density(b, from=0, to=1, n=512) plot(f$x, f$y, type=l, col=blue, xlim=c(0,1), ylim=c(0,1.5)) ## Plot the density estimate x0 - 0.5## Target value of x i0 - max(which(f$x = x0)) i1 - min(which(f$x x0)) ##u0 - f$x[i0] ; v0 - f$y[i1]## WRONG! u0 - f$x[i0] ; v0 - f$y[i0]## Correction u1 - f$x[i1] ; v1 - f$y[i1] y0 - v0 + (v1-v0)*(x0-u0)/(u1-u0) ## Linear interpolation points(x0, y0, pch=+, col=red) ## Add interpolated point Ted. E-Mail: (Ted Harding) ted.hard...@manchester.ac.uk Fax-to-email: +44 (0)870 094 0861 Date: 23-Sep-10 Time: 18:05:41 -- XFMail -- E-Mail: (Ted Harding) ted.hard...@manchester.ac.uk Fax-to-email: +44 (0)870 094 0861 Date: 23-Sep-10 Time: 20:38:53 -- XFMail -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
