Re: [R] Superscript and subscrib R for legend x-axis and y-axis and colour different subjects in longitudinal data with different colours
Dear all and Pikal (in particular :)) Pikal, I’m sorry. It works!!! Thank you very much :) Best, Rosa Oliveira > On 1 Aug 2017, at 16:31, PIKAL Petr <petr.pi...@precheza.cz> wrote: > > Hi > > Keep your messages coppied to R helplist, others could give you answers too. > > See in line > > From: Rosa Oliveira [mailto:rosit...@gmail.com] > Sent: Tuesday, August 1, 2017 4:38 PM > To: PIKAL Petr <petr.pi...@precheza.cz> > Subject: Re: [R] Superscript and subscrib R for legend x-axis and y-axis and > colour different subjects in longitudinal data with different colours > > Hi Pikal, > > I looked your message, but I think you do not answer the question, right? > > > Why do you think so. did you try it? > > I basically suggested > > plot(c(1:5), type = "n", > xlab=expression("t"[i]), > ylab=expression("y"^ij)) > mtext(expression(lambda^2)) > > which by my humble opinion is precisely what you wanted. > > I want to write t_i instead of "Day in ICU” [i subscript for t] and y_ij > instead of "CRP (mg/dL)” [ij superscript for y]. > > If it is not you need to express yourself more clearly. My crystal ball is > broken. > > Cheers > Petr > > On 31 Jul 2017, at 13:44, PIKAL Petr <petr.pi...@precheza.cz > <mailto:petr.pi...@precheza.cz>> wrote: > > Hi > > > From: Rosa Oliveira [mailto:rosit...@gmail.com <mailto:rosit...@gmail.com>] > Sent: Monday, July 31, 2017 11:47 AM > To: Martin Maechler <maech...@stat.math.ethz.ch > <mailto:maech...@stat.math.ethz.ch>> > Cc: PIKAL Petr <petr.pi...@precheza.cz <mailto:petr.pi...@precheza.cz>>; > r-help mailing list <r-help@r-project.org <mailto:r-help@r-project.org>> > Subject: Re: [R] Superscript and subscrib R for legend x-axis and y-axis and > colour different subjects in longitudinal data with different colours > > Hi, everyone, > > Before everything, thanks. Lots of thanks ;) > > I don’t think you understood everything I need to do. > > Most probably because you did not tell precisely what you really want. > > I want to write t_i instead of "Day in ICU” [i subscript for t] and y_ij > instead of "CRP (mg/dL)” [ij superscript for y]. The label of the axis… :( > > So something like that. > plot(c(1:5), CRP7raw[1,], type = "n", xlim=c(1,5), > ylim=c(-10,5) , > xlab=expression("t"[i]), > ylab=expression("y"^ij)) > mtext(expression(lambda^2)) > > Cheers > Petr > > Can you help me on that task? > > Thanks! > > Best, > Rosa Oliveira > > > On 31 Jul 2017, at 10:28, Martin Maechler <maech...@stat.math.ethz.ch > <mailto:maech...@stat.math.ethz.ch>> wrote: > > PIKAL Petr <petr.pi...@precheza.cz <mailto:petr.pi...@precheza.cz>> >on Mon, 31 Jul 2017 09:11:18 + writes: > > Hi Martin see in line > > -Original Message- From: Martin Maechler > [mailto:maech...@stat.math.ethz.ch <mailto:maech...@stat.math.ethz.ch>] Sent: > Monday, July > 31, 2017 10:52 AM To: PIKAL Petr <petr.pi...@precheza.cz > <mailto:petr.pi...@precheza.cz>> > Cc: Rosa Oliveira <rosit...@gmail.com <mailto:rosit...@gmail.com>>; r-help > mailing > list <r-help@r- >project.org <http://project.org/>> > > Subject: Re: [R] Superscript and subscrib R for legend > x-axis and y-axis and colour different subjects in > longitudinal data with different colours > > > > Hi Rosa > something like > > plot(1,1, sub=expression(lambda^"2")) > > So with your example, do you want something like > > plot(c(1:5), CRP7raw[1,], type = "n", xlim=c(1,5), > ylim=c(-10,5) , > xlab="Day in ICU", > ylab="CRP > (mg/dL)", > sub = mtext(expression(lambda^2))) > > OOps! Either plot( ..., sub = *) or plot( ... ) ; > mtext(*) > > but not both! > > You are right, I used a code from OP and did not much > think about it. Strangely enough, the code worked without > any complain. Probably mtext is "stronger" than sub and > overrides it. > > Well, well, "the magic of R" . > Not quite: mtext(..) is a valid function call that is evaluated > before plot() finishes; then it returns NULL invisibly, and > 'plot(*, sub=NULL)' does nothing "coincidentally" (but for good > reasons) > > Martin > > > Cheers Petr > > > CRP7graph <- apply(CRP7, 1, lines, col="gray") > > Cheers > Petr > > > > -Orig
Re: [R] Superscript and subscrib R for legend x-axis and y-axis and colour different subjects in longitudinal data with different colours
Hi, everyone, Before everything, thanks. Lots of thanks ;) I don’t think you understood everything I need to do. I want to write t_i instead of "Day in ICU” [i subscript for t] and y_ij instead of "CRP (mg/dL)” [ij superscript for y]. The label of the axis… :( Can you help me on that task? Thanks! Best, Rosa Oliveira > On 31 Jul 2017, at 10:28, Martin Maechler <maech...@stat.math.ethz.ch> wrote: > >>>>>> PIKAL Petr <petr.pi...@precheza.cz <mailto:petr.pi...@precheza.cz>> >>>>>>on Mon, 31 Jul 2017 09:11:18 + writes: > >> Hi Martin see in line > >>> -Original Message- From: Martin Maechler >>> [mailto:maech...@stat.math.ethz.ch] Sent: Monday, July >>> 31, 2017 10:52 AM To: PIKAL Petr <petr.pi...@precheza.cz> >>> Cc: Rosa Oliveira <rosit...@gmail.com>; r-help mailing >>> list <r-help@r- >project.org> >>> Subject: Re: [R] Superscript and subscrib R for legend >>> x-axis and y-axis and colour different subjects in >>> longitudinal data with different colours >>> >>> >>>> Hi Rosa > something like >>> >>>> plot(1,1, sub=expression(lambda^"2")) >>> >>>> So with your example, do you want something like >>> >>>> plot(c(1:5), CRP7raw[1,], type = "n", xlim=c(1,5), >>> ylim=c(-10,5) , > xlab="Day in ICU", > ylab="CRP >>> (mg/dL)", > sub = mtext(expression(lambda^2))) >>> >>> OOps! Either plot( ..., sub = *) or plot( ... ) ; >>> mtext(*) >>> >>> but not both! > >> You are right, I used a code from OP and did not much >> think about it. Strangely enough, the code worked without >> any complain. Probably mtext is "stronger" than sub and >> overrides it. > > Well, well, "the magic of R" . > Not quite: mtext(..) is a valid function call that is evaluated > before plot() finishes; then it returns NULL invisibly, and > 'plot(*, sub=NULL)' does nothing "coincidentally" (but for good > reasons) > > Martin > >> Cheers Petr > >>> >>>> CRP7graph <- apply(CRP7, 1, lines, col="gray") >>> >>>> Cheers > Petr >>> >>> >>>>> -Original Message- > > From: R-help >>> [mailto:r-help-boun...@r-project.org] On Behalf Of Rosa > >>>> Oliveira > > Sent: Friday, July 28, 2017 5:07 PM > > >>> To: r-help mailing list <r-help@r-project.org>; >>> R-help@r-project.org > > Subject: [R] Superscript and >>> subscrib R for legend x-axis and y-axis > > and colour >>> different subjects in longitudinal data with different > >>>> colours >>>>> >>>>> I am trying to make a x-axis and y-axis titles with >>> both a special > > character and a subscript. I am not >>> being able to do this. I think > > its just a placing of >>> my parenthesis, but I've tried (seemingly) everything. >>>>> >>>>> Even more, when I try the blog users code it works. >>>>> >>>>> >>>>> >>>>> Is it because I’m using longitudinal data? >>>>> >>>>> >>>>> >>>>> Even more. Is it possible to colour each one of the >>> 15 lines with a > > different colour? >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> library(ggplot2) > > library(reshape) > > >>> library(lattice) > > library(gtable) > > library(grid) >>>>> >>>>> attach(mtcars) >>>>> >>>>> beta0 = rnorm (15, 1, .5) > > beta1 = rnorm (15, -1, >>> .5) >>>>> >>>>> tempo = seq(1:5) >>>>> >>>>> CRP7raw = matrix(NA, 15, 5) > > CRP7 = matrix(NA, 15, >>> 5) > > CRP98raw = matrix(NA, 15, 5) > > CRP98 = >>> matrix(NA, 15, 5) >>>>> >>>>> crp <- for (i in 1:15) { > > CRP7raw[i,] = beta0[i] + >>> beta1[i] * tempo > > CRP7[i,] = CRP7raw[i,] + rnorm(5, 0, >>> 2.14) >>>>> >>>>> CRP98raw[i,] = beta0[i] + beta1[i] * tempo > > >>> CRP98[i,] = CRP98raw[i,] + rnorm(5, 0, .1) > > } >>>>> >>>>> >>>>> # plot(c(1:5), CRP7raw[
[R] Superscript and subscrib R for legend x-axis and y-axis and colour different subjects in longitudinal data with different colours
I am trying to make a x-axis and y-axis titles with both a special character and a subscript. I am not being able to do this. I think its just a placing of my parenthesis, but I've tried (seemingly) everything. Even more, when I try the blog users code it works. Is it because I’m using longitudinal data? Even more. Is it possible to colour each one of the 15 lines with a different colour? library(ggplot2) library(reshape) library(lattice) library(gtable) library(grid) attach(mtcars) beta0 = rnorm (15, 1, .5) beta1 = rnorm (15, -1, .5) tempo = seq(1:5) CRP7raw = matrix(NA, 15, 5) CRP7 = matrix(NA, 15, 5) CRP98raw = matrix(NA, 15, 5) CRP98 = matrix(NA, 15, 5) crp <- for (i in 1:15) { CRP7raw[i,] = beta0[i] + beta1[i] * tempo CRP7[i,] =CRP7raw[i,] + rnorm(5, 0, 2.14) CRP98raw[i,] = beta0[i] + beta1[i] * tempo CRP98[i,] =CRP98raw[i,] + rnorm(5, 0, .1) } # plot(c(1:5), CRP7raw[1,], type = "n", xlim=c(1,5), ylim=c(-10,5) , # xlab="Day in ICU", # ylab="CRP (mg/dL)", # sub = mtext(expression(paste(lambda))) # # CRP7graph <- apply(CRP7, 1, lines, col="gray") # plot(c(1:5), CRP98raw[1,], type = "n", xlim=c(1,5), ylim=c(-10,5), # xlab="Day in ICU", # ylab="CRP (mg/dL)") # CRP98graph <- apply(CRP98, 1, lines, col="gray") par(mfrow=c(1,2)) plot(c(1:5), CRP7raw[1,], type = "n", xlim=c(1,5), ylim=c(-10,5) , xlab="t_i", ylab="y_ij", sub = "lambda = 0.7") CRP7graph <- apply(CRP7, 1, lines, col="gray") plot(c(1:5), CRP98raw[1,], type = "n", xlim=c(1,5), ylim=c(-10,5), xlab="Day in ICU", ylab="CRP (mg/dL", sub = "lambda = 0.98") CRP98graph <- apply(CRP98, 1, lines, col="gray") [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] spaghetti plot - urgent
Thanks for the tip! Ulrik, I've solved the problem with a different code Best ;) Ulrik Stervbo <ulrik.ster...@gmail.com> escreveu em qua, 19/07/2017 às 20:28 : > Hi Rosa, > > You pass a vector to ggplot, which expects a data.frame. I am sure you > meant to do this: > > point7$y_point7 <- point7$beta0_7 + point7$beta1_7*point7$time + point7 > $epsilon_7 > > ggplot(point7, aes(time, y_point7)) + geom_line() > > HTH > Ulrik > > > On Wed, 19 Jul 2017 at 20:37 Rosa Oliveira <rosit...@gmail.com> wrote: > >> Hi everyone, >> >> I’m trying to do a spaghetti plot and I know I’m doing all wrong, It must >> be. >> >> What I need: >> >> 15 subjects, each with measurements over 5 different times (t1, ..., t5), >> and the variable that I need to represent in the spaguetti plot is given by: >> >> PCR = b0 + b1 * ti + epsilon >> >> B0, - baseline of each subject >> B1 - trajectory of each subject over time (so multiply by t) >> Epsilon - error associated with each subject >> >> Regression model with mixed effects. >> >> Thus, I generated b0, b1, epsilon and time created sequence. >> >> But I need to do spaguetti plot of the outcome and I can not understand >> how much I search the publications. >> >> Sorry for the stupidity, but I do not even know how to do it and I need >> it with the utmost urgency to finish a publication proposal :( >> >> Follows what I tried to do :( :( :( >> >> >> library(ggplot2) >> library(reshape) >> library(lattice) >> library(gtable) >> library(grid) >> >> >> set.seed(9027) >> >> n.longitudinal.observations = 5 # number of PCR >> measures (per subject) in the hospital period >> subjects = 15 # Number >> of simulations (1 per subject in the study) >> >> beta0_7_gerar = rnorm(subjects, mean = 1, sd = .5) >> beta0_7= >> as.data.frame(matrix(beta0_7_gerar,nrow=subjects,ncol=1)) # beta 0 - >> input variable used to calculate PCR (the outcome) >> beta1_7_gerar = rnorm(subjects, mean = -1, sd = .5) >> beta1_7 = >> as.data.frame(matrix(beta1_7_gerar,nrow=subjects,ncol=1) ) # beta 1 - >> input variable used to calculate PCR (the outcome) >> >> tj_gerar= seq.int(1, >> n.longitudinal.observations, 1) >> epsilon_7_gerar = rnorm(5*subjects, mean = 0, sd = .1) >> epsilon_7 = >> as.data.frame(matrix(epsilon_7_gerar,nrow=subjects,ncol=1) ) # epsilon_7 >> - input variable used to calculate PCR (the outcome) - associated with each >> subject >> >> tj = >> as.data.frame(matrix(tj_gerar,nrow=subjects,ncol=1) ) # >> time >> >> point7 <- cbind(beta0_7, beta1_7, tj, epsilon_7) >> point7 >> point7 <- as.data.frame(point7) >> >> colnames(point7) = c("beta0_7","beta1_7","time", "epsilon_7") >> >> >> y_point7 <- point7$beta0_7 + point7$beta1_7*point7$time + point7 >> $epsilon_7 (the outcome of the study - PCR) >> y_point7 >> >> require(ggplot2) >> >> png('test.png') >> p = ggplot(y_point7, aes(time, y_point7)) + geom_line() >> print(p) >> dev.off() >> savehistory() >> >> >> >> >> >> >> OR: >> >> In the last part I also tried: >> >> >> ID = rep(1:3, each = 5) >> >> >> point7 <- cbind(ID,beta0_7, beta1_7, tj, epsilon_7) >> point7 >> point7 <- as.data.frame(point7) >> >> colnames(point7) = c("ID","beta0_7","beta1_7","time", "epsilon_7") >> >> >> >> >> >> y_point7 <- point7$beta0_7 + point7$beta1_7*point7$time + point7 >> $epsilon_7 >> y_point7 >> >> crp7 <- y_point7 >> >> head(point7, n = 15) >> >> >> ggplot(aes(x = tj_gerar, y = crp7), data = point7) + >> geom_line(aes(group = ID), color = "gray") + >> geom_smooth(aes(group = 1), method = "lm", size = 3, color = "red", se >> = FALSE) + >> theme_bw() >> >> But none of these worked :( >> >> I was looking to have something like: >> >> >> Being the outcome PCR and the year the times (1, 2, 3, 4, 5). >> >> Can someone help me please? >> >> >> Thanks, >> >> Best Rosa >> >> >> >> __ >> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. > > -- Enviado do Gmail Mobile [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] spaghetti plot - urgent
Hi everyone, I’m trying to do a spaghetti plot and I know I’m doing all wrong, It must be. What I need: 15 subjects, each with measurements over 5 different times (t1, ..., t5), and the variable that I need to represent in the spaguetti plot is given by: PCR = b0 + b1 * ti + epsilon B0, - baseline of each subject B1 - trajectory of each subject over time (so multiply by t) Epsilon - error associated with each subject Regression model with mixed effects. Thus, I generated b0, b1, epsilon and time created sequence. But I need to do spaguetti plot of the outcome and I can not understand how much I search the publications. Sorry for the stupidity, but I do not even know how to do it and I need it with the utmost urgency to finish a publication proposal :( Follows what I tried to do :( :( :( library(ggplot2) library(reshape) library(lattice) library(gtable) library(grid) set.seed(9027) n.longitudinal.observations = 5 # number of PCR measures (per subject) in the hospital period subjects = 15 # Number of simulations (1 per subject in the study) beta0_7_gerar = rnorm(subjects, mean = 1, sd = .5) beta0_7= as.data.frame(matrix(beta0_7_gerar,nrow=subjects,ncol=1)) # beta 0 - input variable used to calculate PCR (the outcome) beta1_7_gerar = rnorm(subjects, mean = -1, sd = .5) beta1_7 = as.data.frame(matrix(beta1_7_gerar,nrow=subjects,ncol=1) ) # beta 1 - input variable used to calculate PCR (the outcome) tj_gerar= seq.int(1, n.longitudinal.observations, 1) epsilon_7_gerar = rnorm(5*subjects, mean = 0, sd = .1) epsilon_7 = as.data.frame(matrix(epsilon_7_gerar,nrow=subjects,ncol=1) ) # epsilon_7 - input variable used to calculate PCR (the outcome) - associated with each subject tj = as.data.frame(matrix(tj_gerar,nrow=subjects,ncol=1) ) # time point7 <- cbind(beta0_7, beta1_7, tj, epsilon_7) point7 point7 <- as.data.frame(point7) colnames(point7) = c("beta0_7","beta1_7","time", "epsilon_7") y_point7 <- point7$beta0_7 + point7$beta1_7*point7$time + point7 $epsilon_7 (the outcome of the study - PCR) y_point7 require(ggplot2) png('test.png') p = ggplot(y_point7, aes(time, y_point7)) + geom_line() print(p) dev.off() savehistory() OR: In the last part I also tried: ID = rep(1:3, each = 5) point7 <- cbind(ID,beta0_7, beta1_7, tj, epsilon_7) point7 point7 <- as.data.frame(point7) colnames(point7) = c("ID","beta0_7","beta1_7","time", "epsilon_7") y_point7 <- point7$beta0_7 + point7$beta1_7*point7$time + point7 $epsilon_7 y_point7 crp7 <- y_point7 head(point7, n = 15) ggplot(aes(x = tj_gerar, y = crp7), data = point7) + geom_line(aes(group = ID), color = "gray") + geom_smooth(aes(group = 1), method = "lm", size = 3, color = "red", se = FALSE) + theme_bw() But none of these worked :( I was looking to have something like: Being the outcome PCR and the year the times (1, 2, 3, 4, 5). Can someone help me please? Thanks, Best Rosa __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] outliers in Box Plot
Dear all, I have searched all over and didn´t found an answer :( I need urgently to "extract" de identification of the weight outliers of the participants of a study. So, I have a data base with several variables: id weight are 2 off them. So, I've done a boxplot and found the weight have outliers. Now I need to identify the id's of those participants. I can get the outliers values, nonetheless lots os them are not correct for some vriables. Can you please help me? boxplot(WEIGHT~AGE,baseR,range=3) #print WEIGHT boxplot boxplot(WEIGHT~AGE==1,,baseR,range=3, plot=FALSE)$out # find outliers values to age 1 (example) I attach the data. Best, Rosa Oliveira __ Antes de imprimir este e-mail pense bem se tem mesmo que o fazer. Há cada vez menos árvores. Não imprima, pense na sua responsabilidade e compromisso com o MEIO AMBIENTE! <http://pt.dreamstime.com/cora-ccedil-atildeo-criado-das-folhas-de-aacutervores-diferentes-thumb12275776.jpg> <http://pt.dreamstime.com/cora-ccedil-atildeo-criado-das-folhas-de-aacutervores-diferentes-thumb12275776.jpg> __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] substring simultaneous conditions
Dear all, I have searched all over and didn´t found an answer :( Sorry, I'm new. I need urgently to "not analyse"the weight the ID's that have 07 in the position 3 and 4 respectively, 01 or 11 in positions 11 and 12 of ID variable. . I used the following code: base<--baseR[substr(baseR$'ID',3,4)!='03' (baseR$'ID',11,12)!='01' (baseR$'ID',11,12)!='11',] But, instead of removing just the id's that respect the 3 conditions simultaneously, base don't have all the id's that have 03 in the 3 and 4 positions os id variable, neither 01 in positions 11 and 12 of ID, neither 11 in positions 11 and 12 of ID variable.variable. So, it seems, that the code exclude all the conditions, as it was a OR (|) condition in spite of AND (&) condition. Can anyone help me please? I attach the data. Best, Rosa Oliveira __ Antes de imprimir este e-mail pense bem se tem mesmo que o fazer. Há cada vez menos árvores. Não imprima, pense na sua responsabilidade e compromisso com o MEIO AMBIENTE! <http://pt.dreamstime.com/cora-ccedil-atildeo-criado-das-folhas-de-aacutervores-diferentes-thumb12275776.jpg> <http://pt.dreamstime.com/cora-ccedil-atildeo-criado-das-folhas-de-aacutervores-diferentes-thumb12275776.jpg> __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] boostrap - not converging cycles
Dear all, I have a problem, rather serious and I am not being able to resolve it. I started working with R a few time ago :( I’m running a bootstrap analyses for 2 different methods and It was expected that I get the same coefficients and standard errors for both, consequently the same coverage - but I’m not getting the same values. With the estimates obtained in bootstrap we don’t get that, which means that there are resampling cycles that are different, even more, as the coverage is different, it leads me to think that there are bootstrap cycles that do not converge and therefore left out of coverage. A. Is there a way of knowing which were the cycles that did not converge and correct the estimates? B. For this I have the "results” saved - save.image("~/Documents/phd_april_v16.RData”) will this be useful for something? Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com <mailto:rosit...@gmail.com> Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira <https://pt.linkedin.com/in/rosacsoliveira> "Many admire, few know" Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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 with OR confidence interval using probit link
Howdy everyone I’m trying to get Odds ratio and OR confidence intervals using a probit model, but I'm not getting. Do you think you can help me? I’m new with R L naive = summary(glm(pcr.data[,7]~boldBeta_individual+pcr.data$age,family=binomial(link=probit))) naive_answer= c(naive$coefficients[,1:3]) #naive estimates for #alpha (first 4 collumns: intercept; beta_intercept, beta_slope and age) and #and SE(last 4 collumns: intercept; beta_intercept, beta_slope and age) OR.naive = exp(1.6*coef(naive)) (till here works, the problem is with the confidence interval) I tried to get the Standard error from the variance, but I’m not sure if this can be done as I’ve done. Var_coef <- 1.6^2*var(coef(naive)) SE_coef <- Var_coef/sqrt(nsample)## I thi k this is correct OR.naive.inf <- exp(OR.naive - (1.96 * SE_coef)) OR.naive.sup <- exp(OR.naive + (1.96 * SE_coef)) if I used logit link I would get the CI with confint(naïve) command, but with probit I don't think so. Is there a way? What should I do? Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com <mailto:rosit...@gmail.com> Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira <https://pt.linkedin.com/in/rosacsoliveira> "Many admire, few know" Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] Odds Ratio and OR CI
Howdy everyone I’m trying to get Odds ratio and OR confidence intervals using a probit model, but I'm not getting. Do you think you can help me? I’m new with R L naive = summary(glm(pcr.data[,7]~boldBeta_individual+pcr.data$age,family=binomial(link=probit))) naive_answer= c(naive$coefficients[,1:3]) #naive estimates for #alpha (first 4 collumns: intercept; beta_intercept, beta_slope and age) and #and SE(last 4 collumns: intercept; beta_intercept, beta_slope and age) OR.naive = exp(1.6*coef(naive)) (till here works, the problem is with the confidence interval) I tried to get the Standard error from the variance, but I’m not sure if this can be done as I’ve done. Var_coef <- 1.6^2*var(coef(naive)) SE_coef <- Var_coef/sqrt(nsample)## I thi k this is correct OR.naive.inf <- exp(OR.naive - (1.96 * SE_coef)) OR.naive.sup <- exp(OR.naive + (1.96 * SE_coef)) if I used logit link I would get the CI with confint(naïve) command, but with probit I don't think so. Is there a way? What should I do? Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com <mailto:rosit...@gmail.com> Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira <https://pt.linkedin.com/in/rosacsoliveira> "Many admire, few know" Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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 meta-analysis (URGENT please)
Dear all, I’m conducting a met analysis and I usually use Revman, bur as I’m trying to use R more and more, I would like to conduct the met analysis here, in R (R-studio). One off my problems, I think, is that: 1st. it’s the first time :) 2. I only have data for 1 arm as you can see on the data that follows. ARTIGO qt tt qc tc Personal Notes qt2 Giuliani M. (2014) - - 15151862only MSM. 347 Diaz A. (2015) - - - - only MSM (n=3081) 2499 Niedźwiedzka-Stadnik M. (2015) 828 1098- - 326 Hoenig M. (2015)- - 8506- MSM (n=8925)419 Wu H. (2015)58 145 - 16713 n=16892 87 Pan X. (2015) - - only MSM (n=1316) - Ma Q. (2015)- - only MSM (n=424)- Op de Coul E. (2015)8596- HIV-infected patients (n=20965) 12369 Liu G. (2015) - - 1003 (?)- only MSM (n=1041) - some converted to HIV+ during the study - Hoenigl M. (2015) - - only MSM (n=8935) analysis HIV tests repetitions- Moller L. M. (2015) - - 469 - only MSM (N=561) 92 watkins (2015) - - - - only MSM (n=1154) only analysis believes concerning the risk- den Dass C. (2015) - - 2408- only MSM (n=3787) 589 Jia Z. (2015) - - 5314- only MSM (n=5800) 486 solomon S. (2015) - - 10875 - only MSM (n=12022) 1147 Diez M. (2014) - 3599- - n=145337 legend: qt the number of hiv subjetcs who are not sms. qt2 the number of hiv subjetcs who are sms. tt the total number of hiv subjects. qc the number of subjetcs who are sms without being hiv. tc the total number of subject not hiv. Is it possible to conduct a met analysis concerning the risk of hiv among MSM relative to the ones that are not MSM? or simply concerning the risk of hiv among MSM If yes, How? Metafor? I’ve tried, but wasn’t succeed :( Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com <mailto:rosit...@gmail.com> Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira <https://pt.linkedin.com/in/rosacsoliveira> "Many admire, few know" Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] doubt with Odds ratio - URGENT HELP NEEDED
Dear Michael (and all :)) Thank you very much. I fixed my problem, I think ;) Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira "Many admire, few know" Hippocrates > On 24 Sep 2015, at 08:51, Michael Dewey <li...@dewey.myzen.co.uk> wrote: > > Dear Rosa > > Please keep the list on the recipients as others may be able to help. > > See inline > > On 23/09/2015 19:19, Rosa Oliveira wrote: >> Dear Michael, >> >> *New cleaned code :)(I think :))* >> >> casedata <-read.spss("tas_05112008.sav") >> tas.data<-data.frame(casedata) >> >> #Delete patients that were not discharged >> tas.data <- tas.data[ tas.data$hosp!="si ",] >> tas.data$resultado.hosp <- ifelse(tas.data$hosp=="l", 0, 1) >> >> tas.data$tas_d2 <- >> log(ifelse(tas.data$tas_d2==8|tas.data$tas_d2==9, NA, >> tas.data$tas_d2)) >> tas.data$tas_d3 <- >> log(ifelse(tas.data$tas_d3==8|tas.data$tas_d3==9, NA, >> tas.data$tas_d3)) >> tas.data$tas_d4 <- >> log(ifelse(tas.data$tas_d4==8|tas.data$tas_d4==9, NA, >> tas.data$tas_d4)) >> tas.data$tas_d5 <- >> log(ifelse(tas.data$tas_d5==8|tas.data$tas_d5==9, NA, >> tas.data$tas_d5)) >> tas.data$tas_d6 <- >> log(ifelse(tas.data$tas_d6==8|tas.data$tas_d6==9, NA, >> tas.data$tas_d6)) >> >> tas.data$age <- >> ifelse(tas.data$age==8|tas.data$age==9, NA, tas.data$age) >> >> >> tas.data <- data.frame(tas1 = >> tas.data$tas_d2, tas2 = tas.data$tas_d3, >> tas3 = tas.data$tas_d4, >> tas4 = tas.data$tas_d5, >> tas5 = tas.data$tas_d6, >> age = tas.data$age, >> discharge = >> tas.data$resultado.hosp, id.pat=tas.data$id) >> >> #tas.data$discharge <- factor( tas.data$discharge , >> levels=c(0,1), labels = c("dead", "alive")) >> >> #select only cases that have more than 3 tas >> tas.data <- tas.data[apply(tas.data[,-8:-6], >> 1, function(x) sum(!is.na(x)))>2,] >> >> >> nsample <- n.obs <- dim(tas.data)[1] #nr of patients >> with more than 2 tas measurements >> >> tas.data.long <- data.frame( >> tas=c(unlist(tas.data[,-8:-6])), time=rep(c(0:4), each=n.obs), >> age=rep(tas.data$age, 5), discharge=rep(tas.data$discharge, 5), >>id=rep(c(1:n.obs), 5)) >> tas.data.long <- tas.data.long >> [order(tas.data.long$id),] >> >> age=tas.data$age >> >> >> >> library(verification) > > What does that do? > >> prevOR1 <- >> summary(glm(tas.data[,7]~tas.data[,4],family=binomial(link=probit))) >> ORmodel1 <- exp(prevOR1$coeff[,1])#computes OR? >> ORmodel1 >> >> prevOR2 <- >> summary(glm(tas.data[,7]~tas.data[,4]+tas.data[,6],family=binomial(link=probit))) >> ORmodel2 <- exp(prevOR2$coeff[,1])#computes OR? >> ORmodel2 >> > > So are you happy that those are odds ratios but you need the confidence > intervals now? > > ?confint > > may help you >> >> Nonetheless I can’t get OR confidence intervals :( and i’m not sure if I >> have it right :( >> >> Best, >> RO >> >> >> >> Atenciosamente, >> Rosa Oliveira >> >> -- >> >> >> >> Rosa Celeste dos Santos Oliveira, >> >> E-mail:rosit...@gmail.com <mailto:rosit...@gmail.com> >> Tlm: +351 939355143 >> Linkedin: https://pt.linkedin.com/in/rosacsoliveira >> >> "Many admire, few know" >> Hippocrates >> >>> On 23 Sep 2015, at 16:29, Michael Dewey <li...@dewey.myzen.co.uk >>> <mailto:li...@dewey.myzen.co.uk>> wrote: >>> >>> Dear Rosa >>> >>> Can you remove all the code which is no
Re: [R] doubt with Odds ratio - URGENT HELP NEEDED
emented for type 'list' > > > > modelo<-function (dataainit) + + { + + #dataa<-tas.data + dataa<-dataainit + + dataa$ident<-seq(1:90) + tas.6days<-cbind(id=rep(dataa$id,5),tas=c(dataa$tas_d2, dataa$tas_d3, + dataa$tas_d4, dataa$tas_d5, dataa$tas_d6), + time=rep(c(2:6)-2, each=90), out.come=rep(dataa$hosp,5), ident=rep(dataa$ident,5)) + + tas.6days<-data.frame(tas.6days[order(tas.6days[,1]),]) + tas.6days$tas[tas.6days$tas==8|tas.6days$tas==9 ]<-NA + + #mixed model for the longitudinal tas + lme.1 <- lme(tas~ time+1, random = ~ time+1 |ident, data=tas.6days, na.action=na.exclude ) + + #Random intercept and slopes + pred.lme<-predict(lme.1) + lme.intercept<-lme.1$coef$random$ident[,1]+lme.1$coef$fixed[1] + lme.slope<- lme.1$coef$random$ident[,2]+lme.1$coef$fixed[2] + selector<-as.numeric(names(lme.intercept)) #to select not NA rows. Apply to the vector in the dataset + + test(dataa$intercept[resultado.hosp==1], dataa$intercept[resultado.hosp==0]) + + print(summary(model.surv1)) + return(model.surv1$coef) + + } > I can’t get the OR and OR CI :( DATA: Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira "Many admire, few know" Hippocrates > On 23 Sep 2015, at 12:02, Michael Dewey <li...@dewey.myzen.co.uk> wrote: > > Dear Rosa > > It would help if you posted the error messages where they occur so that we > can see which of your commands caused which error. However see comment inline > below. > > On 22/09/2015 22:17, Rosa Oliveira wrote: >> Dear all, >> >> >> I’m trying to compute Odds ratio and OR confidence interval. >> >> I’m really naive, sorry for that. >> >> >> I attach my data and my code. >> >> I’m having lots of errors: >> >> 1. Error in data.frame(tas1 = tas.data$tas_d2, tas2 = tas.data$tas_d3, tas3 >> = tas.data$tas_d4, : >> arguments imply differing number of rows: 90, 0 >> > > At least one of tas_d2, tas_d3, tas_d4 does not exist > > I suggest fixing that one and hoping the rest go away. > >> 2. Error in data.frame(tas = c(unlist(tas.data[, -8:-6])), time = >> rep(c(0:4), : >> arguments imply differing number of rows: 630, 450, 0 >> >> 3. Error: object 'tas.data.long' not found >> >> 4. Error in data.frame(media = c(mean.dead, mean.alive), standarderror = >> c(se.dead, : >> arguments imply differing number of rows: 14, 10 >> >> 5. Error in ggplot(summarytas, aes(x = c(c(1:5), c(1:5)), y = mean, colour = >> discharge)) : >> object 'summarytas' not found >> >> 6. Error in summary(glm(tas.data[, 6] ~ tas.data[, 4], family = >> binomial(link = probit))) : >> error in evaluating the argument 'object' in selecting a method for >> function 'summary': Error in eval(expr, envir, enclos) : y values must be 0 >> <= y <= 1 >> >> 7. Error in wilcox.test.default(pred[obs == 1], pred[obs == 0], alternative >> = "great") : >> not enough (finite) 'x' observations >> In addition: Warning message: >> In is.finite(x) & apply(pred, 1, f) : >> longer object length is not a multiple of shorter object length >> >> >> and off course I’m not getting OR. >> >> Nonetheless all this errors, I think I have not been able to compute de code >> to get OR and OR confidence interval. >> >> >> Can anyone help me please. It’s really urgent. >> >> PLEASE >> >> THE CODE: >> >> the hospital outcome is discharge. >> >> require(gdata) >> library(foreign) >> library(nlme) >> library(lme4) >> library(boot) >> library(MASS) >> library(Hmisc) >> library(plotrix) >> library(verification) >> library(mvtnorm) >> library(statmod) >> library(epiR) >> >> # >> # Data preparation >># >> # >
[R] doubt with Odds ratio - URGENT HELP NEEDED
Dear all, I’m trying to compute Odds ratio and OR confidence interval. I’m really naive, sorry for that. I attach my data and my code. I’m having lots of errors: 1. Error in data.frame(tas1 = tas.data$tas_d2, tas2 = tas.data$tas_d3, tas3 = tas.data$tas_d4, : arguments imply differing number of rows: 90, 0 2. Error in data.frame(tas = c(unlist(tas.data[, -8:-6])), time = rep(c(0:4), : arguments imply differing number of rows: 630, 450, 0 3. Error: object 'tas.data.long' not found 4. Error in data.frame(media = c(mean.dead, mean.alive), standarderror = c(se.dead, : arguments imply differing number of rows: 14, 10 5. Error in ggplot(summarytas, aes(x = c(c(1:5), c(1:5)), y = mean, colour = discharge)) : object 'summarytas' not found 6. Error in summary(glm(tas.data[, 6] ~ tas.data[, 4], family = binomial(link = probit))) : error in evaluating the argument 'object' in selecting a method for function 'summary': Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1 7. Error in wilcox.test.default(pred[obs == 1], pred[obs == 0], alternative = "great") : not enough (finite) 'x' observations In addition: Warning message: In is.finite(x) & apply(pred, 1, f) : longer object length is not a multiple of shorter object length and off course I’m not getting OR. Nonetheless all this errors, I think I have not been able to compute de code to get OR and OR confidence interval. Can anyone help me please. It’s really urgent. PLEASE THE CODE: the hospital outcome is discharge. require(gdata) library(foreign) library(nlme) library(lme4) library(boot) library(MASS) library(Hmisc) library(plotrix) library(verification) library(mvtnorm) library(statmod) library(epiR) # # Data preparation # # setwd("/Users/RO/Desktop") casedata <-read.spss("tas_05112008.sav") tas.data<-data.frame(casedata) #Delete patients that were not discharged tas.data <- tas.data[ tas.data$hosp!="si ",] tas.data$resultado.hosp <- ifelse(tas.data$hosp=="l", 0, 1) tas.data$tas_d2 <- log(ifelse(tas.data$tas_d2==8|tas.data$tas_d2==9, NA, tas.data$tas_d2)) tas.data$tas_d3 <- log(ifelse(tas.data$tas_d3==8|tas.data$tas_d3==9, NA, tas.data$tas_d3)) tas.data$tas_d4 <- log(ifelse(tas.data$tas_d4==8|tas.data$tas_d4==9, NA, tas.data$tas_d4)) tas.data$tas_d5 <- log(ifelse(tas.data$tas_d5==8|tas.data$tas_d5==9, NA, tas.data$tas_d5)) tas.data$tas_d6 <- log(ifelse(tas.data$tas_d6==8|tas.data$tas_d6==9, NA, tas.data$tas_d6)) tas.data$age <- ifelse(tas.data$age==8|tas.data$age==9, NA, tas.data$age) tas.data <- data.frame(tas1 = tas.data$tas_d2, tas2 = tas.data$tas_d3, tas3 = tas.data$tas_d4, tas4 = tas.data$tas_d5, tas5 = tas.data$tas_d6, age = tas.data$age, discharge = tas.data$resultado.hosp, id.pat=tas.data$ID) #tas.data$discharge <- factor( tas.data$discharge , levels=c(0,1), labels = c("dead", "alive")) #select only cases that have more than 3 tas tas.data <- tas.data[apply(tas.data[,-8:-6], 1, function(x) sum(!is.na(x)))>2,] nsample <- n.obs <- dim(tas.data)[1] #nr of patients with more than 2 tas measurements tas.data.long <- data.frame( tas=c(unlist(tas.data[,-8:-6])), time=rep(c(0:4), each=n.obs), age=rep(tas.data$age, 5), discharge=rep(tas.data$discharge, 5), id=rep(c(1:n.obs), 5)) tas.data.long <- tas.data.long [order(tas.data.long$id),] age=tas.data$age ## #PLOT EMPIRICAL MEANS OF CRP FOR ALIVE DEATh ## mean.alive <- apply(tas.data[tas.data$discharge==0, -8:-6], 2, mean, na.rm=T) mean.dead <- apply(tas.data[tas.data$discharge==1, -8:-6], 2, mean, na.rm=T) stderr <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x))) se.alive<- apply(tas.data[tas.data$discharge==0, -8:-6], 2, stderr) se.dead <- apply(tas.data[tas.data$discharge==1, -8:-6], 2, stderr) summarytas <- data.frame (media = c(mean.dead, mean.alive),
Re: [R] HELP IN GRAPHS - slip screen (URGENT)
mple==1000 matplot(x=SE.alpha2$lambda[n1000],y=SE.alpha2[n1000,3:5], type="l",pch=1:3,col=c(4,2,3),yaxt="n",ylim=c(0, 1.1),ylab="") abline(h = -1, col = "gray60") mtext(expression(paste(lambda)),side=1,line=2, , cex.main=1.5) screen(2) par(mar=c(0,0,0,0)) # plot an empty plot to get the coordinates plot(0:1,0:1,type="n",axes=FALSE) legend(0,0.6,c("OLS", "GLS", "Reg. Cal.", "0"),bty = "n", lty=1:3,col=c(4,2,3,"gray60"),xpd=TRUE) close.screen(all=TRUE) Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira "Many admire, few know" Hippocrates > On 18 Sep 2015, at 10:38, Jim Lemon <drjimle...@gmail.com> wrote: > > Hi Rosa, > I have had a moment to look at your code. First I think you should start your > device as: > > quartz(width=12,height=5) > > The split.screen code that I sent seems to work for me, giving the > > 3456 >2 > 78910 > > layout of screens. To get the aspect ratio of the plots more similar, try > this: > > # do the first split, to get the rightmost screen for the legend > split.screen(figs=matrix(c(0,0.84,0,1,0.84,1,0,1),nrow=2,byrow=TRUE)) > # now split the first screen to get your eight screens (numbered 3 to 10) for > the plots > split.screen(figs=matrix(c(0,0.31,0.5,1, >0.31,0.54,0.5,1, >0.54,0.77,0.5,1, >0.77,1,0.5,1, >0,0.31,0,0.5, >0.31,0.54,0,0.5, >0.54,0.77,0,0.5, >0.77,1,0,0.5), > ncol=4,byrow=TRUE),screen=1) > > I'm not sure of which plots should go on the top line and which on the > bottom, but I think you want margins like this: > > screen(3) > par(mar=c(0,3.5,3,0)) > screen(4) > par(mar=c(0,0,3,0)) > screen(5) > par(mar=c(0,0,3,0)) > screen(6) > par(mar=c(0,0,3,0)) > screen(7) > par(mar=c(3,3.5,0,0)) > screen(8) > par(mar=c(3,0,3,0)) > screen(9) > par(mar=c(3,0,3,0)) > screen(10) > par(mar=c(3,0,3,0)) > > Perhaps this will help. > > Jim > > > On Fri, Sep 18, 2015 at 6:14 AM, Jim Lemon <drjimle...@gmail.com > <mailto:drjimle...@gmail.com>> wrote: > Hi Rosa, > I don't think the problem is with the split.screen command, for you are > getting the eight plots and the screen at the right as you requested. It > looks like your margins for each plot need adjusting, and I also think you > should have about a 2.2 to 1 width to height ratio in the graphics device. I > can't analyze the rest of the code at the moment, but perhaps tomorrow if you > can't work it out I can provide some suggestions. > > Jim > > > On Fri, Sep 18, 2015 at 1:16 AM, Rosa Oliveira <rosit...@gmail.com > <mailto:rosit...@gmail.com>> wrote: > Dear Jim, > > It works, nonetheless, it doesn't slip the screen correctly :( > > Do you have any idea? > > > I used the code: > > > #setwd("/Users/RO/Dropbox/LMER - 3rdproblem/R/latest_version/graphs/data") > setwd("~/Dropbox/LMER - 3rdproblem/R/latest_version/graphs/data") > > > library(ggplot2) > library(reshape) > library(lattice) > > > # read in what looks like half of the data > > bias.alpha2<-read.csv("graphs_bias_alpha2.csv") > SE.alpha2<-read.csv("graphs_SE_alpha2.csv") > bias.alpha1<-read.csv("graphs_bias_alpha1.csv") > SE.alpha1<-read.csv("graphs_SE_alpha1.csv") > > > > quartz(width=10,height=6) > > # do the first split, to get the rightmost screen for the legend > split.screen(figs=matrix(c(0,0.84,0,1,0.84,1,0,1),nrow=2,byrow=TRUE)) > # now split the first screen to get your eight screens (numbered 3 to 10) for > the plots > split.screen(figs=matrix(c(0,0.25,0.5,1, >0.25,0.5,0.5,1, >0.5,0.75,0.5,1, >0.75,1,0.5,1, >0,0.25,0,0.5, >0.25,0.5,0,0.5, >0.5,0.75,0,0.5, >0.75,1,0,0.5), > ncol=4,byrow=TRUE),screen=1) > > > > #split.screen(figs=matrix(c(0,0.5,0.5,1,#primeira linha primeira coluna > #
[R] HELP IN GRAPHS - slip screen
Dear all, I’m trying to do a graph, 3 rows, 5 columns, with the design: # 3 4 5 6 #2 # 7 8 9 10 I had a code for 3 rows, 3 columns, with the design:: # 3 4 #2 # 7 8 and I tried to modify it, but I had no success :( I suppose the problem is in the slip.screen code (red part of the code). I attach my code, can anyone please help me? Best, RO setwd("/Users/RO/Dropbox/LMER - 3rdproblem/R/latest_version/graphs/data") library(ggplot2) library(reshape) library(lattice) # read in what looks like half of the data bias.alpha2<-read.csv("graphs_bias_alpha2.csv") SE.alpha2<-read.csv("graphs_SE_alpha2.csv") bias.alpha1<-read.csv("graphs_bias_alpha1.csv") SE.alpha1<-read.csv("graphs_SE_alpha1.csv") quartz(width=10,height=6) # do the first split, to get the rightmost screen for the legend split.screen(figs=matrix(c(0,0.8,0,1,0.8,1,0,1),nrow=2,byrow=TRUE)) # now split the first screen to get your six screens for the plots split.screen(figs=matrix(c(0,0.5,0.5,1,#primeira linha primeira coluna 0.5,1,0.5,1,#primeira linha segunda coluna 0,0.5,0,0.5,#segunda linha primeira coluna 0.5,1,0,0.5),#segunda linha segunda coluna ncol=4,byrow=TRUE),screen=1) # this produces seven screens numbered like this: # 3 4 5 6 #2 # 7 8 9 10 # select the upper left screen screen(3) par(mar=c(0,3.5,3,0)) # now the second set n250<-bias.alpha1$nsample==250 matplot(x=bias.alpha1$lambda[n250],y=bias.alpha1[n250,3:5], type="l",pch=1:3,col=c(4,2,3),xaxt="n",ylim=c(-.1, .6),main="nsample=250",ylab="", cex.main=1) abline(h = 0, col = "gray60") mtext(expression(paste("Bias av. for ",alpha[1])),side=2,line=2, cex.main=1) screen(4) par(mar=c(0,0,3,0)) # now the second set n1000<-bias.alpha1$nsample==1000 matplot(x=bias.alpha1$lambda[n1000],y=bias.alpha1[n1000,3:5], type="l",pch=1:3,col=c(4,2,3),xaxt="n",yaxt="n",ylim=c(-.1, .6),main="nsample=1000",ylab="") abline(h = 0, col = "gray60") screen(5) par(mar=c(0,3.5,3,0)) # now the second set par(mar=c(3,3.5,0,0)) # now the second set n250<-bias.alpha2$nsample==250 matplot(x=bias.alpha2$lambda[n250],y=bias.alpha2[n250,3:5], type="l",pch=1:3,col=c(4,2,3),ylim=c(-.1, .6),ylab="") abline(h = 0, col = "gray60") mtext(expression(paste("Bias av. for ",alpha[2])),side=2,line=2, cex.main=1.5) screen(6) par(mar=c(3,0,0,0)) # now the second set n1000<-bias.alpha2$nsample==1000 matplot(x=bias.alpha2$lambda[n1000],y=bias.alpha2[n1000,3:5], type="l",pch=1:3,col=c(4,2,3),yaxt="n",ylim=c(-.1, .6)) abline(h = 0, col = "gray60") screen(7) par(mar=c(0,3.5,3,0)) # now the second set n250<-SE.alpha1$nsample==250 matplot(x=SE.alpha1$lambda[n250],y=SE.alpha1[n250,3:5], type="l",pch=1:3,col=c(4,2,3),xaxt="n",ylim=c(0, 1.1),main="nsample=250",ylab="", cex.main=1) abline(h = -1, col = "gray60") mtext(expression(paste("SE av. for ",alpha[1])),side=2,line=2, cex.main=1) mtext(expression(paste(lambda)),side=1,line=2, cex.main=1.5) screen(8) par(mar=c(0,0,3,0)) # now the second set n1000<-SE.alpha1$nsample==1000 matplot(x=SE.alpha1$lambda[n1000],y=SE.alpha1[n1000,3:5], type="l",pch=1:3,col=c(4,2,3),xaxt="n",yaxt="n",ylim=c(0, 1.1),main="nsample=1000",ylab="") abline(h = -1, col = "gray60") screen(9) par(mar=c(3,3.5,0,0)) # now the second set n250<-SE.alpha2$nsample==250 matplot(x=SE.alpha2$lambda[n250],y=SE.alpha2[n250,3:5], type="l",pch=1:3,col=c(4,2,3),ylim=c(0, 1.1),ylab="") abline(h = -.5, col = "gray60") mtext(expression(paste("SE av. for ",alpha[2])),side=2,line=2, cex.main=1.5) mtext(expression(paste(lambda)),side=1,line=2, cex.main=1.5) screen(10) par(mar=c(3,0,0,0)) # now the second set n1000<-SE.alpha2$nsample==1000 matplot(x=SE.alpha2$lambda[n1000],y=SE.alpha2[n1000,3:5], type="l",pch=1:3,col=c(4,2,3),yaxt="n",ylim=c(0, 1.1)) abline(h = -.5, col = "gray60") mtext(expression(paste(lambda)),side=1,line=2, , cex.main=1.5) screen(2) par(mar=c(0,0,0,0)) # plot an empty plot to get the coordinates plot(0:1,0:1,type="n",axes=FALSE) legend(0,0.6,c("OLS", "GLS", "Reg. Cal.", "0"),bty = "n", lty=1:3,col=c(4,2,3,"gray60"),xpd=TRUE) close.screen(all=TRUE) Best, RO Atenciosamente, Rosa Oliveira -- __
[R] split.screen to draw graphs - ggplot2 and lattice (can't slip in 4 cells)
Dear all, I need your urgent help J I’m naïve, and I’m pretty sure my doubt is very simple to solve, but I’m not getting it. I used the following code to produce my research graphs, nonetheless, is this problem, I do not have 6 graphs (1 – 6), # 3 4 5 #2 # 6 7 8 but only 4,instead. So, I need to reformulate the code, just to # 3 4 #2 # 6 7 Can you please help me reformulating? I suppose I have to change something in the split.screen code, because nowadays, my “third” column is blank J I attach the code: library(ggplot2) library(reshape) library(lattice) mean.alpha1<-read.csv("graphs_mean_alpha1.csv") mean.alpha2<-read.csv("graphs_mean_alpha2.csv") quartz(width=10,height=5) split.screen(figs=matrix(c(0,0.4,0.5,1, 0.4,0.7,0.5,1, 0.7,1,0.5,1, 0,0.4,0,0.5, 0.4,0.7,0,0.5, 0.7,1,0,0.5),nrow=6,byrow=TRUE), screen=1) screen(3) par(mar=c(0,3.5,3,0)) # now the second set n250<-mean.alpha1$nsample==250 matplot(x=mean.alpha1$lambda[n250],y=mean.alpha1[n250,3:5], type="l",pch=1:3,col=c(4,2,3),xaxt="n",ylim=c(-1.2, -0.25),main="nsample=250",ylab="", cex.main=1) abline(h = -1, col = "gray60") mtext(expression(paste("mean av. for ",alpha[1])),side=2,line=2, cex.main=1) screen(4) par(mar=c(0,0,3,0)) # now the second set n250<-mean.alpha1$nsample==1000 matplot(x=mean.alpha1$lambda[n1000],y=mean.alpha2[n1000,3:5], type="l",pch=1:3,col=c(4,2,3),xaxt="n",yaxt="n",ylim=c(-1.2, -0.25),main="nsample=1000", cex.main=1) abline(h = -1, col = "gray60") screen(6) par(mar=c(3,3.5,0,0)) # now the second set n250<-mean.alpha2$nsample==250 matplot(x=mean.alpha2$lambda[n250],y=mean.alpha2[n250,3:5], type="l",pch=1:3,col=c(4,2,3),ylim=c(-.6, -.35),ylab="") abline(h = -.5, col = "gray60") mtext(expression(paste(lambda)),side=1,line=2, cex.main=1.5) mtext(expression(paste("mean av. for ",alpha[2])),side=2,line=2, cex.main=1.5) screen(7) par(mar=c(3,0,0,0)) # now the second set n1000<-mean.alpha2$nsample==1000 matplot(x=mean.alpha2$lambda[n1000],y=mean.alpha2[n1000,3:5], type="l",pch=1:3,col=c(4,2,3),yaxt="n",ylim=c(-.6, -.35)) abline(h = -.5, col = "gray60") mtext(expression(paste(lambda)),side=1,line=2, , cex.main=1.5) screen(2) par(mar=c(0,0,0,0)) # plot an empty plot to get the coordinates plot(0:1,0:1,type="n",axes=FALSE) legend(0,0.6,c("OLS", "GLS", "Reg. Cal.", "true coefficient"),bty = "n", lty=1:3,col=c(4,2,3,"gray60"),xpd=TRUE) close.screen(all=TRUE) BEST, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira "Many admire, few know" Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] LMER - generate data and define model (2 fixed effects and 1 random effect or 1 fixed effect and 2 random effects)
Dear Thierry id.pat = rep(seq(1:nsample), each =n.longitudinal.observations) time = rep(seq(1:n.longitudinal.observations)-1, nsample) age = rnorm(nsample, mean = 36, sd = .8) f.age= table(cut(age, breaks = 8)) id.age = rep(f.age, each =5) visit = rnorm(nsample, mean = 2, sd = .03) f.visit = table(cut(visit, breaks = 3)) id.visit = rep(f.visit, each =5) 1. I was able to solve the age problem and the number of visits problem. Now the problem is another: yy~1+id.age+time+(time|id.pat) or yy~1+id.age + time+(time|id.pat) + (1|id.visit) there is no problem, but when I perform LMER, for example: lmer(yy~1+id.age+time+(time|id.pat)), there are errors: Error in model.frame.default(drop.unused.levels = TRUE, formula = yy ~ : variable lengths differ (found for 'time’) or Error in model.frame.default(drop.unused.levels = TRUE, formula = yy ~ : variable lengths differ (found for 'id.age') 2. I wasn’t trying to add the same variable to both random and fixed effects. I’m simulating a model and I’m trying to figure if it makes sense to add some variables into the model. Another doubt, 3. LMER only accepts categorical variables? Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 11 Aug 2015, at 09:00, Thierry Onkelinx thierry.onkel...@inbo.be wrote: Dear Rosa, 1) use cut() to convert a continuous variable into a factor. See ?cut for the details. 2) The syntax for factors is the same as for continuous variables. Just add the name of the factor variable to the formula fAge - cut(age) yy~1+fAge+time+(time|id.pat) 3) Add + (1|fAge) to the formula. Note that adding fAge to both the fixed and the random effect doesn't make sense. yy~1+time+(time|id.pat) + (1|fAge) Best regards, ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest team Biometrie Kwaliteitszorg / team Biometrics Quality Assurance Kliniekstraat 25 1070 Anderlecht Belgium 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 2015-08-11 1:44 GMT+02:00 Rosa Oliveira rosit...@gmail.com mailto:rosit...@gmail.com: ### # Clear memory and set the working directory and the seed ### rm(list = ls()) setwd(/Dropbox/LMER/R ) set.seed(7010) ### # Load up needed packages and do a require ### # install.packages(gdata) library(nlme) library(lme4) nsample = 1000# Number of subjects n.longitudinal.observations = 5 # number of observations per subject ### # Set the other parameters ### id.pat = rep(seq(1:nsample), each =n.longitudinal.observations) time = rep(seq(1:n.longitudinal.observations)-1, nsample) age = rnorm(nsample, mean = 36, sd = .8) id.age = rep(seq(1: n.longitudinal.observations), each =age) ### MODEL WITHOUT AGE boldBeta_individual_blup = coef(lmer(yy~1+time+(time|id.pat) ))$id.pat #mixed model ### MODEL WITH AGE boldBeta_individual_blup = coef(lmer(yy~1+age+time+(time|id.pat) ))$id.pat #mixed model Dear all, I’m trying to use LMER in my simulation problem, and I’m having problems ate the very begging L I’m new in LMER. Can you please help me? 1st problem: how do I generate age so I can use it as a fixed factor? 2nd problem: how do I
Re: [R] LMER - generate data and define model (2 fixed effects and 1 random effect or 1 fixed effect and 2 random effects)
Dear Thierry, even using the code: id.pat = rep(seq(1:nsample), each =n.longitudinal.observations) time = rep(seq(1:n.longitudinal.observations)-1, nsample) age = rnorm(nsample, mean = 36, sd = .8) f.age= cut(age, breaks = 8) id.age = rep(f.age, each =5) visit = rnorm(nsample, mean = 2, sd = .03) f.visit = cut(visit, breaks = 3) id.visit = rep(f.visit, each =5) 1st lmer(yy~1+id.age+time+(time|id.pat)) 2nd lmer(yy~1+time+(time|id.pat)+(1|id.visit)) the problem remains: variable lengths differ - for the first fixed effect in the model id.age in the first or time in the second Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 11 Aug 2015, at 09:00, Thierry Onkelinx thierry.onkel...@inbo.be wrote: Dear Rosa, 1) use cut() to convert a continuous variable into a factor. See ?cut for the details. 2) The syntax for factors is the same as for continuous variables. Just add the name of the factor variable to the formula fAge - cut(age) yy~1+fAge+time+(time|id.pat) 3) Add + (1|fAge) to the formula. Note that adding fAge to both the fixed and the random effect doesn't make sense. yy~1+time+(time|id.pat) + (1|fAge) Best regards, ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest team Biometrie Kwaliteitszorg / team Biometrics Quality Assurance Kliniekstraat 25 1070 Anderlecht Belgium 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 2015-08-11 1:44 GMT+02:00 Rosa Oliveira rosit...@gmail.com mailto:rosit...@gmail.com: ### # Clear memory and set the working directory and the seed ### rm(list = ls()) setwd(/Dropbox/LMER/R ) set.seed(7010) ### # Load up needed packages and do a require ### # install.packages(gdata) library(nlme) library(lme4) nsample = 1000# Number of subjects n.longitudinal.observations = 5 # number of observations per subject ### # Set the other parameters ### id.pat = rep(seq(1:nsample), each =n.longitudinal.observations) time = rep(seq(1:n.longitudinal.observations)-1, nsample) age = rnorm(nsample, mean = 36, sd = .8) id.age = rep(seq(1: n.longitudinal.observations), each =age) ### MODEL WITHOUT AGE boldBeta_individual_blup = coef(lmer(yy~1+time+(time|id.pat) ))$id.pat #mixed model ### MODEL WITH AGE boldBeta_individual_blup = coef(lmer(yy~1+age+time+(time|id.pat) ))$id.pat #mixed model Dear all, I’m trying to use LMER in my simulation problem, and I’m having problems ate the very begging L I’m new in LMER. Can you please help me? 1st problem: how do I generate age so I can use it as a fixed factor? 2nd problem: how do I insert age as a fixed factor? 3rd problem: what if I wanted to insert a 2nd random effect based on age? Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com mailto:rosit...@gmail.com Tlm: +351 939355143 tel:%2B351%20939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira https://pt.linkedin.com/in/rosacsoliveira
[R] LMER - generate data and define model (2 fixed effects and 1 random effect or 1 fixed effect and 2 random effects)
### # Clear memory and set the working directory and the seed ### rm(list = ls()) setwd(/Dropbox/LMER/R ) set.seed(7010) ### # Load up needed packages and do a require ### # install.packages(gdata) library(nlme) library(lme4) nsample = 1000# Number of subjects n.longitudinal.observations = 5 # number of observations per subject ### # Set the other parameters ### id.pat = rep(seq(1:nsample), each =n.longitudinal.observations) time = rep(seq(1:n.longitudinal.observations)-1, nsample) age = rnorm(nsample, mean = 36, sd = .8) id.age = rep(seq(1: n.longitudinal.observations), each =age) ### MODEL WITHOUT AGE boldBeta_individual_blup = coef(lmer(yy~1+time+(time|id.pat) ))$id.pat #mixed model ### MODEL WITH AGE boldBeta_individual_blup = coef(lmer(yy~1+age+time+(time|id.pat) ))$id.pat #mixed model Dear all, I’m trying to use LMER in my simulation problem, and I’m having problems ate the very begging L I’m new in LMER. Can you please help me? 1st problem: how do I generate age so I can use it as a fixed factor? 2nd problem: how do I insert age as a fixed factor? 3rd problem: what if I wanted to insert a 2nd random effect based on age? Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] multiple graphs with a single legend and trellis graph
Dear Jim, first of all, thank you very much :) can you please explain me how to use split.screen? I’m felling so silly, I could not run your example because of x11(width=10,height=4). I already installed package XQuartz because X11 library was missing , nonetheless, after I have installed it, “ Error in .External2(C_X11, d$display, d$width, d$height, d$pointsize, : unable to start device X11 In addition: Warning message: In x11(width = 10, height = 4) : unable to open connection to X11 display ‘' I didn’t started the x11 server on my mac ;) Now I got the graphs ;) I attach my previous graphs and my data, so you can see :) I’m very naive and new in R :( Thanks again for your help ;) Atenciosamente, Rosa Oliveira Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 08 Jul 2015, at 11:20, Jim Lemon drjimle...@gmail.com wrote: Hi Rosa, As you are using base graphics, here is an example that might be of use. As we don't have access to your data, I have used something similar to the toy data in the example for the panes function. This could be done better using the split.screen function, so let me know if you would like an example using that. library(plotrix) # start a wide plotting device x11(width=10,height=4) y-runif(100) oldpar-panes(matrix(1:6,nrow=2,byrow=TRUE),widths=c(1,1,1.7)) par(mar=c(0,2,1.8,0)) boxplot(y,axes=FALSE) axis(2) box() par(mar=c(0,0,1.8,0)) tab.title(Boxplot of y,tab.col=#88dd88,cex=1) y_hist-hist(y,axes=FALSE,breaks=seq(0,1,length.out=5)) box() tab.title(Histogram of y,tab.col=#dd8800,cex=1) par(mar=c(0,0,1.8,12)) pie(y_hist$counts,col=2:9) tab.title(Pie chart of y categories,tab.col=#dd,cex=1) box() par(mar=c(2,2,1.8,0)) plot(y,xaxs=i,xlim=c(0,101),axes=FALSE,col=2:9) axis(2) box() tab.title(Scatterplot of y,tab.col=#aabbcc,cex=1) par(mar=c(2,0,1.8,0)) plot(sort(y),xaxs=i,xlim=c(0,101),axes=FALSE,col=2:9) box() tab.title(Scatterplot of y sorted,tab.col=#ddbc44,cex=1) # center the title at the middle of the fifth plot mtext(Overall title of plot,side=1,line=0.8,cex=1.5) par(mar=c(2,0,1.8,12)) plot(diff(y),xaxs=i,xlim=c(0,100),axes=FALSE,col=2:9) axis(4) box() tab.title(Scatterplot of diff(y),tab.col=#ff33cc,cex=1) legend(115,1.8, c(Boxplot,Histogram,Pie chart,Scatterplot,Sort,Diff), fill=c(#88dd88,#dd8800,#dd,#aabbcc,#ddbc44,#ff33cc), xpd=NA) Jim On Wed, Jul 8, 2015 at 1:05 PM, David Winsemius dwinsem...@comcast.net wrote: On Jul 7, 2015, at 2:45 PM, Rosa Oliveira wrote: Iam trying to plot 6 graphs in one single plot and I was able to, nonetheless I wanted that all graphs had just 1 common legend, as the legend is the same for all the 6 graphs and there is no sense in repeating it 6 times and even more, the legends in each graph sometimes don’t fit the graph. Is there a way to put just one legend for all the 6 graphs ate the same time? I was told to use a trellis graph, but after days of trying to do that I wasn’t able to. Can anyone help me? library(ggplot2) library(reshape) library(lattice) Why did you load those packages above? As far as I can see you did not use any lattice or ggplot2 functions. (Also see no reshape or reshape2 functions in use.) par(mfrow=c(2,3)) mse.alpha1 -read.csv(file=graphs_mse_alpha1.csv,head=TRUE,sep=,) And there you lose us. We are unable to see your data. The `legend` function can put the legend anywhere. You may need to set par(xpd=TRUE) if the location is outside the current plot area. If you wnated just one legend then you mus ask yourself why you are issuing multiple legend calls. I see the token-`legend` a total of 12 times in the code below. attach(mse.alpha1) names(mse1000.alpha1) mse.alpha2 -read.csv(file=graphs_mse_alpha2.csv,head=TRUE,sep=,) attach(mse.alpha2) names(mse.alpha2) nsample==50 plot(mse.alpha1$lambda[mse.alpha1$nsample==50], mse.alpha1$mse.naive[mse.alpha1$nsample==50], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,1), cex.lab=1.5 ) lines(mse.alpha1$lambda[mse.alpha1$nsample==50],mse.alpha1$mse.RegCal[mse.alpha1$nsample==50],col=2,lty=2) lines(mse.alpha1$lambda[mse.alpha1$nsample==50],mse.alpha1$mse.PL[mse.alpha1$nsample==50],col=3,lty=3) title ( expression (paste (Mean squared error for , alpha[1])), cex.main=1.5) title(\n\n sample size=50) legend(.7,1, legend= c(Naive, Regression Calibration, Pseudo Likelihood), bty = n,col=c(4,2,3),lty=c(4,2,3)) plot(mse.alpha1$lambda[mse.alpha1$nsample==250], mse.alpha1$mse.naive[mse.alpha1$nsample==250], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4
Re: [R] multiple graphs with a single legend and trellis graph
Dear Jim, first of all, thank you very much :) can you please explain me how to use split.screen? I attach my previous graphs and my data, so you can see :) I’m very naive and new in R :( I really tried: library(plotrix) # start a wide plotting device x11(width=10,height=4) y-runif(100) oldpar-panes(matrix(1:6,nrow=2,byrow=TRUE),widths=c(1,1,1.7)) par(mar=c(0,2,1.8,0)) mse - plot(mse.alpha1$lambda[mse.alpha1$nsample==50], mse.alpha1$mse.naive[mse.alpha1$nsample==50], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,1), cex.lab=1.5 #,main=yaxs default ) lines(mse.alpha1$lambda[mse.alpha1$nsample==50],mse.alpha1$mse.RegCal[mse.alpha1$nsample==50],col=2,lty=2) lines(mse.alpha1$lambda[mse.alpha1$nsample==50],mse.alpha1$mse.PL[mse.alpha1$nsample==50],col=3,lty=3) tab.title(alpha 1 N sample=50,tab.col=#88dd88,cex=1) # tab.title(\n\n sample size=50) # problem: when I run: tab.title(Mean squared error for , paste(alpha[1])),tab.col=#88dd88,cex=1) # I get an error message: unexpected string constant in tab.title(Mean squared error for , paste(alpha[1])),tab.col= # I searched but wasn't able to fix this one, neither the other subtitle: # tab.title(\n\n sample size=50) box() par(mar=c(0,0,1.8,0)) plot(mse.alpha1$lambda[mse.alpha1$nsample==250], mse.alpha1$mse.naive[mse.alpha1$nsample==250], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,1), cex.lab=1.5 #,main=yaxs default ) lines(mse.alpha1$lambda[mse.alpha1$nsample==250],mse.alpha1$mse.RegCal[mse.alpha1$nsample==250],col=2,lty=2) lines(mse.alpha1$lambda[mse.alpha1$nsample==250],mse.alpha1$mse.PL[mse.alpha1$nsample==250],col=3,lty=3) tab.title ( alpha 1 N sample=250,tab.col=#dd8800,cex=1) box() plot(mse.alpha1$lambda[mse.alpha1$nsample==1000], mse.alpha1$mse.naive[mse.alpha1$nsample== 1000], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,1), cex.lab=1.5 #,main=yaxs default ) lines(mse.alpha1$lambda[mse.alpha1$nsample== 1000],mse.alpha1$mse.RegCal[mse.alpha1$nsample== 1000],col=2,lty=2) lines(mse.alpha1$lambda[mse.alpha1$nsample== 1000],mse.alpha1$mse.PL[mse.alpha1$nsample== 1000],col=3,lty=3) tab.title(alpha 1 N sample=1000,tab.col=#dd,cex=1) box() par(mar=c(2,2,1.8,0)) plot(mse.alpha2$lambda[mse.alpha2$nsample==50], mse.alpha2$mse.naive[mse.alpha2$nsample==50], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,.17), cex.lab=1.5 #,main=yaxs default ) lines(mse.alpha2$lambda[mse.alpha2$nsample==50],mse.alpha2$mse.RegCal[mse.alpha2$nsample==50],col=2,lty=2) lines(mse.alpha2$lambda[mse.alpha2$nsample==50],mse.alpha2$mse.PL[mse.alpha2$nsample==50],col=3,lty=3) box() tab.title(alpha 2 N sample=50,tab.col=#aabbcc,cex=1) par(mar=c(2,0,1.8,0)) plot(mse.alpha2$lambda[mse.alpha2$nsample==250], mse.alpha2$mse.naive[mse.alpha2$nsample==250], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,.17), cex.lab=1.5 #,main=yaxs default ) lines(mse.alpha2$lambda[mse.alpha2$nsample==250],mse.alpha2$mse.RegCal[mse.alpha2$nsample==250],col=2,lty=2) lines(mse.alpha2$lambda[mse.alpha2$nsample==250],mse.alpha2$mse.PL[mse.alpha2$nsample==250],col=3,lty=3) box() tab.title(alpha 2 N sample=250,tab.col=#ddbc44,cex=1) # center the title at the middle of the fifth plot mtext(Mean Squared Error,side=1,line=0.8,cex=1.5) par(mar=c(2,0,1.8,12)) plot(mse.alpha2$lambda[mse.alpha2$nsample==1000], mse.alpha2$mse.naive[mse.alpha2$nsample== 1000], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,.17), cex.lab=1.5 #,main=yaxs default ) lines(mse.alpha2$lambda[mse.alpha2$nsample== 1000],mse.alpha2$mse.RegCal[mse.alpha2$nsample==250],col=2,lty=2) lines(mse.alpha2$lambda[mse.alpha2$nsample== 1000],mse.alpha2$mse.PL[mse.alpha2$nsample== 1000],col=3,lty=3) box() tab.title(alpha 2 N sample=1000,tab.col=#ff33cc,cex=1) legend(115,1.8, c(alpha 1 N sample=50,alpha 1 N sample=250,alpha 1 N sample=1000, alpha 2 N sample=50,alpha 2 N sample=250,alpha 2 N sample=1000), fill=c(#88dd88,#dd8800,#dd,#aabbcc,#ddbc44,#ff33cc), xpd=NA) legend(115,1.8, c(Naive, Regression Calibration, Pseudo Likelihood), bty = n,col=c(4,2,3),lty=c(4,2,3), xpd=NA) and I got: Thanks again for your help ;) Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 08 Jul 2015, at 11:35, Dennis Murphy djmu...@gmail.com wrote: Hi: The general process for doing this kind of thing in either ggplot2 or lattice
Re: [R] multiple graphs with a single legend and trellis graph
Dear Jim, first of all, thank you very much :) when I run the code: myDF - rbind(mse.alpha1, mse.alpha2) # assumes both data frames have the same variables in the same order myDF$ID - factor(rep(c(alpha1, alpha2), times = c(nrow(mse.alpha1), nrow(mse.alpha2))) ) library(reshape2) myDF.melt - melt(myDF, id = c(ID, nsample, lambda), measure = c(mse.naive, mse.RegCal, mse.PL)) #Melt the three columns to be plotted, something like library(ggplot2) ggplot(myDF.melt, aes(x = lambda, y = value, color = variable)) + geom_point() + geom_line() + facet_grid(ID ~ nsample) + labs(x = expression(paste(lambda)), y = MSE, color = Type”) I get the attached graph. I also attach my data, so you can see :) I’m I able to resize the graphs differently? I.e. for alpha1: ylim=ylim=c(0,.6) and for alpha2: ylim=c(0,.17)? I reshaped, a new variable in sample was created, NA, Was it me that made something wrong? I’m very naive and new in R :( Thanks again ;) Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] multiple graphs with a single legend and trellis graph
Iam trying to plot 6 graphs in one single plot and I was able to, nonetheless I wanted that all graphs had just 1 common legend, as the legend is the same for all the 6 graphs and there is no sense in repeating it 6 times and even more, the legends in each graph sometimes don’t fit the graph. Is there a way to put just one legend for all the 6 graphs ate the same time? I was told to use a trellis graph, but after days of trying to do that I wasn’t able to. Can anyone help me? library(ggplot2) library(reshape) library(lattice) par(mfrow=c(2,3)) mse.alpha1 -read.csv(file=graphs_mse_alpha1.csv,head=TRUE,sep=,) attach(mse.alpha1) names(mse1000.alpha1) mse.alpha2 -read.csv(file=graphs_mse_alpha2.csv,head=TRUE,sep=,) attach(mse.alpha2) names(mse.alpha2) nsample==50 plot(mse.alpha1$lambda[mse.alpha1$nsample==50], mse.alpha1$mse.naive[mse.alpha1$nsample==50], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,1), cex.lab=1.5 ) lines(mse.alpha1$lambda[mse.alpha1$nsample==50],mse.alpha1$mse.RegCal[mse.alpha1$nsample==50],col=2,lty=2) lines(mse.alpha1$lambda[mse.alpha1$nsample==50],mse.alpha1$mse.PL[mse.alpha1$nsample==50],col=3,lty=3) title ( expression (paste (Mean squared error for , alpha[1])), cex.main=1.5) title(\n\n sample size=50) legend(.7,1, legend= c(Naive, Regression Calibration, Pseudo Likelihood), bty = n,col=c(4,2,3),lty=c(4,2,3)) plot(mse.alpha1$lambda[mse.alpha1$nsample==250], mse.alpha1$mse.naive[mse.alpha1$nsample==250], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,1), cex.lab=1.5 ) lines(mse.alpha1$lambda[mse.alpha1$nsample==250],mse.alpha1$mse.RegCal[mse.alpha1$nsample==250],col=2,lty=2) lines(mse.alpha1$lambda[mse.alpha1$nsample==250],mse.alpha1$mse.PL[mse.alpha1$nsample==250],col=3,lty=3) title ( expression (paste (Mean squared error for , alpha[1])), cex.main=1.5) title(\n\n sample size=250) legend(.7,1, legend= c(Naive, Regression Calibration, Pseudo Likelihood), bty = n,col=c(4,2,3),lty=c(4,2,3)) plot(mse.alpha1$lambda[mse.alpha1$nsample==1000], mse.alpha1$mse.naive[mse.alpha1$nsample== 1000], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,1), cex.lab=1.5 ) lines(mse.alpha1$lambda[mse.alpha1$nsample== 1000],mse.alpha1$mse.RegCal[mse.alpha1$nsample== 1000],col=2,lty=2) lines(mse.alpha1$lambda[mse.alpha1$nsample== 1000],mse.alpha1$mse.PL[mse.alpha1$nsample== 1000],col=3,lty=3) title ( expression (paste (Mean squared error for , alpha[1])), cex.main=1.5) title(\n\n sample size=1000) legend(.7,1, legend= c(Naive, Regression Calibration, Pseudo Likelihood), bty = n,col=c(4,2,3),lty=c(4,2,3)) plot(mse.alpha2$lambda[mse.alpha2$nsample==50], mse.alpha2$mse.naive[mse.alpha2$nsample==50], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,.17), cex.lab=1.5 ) lines(mse.alpha2$lambda[mse.alpha2$nsample==50],mse.alpha2$mse.RegCal[mse.alpha2$nsample==50],col=2,lty=2) lines(mse.alpha2$lambda[mse.alpha2$nsample==50],mse.alpha2$mse.PL[mse.alpha2$nsample==50],col=3,lty=3) title ( expression (paste (Mean squared error for , alpha[2])), cex.main=1.5) title(\n\n sample size=50) legend(.7,.17, legend= c(Naive, Regression Calibration, Pseudo Likelihood), bty = n,col=c(4,2,3),lty=c(4,2,3)) plot(mse.alpha2$lambda[mse.alpha2$nsample==250], mse.alpha2$mse.naive[mse.alpha2$nsample==250], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,.17), cex.lab=1.5 ) lines(mse.alpha2$lambda[mse.alpha2$nsample==250],mse.alpha2$mse.RegCal[mse.alpha2$nsample==250],col=2,lty=2) lines(mse.alpha2$lambda[mse.alpha2$nsample==250],mse.alpha2$mse.PL[mse.alpha2$nsample==250],col=3,lty=3) title ( expression (paste (Mean squared error for , alpha[2])), cex.main=1.5) title(\n\n sample size=250) legend(.7,.17, legend= c(Naive, Regression Calibration, Pseudo Likelihood), bty = n,col=c(4,2,3),lty=c(4,2,3)) plot(mse.alpha2$lambda[mse.alpha2$nsample==1000], mse.alpha2$mse.naive[mse.alpha2$nsample== 1000], xlab=expression(paste(lambda)),ylab=MSE,type=l,col=4,lty=4, xlim=c(.6,1), ylim=c(0,.17), cex.lab=1.5 ) lines(mse.alpha2$lambda[mse.alpha2$nsample== 1000],mse.alpha2$mse.RegCal[mse.alpha2$nsample==250],col=2,lty=2) lines(mse.alpha2$lambda[mse.alpha2$nsample== 1000],mse.alpha2$mse.PL[mse.alpha2$nsample== 1000],col=3,lty=3) title ( expression (paste (Mean squared error for , alpha[2])), cex.main=1.5) title(\n\n sample size=1000) legend(.7,.17, legend= c(Naive, Regression Calibration, Pseudo Likelihood), bty = n,col=c(4,2,3),lty=c(4,2,3)) Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates
Re: [R] graphs, need urgent help (deadline :( )
Dear Don, thank you very much. I really wasn’t being able to figure the problem. You were a big (huge) help. Seeing the graphs, I think I’ll try to put the 3 settings (sample size) in different graphs. I’ll try to use trellis graphs :) using sample size as the “factor” Thank you very much ;) Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 10 Jun 2015, at 20:07, Don McKenzie d...@u.washington.edu wrote: Here is code that IS tested. I am sending Rosa the (ugly) output in a separate file. Crazy problems with argument order; I never figured out exactly what was wrong. # therapy plot plot(therapy.df$Region[therapy.df$sample==50],therapy.df$factor.a[therapy.df$sample==50],xlab=Region,ylab=factor,type=l,col=4,ylim=c(0,1.5)) lines(therapy.df$Region[therapy.df$sample==50],therapy.df$factor.b[therapy.df$sample==50],col=2) lines(therapy.df$Region[therapy.df$sample==50],therapy.df$factor.c[therapy.df$sample==50],col=3) lines(therapy.df$Region[therapy.df$sample==250],therapy.df$factor.a[therapy.df$sample==250],col=4,lty=2) lines(therapy.df$Region[therapy.df$sample==250],therapy.df$factor.b[therapy.df$sample==250],col=2,lty=2) lines(therapy.df$Region[therapy.df$sample==250],therapy.df$factor.c[therapy.df$sample==250],col=3,lty=2) lines(therapy.df$Region[therapy.df$sample==1000],therapy.df$factor.a[therapy.df$sample==1000],col=4,lty=3) lines(therapy.df$Region[therapy.df$sample==1000],therapy.df$factor.b[therapy.df$sample==1000],col=2,lty=3) lines(therapy.df$Region[therapy.df$sample==1000],therapy.df$factor.c[therapy.df$sample==1000],col=3,lty=3) legend(7,1.4,c(factor.a,factor.b,factor.c),col=c(4,2,3),lty=1) On Jun 10, 2015, at 11:03 AM, Rosa Oliveira rosit...@gmail.com mailto:rosit...@gmail.com wrote: Sorry, I taught I attached the cvs file :) therapy.csv Don, I tried, but I got an error: my.data$Region [1] 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 my.data$sample [1] 50 50 50 50 50 50 50 50 50 50 250 250 250 250 250 250 250 250 250 250 1000 1000 1000 1000 1000 1000 1000 1000 [29] 1000 1000 my.data$factor.a [1] 0.895 0.811 0.685 0.777 0.600 0.466 0.446 0.392 0.256 0.198 0.136 0.121 0.875 0.777 0.685 0.626 0.550 0.466 0.384 0.330 0.060 0.138 0.065 [24] 0.034 0.931 0.124 0.060 0.028 0.017 0.014 plot(my.data$Region[my.data$sample==50],my.data$factor.a[my.data$sample==50],col=4,type=“l”,xlab=“Region”,ylab=“factor) Error: unexpected input in plot(my.data$Region[my.data$sample==50],my.data$factor.a[my.data$sample==50],col=4,type=�” I’m really naive, right? Best, RO Atenciosamente, Rosa Oliveira -- smile.jpg Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com mailto:rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 10 Jun 2015, at 18:10, Don McKenzie d...@u.washington.edu wrote: For a legend, try (untested) legend(0.15,0.9,c(factora,factorb,factorc),col=c(4,2,3),lty=1) If it overlaps data points move the first two arguments (0.15 and 0.9) around, or change the “ylim” argument in the plot() to ~1.2. to avoid clutter, put the line-types information in the figure caption (IMO) On Jun 10, 2015, at 10:03 AM, Don McKenzie d...@u.washington.edu wrote: On Jun 10, 2015, at 9:08 AM, Rosa Oliveira rosit...@gmail.com wrote: Dear All, I attach my data. Dear Jim, when I run your code (even the one you send me, not in my data), I get: Don't know how to automatically pick scale for object of type function. Defaulting to continuous Error in data.frame(x = c(0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, : arguments imply differing number of rows: 24, 0 Dear Don, It’s meant that I will have 12 lines: 3 factors - lines colors with 3 different values of “sample” for each - line types [Three colors, one for each factor, and three line types (lty=1,2,3), one for eachvalue of “sample - preferable dash, thin and thick). in the X - I should have region (because I have 10 regions) for each region I have the outcome of 3 different treatments (factor) for each region and each treatment I have 3 different sample size. But in your original post you had 4 sample sizes: 10,20,30,40. I need to “see” the the influence of the region in the treatment outcome for each sample size. So, at the end I should
Re: [R] graphs, need urgent help (deadline :( )
Dear Don and all, I’ve read the tutorial and tried several codes before posting :) I’m really naive. what I was trying to : is something like the graph in the picture I drawee. Is it more clear now? Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com mailto:rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 09 Jun 2015, at 19:23, Don McKenzie d...@u.washington.edu mailto:d...@u.washington.edu wrote: The answer lies in learning to use the help (and knowing where to start). Did you look at the tutorial that comes with the R installation? ?plot ?lines ?par In the last, look for the descriptions of “col” and “lty”. Using plot() and lines(), and subsetting the four unique values of “sample”, you can create your lines. Here is a crude start, assuming your columns are part of a data frame called “my.data”. Untested... plot(my.data$region[my.data$sample==10],my.data$factora[my.data$sample==10],col=4) # blue line, not dashed . . . lines(my.data$region[my.data$sample==20],my.data$factorb[my.data$sample==20],col=2,lty=2) # red dashed line On Jun 9, 2015, at 10:36 AM, Rosa Oliveira rosit...@gmail.com mailto:rosit...@gmail.com wrote: Hi, another naive question (i’m pretty sure :( ) I’m trying to plot a multiple line graph: region sample factora factorb factorc 0.1 10 0.895 0.903 0.378 0.2 10 0.811 0.865 0.688 0.1 20 0.735 0.966 0.611 0.2 20 0.777 0.732 0.653 0.1 30 0.600 0.778 0.694 0.2 30 0.466 174.592 0.461 0.1 40 0.446 0.432 0.693 0.2 40 0.392 0.294 0.686 The first column should be the independent variable, the second should compute a bold line for sample(10) and dash line for sample 20. What about the other two values of “sample”? The others variables are outcomes for each of the first scenarios, and so it should: the 3rd, 4th and 5th columns should be blue, red and green respectively. Resume :) I should have a graph, in the x-axe should have the region and in the y axe, the factor. Lines: 1 - blue and bold for region 0.1, sample 10 and factor a 2 - blue and dash for region 0.2, sample 10 and factor a 3 - red and bold for region 0.1, sample 10 and factor b 4 - red and dash for region 0.2, sample 10 and factor b 5 - green and bold for region 0.1, sample 10 and factor c 6 - green and dash for region 0.2, sample 10 and factor c Not consistent with what you said above. These are no longer lines, but points. nonetheless the independent variable is nominal, I should plot a line graph. Can anyone help me please? I have my file as a cvs file, so I first read that file (that I know how to do :)). But I have it in that format. Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com mailto:rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates [[alternative HTML version deleted]] __ R-help@r-project.org mailto:R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html http://www.r-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. PastedGraphic-1.tiff __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] graphs, need urgent help [from Rosa Oliveira]
Dear Don, I done the plot and the lines, and it’s fine. I’ll have 10 values on sample. It’s generating (on simulation), that’s why that huge outlier, and the other missing points. The graph I’ve done, is just an example, just to illustrate what I have to get, but off course with 10 points in sample, and all the other specificityies. Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 10 Jun 2015, at 02:41, Don McKenzie d...@u.washington.edu wrote: The R function plot() will draw the first line and the two axes. You need to tell it which subsample of your data to plot, as in my example below. So start with those two observations for which “sample” = 10. But if you want separate lines for each unique value of “sample”, your lines will connect only two data points, because you have only two instances of each of those unique values, unlike the lines in your hand-drawn graph. Another issue is that you have a huge outlier (value very much larger than the others) in the 6th row of “factorb”. Is this an error? If not, your other lines will be indistinguishable when you try to plot everything. Part of the reason no one else has responded may be that it appears that you are confused about your own data in a way that makes it very difficult for us to help you. Can you get some basic advice from someone local? I or someone else on the list could give you the code line-by-line that we THINK you need, but it could be wrong, given the inconsistencies in what you have shown us, and that would make everything worse. plot(my.data$region[my.data$sample==10],my.data$factora[my.data$sample==10],col=4) # blue line, not dashed Did you try plotting just this line? What happened? On Jun 9, 2015, at 5:53 PM, Rosa Oliveira rosit...@gmail.com mailto:rosit...@gmail.com wrote: Dear Don and all, I’ve read the tutorial and tried several codes before posting :) I’m really naive. what I was trying to : is something like the graph in the picture I drawee. FullSizeRender.jpg Is it more clear now? Atenciosamente, Rosa Oliveira -- smile.jpg Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com mailto:rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 09 Jun 2015, at 19:23, Don McKenzie d...@u.washington.edu mailto:d...@u.washington.edu wrote: The answer lies in learning to use the help (and knowing where to start). Did you look at the tutorial that comes with the R installation? ?plot ?lines ?par In the last, look for the descriptions of “col” and “lty”. Using plot() and lines(), and subsetting the four unique values of “sample”, you can create your lines. Here is a crude start, assuming your columns are part of a data frame called “my.data”. Untested... plot(my.data$region[my.data$sample==10],my.data$factora[my.data$sample==10],col=4) # blue line, not dashed Did you try plotting just this line? What happened? . . . lines(my.data$region[my.data$sample==20],my.data$factorb[my.data$sample==20],col=2,lty=2) # red dashed line On Jun 9, 2015, at 10:36 AM, Rosa Oliveira rosit...@gmail.com mailto:rosit...@gmail.com wrote: Hi, another naive question (i’m pretty sure :( ) I’m trying to plot a multiple line graph: regionsample factora factorb factorc 0.110 0.895 0.903 0.378 0.210 0.811 0.865 0.688 0.120 0.735 0.966 0.611 0.220 0.777 0.732 0.653 0.130 0.600 0.778 0.694 0.230 0.466 174.592 0.461 0.140 0.446 0.432 0.693 0.240 0.392 0.294 0.686 The first column should be the independent variable, the second should compute a bold line for sample(10) and dash line for sample 20. What about the other two values of “sample”? The others variables are outcomes for each of the first scenarios, and so it should: the 3rd, 4th and 5th columns should be blue, red and green respectively. Resume :) I should have a graph, in the x-axe should have the region and in the y axe, the factor. Lines: 1 - blue and bold for region 0.1, sample 10 and factor a 2 - blue and dash for region 0.2, sample 10 and factor a 3 - red and bold for region 0.1, sample 10 and factor
Re: [R] graphs, need urgent help (deadline :( )
Dear All, I attach my data. Dear Jim, when I run your code (even the one you send me, not in my data), I get: Don't know how to automatically pick scale for object of type function. Defaulting to continuous Error in data.frame(x = c(0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, : arguments imply differing number of rows: 24, 0 Dear Don, It’s meant that I will have 12 lines: 3 factors - lines colors with 3 different values of “sample” for each - line types [Three colors, one for each factor, and three line types (lty=1,2,3), one for eachvalue of “sample - preferable dash, thin and thick). in the X - I should have region (because I have 10 regions) for each region I have the outcome of 3 different treatments (factor) for each region and each treatment I have 3 different sample size. I need to “see” the the influence of the region in the treatment outcome for each sample size. So, at the end I should have 9 lines 3 red (1 dash, 1 thin, 1 thick) - concerning factor a (dash for sample size 50, thin for sample size 250 and thick for sample size 1000) 3 blue (1 dash, 1 thin, 1 thick) - concerning factor b (dash for sample size 50, thin for sample size 250 and thick for sample size 1000) 3 green (1 dash, 1 thin, 1 thick) - concerning factor c (dash for sample size 50, thin for sample size 250 and thick for sample size 1000) Hope this time is clear. I also though about doing 3 different graphs, each one for 1 different sample size, and in that case I should have 3 graphs each one with 3 lines 1 red to factor a, 1 blue to factor b and 1 green to factor c. Do you all think is better? Nonetheless I can’t do it :( best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 10 Jun 2015, at 14:13, John Kane jrkrid...@inbox.com wrote: Hi Jim, I was looking at that last night and had the same problem of visualizing what Rosa needed. Hi Rosa This is nothing like what you wanted and I really don't understand your data but would something like this work as a substitute or am I completely lost? dat1 - structure(list(region = c(0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0.2), sample = c(10L, 10L, 20L, 20L, 30L, 30L, 40L, 40L), factora = c(0.895, 0.811, 0.735, 0.777, 0.6, 0.466, 0.446, 0.392), factorb = c(0.903, 0.865, 0.966, 0.732, 0.778, 0.592, 0.432, 0.294), factorc = c(0.37, 0.688, 0.611, 0.653, 0.694, 0.461, 0.693, 0.686)), .Names = c(region, sample, factora, factorb, factorc), class = data.frame, row.names = c(NA, -8L)) mdat1 - melt(dat1, id.var = c(region, sample), variable.name = factor, value.name = value) str(mdat1) ggplot(mdat1, aes(region, value, colour = factor)) + geom_line() + facet_grid(sample ~ .) John Kane Kingston ON Canada -Original Message- From: drjimle...@gmail.com Sent: Wed, 10 Jun 2015 20:51:52 +1000 To: rosit...@gmail.com Subject: Re: [R] graphs, need urgent help (deadline :( ) Hi Rosa, Like Don, I can't work out what you want and I don't even have the picture. For example, your specification of color and line type leaves only one point for each color and line type, and the line from one point to the same point is not going to show up. Here is a possibility that may lead (eventually) to a solution. library(plotrix) par(tcl=-0.1) gap.plot(x=rep(seq(10,45,by=5),3), y=unlist(my.data[,c(factora,factorb,factorc)]), main=A plot of factorial mystery, gap=c(1.1,174),ylim=c(0,175),ylab=factor score,xlab=Group, xticlab=c( \n0.1\n10, \n0.2\n10, \n0.1\n20, \n0.2\n20, \n0.1\n30, \n0.2\n30, \n0.1\n40, \n0.2\n40), ytics=c(0,0.5,1,174.59),pch=rep(1:3,each=8),col=rep(c(4,2,3),each=8)) mtext(c(Region,Sample),side=1,at=6,line=c(0,1)) lines(seq(10,45,by=5),my.data$factora,col=4) lines(seq(10,45,by=5),my.data$factorb[c(1:5,NA,7,8)],col=2) lines(seq(10,45,by=5),my.data$factorc,col=3) Jim On Wed, Jun 10, 2015 at 10:53 AM, Rosa Oliveira rosit...@gmail.com wrote: Dear Don and all, I’ve read the tutorial and tried several codes before posting :) I’m really naive. what I was trying to : is something like the graph in the picture I drawee. Is it more clear now? Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com mailto:rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 09 Jun
Re: [R] graphs, need urgent help (deadline :( )
Sorry, I taught I attached the cvs file :) Don, I tried, but I got an error: my.data$Region [1] 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 my.data$sample [1] 50 50 50 50 50 50 50 50 50 50 250 250 250 250 250 250 250 250 250 250 1000 1000 1000 1000 1000 1000 1000 1000 [29] 1000 1000 my.data$factor.a [1] 0.895 0.811 0.685 0.777 0.600 0.466 0.446 0.392 0.256 0.198 0.136 0.121 0.875 0.777 0.685 0.626 0.550 0.466 0.384 0.330 0.060 0.138 0.065 [24] 0.034 0.931 0.124 0.060 0.028 0.017 0.014 plot(my.data$Region[my.data$sample==50],my.data$factor.a[my.data$sample==50],col=4,type=“l”,xlab=“Region”,ylab=“factor) Error: unexpected input in plot(my.data$Region[my.data$sample==50],my.data$factor.a[my.data$sample==50],col=4,type=�” I’m really naive, right? Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 10 Jun 2015, at 18:10, Don McKenzie d...@u.washington.edu wrote: For a legend, try (untested) legend(0.15,0.9,c(factora,factorb,factorc),col=c(4,2,3),lty=1) If it overlaps data points move the first two arguments (0.15 and 0.9) around, or change the “ylim” argument in the plot() to ~1.2. to avoid clutter, put the line-types information in the figure caption (IMO) On Jun 10, 2015, at 10:03 AM, Don McKenzie d...@u.washington.edu mailto:d...@u.washington.edu wrote: On Jun 10, 2015, at 9:08 AM, Rosa Oliveira rosit...@gmail.com mailto:rosit...@gmail.com wrote: Dear All, I attach my data. Dear Jim, when I run your code (even the one you send me, not in my data), I get: Don't know how to automatically pick scale for object of type function. Defaulting to continuous Error in data.frame(x = c(0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, : arguments imply differing number of rows: 24, 0 Dear Don, It’s meant that I will have 12 lines: 3 factors - lines colors with 3 different values of “sample” for each - line types [Three colors, one for each factor, and three line types (lty=1,2,3), one for eachvalue of “sample - preferable dash, thin and thick). in the X - I should have region (because I have 10 regions) for each region I have the outcome of 3 different treatments (factor) for each region and each treatment I have 3 different sample size. But in your original post you had 4 sample sizes: 10,20,30,40. I need to “see” the the influence of the region in the treatment outcome for each sample size. So, at the end I should have 9 lines 3 red (1 dash, 1 thin, 1 thick) - concerning factor a (dash for sample size 50, thin for sample size 250 and thick for sample size 1000) 3 blue (1 dash, 1 thin, 1 thick) - concerning factor b (dash for sample size 50, thin for sample size 250 and thick for sample size 1000) 3 green (1 dash, 1 thin, 1 thick) - concerning factor c (dash for sample size 50, thin for sample size 250 and thick for sample size 1000) Hope this time is clear. I also though about doing 3 different graphs, each one for 1 different sample size, and in that case I should have 3 graphs each one with 3 lines 1 red to factor a, 1 blue to factor b and 1 green to factor c. Do you all think is better? A matter of style perhaps but I would use dotplots because you have only two data points for each “line”. The lines will be misleading. You also could use panel plots, but given your skill set (unless someone wants to spend a fair bit of time with you), it’s probably best to stay as simple as possible. But given your original post (cleaned up) # untested: apologies for any typos region sample factora factorb factorc 0.1 10 0.895 0.903 0.378 0.2 10 0.8110.865 0.688 0.1 20 0.735 0.966 0.611 0.2 20 0.777 0.732 0.653 0.1 30 0.600 0.778 0.694 0.2 30 0.466 174.592 0.461 0.1 40 0.446 0.432 0.693 0.2 40 0.392 0.294 0.686 plot(my.data$region[my.data$sample==10],my.data$factora[my.data$sample==10],col=4,type=“l”,ylim=c(0,1),xlab=“region”,ylab=“factor) lines(my.data$region[my.data
[R] graphs, need urgent help (deadline :( )
Hi, another naive question (i’m pretty sure :( ) I’m trying to plot a multiple line graph: regionsample factora factorbfactorc 0.1 10 0.895 0.903 0.378 0.2 10 0.811 0.865 0.688 0.1 20 0.735 0.966 0.611 0.2 20 0.777 0.732 0.653 0.1 30 0.600 0.778 0.694 0.2 30 0.466 174.592 0.461 0.1 40 0.446 0.432 0.693 0.2 40 0.392 0.294 0.686 The first column should be the independent variable, the second should compute a bold line for sample(10) and dash line for sample 20. The others variables are outcomes for each of the first scenarios, and so it should: the 3rd, 4th and 5th columns should be blue, red and green respectively. Resume :) I should have a graph, in the x-axe should have the region and in the y axe, the factor. Lines: 1 - blue and bold for region 0.1, sample 10 and factor a 2 - blue and dash for region 0.2, sample 10 and factor a 3 - red and bold for region 0.1, sample 10 and factor b 4 - red and dash for region 0.2, sample 10 and factor b 5 - green and bold for region 0.1, sample 10 and factor c 6 - green and dash for region 0.2, sample 10 and factor c nonetheless the independent variable is nominal, I should plot a line graph. Can anyone help me please? I have my file as a cvs file, so I first read that file (that I know how to do :)). But I have it in that format. Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] problems editing R console
Thanks all off you ;) I think I got it. I was saving the workplace and loading it, but after that I wasn’t calling my data ;) really naive. Thanks very much. best RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 08 Jun 2015, at 01:30, Mark Sharp msh...@txbiomed.org wrote: Rosa, See save() and load() functions for background. However, I suspect you will want to do something as described in the article in this link http://www.fromthebottomoftheheap.net/2012/04/01/saving-and-loading-r-objects/ Mark R. Mark Sharp, Ph.D. Director of Primate Records Database Southwest National Primate Research Center Texas Biomedical Research Institute P.O. Box 760549 San Antonio, TX 78245-0549 Telephone: (210)258-9476 e-mail: msh...@txbiomed.org On Jun 7, 2015, at 5:58 PM, Rosa Oliveira rosit...@gmail.com wrote: Dear Mark, I’ll try to explain better. Imagine I write: library(foreign) library(nlme) set.seed(1000) n.sample-1 #sample size M - 5 DP_x - 2 x - rnorm(n.sample,M,DP_x) p - pnorm(-3+x) y - rbinom(n.sample,1,p) dp_erro - 0.01 erro - rnorm(n.sample,0,dp_erro) x.erro - x+erro but with a function, with 2000 simulations. I save my “output” and I get X.erro in a .txt file. (text edit file). I do another setting with DP_x=3 and save, and so on. For some reason I realize I’ve done my simulation the wrong way and I have to apply a correction, for example: x.erro = 1.4X+erro, i.e. in the truth I could use my first X and erro values in each setting, but as it is in a .txt file I can’t use them any more. Is there a way to save the results in a format that I can use the values? Just apply my corrections and don’t have to do the 2000 simulations for each setting again? My problem is that the function I use takes 3 days running, and just 500 simulations :( Best, RO Atenciosamente, Rosa Oliveira -- smile.jpg Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 07 Jun 2015, at 23:03, Mark Sharp msh...@txbiomed.org wrote: I cannot understand your request as stated. Can you provide a small example? Mark R. Mark Sharp, Ph.D. msh...@txbiomed.org On Jun 7, 2015, at 2:49 PM, Rosa Oliveira rosit...@gmail.com wrote: Dear all, I’m doing simulations on R, and as my code is being changed and improved I need to, sometimes, work in finished simulations, i.e, After my simulation is over I need to settle another setting. The problem is that I need to get back to the previous result. When I save the result it saves as txt, so I can’t edit that result any more. Imagine I save a setting and save the mean, nonetheless, in another setting the mean as problems, so I have to ask the median. As I have to have the same statistics to all settings, nowadays I have to run my first setting again. My advisor told me that I could save another way so I can “edit” my first result. Is it possible? I tried to save as save my workplace, … but after I don’t know what to do with it. Can you please help me? I know is a naive question, but I have to go through this every 3 days (time each simulation takes long). And my work is being delayed :( Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] problems editing R console
Dear all, I’m doing simulations on R, and as my code is being changed and improved I need to, sometimes, work in finished simulations, i.e, After my simulation is over I need to settle another setting. The problem is that I need to get back to the previous result. When I save the result it saves as txt, so I can’t edit that result any more. Imagine I save a setting and save the mean, nonetheless, in another setting the mean as problems, so I have to ask the median. As I have to have the same statistics to all settings, nowadays I have to run my first setting again. My advisor told me that I could save another way so I can “edit” my first result. Is it possible? I tried to save as save my workplace, … but after I don’t know what to do with it. Can you please help me? I know is a naive question, but I have to go through this every 3 days (time each simulation takes long). And my work is being delayed :( Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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] problems editing R console
Dear Mark, I’ll try to explain better. Imagine I write: library(foreign) library(nlme) set.seed(1000) n.sample-1 #sample size M - 5 DP_x - 2 x - rnorm(n.sample,M,DP_x) p - pnorm(-3+x) y - rbinom(n.sample,1,p) dp_erro - 0.01 erro - rnorm(n.sample,0,dp_erro) x.erro - x+erro but with a function, with 2000 simulations. I save my “output” and I get X.erro in a .txt file. (text edit file). I do another setting with DP_x=3 and save, and so on. For some reason I realize I’ve done my simulation the wrong way and I have to apply a correction, for example: x.erro = 1.4X+erro, i.e. in the truth I could use my first X and erro values in each setting, but as it is in a .txt file I can’t use them any more. Is there a way to save the results in a format that I can use the values? Just apply my corrections and don’t have to do the 2000 simulations for each setting again? My problem is that the function I use takes 3 days running, and just 500 simulations :( Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates On 07 Jun 2015, at 23:03, Mark Sharp msh...@txbiomed.org wrote: I cannot understand your request as stated. Can you provide a small example? Mark R. Mark Sharp, Ph.D. msh...@txbiomed.org On Jun 7, 2015, at 2:49 PM, Rosa Oliveira rosit...@gmail.com wrote: Dear all, I’m doing simulations on R, and as my code is being changed and improved I need to, sometimes, work in finished simulations, i.e, After my simulation is over I need to settle another setting. The problem is that I need to get back to the previous result. When I save the result it saves as txt, so I can’t edit that result any more. Imagine I save a setting and save the mean, nonetheless, in another setting the mean as problems, so I have to ask the median. As I have to have the same statistics to all settings, nowadays I have to run my first setting again. My advisor told me that I could save another way so I can “edit” my first result. Is it possible? I tried to save as save my workplace, … but after I don’t know what to do with it. Can you please help me? I know is a naive question, but I have to go through this every 3 days (time each simulation takes long). And my work is being delayed :( Best, RO Atenciosamente, Rosa Oliveira -- Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira Many admire, few know Hippocrates __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.