Re: [R] Visualising multiple response contingency tables
Not sure I understand your question (or if there is one) and I am not
familiar with vcd::mosaic. But if you are asking is there a simpler
way ? than yes:
1. work with ?array and ?aperm
2. create the array directly in R from the original data - not excel
3. ?mosaicplot (no package required - it's in grid)
Here is what I mean based on your f.tbl:
>> f.tbl = structure(c(10, 15, 25, 45, 30, 50), .Dim = 2:3, .Dimnames =
>> structure(list(Sex = c("F", "M"), Responses = c("A", "B", "total
>> subjects")), .Names = c("Sex", "Responses")), class = "table")
# Calculate the No-A No-B columns:
(ff.tbl <- rbind(f.tbl[,1:2],f.tbl[,3]-f.tbl[,1:2]))
# rearrange to a CxRxB (in this case 2x2x2) array:
dim(ff.tbl) <- c(2,2,2)
# give some names
dimnames(ff.tbl) <- list(Sex=c('F','M'),c('yes','no'),Response=c('A','B'))
ff.tbl
# plot
mosaicplot(ff.tbl)
# or plot
mosaicplot(aperm(ff.tbl,3:1))
# or test
apply(ff.tbl, 3 , chisq.test) # and sum the result
Hope this helps get you started
> f.tbl Responses
> Sex A B total subjects
> F 10 25 30
> M 15 45 50
>
>
> The answer I have is to adjust my data and then use the mosaic() function
> in package:vcd; however, I'm not sure that's the best way forward and I
> don't have a very efficient way of getting there. I will present my
> solution so you guys can take a look.
>
> The fundamental problem is that because of the multiple response data, you
> can't simply apply a normal Chi-square test to the contingency table.
> There's a raft of approaches, but I've decided to use a simple technique
> introduced by (A. Agresti, I. Liu, Modeling a categorical variable allowing
> arbitrarily many category choices, Biometrics 55 (1999) 936-43.) and
> refined by Thomas and Decady and Bilder and Loughin. In summary, the test
> statistic (a modified Chi square statistic) is calculated by summing up the
> individual chi-square statistics for each of the c marginal r в 2 tables
> relating the single response variable to the multiple response variable
> with df = c(r - 1)). Note, that instead of using the row totals (total
> number of responses) the test statistic is calculated with the total number
> of subjects per row.
>
> (phew, I hope that made sense :) ) Unfortunately, my google-research has
> not revealed an easy way to transform my one data table into c x r x 2
> tables for analysis. So I end up having to create the two different tables
> myself, shown below (note that the Not-A/B columns are calculated as the
> difference between the main data column (A/B) and the total number of
> subjects listed above.
>
>> g.mtrx=matrix(c(10,15,20,35),nrow=2)> g.tbl=as.table(g.mtrx)>
>> dimnames(g.tbl)=list(Sex=c("F","M"),Responses=c("A","Not-A"))> g.tbl
>> Responses
> Sex A Not-A
> F 10 20
> M 15 35
>
>> h.tbl=as.table(h.mtrx)> h.mtrx=matrix(c(25,45,5,5),nrow=2)>
>> h.tbl=as.table(h.mtrx)>
>> dimnames(h.tbl)=list(Sex=c("F","M"),Responses=c("B","Not-B"))> h.tbl
>> Responses
> Sex B Not-B
> F 25 5
> M 45 5
>
>
> If I then preform the normal Chi-square test on each of the two tables
> (chisq.test()) and then sum up the results, I get the answer I want.
> Clearly this is cumbersome, which is why I do it in Excel at the moment (I
> know shame on me). However, I really want to take advantage of the mosaic
> function in vcd. So what I have to do at the moment is create the tables
> above and use abind() (package:abind) to bring my two matrices together to
> form a multidimensional matrix. Example:
>
>> gh.abind = abind(g.mtrx,h.mtrx,along=3)>
>> dimnames(gh.abind)=list(Sex=c("F","M"),Responses=c("Yes","No"),Factors=c("A","B"))>
>> gh.abind, , Factors = A
>
> Responses
> Sex Yes No
> F 10 20
> M 15 35
>
> , , Factors = B
>
> Responses
> Sex Yes No
> F 25 5
> M 45 5
>
> Now I can use the simple mosaic function to plot the combined matrix
>
>> mosaic(gh.abind)
>
> So that's it. I don't use any pearson-r shading in mosaic since I
> don't think it would be appropriate to try and model my weird multiple
> response tables (at the moment), but what I will do is look at the
> odds-ratio table and then manually colour the mosaic cells with high
> odds-ratios (greater than 2).
>
> I am literally having to type all this by hand into R, and as you can
> imagine, it gets cumbersome with large multi column tables (which I
> have). Does any body have any thoughts on my approach of using mosaic
> for this sort of data? And if so, any insight on how I can be a bit
> slicker with my R code?
>
> All help is appreciated and I hope that this question wasn't too long
> to read through.Not sure I uderstand your question (or if there is one) and I
> am not familiar with vcd::mosaic. But if you are asking is there a simpler
> way ? than yes:
1. work with ?array and ?aperm not tables
2. create the array directly in R from the original data
3. ?mosaicplot (no package required - it's in grid)
Here is what I mean based on your f.tbl:
>> f.tbl = structu
Re: [R] 3D Black-Scholes Graph Help!
On Tue, Mar 13, 2012 at 3:34 PM, David Winsemius wrote: > > When I got around to running it I was hampered by a lack of knowledge about > what sort of data-object "price" might have been. I tried putting in a > single number on hte theory that it would saitisfy the seq() call, and also > got the error you report. More input is needed from the OP about the > problem specification, and hopefully she will provide a test dataset. It is just a number. The OP function works fine, the error was generated by wireframe because of the "partial" data passed to formula - OptionPrice is a matrix, not a column in the grid data.frame which only holds the scales. Anna, Replace call to wireframe in your function with wireframe(OptionPrice,main="3D Option", drape=T, col.regions=heat.colors(100), scales = list(arrows=F,x=list(at=1:length(t),labels=t),y=list(at=1:length(s),labels=s))) Then: plotbs(16) Note, the first value is always NaN but that's from your calculations - look at the OptionPrice matrix. Cheers Elai > -- > David. > > >> >> Berend >> > > David Winsemius, MD > West Hartford, CT > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] barplot and NA
d2 <- as.matrix(c(2,NA,4)) barplot(d2,beside=T) barplot(c(d2)) barplot(na.omit(d2)) d2[2,] <- 0 barplot(d2) # So barplot is not "stopping" at the first NA (first 2 plots). But what does stacking even mean when you have a missing group in the middle ? you can't expect barplot to know... if you think it means 0 and the rest can just be stacked on top - define it that way. Cheers On Sun, Mar 11, 2012 at 10:22 PM, John D. Muccigrosso wrote: > Am I wrong that barplot is supposed to just skip NAs, and continue with the > rest of the data in a matrix column? That's how I read various posts on the > subject. > > But that's not what happens for me with R64.app (on a Mac, obviously). For > example: > > d0 <- as.matrix(c(2,3,4)) > d1 <- as.matrix(c(2,3,NA)) > d2 <- as.matrix(c(2,NA,4)) > d3 <- as.matrix(c(NA,3,4)) > barplot(d0) > barplot(d1) > barplot(d2) > barplot(d3) > > generates four bar plots. The first has one bar with three visible bands, as > expected. The second has two bands; still OK. But the third has only one band > (at 2) and the fourth has none. > > So it appears that barplot is barfing on those NAs and stopping its plot at > those points. > > Is that the expected behavior? > > Thanks. > > John > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Use different panel functions with lattice
Inline
On Sat, Mar 10, 2012 at 1:47 PM, Balaitous wrote:
> Le samedi 10 mars 2012 à 12:25 -0700, ilai a écrit :
>> On Sat, Mar 10, 2012 at 9:33 AM, Balaitous wrote:
>
> Var1 and Var2 are 2 two different observed variables (with different scales)
You might want to consider scales=list(y=list(relation='free')) in ?xyplot
> Var3 is the time
> Var4 is the point of observation
>
> I have also a Var5 for groups, but I just want groups for the Var1.
>
> But I don't know how to make the test
> if(Varx)
> in the function panel.mypanel, because I need
>
> Var1 -> panel.superpose (It's OK)
> Var2 -> panel.lines (I don't want groups for this)
>
> (And I will have others variables with other panel functions to use)
>
Since outer=T (i.e. Var1 and Var2 are in different panels), at the
beginning of the panel or panel.groups function, try
if(packet.number() %in% 1:3) {
panel.rect(x,y,groups,...) # or whatever for panels 1:3
}
else{
panel.rect(x,y,groups,col=constant,...) # or some other stuff for panels 4:6
}
Hope that works better.
>> >
>> > Something like :
>> >
>> > panel.mypanel = function(x, y, ...) {
>> > if (Var1) panel.Var1Panel(x, y, ...)
>> > else panel.Var2Panel(x, y, ...)
>> > }
>> > xyplot(Var1+Var2~Var3|Var4, data=df, panel=panel.mypanel)
>> >
>> > (I have search with google, but I found nothing)
>> >
>> > Thanks
>> >
>> > __
>> > R-help@r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide
>> > http://www.R-project.org/posting-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Use different panel functions with lattice
On Sat, Mar 10, 2012 at 9:33 AM, Balaitous wrote:
> Hi,
>
> I have a data.frame df with
> names(df) = c("Var1", "Var2", "Var3", "Var4")
>
> and I plot data with
>
> xyplot(Var1+Var2~Var3|Var4, data=df)
>
> I want to use different panel functions for Var1 and Var2.
> How can I do ?
You didn't specify which "different panel functions" you want. Is
something like this what you're looking for?
xyplot(Var1+Var2~Var3|Var4, data=df, panel=panel.superpose,
panel.groups=function(x , y , group.number , ...){
panel.xyplot(x , y[group.number==1] , ...)
panel.lines(x , y[group.number==2] , lwd=2 , col=1)
})
>
> Something like :
>
> panel.mypanel = function(x, y, ...) {
> if (Var1) panel.Var1Panel(x, y, ...)
> else panel.Var2Panel(x, y, ...)
> }
> xyplot(Var1+Var2~Var3|Var4, data=df, panel=panel.mypanel)
>
> (I have search with google, but I found nothing)
>
> Thanks
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] layer plots.
It's hard to help if you keep changing the framework of your problem,
first two matrices - now it's a data.frame and a list of subset row
names in a plotting method from whatever package "suprow" comes from.
Regardless, Michael's original answer already gave you a solution:
plot(table1,type='l',lwd=2)
points(table1[list1$node,],col=2,pch=19)
The points are overlayed on the line plot, so they are not "obscured"
if you have 20 or 20M values.
On Fri, Mar 9, 2012 at 11:09 AM, aaral singh wrote:
> The response much appreciated. They do match up, one is a small subset of
> the other.
>
> I have this:
>> dput(table1)
> structure(list(var1 = c(2L, 4L, 4L, 1L, 423L), var2 = c(3L, 5L,
> 6L, 342L, 3L)), .Names = c("var1", "var2"), class = "data.frame", row.names
> = c("node1",
> "node2", "node3", "node4", "node5"))
>
>> dput(list1)
> structure(list(node = c("node1", "node2")), .Names = "node")
>
> so one table is a 2 X 5 matrix (called table1) and one table is 1 X 2 table
> (called list1).
>
>
> i then type this:
>
>> plot1 <-plot(table,suprow=c(list1$node),"passive")
> to give me a plot of list1
>
> and this:
>
>> plot2 <-plot(table,suprow=c(list1$node),"active")
> to give me a plot of table1
>
> i want to combine plot 1 and 2.
>
> BUT i know i can do this:
>> plot2 <-plot(table,suprow=c(list1$node),"all") to plot both on the same
> graph,
> but in my actual dataset, the points in list1 are obscured from sight by
> table1, because in reality table 1 may contain 20,000 points and list1 may
> contain 10 points, so i cannot see where my 10 specific nodes of interest
> are on the graph. So i want to plot the graph so that any nodes in list1
> are seen on top of the plot of table 1.
>
>
>
>
> On Fri, Mar 9, 2012 at 5:36 PM, Michael Weylandt [via R] <
> ml-node+s789695n4460118...@n4.nabble.com> wrote:
>
>> Do your matrices "match up" with each other in any meaningful way or
>> do you just want two independent plots on a single page?
>>
>> You should probably provide the dput() output of each table object so
>> we can see what you've got.
>>
>> Michael
>>
>> On Fri, Mar 9, 2012 at 11:07 AM, aoife doherty <[hidden
>> email]>
>> wrote:
>>
>> > Many thanks for reply.
>> > I have trouble understanding how to use response, i am sorry.
>> > My question is i have two matrices. I then plot two matrices. Then I
>> have 2
>> > seperate plots. I can color the nodes in the plots in two different
>> colors.
>> > Then, how do i merge the two plots to view one overlapping the other?
>> i.e.
>> > to view two sets of data in one 2D space?
>> >
>> > Many thanks
>> >
>> >
>> > On Fri, Mar 9, 2012 at 3:51 PM, R. Michael Weylandt
>> > <[hidden email] >
>> wrote:
>> >>
>> >> No idea what table1, table2 are
>> >>
>> >> plot(1:5, type = "l")
>> >> points(5:1, col = 2)
>> >>
>> >> should get you started.
>> >>
>> >> Michael
>> >>
>> >> On Fri, Mar 9, 2012 at 10:17 AM, aaral singh <[hidden
>> >> email]>
>>
>> >> wrote:
>> >> > Hello.
>> >> >
>> >> > I have 2 plots.
>> >> >
>> >> >> plot1 <-plot(table1)
>> >> >> plot2 <-plot(table2)
>> >> >
>> >> > How may i plot these both on the same graph, i.e. layer one graph on
>> top
>> >> > of
>> >> > the other one.
>> >> > The result should look similar to this the image below, where the
>> black
>> >> > lines indicate one plot, and the red dots indicate the second plot.
>> >> >
>> >> > http://r.789695.n4.nabble.com/file/n4459732/R_screen_shot.png
>> >> >
>> >> > Aaral.
>> >> >
>> >> >
>> >> > --
>> >> > View this message in context:
>> >> > http://r.789695.n4.nabble.com/layer-plots-tp4459732p4459732.html
>> >> > Sent from the R help mailing list archive at Nabble.com.
>> >> >
>> >> > __
>> >> > [hidden email]
>> >> > mailing list
>> >> > https://stat.ethz.ch/mailman/listinfo/r-help
>> >> > PLEASE do read the posting guide
>> >> > http://www.R-project.org/posting-guide.html
>> >> > and provide commented, minimal, self-contained, reproducible code.
>> >>
>> >> __
>> >> [hidden email]
>> >> mailing list
>> >> https://stat.ethz.ch/mailman/listinfo/r-help
>> >> PLEASE do read the posting guide
>> >> http://www.R-project.org/posting-guide.html
>> >> and provide commented, minimal, self-contained, reproducible code.
>> >
>> >
>>
>> __
>> [hidden email] mailing
>> list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>> --
>> If you reply to this email, your message wil
Re: [R] How to eliminate for next loops in this script
?by ?aggregate
On Tue, Mar 6, 2012 at 4:14 PM, Walter Anderson wrote:
> I needed to compute a complicated cross tabulation to show weighted means
> and standard deviations and the only method I could get that worked uses a
> series of nested for next loops. I know that there must be a better way to
> do so, but could use some assistance pointing the way.
>
> Here is my working, but inefficient script:
>
> library(Hmisc)
> rm(list=ls())
> load('NHTS.Rdata')
> day.wt <- day[c("HOUSEID","WTTRDFIN", "TRIPPURP","TRVLCMIN")]
> hh.wt <- hh[c("HOUSEID","WTHHFIN","HHSIZE","HHVEHCNT","HOMETYPE")]
> hh.wt$HHBIN <- with(hh.wt,{ cut(HHSIZE,
> breaks=c(0,1,2,3,4,max(HHSIZE)),labels=c("1","2","3","4","5+"),ordered_result=TRUE)})
> hh.wt$VEHBIN <- with(hh.wt,{ cut(HHVEHCNT,
> breaks=c(-1,0,1,2,max(HHVEHCNT)),labels=c("0","1","2","3+"),ordered_result=TRUE)})
> hh.wt$DUTYPE <- factor(hh.wt$HOMETYPE, exclude=c("-7","-8","-9"))
> levels(hh.wt$DUTYPE) <- c("1","1","2","2","2","2","2") # Convert home
> types to 1=SF and 2=MF
> trips <- merge(day.wt,hh.wt,by="HOUSEID")
> trips$PURPOSE <- factor(trips$TRIPPURP, exclude=c("-9"))
> trips <- trips[c("HOUSEID","HHBIN","VEHBIN","DUTYPE","PURPOSE","WTTRDFIN",
> "WTHHFIN")]
> sink('TripRates.csv')
> cat("HHBIN,VEHBIN,DUTYPE,TRIPPURP,COUNT,MEAN,STD.DEV\n")
> for (per in levels(trips$HHBIN))
> for (veh in levels(trips$VEHBIN))
> for (hom in levels(trips$DUTYPE))
> for (pur in levels(trips$PURPOSE))
> {
> cat(per)
> cat(",")
> cat(veh)
> cat(",")
> cat(hom)
> cat(",")
> cat(pur)
> cat(",")
> tmp <- subset(trips,
> trips$HHBIN == per & trips$VEHBIN == veh &
> trips$DUTYPE == hom & trips$PURPOSE == pur,
> select=c("HOUSEID","WTTRDFIN", "WTHHFIN"))
> cat(length(tmp$HOUSEID))
> cat(",")
> tmp$tr <- (tmp$WTTRDFIN/365)/tmp$WTHHFIN
> w.m <- wtd.mean(tmp$tr,weights=tmp$WTHHFIN)
> w.sd <- sqrt(wtd.var(tmp$tr, weights=tmp$WTHHFIN))
> cat(w.m)
> cat(",")
> cat(w.sd)
> cat("\n")
> }
> sink()
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] confidence intervals in dotplots in a for loop
On Tue, Mar 6, 2012 at 8:55 PM, Byerly, Mike M (DFG)
wrote:
>
> estimates <-
c(67.42,30.49,32.95,23.53,10.26,6.03,23.53,0.93,50.72,24.2,25.84,18.54,
7.16,3.6,9.35,0.33,87.28,37.25,40.16,28.59,13.77,8.92,40.74,1.68,48.28,23.09,
24.49,17.7,6.63,3.28,7.79,0.26,91.63,38.74,41.6,29.74,14.49,9.51,44.16,1.88)
> estimates.m<- matrix(estimates, 8,5)
> colnames(estimates.m) <- list("est","lci90","uci90","lci95","uci95")
> id <- c(1,1,2,2,3,3,4,4)
> size <- c("All","All","Large","Large","Medium","Medium","Small","Small")
> gear <- rep(c("Camera","Dredge"),4)
> cds.est <- data.frame(id,size,gear,estimates.m)
# From this point try this instead:
library(lattice)
library(gridExtra)
dat <- by(cds.est,cds.est$id,function(x) x)
plots <- lapply(dat,function(dat){
dotplot(gear ~ est , data=dat, xlim = c(min(dat$lci95-10),
max(dat$uci95 +10)), xlab='Density (num / Nmi^2)',
main = as.vector(unique(dat$size)),
panel=function(x,y) {
panel.xyplot(x,y,pch=16,cex=1)
panel.segments(dat$lci95,as.numeric(y),dat$uci95,as.numeric(y),
lty=1, col=1)
panel.segments(dat$lci90,as.numeric(y),dat$uci90,as.numeric(y),
lty=1,lwd=4, col='grey60')
panel.xyplot(x,y,pch=16,cex=1.2,col='white')
panel.xyplot(x,y,pch=1,cex=1.1, col='black')
}
)
})
do.call(grid.arrange, plots)
> # cds.est data frame and stores the plots in a list. Since dotplot() is
> part of the # lattice package, I used grid.arrange to tile the plots.
In my opinion you are going about it the wrong way. You are not using
any of the functionality provided by lattice or dotplot for scaling
and panel arrangement. All you use are panel.xyplot and
panel.segments, manually calculating limits, labels and arranging the
plots yourself - you could have just done this in base graphics (and a
loop). But if that's how you want it to look...
HTH
>
>
>
> library(grid)
>
> library(gridExtra)
>
>
>
> id2 <- 1:max(cds.est$id)
>
> plots <- vector("list", length(id2))
>
> j <- 0
>
> for (i in id2) {
>
> dat <- cds.est[cds.est$id == i,]
>
> plots[[ j <- j+1]] <-
>
> dotplot(gear ~ est , data=dat, xlim = c(min(dat$lci95
> -10), max(dat$uci95 +10)), xlab='Density (num / Nmi^2)',
>
> main = as.vector(unique(dat$size)),
>
> panel=function(x,y) {
>
> panel.xyplot(x,y,pch=16,cex=1)
>
>
> panel.segments(dat$lci95,as.numeric(y),dat$uci95,as.numeric(y), lty=1,
> col=1)
>
>
> panel.segments(dat$lci90,as.numeric(y),dat$uci90,as.numeric(y), lty=1,
> lwd=4, col='grey60')
>
>
> panel.xyplot(x,y,pch=16,cex=1.2,col='white')
>
> panel.xyplot(x,y,pch=1,cex=1.1,
> col='black')
>
> })
>
>
>
> }
>
> do.call(grid.arrange, plots)
>
>
>
> # The script runs and produces all the plots with the correct estimates,
> but the CI's are not
>
> # plotting correctly. Does anyone have suggestions on what is causing
> this?
>
>
>
> Thanks, Mike
>
>
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] GLMing with all possible combinations (automatically)
Hi try glm(Response~ .^2, data=yourdata.frame) For all predictors (.) and 2-way interactions (^2). You might also want to see ?drop.terms and ?formula for automating the construction of all model combinations. Side note: R is not SAS (fortunately). Interaction is denoted ":", X*Y is shorthand for X+Y+X:Y. Your "X+Y+X*Y" is redundant. HTH On Wed, Mar 7, 2012 at 4:16 AM, Christofer Bogaso wrote: > Dear all, I was working with a GLM model and wondering whether is it > possible to have a model description with all possible main effects as > well as all interactions. Currently I do it manually however this will > become cumbersome if I have large number of variables. > > Means suppose I have 3 explanatory variables and want to feed > automatically with all possible combinations like: > > > glm(Response~1) > > glm(Response~X) > glm(Response~Y) > glm(Response~Z) > > glm(Response~X+Y) > glm(Response~X+Z) > glm(Response~Z+Y) > > glm(Response~X+Y+X*Y) > glm(Response~X+Z+X*Y) > > all combinations > which are possible > > Thanks and regards, > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Grouped barchart confidence intervals in lattice
packet.number() will retrieve the panel you are working in. With a
third conditioning variable, you now have say a list of *.err$x, one
for each panel/level of 3rd var. In that case you need
dsl[[packet.number()]] to get at the right place in the list (or
[,packet.number()] if it's a matrix etc.). Regarding the calculations,
you could type panel.barchart in the console to see Deepayan's
original. Sorry but I just don't remember what my logic was when
figuring it out (it has been a short while since my last barchart),
was probably just trial and error until it worked.
Best,
On Sun, Mar 4, 2012 at 1:15 PM, Nathan Lemoine wrote:
> Also, can you tell me where to add the [[packet.number()]] command for
> condition? I've been trying several things with no luck.
>
> Thanks for helping, sorry if I'm a nuisance, but I'm still learning my way
> around the more advanced functions and there's no one else at my university
> that uses R as extensively as I do.
>
> Nathan
On Sun, Mar 4, 2012 at 11:45 AM, Nathan Lemoine
wrote:
> Thanks! I really appreciate the help, I don't think I ever would have figured
> that out.
>
> If you don't mind, I like to understand exactly whats going on in the code.
>
> I see groupS identifies the start center of each group and identifies which
> group each row belongs to
> w calculates the half-width of each bar with respect to the whole group of
> bars (i.e. a bar with a width of 1in a group of three bars has a half width
> of 1/6)
> groupS - (nv+1)/2 gives each bar's relative position to the central bar
>
> Am I interpreting these correctly? One thing I don't get is where 3/4 came
> from.
>
> Anyway, thanks for all of the help. Wouldn't have gotten it on my own
>
> Nathan
>
>
> On Mar 3, 2012, at 9:46 PM, ilai wrote:
>
>> On Sat, Mar 3, 2012 at 5:48 PM, Nathan Lemoine
>> wrote:
>>> It appears that the subscripts are only passing two values, the center of
>>> each group. There should be six values, one for the center of each bar
>>> (correct?),
>>
>> No. That's also why your code doesn't work. x[subscripts] are not the
>> centroids. To my knowledge there isn't a barchart equivalent to
>> boxplot.stats which means for placement you need to recalculate the
>> location of the bars, as is done internally in panel.barchart. This
>> panel function may need some tweaking depending on complexity of the
>> real data (for conditioning add [[packet.number()]] / for limits and
>> scales - a little more work) but hopefully help get you started:
>>
>> growth <- c(6.6,7.2,6.9,8.3,7.9,9.2,8.3,8.7,8.1,8.5,9.1,9.0)
>> diet <- as.factor(rep(c("A","B","C"),2,each=2))
>> coat <- as.factor(rep(c("light","dark"),each=6))
>> growth.means <- aggregate(growth,list(coat,diet),mean)
>> growth.errs <- aggregate(growth,list(coat,diet),sd)
>> require(lattice)
>> panel.ci <- function(x, y, sdl, groups, subscripts,box.ratio=box.ratio, ...){
>> groupS <- as.numeric(groups[subscripts])
>> nv <- nlevels(groups)
>> w <- box.ratio/(nv*(1 + box.ratio))
>> yy <- as.numeric(y)+sdl
>> xx <- as.numeric(x)+(3/4)*w*(groupS - (nv + 1)/2)
>> panel.barchart(x, y,groups=groups,subscripts=subscripts,...)
>> panel.segments(xx,as.numeric(y),xx,yy,lwd=2)
>> }
>> barchart(x~Group.1, groups=Group.2,data=growth.means,sdl=growth.errs$x,
>> horiz=F,stack=F,
>> ylim=c(0,max(growth.errs$x+growth.means$x)+.2), panel=panel.ci)
>>
>> Cheers
>> Elai
>>
>>
>> but I have no idea how to fix that or what I've done wrong. Here's the
>> (corrected) code I've got so far...
>>>
>>> growth <- c(6.6,7.2,6.9,8.3,7.9,9.2,8.3,8.7,8.1,8.5,9.1,9.0)
>>> diet <- as.factor(rep(c("A","B","C"),2,each=2))
>>> coat <- as.factor(rep(c("light","dark"),each=6))
>>>
>>> growth.means <- aggregate(growth,list(coat,diet),mean)
>>> colnames(growth.means)[3] <- "growth"
>>>
>>> library(lattice)
>>>
>>> panel.ci <- function(x, y, subscripts, ...){
>>> panel.barchart(x, y, horiz=F, subscripts=subscripts,
>>> groups=growth.means$Group.2, ...)
>>> panel.segments(x[subscripts], y, x[subscripts], y+0.5)
>>> print(x[subscripts])
>>> }
>>>
>>> barchart(growth~Group.1, groups=Group.2, data=growth.means,
>>> col=c(1,2,3),
>>> panel=panel.superpo
Re: [R] Grouped barchart confidence intervals in lattice
On Sat, Mar 3, 2012 at 5:48 PM, Nathan Lemoine wrote:
> It appears that the subscripts are only passing two values, the center of
> each group. There should be six values, one for the center of each bar
> (correct?),
No. That's also why your code doesn't work. x[subscripts] are not the
centroids. To my knowledge there isn't a barchart equivalent to
boxplot.stats which means for placement you need to recalculate the
location of the bars, as is done internally in panel.barchart. This
panel function may need some tweaking depending on complexity of the
real data (for conditioning add [[packet.number()]] / for limits and
scales - a little more work) but hopefully help get you started:
growth <- c(6.6,7.2,6.9,8.3,7.9,9.2,8.3,8.7,8.1,8.5,9.1,9.0)
diet <- as.factor(rep(c("A","B","C"),2,each=2))
coat <- as.factor(rep(c("light","dark"),each=6))
growth.means <- aggregate(growth,list(coat,diet),mean)
growth.errs <- aggregate(growth,list(coat,diet),sd)
require(lattice)
panel.ci <- function(x, y, sdl, groups, subscripts,box.ratio=box.ratio, ...){
groupS <- as.numeric(groups[subscripts])
nv <- nlevels(groups)
w <- box.ratio/(nv*(1 + box.ratio))
yy <- as.numeric(y)+sdl
xx <- as.numeric(x)+(3/4)*w*(groupS - (nv + 1)/2)
panel.barchart(x, y,groups=groups,subscripts=subscripts,...)
panel.segments(xx,as.numeric(y),xx,yy,lwd=2)
}
barchart(x~Group.1, groups=Group.2,data=growth.means,sdl=growth.errs$x,
horiz=F,stack=F,
ylim=c(0,max(growth.errs$x+growth.means$x)+.2), panel=panel.ci)
Cheers
Elai
but I have no idea how to fix that or what I've done wrong. Here's the
(corrected) code I've got so far...
>
> growth <- c(6.6,7.2,6.9,8.3,7.9,9.2,8.3,8.7,8.1,8.5,9.1,9.0)
> diet <- as.factor(rep(c("A","B","C"),2,each=2))
> coat <- as.factor(rep(c("light","dark"),each=6))
>
> growth.means <- aggregate(growth,list(coat,diet),mean)
> colnames(growth.means)[3] <- "growth"
>
> library(lattice)
>
> panel.ci <- function(x, y, subscripts, ...){
> panel.barchart(x, y, horiz=F, subscripts=subscripts,
> groups=growth.means$Group.2, ...)
> panel.segments(x[subscripts], y, x[subscripts], y+0.5)
> print(x[subscripts])
> }
>
> barchart(growth~Group.1, groups=Group.2, data=growth.means,
> col=c(1,2,3),
> panel=panel.superpose,
> panel.groups=panel.ci
> )
>
> Begin forwarded message:
>
>> From: Nathan Lemoine
>> Date: March 2, 2012 11:53:24 PM EST
>> To: r-help@r-project.org
>> Subject: Grouped barchart confidence intervals in lattice
>>
>> Hi everyone,
>>
>> I'm having trouble adding error bars to a grouped barchart in lattice. I
>> know that this topic has been addressed quite a bit, as I've been searching
>> the internet for a while to try to troubleshoot the issue, but I've not been
>> able to find any solution that I could get working on my data. I was
>> wondering if someone could look at my code and tell me what I'm doing wrong.
>> I was hoping somebody's found a way to do this (I'm sure they have) and can
>> tell me how to fix my code.
>>
>> # Input example data
>>
>> growth <- c(6.6,7.2,6.9,8.3,7.9,9.2,8.3,8.7,8.1,8.5,9.1,9.0)
>> diet <- as.factor(rep(c("A","B","C"),2,each=2))
>> coat <- as.factor(rep(c("light","dark"),each=6))
>>
>> growth.means <- aggregate(growth,list(coat,diet),mean)
>>
>> library(plotrix)
>>
>> growth.errs <- aggregate(growth,list(coat,diet),std.error)
>>
>> # Try using the superpose call with panel.groups results in an error
>>
>> panel.ci <- function(x, y, subscripts, groups...){
>> panel.barchart(x, y, groups=groups, subscripts=subscripts,
>> horiz=F,...)
>> panel.segments(x[subscripts], y, x[subscripts], y+growth.errs$x, col
>> = 'black')
>> }
>>
>> barchart(growth~Group.1, groups=Group.2, data=growth.means,
>> panel=panel.superpose,
>> panel.groups=panel.ci
>> )
>>
>> # Try using the generic plot.barchart command gives three error bars, all
>> are the appropriate sizes, but all are centered in each group and not on the
>> grouped bars
>>
>> barchart(x~Group.1, groups=Group.2, data=growth.means,
>> panel=function(x,y,subscripts, groups){
>> panel.barchart(x,y,horiz=F,groups=groups, subscripts=subscripts)
>>
>> panel.segments(as.numeric(x)[subscripts],y,as.numeric(x)[subscripts],y+growth.errs$x)
>> }
>> )
>>
>> What am I doing wrong?
>>
>> Thanks,
>>
>> Nathan Lemoine
>>
>>
>>
>>
>>
>>
>>
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.htm
Re: [R] Computing line= for mtext
On Fri, Mar 2, 2012 at 1:17 PM, Frank Harrell wrote:
> Hi Rich and Peter,
>
> What I am trying to do is the right-justify a vector of numbers to the right
> of the y-axis so that the leftmost digit of all of the numbers is one
> character to the right of the axis line. axis() plots tick marks and
No it doesn't (not always)
> left-justifies the numbers.
For axis(4). If you really want right justify on the right hand side,
you could use axis(2) with negative line numbers (device size
dependent). Or you could try something like
plot(1:20)
axis(2,at=seq(1,20,4),labels=T,tick=F,las=T,pos=c(22.25,1))
Hope this is getting there (after only 6 messages...)
Elai
>
> Peter's idea:
> -
> Since you're setting your right margin to 5, why
> not just
>
> mtext(s, side=4, las=1, at=5, adj=1, line = 5, cex=1)
> mtext(s, side=4, las=1, at=7, adj=1, line = 5, cex=2)
>
> i.e. set the line argument to the right margin?
> -
> comes close to what I'm trying to do but I haven't found out how to compute
> "line" from the maximum width over all the strings in the vector being
> plotted vertically, and I'm not sure it scales properly for different cex=.
> strwidth(s, units='inches')/par('cin')[1] doesn't seem to be a complete
> solution.
>
> Thanks
> Frank
>
>
> Richard M. Heiberger wrote
>>
>> Frank,
>>
>> this is it. It uses Peter's idea.
>>
>> plot(1:10)
>> axis(side=2, 1:10, las=1, line=-31.5, lwd=0)
>> axis(side=4, 1:10, las=1, labels=FALSE)
>>
>> Rich
>>
>>
>> On Thu, Mar 1, 2012 at 6:52 PM, Frank Harrell wrote:
>>
>>> Rich's pointers deals with lattice/grid graphics. Does anyone have a
>>> solution for base graphics?
>>> Thanks
>>> Frank
>>>
>>> Richard M. Heiberger wrote
>>> >
>>> > Frank,
>>> >
>>> > This can be done directly with a variant of the panel.axis function.
>>> > See function panel.axis.right in the HH package. This was provided for
>>> me
>>> > by David Winsemius in response to my query on this list in October 2011
>>> > https://stat.ethz.ch/pipermail/r-help/2011-October/292806.html
>>> >
>>> > The email thread also includes comments by Deepayan Sarkar and Paul
>>> > Murrell.
>>> >
>>> > Rich
>>> >
>>> > On Wed, Feb 29, 2012 at 8:48 AM, Frank Harrell wrote:
>>> >
>>> >> I want to right-justify a vector of numbers in the right margin of a
>>> >> low-level plot. For this I need to compute the line parameter to give
>>> to
>>> >> mtext. Is this the correct scalable calculation?
>>> >>
>>> >> par(mar=c(4,3,1,5)); plot(1:20)
>>> >> s <- 'abcde'; w=strwidth(s, units='inches')/par('cin')[1]
>>> >> mtext(s, side=4, las=1, at=5, adj=1, line=w-.5, cex=1)
>>> >> mtext(s, side=4, las=1, at=7, adj=1, line=2*(w-.5), cex=2)
>>> >>
>>> >> Thanks
>>> >> Frank
>>> >>
>>> >> -
>>> >> Frank Harrell
>>> >> Department of Biostatistics, Vanderbilt University
>>> >> --
>>> >> View this message in context:
>>> >>
>>> http://r.789695.n4.nabble.com/Computing-line-for-mtext-tp4431554p4431554.html
>>> >> Sent from the R help mailing list archive at Nabble.com.
>>> >>
>>> >> __
>>> >> R-help@ mailing list
>>> >> https://stat.ethz.ch/mailman/listinfo/r-help
>>> >> PLEASE do read the posting guide
>>> >>
>>> http://www.R-project.org/posting-guide.html;
>>> ;
>>> >> and provide commented, minimal, self-contained, reproducible code.
>>> >>
>>> >
>>> > [[alternative HTML version deleted]]
>>> >
>>> > __
>>> > R-help@ mailing list
>>> > https://stat.ethz.ch/mailman/listinfo/r-help
>>> > PLEASE do read the posting guide
>>> >
>>> http://www.R-project.org/posting-guide.html;
>>> > and provide commented, minimal, self-contained, reproducible code.
>>> >
>>>
>>>
>>> -
>>> Frank Harrell
>>> Department of Biostatistics, Vanderbilt University
>>> --
>>> View this message in context:
>>> http://r.789695.n4.nabble.com/Computing-line-for-mtext-tp4431554p4436923.html
>>> Sent from the R help mailing list archive at Nabble.com.
>>>
>>> __
>>> R-help@ mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html;
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@ mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
> -
> Frank Harrell
> Department of Biostatistics, Vanderbilt University
> --
> View this message in context:
> http://r.789695.n4.nabbl
Re: [R] Parameterization of Inverse Wishart distribution available in MCMCpack and bayesm libraries
On Fri, Mar 2, 2012 at 1:22 AM, peter dalgaard wrote:
>
> Er, yes (scalar does not imply integer)
Dough! awkward... Sorry Shantanu.
I've added
cat('###\n # ',substr(fortunes::fortune(90)$quote,1,146),'\n ### \n')
To .First in my Rhelp directory.
Hope that helps (me).
>
> As a general matter:
>
> 1. This is the Open Source world, you can read the actual function and see
> what it does. It might even say so in the comments.
>
> 2. You can investigate empirically -- the moments are known as a function of
> the parameters (check Wikipedia), so how about simulating a few thousand
> matrices and looking at the means and variances.
>
> (I don't do IW on a daily basis, but AFAICT, the two parametrizations have
> roughly the same mean, but a factor of two between variances, so it should be
> fairly easy to spot whether it is one or the other.)
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd@cbs.dk Priv: pda...@gmail.com
>
>
>
>
>
>
>
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Parameterization of Inverse Wishart distribution available in MCMCpack and bayesm libraries
What do you make of the following from ?riwish " riwish(v, S) v: Degrees of freedom (scalar). " does a m/2 parameterization yield a scalar for, say, 3 dof ? On Thu, Mar 1, 2012 at 9:55 AM, Shantanu MULLICK wrote: > Hello Everyone > > Both the MCMCpack and the bayesm libraries allow us to make draws from the > Inverse Wishart distribution. > > But I wanted to find out how exactly is the Inverse Wishart distribution > parameterized in these libraries. > > The reason I ask is the following: > > Now its generally standard to express Inverse Wishart as IW(0.5 * DOF,0.5* > Scale). (DOF-> Degree of freedom, Scale -> Scale parameter). > > If we follow standard usage when we refer to the Degree of Freedom of the > above IW distribution it is = DOF (and not 0.5* DOF). > > Similarly the Scale parameters of the above IW it is= Scale (and not > 0.5*Scale). > > For the MCMCpack the IW draws are made by riwish(v,S). > > *Question:* Does this compute IW( v,S) or IW(0.5*v,0.5*S) ? > > This is the reason I want to find out the way these libraries parameterize > the Inverted Wishart distribution. > > Best > Shantanu > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to colorize the panel backgrounds of pairs()?
par('bg') is not what you are looking for - it will set the bg of the
whole graphic device, not panels. I think you want:
count <- 0
mypanel <- function(x, y, ...){
count <<- count+1
ll<- par('usr')
if(count %in% c(1,4,9,12)) bg<- "#FDFF65"
else bg<- 'transparent'
rect(ll[1],ll[3],ll[2],ll[4],col=bg)
points(x, y, cex=0.5)
}
Cheers
On Thu, Mar 1, 2012 at 4:49 PM, Marius Hofert
wrote:
> Dear expeRts,
>
> I would like to colorize the backgrounds of a pairs plot according to the
> respective panel number. Here is what I tried (without success):
>
> count <- 0
> mypanel <- function(x, y, ...){
> count <<- count+1
> bg. <- if(count %in% c(1,4,9,12)) "#FDFF65" else NA
> points(x, y, cex=0.5, bg=bg)
> }
>
> U <- matrix(runif(4*500), ncol=4)
> pairs(U, panel=mypanel)
>
> I also tried to set par(bg=bg.) before the call to points(), but that didn't
> work either. The only thing I found is that "bg=" can be used to fill certain
> plot symbols, but I would like to colorize the background of each panel, not
> the drawn circles.
>
> Cheers,
>
> Marius
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] barplots of several variables with different number of categories
Jon, You could create new variables with the combined levels just for the purpose of plotting. Assume I have data.frame bpt str(bpt) 'data.frame': 12 obs. of 2 variables: $ V1: Factor w/ 3 levels "low","med","high": 1 1 1 1 2 2 2 2 3 3 ... $ V2: Factor w/ 6 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ... bpt$V3 <- factor(bpt$V1,c(levels(bpt$V1),levels(bpt$V2))) bpt$V4 <- factor(bpt$V2,c(levels(bpt$V1),levels(bpt$V2))) with(bpt,barplot(cbind(table(V3),table(V4 Hope this helps Elai On Thu, Mar 1, 2012 at 7:19 AM, jon waterhouse wrote: > If I have two factors, v1 and v2 and I want to have a stacked bar graph of > the two variables side by side I could do > > barplot(cbind(table(v1),table(v2))) > > if v1 and v2 have the same number of categories. > > If they don't have the same number of categories this won't work. > > I'm sure there's a simple solution? > > Thanks, > > Jon > > -- > View this message in context: > http://r.789695.n4.nabble.com/barplots-of-several-variables-with-different-number-of-categories-tp4435092p4435092.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Quantile scores as dependent variables.. an R and general method question
On Wed, Feb 29, 2012 at 1:23 PM, Doran, Harold wrote: > > The OP is looking for a way to deal with outcomes scores that are quantiles, > not a method that models different quantiles of the conditional distribution > where an outcome is a continuous variable. So, I don't think QR methods is > what is needed. Huhh? First, "deal with" is rather vague, does it mean summary? plot? delete? as the OP is worried about independence, some modelling exercise is more likely. Second, where do you think the OP quantiles come from if not from the conditional distribution of some continuous RV ? hint: not France. So for "...leads to theory, texts or R code.." on the distributional properties of quantiles/order statistics, do you have a better suggestion for a starting point than QR methods? Cheers > > From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] On Behalf > Of ilai [ke...@math.montana.edu] > Sent: Wednesday, February 29, 2012 1:30 PM > To: Rob James > Cc: r-help@r-project.org > Subject: Re: [R] Quantile scores as dependent variables.. an R and general > method question > > On Tue, Feb 28, 2012 at 3:54 PM, Rob James wrote: >> I have a dataset that does not include native scores, but only serial >> quantile rankings for a set of units. >> >> Clearly these observations are dependent (in that you can't alter one >> observation without also altering others). >> >> Are there methods for dealing with quantile dependent variables. My atempt >> to find such methods has not bee successful. >> > > Really? because google found 227k hits for "R quantile regression" - > none of them lead anywhere ? > > >> Any leads to theory, texts or R code would be most appeciated. >> >> [[alternative HTML version deleted]] >> >> __ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Quantile scores as dependent variables.. an R and general method question
On Tue, Feb 28, 2012 at 3:54 PM, Rob James wrote: > I have a dataset that does not include native scores, but only serial > quantile rankings for a set of units. > > Clearly these observations are dependent (in that you can't alter one > observation without also altering others). > > Are there methods for dealing with quantile dependent variables. My atempt > to find such methods has not bee successful. > Really? because google found 227k hits for "R quantile regression" - none of them lead anywhere ? > Any leads to theory, texts or R code would be most appeciated. > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] ts.plot and x axes customization
On Wed, Feb 29, 2012 at 1:56 AM, Jochem Schuster wrote:
> Hello,
>
> thank you very much for your answer. In the following, I will provide my
> recent code and try to explain again:
>
> series1 = ts(x$france start=c(2000,1), frequency=4)
> series2 = ts(x$germany, start=c(2000,1), frequency=4)
> time = c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011)
> ts.plot(series1, series2, col=c("blue", "red"), lwd=3, gpars=list(xaxt="n")
>
Thank you for an almost reproducible example (we don't have x and you are
missing the closing ')' in plot. No worries, the problem is here:
> axis(1, at=1:12, tick=1:48, labels=time, las=2)
>
> # series* is already of class 'ts', so in the plot the x-axis ranges >
par('usr')[1:2] and "at=1:12" is out of range. I think you wanted:
# x <- matrix(rt(12*4*2,3),nc=2,dimnames=list(NULL,c('france','germany')))
series1 = ts(x[,'france'], start=c(2000,1), frequency=4)
series2 = ts(x[,'germany'], start=c(2000,1), frequency=4)
ts.plot(series1, series2, col=c("blue", "red"), lwd=3, gpars=list(xaxt="n"))
axis(1, at=seq(2000-.25,2012,.25),labels=NA)
axis(1, at=2000:2012,labels=2000:2012,tcl= -.8,lwd.ticks=1.2,las=2)
Hope that did it.
Elai
> There is no error with entering these codes. After the ts.plot command,
> the x axes in the plot is deleted.
> But entering the axis command, the new axes - which should display annual
> labels and quarterly ticks - is not plotted.
>
> Thanks again for your help!
>
> Jochem
>
>
>
> Ihr WEB.DE Postfach immer dabei: die kostenlose WEB.DE Mail App für
> iPhone und Android.
> *https://produkte.web.de/freemail_mobile_startseite/*
Re: [R] ts.plot and x axes customization
On Tue, Feb 28, 2012 at 12:57 PM, Jochem Schuster
wrote:
>
ts.plot(series1, series2, main=main, xlab=xlab, ylab=ylab, col=c("green",
> "red", "blue"), lwd=2)
What I've tried before is deleting the X axes via gpars=list(xaxt="n") in
> the ts.plot-code. But after that I was not aible to add the customized axes
> via axis()...
So what was the problem? A reproducible example with your attempt at
annotating the ?axis would help.
> In advance, thank you very much for your help and hints!
> Regards,
> Jochem
>
>
> Ihr WEB.DE Postfach immer dabei: die kostenlose WEB.DE Mail App für iPhone
> und Android.
> [1]https://produkte.web.de/freemail_mobile_startseite/
>
> References
>
> 1. https://produkte.web.de/freemail_mobile_startseite/
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to replace the values in a column
Hannah,
If Gen is a factor you can simply build the new factor "on top" of it:
dataframe$Gen<- factor( c('Wynda' , 'A_2' , 'B_1' , 'Wynda' , 'Wynda'
, 'OP1_5')[Gen] )
Just make sure the replacement labels are in the same order as levels(Gen).
Cheers
On Tue, Feb 28, 2012 at 8:39 PM, hannahmaohuang
wrote:
> Dear All,
> I've been searching relevant topics about replacing values, none seemed to
> be applicable to me...
>
> I have a file with many many varieties, and want to replace some of them
> into different names.
> I tried various of ways, still don't know how to do that most efficiently..
> Here is part of the example data:
>
>
> Gen Rep
> A_1 1
> A_1 2
> A_2 1
> A_2 2
> B_1 1
> B_1 2
> B_3 1
> B_3 2
> OP1_1 1
> OP1_1 2
> OP1_5 1
> OP1_5 2
>
> For example, I want to replace A_1, B_3, OP1_1 into different name
> "Wynda"
>
> So that the expected file should become:
>
> Gen Rep
> Wynda 1
> Wynda 2
> A_2 1
> A_2 2
> B_1 1
> B_1 2
> Wynda 1
> Wynda 2
> Wynda 1
> Wynda 2
> OP1_5 1
> OP1_5 2
>
>
> I have created a link file, which contains two rows, translating which Gen
> correlating to which Name. Not sure if this file helps or not, example as
> below:
>
> Column1(Gen) Column2(Name)
> A_1 Wynda
> A_2 A_2
> B_1 B_1
> B_3 Wynda
> OP1_1 Wynda
> OP1_5 OP1_5
>
>
> Though I can replace one by one in excel, since there are too many files and
> too many reps, it'll be very time-consuming also easy to make mistakes.
>
> Please give me any guidance or help in terms of finish this with R.
>
> Thanks so much !
>
> Sincerely
> Hannah
>
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/How-to-replace-the-values-in-a-column-tp4430448p4430448.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] group calculations with other columns for the ride
aggregate(val~lvls+nm,data=x,FUN='median')
On Tue, Feb 28, 2012 at 4:43 PM, Ben quant wrote:
> Hello,
>
> I can get the median for each factor, but I'd like another column to go
> with each factor. The nm column is a long name for the lvls column. So
> unique work except for the order can get messed up.
>
> Example:
> x =
> data.frame(val=1:10,lvls=c('cat2',rep("cat1",4),rep("cat2",4),'cat1'),nm=c('longname2',rep("longname1",4),rep("longname2",4),'longname1'))
> x
> val lvls nm
> 1 1 cat2 longname2
> 2 2 cat1 longname1
> 3 3 cat1 longname1
> 4 4 cat1 longname1
> 5 5 cat1 longname1
> 6 6 cat2 longname2
> 7 7 cat2 longname2
> 8 8 cat2 longname2
> 9 9 cat2 longname2
> 10 10 cat1 longname1
>
> unique doesn't work in data.frame:
> mdn = do.call(rbind,lapply(split(x[,1], x[,2]), median))
> data.frame(mdn,ln=as.character(unique(x[,3])))
> mdn ln
> cat1 4 longname2
> cat2 7 longname1
>
> I want:
> mdn ln
> cat1 4 longname1
> cat2 7 longname2
>
> Thank you very much!
>
> PS - looking for simple'ish solutions. I know I can do it with loops and
> merges, but is there an option I am not using here?
>
> Ben
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Need advice on GLM
On Mon, Feb 27, 2012 at 1:44 AM, Christofer Bogaso wrote: > Here I was expecting those 2 approaches should give exactly same result > (i.e. same estimates and same SE), which is not the case. Can somebody > point me what I am missing here? > The vector of weights as described in ?glm which from your example is clear you either didn't look at or completely misinterpreted. summary(glm(YY~XX, binomial(link = "identity"))) summary(glm(Dat~Scores, binomial(link = "identity"))) summary(glm(Proportion~Scores, weights=rowSums(Dat),binomial(link = "identity"))) > Thanks and regards, > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Matrix problem to extract animal associations
set.seed(1)
(DFid <- data.frame(
x = sample(1:20,10),
y = sample(1:20,10),
IDs = sapply(1:10,function(i) paste("ID",i,sep=""
require(spdep)
coordinates(DFid) <- ~x+y
coords <- coordinates(DFid)
dnn4 <- dnearneigh(DFid,0,4)
summary(dnn4)
plot(DFid)
plot(dnn4,coords,add=T,col=2)
nb2mat(dnn4, zero.policy=TRUE)
This just one option from the multitude of spatial packages.
HTH
On Sun, Feb 26, 2012 at 4:55 PM, Ross Dwyer wrote:
> Dear List,
>
> I have been trying to extract associations from a matrix whereby individual
> locations are within a certain distance threshold from one another.
>
> I have been able to extract those individuals where there is 'no interaction'
> (i.e. where these individuals are not within a specified distance threshold
> from another individual) and give these individuals a unique Group ID
> containing that one individual.
>
> i.e.
>
> ID Group
>
> 1 ID1 1
>
> 2 ID3 2
>
> 3 ID4 3
>
> 4 ID5 4
>
> 5 ID7 5
>
> 6 ID8 6
>
> 7 ID9 7
>
>
> What I need assistance with is allocating associations with a unique group id.
> i.e. If we have interactions between "ID2_ID6", "ID6_ID2", "ID6_ID10",
> "ID10_ID6" as in the example code...
>
>
> ID Group
>
> 1 ID1 1
>
> 2 ID3 2
>
> 3 ID4 3
>
> 4 ID5 4
>
> 5 ID7 5
>
> 6 ID8 6
>
> 7 ID9 7
>
> ##
> 8 ID2 8
> 9 ID6 8
> 10 ID10 8
>
> ##
> The code also needs to robust enough to recognize instances where we have an
> interaction in a separate group...
> i.e. "ID11_ID12" should be in a separate group (Group 9) as they don't
> interact with IDs 2, 6, or 10 (not in below code!)
> 11 ID11 9
> 12 ID12 9
>
>
> I've been trying to figure this out but have drawn a blank. My example code
> can be found below.
>
> Very best wishes,
>
> Ross
>
> Dr Ross Dwyer
> Postdoctoral Research Fellow
> University of Queensland
>
>
>
>> ###
>> require(stats)
>> x <- sample(1:20,10)
>> y <- sample(1:20,10)
>> IDs <- sapply(1:10,function(i) paste("ID",i,sep=""))
>> (DFid <- data.frame(x,y))
> x y
> 1 7 20
> 2 5 3
> 3 12 5
> 4 3 12
> 5 18 19
> 6 2 1
> 7 19 15
> 8 20 11
> 9 13 14
> 10 1 2
>>
>>
>> (DMdist <- dist(DFid, method = "euclidean",
> + diag = FALSE, upper = TRUE))
> 1 2 3 4 5 6 7
> 8 9 10
> 1 17.117243 15.811388 8.944272 11.045361 19.646883 13.00
> 15.811388 8.485281 18.973666
> 2 17.117243 7.280110 9.219544 20.615528 3.605551 18.439089
> 17.00 13.601471 4.123106
> 3 15.811388 7.280110 11.401754 15.231546 10.770330 12.206556
> 10.00 9.055385 11.401754
> 4 8.944272 9.219544 11.401754 16.552945 11.045361 16.278821
> 17.029386 10.198039 10.198039
> 5 11.045361 20.615528 15.231546 16.552945 24.083189 4.123106
> 8.246211 7.071068 24.041631
> 6 19.646883 3.605551 10.770330 11.045361 24.083189 22.022716
> 20.591260 17.029386 1.414214
> 7 13.00 18.439089 12.206556 16.278821 4.123106 22.022716
> 4.123106 6.082763 22.203603
> 8 15.811388 17.00 10.00 17.029386 8.246211 20.591260 4.123106
> 7.615773 21.023796
> 9 8.485281 13.601471 9.055385 10.198039 7.071068 17.029386 6.082763
> 7.615773 16.970563
> 10 18.973666 4.123106 11.401754 10.198039 24.041631 1.414214 22.203603
> 21.023796 16.970563
>>
>> #Generate True/False matrix on those individuals < 4 units apart
>> DMTF <- apply(as.matrix(DMdist), c(1,2), function(x) ifelse(x<=4,T,F))
>> diag(DMTF)<- NA #replace diagonal with NA
>> dimnames(DMTF) <- list(IDs, IDs) #add individual's name to matrix
>>
>> DMTF
> ID1 ID2 ID3 ID4 ID5 ID6 ID7 ID8 ID9 ID10
> ID1 NA FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
> ID2 FALSE NA FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
> ID3 FALSE FALSE NA FALSE FALSE FALSE FALSE FALSE FALSE FALSE
> ID4 FALSE FALSE FALSE NA FALSE FALSE FALSE FALSE FALSE FALSE
> ID5 FALSE FALSE FALSE FALSE NA FALSE FALSE FALSE FALSE FALSE
> ID6 FALSE TRUE FALSE FALSE FALSE NA FALSE FALSE FALSE TRUE
> ID7 FALSE FALSE FALSE FALSE FALSE FALSE NA FALSE FALSE FALSE
> ID8 FALSE FALSE FALSE FALSE FALSE FALSE FALSE NA FALSE FALSE
> ID9 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE NA FALSE
> ID10 FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE NA
>>
>> irow <- as.character(gl(length(IDs),length(IDs),labels=IDs))
>> icol <-rep(IDs, length(IDs))
>>
>> AssocMatrix <- matrix(data=paste(irow,"_",icol,"_",sep=""),
> + nrow = length(IDs), ncol = length(IDs),
> + dimnames= list(IDs, IDs))
>>
>> AssocMatrix
>
> ID1 ID2 ID3 ID4 ID5 ID6
> ID7 ID8 ID9 ID10
>
> ID1 "ID1_ID1_" "ID2_ID1_" "ID3_ID1_" "ID4_ID1_" "ID5_ID1_" "ID6_ID1_"
> "ID7_ID1_" "ID8_ID1_" "ID9_ID1_" "ID1
Re: [R] Finding name of variable supplied as function argument
On Sat, Feb 25, 2012 at 12:53 PM, Ted Harding wrote:
> I have defined
> a function med3x3() such that, given vectors X,Y,
> med3x3(X,Y) returns a 3x3 table
...
> But I'd like to simply be able to pick up, within the function,
> the names of the variables that were used as arguments in the
> function call.
>
> The only suggestion .. uses a somewhat rebarbative parsing of the
> call stack.
> Is there a simpler way?
Is this a rhetorical question ? what could be simpler than
deparse(substitute(y)) when setting up the colnames ? If it's not
working for you for some reason then post your function. Or am I
missing something with the definition of "simple" ?
Or should I just use the 4-argument method?
You could, but even more efficient (and better control):
f <- function(x,y,rn=deparse(substitute(x)),cn=deparse(substitute(y))){
cat(rn,'=',x,'\t',cn,'=',y,'\n') }
> f(5,7)
5 = 57 = 7
> g <- 5 ; gs <- 7
> f(g,gs)
g = 5gs = 7
> f(g,gs,'g','gs2')
g = 5gs2 = 7
HTH
Elai
> With thanks,
> Ted.
>
> -
> E-Mail: (Ted Harding)
> Date: 25-Feb-2012 Time: 19:53:49
> This message was sent by XFMail
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Optim() package restriction
On Fri, Feb 24, 2012 at 4:03 PM, nserdar wrote:
> I did it like above but got an error message.
>
>> estimate<- optim(init.par,Linn,gr=NULL,method= "L-BFGS-B",
>> hessian=FALSE,control =
>> list(trace=1),lower=c(0,-Inf,Inf,Inf),upper=c(1,Inf,Inf,Inf))
Your lower bound for parameters 3,4 needs to be -Inf not Inf.
Otherwise this looks right, any other errors are your data/function.
> Error in solve.default(sig[, , 1]) :
> system is computationally singular: reciprocal condition number = 0
>
> "how to derive constrain for only 1 variable"
>
Adopted from the examples in optim:
fr <- function(x) { ## Rosenbrock Banana function
x1 <- x[1]
x2 <- x[2]
100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}
optim(c(-1.2,1), fr, method = "L-BFGS-B")
optim(c(-1.2,1), fr, method = "L-BFGS-B",upper=c(Inf,Inf))
optim(c(-1.2,1), fr, method = "L-BFGS-B",upper=c(Inf,.7))
HTH
> Regards,
> Ser
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Optim-package-restriction-tp4418379p4418865.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] tcl tk command function with arguments ??
On Fri, Feb 24, 2012 at 5:03 PM, peter dalgaard wrote:
>
>
> Now that's just weird... Firstly, it has nothing to do with sapply vs. for
> loops. It just works because you are inserting yet another function
> environment.
Thank you Peter, that makes more sense. As you can probably imagine,
first I made it work, then I convinced myself why. Just another of my
bad habits...
Secondly, if there was ever a point to having the buttons in the
OK.but list, it is not going to work anymore because you are assigning
each element to a different list, one per function environment.
Yes. I was aware of this but what is not going to work ? isn't the
OK.but list simply "storage" ? can it be used in some other way ? The
OP should answer if there is indeed some other application intended
that might fail, but I assumed there really wasn't a reason for this
exercise of creating the same button in a loop except for convenience.
> A more straightforward fix could be to do
>
> command = local({i <- i; function() PressedOK(i)})
Ouch! the thought to use local crossed my mind and was discarded in
favor of other (failed) attempts to "store" i before evaluation. This
is much better. Thank you again. I change my vote from my own solution
to this one.
>
> Other options might use the fact that command arguments can be R call
> objects, so you can insert the value of i directly into the call using
>
> command = substitute(PressedOK(i), list(i=i))
>
> or
>
> command = bquote(PressedOK(.(i)))
>
> [ There's a bug fixed in 2.14.2-RC which prevents this from working in a for
> loop; for 2.14.1 and earlier, you'll need .(force(i)) ]
>
>> Cheers,
>> Elai
>>
>>
>>> Thanks
>>>
>>>
>>> --
>>> View this message in context:
>>> http://r.789695.n4.nabble.com/tcl-tk-command-function-with-arguments-tp4416470p4416470.html
>>> Sent from the R help mailing list archive at Nabble.com.
>>>
>>> __
>>> R-help@r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd@cbs.dk Priv: pda...@gmail.com
>
>
>
>
>
>
>
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optim() package restriction
On Fri, Feb 24, 2012 at 12:18 PM, nserdar wrote: > Hi > > I need a "phi" restriction in my code. That is "0 > How can I do that ? > init.par<-c(1,1,1,1) > estimate<- optim(init.par,Linn,gr=NULL,method= "BFGS", hessian=FALSE,control > = list(trace=1)) > You want method "L-BFGS-B" not "BFGS". See details in ?optim Best > Regards, > res > > -- > View this message in context: > http://r.789695.n4.nabble.com/Optim-package-restriction-tp4418379p4418379.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] data frame manipulation with conditions
Ahh, I see it now.
For some reason your original post popped up on the list again, could
be just my mail server, sorry.
Looks like Uwe gave you the same solution (in two lines for better
clarity) right away. Depending on your level of "noobiness", my advice
would have been to ignore everything after that. Although if the other
approaches worked better for you, cheers.
Again sorry for this double thread.
On Fri, Feb 24, 2012 at 12:31 PM, Arnaud Gaboury
wrote:
> TY Elai for your answer. One solution has been given earlier in this list by
> Sarah Goslee and William Dunlap.
>
> Arnaud Gaboury
>
> A2CT2 Ltd.
>
>
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf Of ilai
> Sent: vendredi 24 février 2012 20:14
> To: A2CT2 Trading
> Cc: r-help@r-project.org
> Subject: Re: [R] data frame manipulation with conditions
>
> On Fri, Feb 24, 2012 at 8:11 AM, A2CT2 Trading wrote:
>> Dear list,
>>
>> n00b question, but still can't find any easy answer.
>>
>> Here is a df:
>>
>>> df<-data.frame(cbind(x=c("AA","BB","CC","AA"),y=1:4))
>
> # No, your y is a factor
> str(df)
> 'data.frame': 4 obs. of 2 variables:
> $ x: Factor w/ 3 levels "AA","BB","CC": 1 2 3 1 $ y: Factor w/ 4 levels
> "1","2","3","4": 1 2 3 4
>
> # You want to remove the cbind
>> df<-data.frame(x=c("AA","BB","CC","AA"),y=1:4)
>> str(df)
> 'data.frame': 4 obs. of 2 variables:
> $ x: Factor w/ 3 levels "AA","BB","CC": 1 2 3 1 $ y: int 1 2 3 4
>
>> I want to modify this df this way :
>> if df$x=="AA" then df$y=df$y*10
>> if df$x=="BB" then df$y=df$y*25
>>
>> and so on with other conditions.
>>
>
> df$y<- df$y * c(10,25,.5)[df$x]
> [1] 10.0 50.0 1.5 40.0
> # 1*10 2*25 3*.5 4*10
>
> HTH
>
> Elai
>
>> TY for any help.
>>
>> Trading
>>
>> A2CT2 Ltd.
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] data frame manipulation with conditions
On Fri, Feb 24, 2012 at 8:11 AM, A2CT2 Trading wrote:
> Dear list,
>
> n00b question, but still can't find any easy answer.
>
> Here is a df:
>
>> df<-data.frame(cbind(x=c("AA","BB","CC","AA"),y=1:4))
# No, your y is a factor
str(df)
'data.frame': 4 obs. of 2 variables:
$ x: Factor w/ 3 levels "AA","BB","CC": 1 2 3 1
$ y: Factor w/ 4 levels "1","2","3","4": 1 2 3 4
# You want to remove the cbind
> df<-data.frame(x=c("AA","BB","CC","AA"),y=1:4)
> str(df)
'data.frame': 4 obs. of 2 variables:
$ x: Factor w/ 3 levels "AA","BB","CC": 1 2 3 1
$ y: int 1 2 3 4
> I want to modify this df this way :
> if df$x=="AA" then df$y=df$y*10
> if df$x=="BB" then df$y=df$y*25
>
> and so on with other conditions.
>
df$y<- df$y * c(10,25,.5)[df$x]
[1] 10.0 50.0 1.5 40.0
# 1*10 2*25 3*.5 4*10
HTH
Elai
> TY for any help.
>
> Trading
>
> A2CT2 Ltd.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] tcl tk command function with arguments ??
On Fri, Feb 24, 2012 at 12:58 AM, Alexander wrote:
> I would like to know if its possible to use a function with arguments as a
> command in tcl tk.
Yes
I think
> this is due to the fact that the PressedOK(3) was the last call of the
> function, but I don't understand why all the other buttons have now a
> different command. Any idea?
>
Because in for loop i is overwritten every time - as you said i=3 is
the last command.
Replacing for with sapply should fix it:
require(tcltk)
OK.but <- NULL
PressedOK <- function(i)
{
tkmessageBox(message=paste("You pressed OK!",i,sep=""))
}
tt <- tktoplevel()
sapply(1:3,function(i){
OK.but[[i]] <- tkbutton(tt,text="OK",command=function()PressedOK(i))
tkgrid(OK.but[[i]])
})
tkfocus(tt)
Cheers,
Elai
> Thanks
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/tcl-tk-command-function-with-arguments-tp4416470p4416470.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] help with winbugs glm
My apologies to R-list members, this discussion diverged and is no
longer R related. Probably more fitting in one of the BUGS forums, but
I keep cc ing for any future interested reader that stumbled upon this
post in r-help.
Now to the point. Thank you so much for the reference, I was not aware
of this book. I will check it out as it seems can be useful for me
too.
However, compare McArthy's model to your own:
He is using Refuge[i] i=1,2,3 BUT ONLY AS A PLACE HOLDER for b[1:3] !!!
"... b[Refuge[i]]*Density[i] " where density is numeric - this is as
you said ANCOVA - fitting different slopes to Density BY levels of
Refuge. Not the same thing as your
"... b[Depth[i]] * Depth[i] "
Which is senseless as already discussed.
Side note, he can enter an explicit intercept "a" because it is the
estimate for the Med level, in the data he has dummy vars only for
Small and Large (e.g. in the first two columns of his design matrix,
lines 6 and 7 are 0 0 ). Thus, HIS model is identifiable.
Brad Efron (1986) said "Bayesian theory requires a great deal of
thought about the given situation to apply sensibly". And maybe this
is a good time to revisit another statement I think attributed to him:
"Recommending that scientists use Bayes’ theorem is like giving the
neighborhood kids the key to your F-16. "
Present company excluded :) as you seem to be more careful in your
application than some.
Cheers,
Elai
On Fri, Feb 24, 2012 at 8:46 AM, Adan Jordan-Garza
wrote:
> Ohh! about the reference code for the categorical predictor,
> from the McCarthy book, Bayesian Methods for Ecology (ISBN978-0-521-85057-5)
> the following ANCOVA model has a Refuge as a categorical predictor with 3
> levels coded as "1", "2" and "3".
>
> This is the cited model:
>
> model
> {
> a ~ dnorm(1.01, 24.826)
>
> # intercept
> b[1] ~ dnorm(0.628, 56.70) # effect of density at resource level 1 (low)
> b[2] ~ dnorm(0.488, 110.83) # effect of density at resource level 2 (medium)
> b[3] ~ dnorm(0.180, 70.144) # effect of density at resource level 3 (high)
> intnS[1] ~ dnorm(0, 1.0E-6) # interaction terms for small plots at 3 refuge
> densities
> intnS[2] ~ dnorm(0, 1.0E-6)
> intnS[3] ~ dnorm(0, 1.0E-6)
> intnL[1] ~ dnorm(0, 1.0E-6) # interaction terms for large plots at 3 refuge
> densities
> intnL[2] ~ dnorm(0, 1.0E-6)
> intnL[3] ~ dnorm(0, 1.0E-6)
> prec ~ dgamma(0.001, 0.001) # precision
> for (i in 1:25) # for each of the 25 plots
> {
> pred[i] <- a + b[Refuge[i]]*Density[i] +
> intnS[Refuge[i]]*Small[i]*Density[i] + intnL[Refuge[i]]*Large[i]*Density[i]
> Mortality[i] ~ dnorm(pred[i], prec)
> }
> }
> Initial values
> list(a=0, b=c(0,0,0), intnS=c(0,0,0), intnL=c(0,0,0), prec=1)
> Data
> Small[] Large[] Density[] Refuge[] Mortality[]
>
> 0 1 0.3125 1 1.6
>
> 1 0 0.75 1 0
>
> 0 1 0.75 1 2.22
>
> 1 0 1 1 1.19
>
> 1 0 1.25 1 5.19
>
> 0 0 1.5 1 1.73
>
> 0 0 1.9375 1 5.15
>
> 1 0 1 2 0.83
>
> 0 1 1.1875 2 1.44
>
> 0 0 1.375 2 1.5
>
> 0 0 1.4375 2 1.73
>
> 1 0 1.75 2 0
>
> 1 0 1.75 2 2.5
>
> 1 0 3.5 2 2.82
>
> 0 0 4 2 1.81
>
> 1 0 1.75 3 1.5
>
> 0 1 1.84375 3 1.33
>
> 1 0 2 3 1.66
>
> 1 0 2.25 3 0.73
>
> 0 0 2.5 3 1.6
>
> 0 1 3 3 1.99
>
> 1 0 3.25 3 1.09
>
> 1 0 4.25 3 1.54
>
> 1 0 5 3 1.6
>
> 1 0 5.5 3 1.73
>
> END
>
>
>
>
>
>
> On Fri, Feb 24, 2012 at 9:28 AM, Adan Jordan-Garza
> wrote:
>>
>> Ok this makes a lot of sense, thank you very much Ilai!
>> Cheers
>> Guillermo
>>
>> On Fri, Feb 24, 2012 at 12:16 AM, ilai wrote:
>>>
>>> On Thu, Feb 23, 2012 at 8:32 PM, Adan Jordan-Garza
>>> wrote:
>>> > Hello Ilai,
>>> > thank you very much for your response,
>>> > can I bother you a little further?
>>> > What do you mean my model is not identifiable?
>>>
>>> You should read up on model identifiably or consult your local
>>> statistician for a complete explanation, especially if you are to be
>>> doing more analysis. An incomplete explanation just regarding your
>>> bugs code: you are estimating "a" (an overall mean) and the three
>>> levels b[1:3]. That's too much, you need to remove "a". Think what you
>>> are asking of bugs to estimate (I'll make up some numbers): for bugs,
>>> the solution
>>> a = 10
>>> b[1] = 0
>>> b[2] = -2
>>> b[3] = -5
>>> is no different than
>>> a = 8
>>> b[1] = 2
>>> b[2] = 0
>>> b[3] = -3
>>> or
>>> a = 5
>>> b[1] = 5
>>> b[2] = 3
>&
Re: [R] Advice on exploration of sub-clusters in hierarchical dendrogram
Inline: On Thu, Feb 23, 2012 at 8:23 PM, R. Michael Weylandt wrote: > Inline: > > On Feb 23, 2012, at 6:20 PM, kosmo7 wrote: > >> Dear Elai, >> thank you very much for your suggestion. I tried cutting the dendrogram >> instead of the hclust tree with: >> clusters<-cut(x, h=1.6) >> >> but then when I try to call/plot cluster 1 for example, with: >> plot(clusters$lower[[1]]) >> >> I get only 2 members that are joined together at distance=0 (cluster 1 for >> instance, consists of several hundred of members). >> So it looks like / plot(clusters$lower[[1]])/ only calls the very first node >> of the tree and not the content of the respective cluster [[1]] at the >> defined cutoff=1.6. The "suggestions" in my original post are just pointers to the fact there are methods for class dendrogram to achieve what you wanted. Since you got as far as x<-as.dendrogram(z) I assumed that's all you needed. Maybe /cut/ instead of /cutree/ doesnt do the work? Or >> maybe I am just doing something wrong?... The examples in ?as.dendrogram and ?dendrapply are self contained, very clear and straight forward. If you haven't done so already I suggest you try them. Most likely the problem is in your data (row.names ? ) or your interpretation of who is "cluster1" or the 1.6 cutoff. >> >> >> >> In another post I read that with /df[value %in% v, ] / I can extract >> specific subsets of a data frame/table. > Seems I missed some back and forth on this post already, so my apologies if this is no longer an issue. Personally I find that because there are many more nodes and info in a tree than rows in the data set (leaf nodes only) much of the "usual" generic R solutions get distorted when it comes to trees. Better to use appropriate methods for the class (dendrapply helps as I've said before). Hope that helps dig you out of the hole. Elai > That was me and there's a slight mistake in that post (corrected by Sarah): > should be > > df[df$value %in% v, ] > > Sorry for any confusion that might have caused > > Michael >> Maybe I could use this to extract >> only the distances of members of a specific cluster as defined by cutree >> from the initial distance matrix? But still, I am afraid I don't get what I >> should use as /value/ and /v/ >> >> -- >> View this message in context: >> http://r.789695.n4.nabble.com/Advice-on-exploration-of-sub-clusters-in-hierarchical-dendrogram-tp4414277p4415589.html >> Sent from the R help mailing list archive at Nabble.com. >> >> __ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] cor() on sets of vectors
On Thu, Feb 23, 2012 at 3:24 PM, Bert Gunter wrote: > Use 1:n as an index. > > e.g. > sapply(1:n, function(i) cor(x[,i],y[,i])) ## sapply is a good solution (the only one I could think of too), but not always worth it: # for 100 x 1000 x <- data.frame(matrix(rnorm(10),nc=1000)) y <- data.frame(matrix(rnorm(10),nc=1000)) system.time(diag(cor(x,y))) # user system elapsed # 0.592 0.008 0.623 system.time(sapply(1:1000,function(i) cor(x[,i],y[,i]))) # user system elapsed # 0.384 0.000 0.412 # Great. but for 10 x 1000 x <- data.frame(matrix(rnorm(1),nc=1000)) y <- data.frame(matrix(rnorm(1),nc=1000)) system.time(diag(cor(x,y))) # user system elapsed # 0.256 0.008 0.279 system.time(sapply(1:1000,function(i) cor(x[,i],y[,i]))) # user system elapsed # 0.376 0.000 0.388 # or 100 x 100 system.time(diag(cor(x,y))) # user system elapsed # 0.016 0.000 0.014 system.time(sapply(1:100,function(i) cor(x[,i],y[,i]))) # user system elapsed # 0.036 0.000 0.036 # Not so great. Bottom line, as always, it depends. Cheers Elai > > -- Bert > > > > On Thu, Feb 23, 2012 at 2:10 PM, Sam Steingold wrote: >> suppose I have two sets of vectors: x1,x2,...,xN and y1,y2,...,yN. >> I want N correlations: cor(x1,y1), cor(x2,y2), ..., cor(xN,yN). >> my sets of vectors are arranged as data frames x & y (vector=column): >> >> x <- data.frame(a=rnorm(10),b=rnorm(10),c=rnorm(10)) >> y <- data.frame(d=rnorm(10),e=rnorm(10),f=rnorm(10)) >> >> cor(x,y) returns a _matrix_ of all pairwise correlations: >> >> cor(x,y) >> d e f >> a 0.2763696 -0.3523757 -0.373518870 >> b 0.5892742 -0.1969161 -0.007159589 >> c 0.3094301 0.997 -0.094970748 >> >> which is _not_ what I want. >> >> I want diag(cor(x,y)) but without the N^2 calculations. >> >> thanks. >> >> -- >> Sam Steingold (http://sds.podval.org/) on Ubuntu 11.10 (oneiric) X >> 11.0.11004000 >> http://www.childpsy.net/ http://iris.org.il http://americancensorship.org >> http://dhimmi.com http://www.PetitionOnline.com/tap12009/ >> http://jihadwatch.org >> Never argue with an idiot: he has more experience with idiotic arguments. >> >> __ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. > > > > -- > > Bert Gunter > Genentech Nonclinical Biostatistics > > Internal Contact Info: > Phone: 467-7374 > Website: > http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help with winbugs glm
Adan,
How many levels does Depth have? my wild guess: 3 and your bugs model
is not identifiable.
Second, I think you may have a critical error in the way you formatted
the data for the bugs model. From your code it looks like you are just
using the factor Depth and not a design matrix of dummy variables.
I may be wrong with respect to WinBugs (I use JAGS), but if Depth is
denoted as, for e.g., "low","med","high" wouldn't your multiply
operation "...*Depth[i] " on line 5 fail ?
More likely Depth is denoted "1","2","3" and WinBugs thinks it's
numerical. Well, in that case clearly coefficients for this model
don't make any sense (you'll only need one b for the slope). You can
use model.matrix(~Depth) to get the proper format for your data.
Hope this solves it. Next time, knowing n.chains n.iter and if they
achieved convergence (with different starting values) can help sort
through these sort of issues.
Cheers
Elai
On Thu, Feb 23, 2012 at 9:57 AM, Adan Jordan-Garza
wrote:
> Hi,
> I am running a model with count data and one categorical predictor (simple
> model for me to understand it fully), I did in R a glm like this:
> glm(Recruitment~Depth, family=poisson). I get the coefficientes and
> confidence intervals and all is ok. But then I want to do the same model
> with Bayesian stats, here is my code:
>
> model
> { for (i in 1:232)
> {
> Recruitment[i]~dpois(lambda[i])
> log(lambda[i])<-a+b[Depth[i]]*Depth[i]
> }
> a~dnorm(0,0.01)
> b[1]~dnorm(0,0.01)
> b[2]~dnorm(0,0.01)
> b[3]~dnorm(0,0.001)
> }
> list(a=0, b=c(0,0,0))
>
> I have two problems: 1) the resulting credible intervals for the
> coefficients (a, b1, b2 and b3) are HUGE don t make any reasonable sense;
> 2) Using OpenBugs and Winbugs I get different results,
>
> if anyone can help me I appreciate a lot your time,
>
> thanks
>
> Guillermo
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Advice on exploration of sub-clusters in hierarchical dendrogram
See inline
On Thu, Feb 23, 2012 at 8:54 AM, kosmo7 wrote:
> Dear R user,
> In other words, I am trying to obtain/read the sub-clusters of a specific
> cluster in the dendrogram, by isolating a specific node and exploring
> locally its lower hierarchy.
To explore or "zoom in" on elements of z you had the first step right:
create x<-as.dendrogram(z) but then you didn't use x anymore (except
for the plot which could have been done on z). Maybe you wanted:
> df=read.table('mydata.txt', head=T, row.names=1) #read file with distance
> matrix
> d=as.dist(df) #format table as distance matrix
> z<-hclust(d,method="complete", members=NULL)
> x<-as.dendrogram(z)
> plot(x, xlab="mydata complete-LINKAGE", ylim=c(0,4)) #visualization of the
> dendrogram
>From this point
clusters<-cut(x, h=1.6) #obtain clusters at cutoff height=1.6
# clusters is now (after cut x not cutree z) a list of two components:
upper and lower. Each is in itself a list of dendrograms: the
structure above 1.6, and the local clusters below:
plot(clusters$upper) # the structure above 1.6
plot(clusters$lower[[1]]) # cluster 1
# To print the details of cluster 1 (this output maybe very long
depending on how many members):
str(clusters$lower[[1]])
To extract specific details from the list and automate for all or some
of the clusters ?dendrapply is your friend.
I'm assuming your attempts at reclustering locally later in your post
are no longer necessary, unless I'm missing something on what exactly
you are trying to do.
Hope this helps
Elai
> ord<-cmdscale(d, k=2) #Multidimensional scaling of the data down to 2
> dimensions
> clusplot(ord,clusters, color=TRUE, shade=TRUE,labels=4, lines=0)
> #visualization of the clusters in 2D map
> var1<-var(clusters==1) #variance of cluster 1
>
> #extract cluster memberships:
> clids = as.data.frame(clusters)
> names(clids) = c("id")
> clids$cdr = row.names(clids)
> row.names(clids) = c(1:dim(clids)[1])
> clstructure = lapply(unique(clids$id), function(x){clids[clids$id ==
> x,'cdr']})
>
> clstructure[[1]] #get memberships of cluster 1
>
>
>
> >From this point, eventually, I could recreate a distance matrix with only
> the members of a specific cluster and then re-apply hierarchical clustering
> and start all over again.
> But this would take me ages to perform individually for hundred of clusters.
> So, I was hoping if anyone could point me to a direction as to how to take
> advantage of the initial dendrogram and focus on specific clusters from
> which to derive the sub-clusters at a new given cutoff height.
>
> I recently found in this page
> http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual
> http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual
>
> the following code:
> clid <- c(1,2)
> ysub <- y[names(mycl[mycl%in%clid]),]
> hrsub <- hclust(as.dist(1-cor(t(ysub), method="pearson")),
> method="complete") # Select sub-cluster number (here: clid=c(1,2)) and
> generate corresponding dendrogram.
>
> Even with this given example I am afraid I can't work my way around.
> So I guess in my case I could grab all the members of a specific cluster
> using my existing code and try to reformat the distance matrix in one that
> only contains the distances of those members:
> cluster1members<-clstructure[[1]]
>
> Then I need to reformat the distance matrix into a new one, say d1, which I
> can feed to a new -local- hierarchical clustering:
> hrsub<-hclust(d1, method="complete")
>
> Any ideas on how I can obtain a new distance matrix with just the distances
> of the members in that clusters, with names contained in vector
> "cluster1members" ?
>
> Apologies if this seems trivial, but I really can't find the correct
> functions to use for this task.
> Thank you very much in advance - as I am really a novice with R, small
> chunks of code as example would be of great help.
>
> Take care all -
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Advice-on-exploration-of-sub-clusters-in-hierarchical-dendrogram-tp4414277p4414277.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Lattice and horizontally stacked density plots
The plot you referred to depends on packages flowViz and flowCore from R-bioconductor. With lattice alone you can easily get all curves on the same level: densityplot(~ val | factor(id2), groups=factor(id1),data=a_df,pch='|') But if that doesn't do it for you, you could write your own panel function. I don't have time to try it but I'm thinking one of these might work 1) create new ylim[1] from current.panel.limits()$ylim / number of groups dens<- density(x) use lines and polygons to draw the curves dens$x,dens$y at each new ylim level. 2) grid may come in handy here, splitting each panel into several viewports ? 3) use bwplot for "setup" but plot polygons of density (steps 2:3 from option 1) Good luck with that. Elai On Wed, Feb 22, 2012 at 3:01 PM, Manish Nag wrote: > Hello, > > I am try to make a density plot where plots are stacked like the one > found here: > > http://dsarkar.fhcrc.org/lattice/book/images/Figure_14_03_stdBW.png > > I am facing problems, however. Using the code example below, I'd like > to generate a separate panel for each val of id2. Within each panel, > I'd like to have individual histograms each on separate lines based on > the value of id1. Note that the code example works fine if I use > "boxplot" instead of "densityplot". Any pointers would be much > appreciated. > > > library(lattice) > val<-rep(rnorm(10),100) > id1<-sample(c(1:5), 100, replace = TRUE) > id2<-rep(c(6:10),100, replace = TRUE) > a_df<-data.frame(cbind(id1, id2, val)) > densityplot(factor(id1) ~ val | factor(id2), data=a_df) > > -Manish > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] tapply for enormous (>2^31 row) matrices
On Tue, Feb 21, 2012 at 4:04 PM, Matthew Keller wrote:
> X <- read.big.matrix("file.loc.X",sep=" ",type="double")
> hap.indices <- bigsplit(X,1:2) #this runs for too long to be useful on
> these matrices
> #I was then going to use foreach loop to sum across the splits
> identified by bigsplit
How about just using foreach earlier in the process ? e.g. split
file.loc.X to (80) sub files and then run
read.big.matrix/bigsplit/sum inside %dopar%
If splitting X beforehand is a problem, you could also use ?scan to
read in different chunks of the file, something like (untested
obviously):
# for X a matrix 800x4
lineind<- seq(1,800,100) # create an index vec for the lines to read
ReducedX<- foreach(i = 1:8) %dopar%{
x <-
scan('file.loc.X',list(double(0),double(0),double(0),double(0)),skip=lineind[i],nlines=100)
... do your thing on x (aggregate/tapply etc.)
}
Hope this helped
Elai.
>
> SO - does anyone have ideas on how to deal with this problem - i.e.,
> how to use a tapply() like function on an enormous matrix? This isn't
> necessarily a bigtabulate question (although if I screwed up using
> bigsplit, let me know). If another package (e.g., an SQL package) can
> do something like this efficiently, I'd like to hear about it and your
> experiences using it.
>
> Thank you in advance,
>
> Matt
>
>
>
> --
> Matthew C Keller
> Asst. Professor of Psychology
> University of Colorado at Boulder
> www.matthewckeller.com
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] line width in legend of interaction.plot
On Wed, Feb 22, 2012 at 11:37 AM, Matthias Gondan
wrote:
> Dear R developers,
>
> The following command produces an interaction plot with lwd=2.
>
> interaction.plot(c(1, 2, 1, 2), c(1, 1, 2, 2), 1:4, lwd=2)
>
> In the legend, however, lwd seems to be 1, which does not seem
> to be intended behavior.
Why not? Your suggestion will cause something like this to be the default
interaction.plot(c(1, 2, 1, 2), c(1, 1, 2, 2), 1:4, lwd=4,legend=F)
legend('topleft',c('1','2'),lty=1:2,lwd=4)
Which seems worse to me.
Cheers
Probably the lwd is not correctly forwarded
> to legend:
>
> from the interaction.plot source:
>
> legend(xleg, yleg, legend = ylabs, col = col, pch = if (type %in%
> c("p", "b"))
> pch, lty = if (type %in% c("l", "b"))
> lty, bty = leg.bty, bg = leg.bg) <- here I would add lwd= the lwd from the ... argument, or perhaps something like leg.lwd>
>
> Best wishes,
>
> Matthias Gondan
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Installing rgl
On Wed, Feb 22, 2012 at 3:23 AM, aoife wrote: > Hey guys,I'm working with R version 2.14.1 (2011-12-22) on a unix machine You may be missing some openGL libraries (mesa ? ) On Ubuntu, this solved my problem of installing rgl: sudo apt-get build-dep r-cran-rgl Can't vouch for any future complications with rendering in other applications...so as always, use at your own risk. Hope this helps Elai and > I'm having a similar problem. I have a matrix and i want to do a simple > correspondance analysis. I tried to install the ca package i get this error: > >> library(ca) > Loading required package: rgl > Warning messages: > 1: In rgl.init(initValue) : RGL: no suitable visual available > 2: In fun(libname, pkgname) : error in rgl_init >> > > so then i downloaded rgl, at the end of the installation process for rgl it > says: > > Warning messages: > 1: In rgl.init(initValue) : RGL: no suitable visual available > 2: In fun(libname, pkgname) : error in rgl_init > > > and when you look in the library for rgl you get: >> library(rgl) >> > > > Would anyone know (in simple language, i'm a biologist, not a statistician > :( ) either how to fix this, or an alternative package that will do a simple > correpondance analysis on a matrix? > > Many thanks > Aoife > > -- > View this message in context: > http://r.789695.n4.nabble.com/Installing-rgl-tp865867p4409845.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Several densityplots in single figure
On Wed, Feb 22, 2012 at 8:49 AM, David Winsemius wrote:
>
> After going back and constructing a proper dataset, you should be passing
> 'groups' into the panel function and picking it up inside panel.abline.
Close, but unfortunately things get more complicated when using groups
in densityplot. A straight up panel function wouldn't work for him out
of the box (and you/list be getting more follow ups). He'll need
something like:
densityplot(~value | g1, groups=g2, data= combinedlongformat,
panel=panel.superpose, panel.groups=function(x,...,group.number){
panel.abline(v=mean(x),lwd=(1:0)[group.number]) # seems like he only wants x1
panel.densityplot(x,...)
})
At this stage in the game
trellis.focus('panel',1,1) ; panel.abline(v=calculatedmean1)
trellis.focus('panel',2,1) ; panel.abline(v=calculatedmean2)
...
trellis.unfocus()
Would probably go down easier...
Just a thought
Cheers,
Elai
You
> can also recover the 'panel.number()' and then use them both to recover the
> means from the global environment when inside panel.abline. I probably would
> have used tapply to create a table with an index (in the enclosing
> enviroment) rather than creating 4 separate values to choose That way you
> could use one object name but construct the 4 indices. But for that to work
> you would have needed to follow my example about binding the 4 datasets into
> one.
>
>
> One way to find previous worked examples in the rhelp Archives is to search
> with your favorite engine on strategies like "panel.number groups Sarkar" or
> "panel.number groups Andrews"
>
> Deepayan Sarkar and Felix Andrews are the two persons from whom I have
> learned the most regarding the fine points of lattice plots.
>
> And will you please learn to post in plain text?
>
> --
> david
>
>
>> }
>> )
>>
>>
>>
>>
>>
>>
>> On 22 February 2012 13:48, David Winsemius wrote:
>>
>> On Feb 22, 2012, at 5:28 AM, josh rosen wrote:
>>
>> Hi,
>>
>> I have created two separate overlapping density plots- see example code
>> below.
>> What I wish now to do is combine them into one figure where they sit side
>> by side.
>> Any help would be great!
>>
>> many thanks in advance, josh.
>>
>> #
>> thedataA <- data.frame(x1=rnorm(100,1,1),x2=rnorm(100,3,1)) #create data
>> thedataA.m<-melt(thedataA)
>>
>> densityplot(~value, thedataA.m, groups=variable,auto.key=list(columns=2),
>>
>> panel = function(x, y, ...) {
>> panel.densityplot(x, ...)
>> panel.abline(v=0)
>> }
>> )
>>
>>
>> The syntax for grouping (which gives theoverlaid but different colors as
>> default output) and "separation" is fairly simple. Use the "|" operator for
>> separated plots and the " .., groups= , .." parameter for overlaying
>> results. It's only going to work easily if they are all in the same dataset.
>>
>> Try:
>>
>> bigset <- cbind( rbind(thedataA.m, thedataB.m), ABgrp=rep(c("datA",
>> "datB"), each=200) )
>> densityplot(~value|ABgrp, data=bigset, groups=variable,
>> auto.key=list(columns=2),
>>
>> panel = function(x, y, ...) {
>> panel.densityplot(x, ...)
>> panel.abline(v=0) } )
>>
>>
>> And please work on whatever practice is producing duplicate postings.
>>
>> --
>> David.
>>
>>
>>
>>
>> thedataB <- data.frame(x1=rnorm(100,2,1),x2=rnorm(100,4,1)) #create data
>>
>> thedataB.m<-melt(thedataA)
>>
>> I assume that is a copy-paste-fail-to-correct error.
>>
>>
>> densityplot(~value, thedataB.m, groups=variable,auto.key=list(columns=2),
>>
>> panel = function(x, y, ...) {
>> panel.densityplot(x, ...)
>> panel.abline(v=0)
>> }
>> )
>> ##
>>
>> [[alternative HTML version deleted]]
>
>
>
> David Winsemius, MD
> West Hartford, CT
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Time taken to process a file after a socket connection was made
Aniruddha, Thank you, this is much more clear (and makes sense). Unfortunately I am out of my element here so I don't think I'll be much help, but see inline On Tue, Feb 21, 2012 at 4:14 AM, Aniruddha Mukherjee < aniruddha.mukher...@tcs.com> wrote: > > > -- I ran the "socketConnection" command and the R session went into a > 'wait' state > -- In another window, I ran the Java code thru which the socket connection > was established and the 'R' prompt was back > -- Now in the meantime, the Java code has written a line onto the socket > whereby a letter, say '*d*', was sent > -- I ran the "readLines" command at 'R' prompt and after waiting for > *61*seconds, it displayed the ' > *d*'. > -- Again later, in a next trial, I sent a matrix of 20,000 values (in > 1,000 rows & 20 cols) from Java and that took almost the same time (*62*sec.) > by the 'readLine' command OK. So that is much clearer. I think, contrary to your statement in the original post, it is not readLines taking up the 60 sec, but something with the communication. Since you get an answer eventually, it must be listening, and if you opened the connection with default, it should not be blocking either (no timeout)... Like I said I don't have much to add but maybe try readLines(connection,n=1000) I wouldn't think it matters, but maybe explicitly declaring the number of lines to read will speed it up because R is waiting for more input and only after 60 sec it times out and reads what's available ? again, the default is blocking=FALSE, but... Second, make sure the connection is established (ready to read/write): socketSelect(list(conname,conname),c(F,T),timeout=0) Bottom line, I have no idea what is going on except 60 sec doesn't sound right to me... Good luck Elai > > Here lies my query : Is there any kind of "*timeout*" factor which is > causing the delay? > Please share your thought. > Also do you think that instead of *socketConnection & readLines*, *make.socket > & read.socket* would perform better? > > Thanks and regards. > Aniruddha. > > > From: ilai To: Aniruddha Mukherjee < > aniruddha.mukher...@tcs.com> Date: 02/21/2012 10:42 AM Subject: Re: [R] > Time taken to process a file after a socket connection was made Sent by: > ila...@gmail.com > -- > > > > That sounds very strange. Are you sure readLines consumed 63 sec for > 2000 vals ? doubt it. When you say "...started by an "R" session..." > does it imply the R session actually doing the work is already up and > running ? a shot in the dark, but maybe the 60 sec are the time for R > to load on your machine (even though that's still 10 times what I'd > expect from my dino)? Did you test the R part of it directly ? i.e. > could you post the result from an open R session or running a batch > script starting with readLines ? > see ?system.time or ?proc.time if you need them. > BTW hope this is just a test and you're not firing up R just to > calculate a mean... > > Cheers > > > On Mon, Feb 20, 2012 at 6:21 AM, Aniruddha Mukherjee > wrote: > > Hello R people, > > > > I have created a '.csv' file of 100 rows by 20 columns whose each cell > > contains a random numbers between 0 & 1, thru a Java program. Once that > is > > created a signal (just a letter) is send to the port of a socket > > connection at "localhost", which was earlier started by an "R" session. > > Now the "R" reads the '.csv' file into a data frame and calculates the > > average of 2000 numbers. This mean value was then written to the socket > > connection and Java received that successfully. The duration of this > > processing at "R" session took 1 min 03 seconds. > > My question : Is the time duration of 63 seconds OK for only this set of > > activities or it should be less? If so how it can be reduced? Point to be > > noted that both the Java and R were running on the same PC and the 62 > > seconds were consumed for the "readLines" command. > > > > Thanks and regards. > > Aniruddha. > > =-=-= > > Notice: The information contained in this e-mail > > message and/or attachments to it may contain > > confidential or privileged information. If you are > > not the intended recipient, any dissemination, use, > > review, distribution, printing or copying of the > > information contained in this e-mail message > > and/or attachments to it are strictly prohibited. If > > you have received this communication in error, > &g
[R] Fwd: bigmemory not really parallel
The default nstart in ?bigkmeans is 1. Try
ans<-bigkmeans(data,k,nstart=8)
Good luck
On Mon, Feb 20, 2012 at 4:22 PM, Lishu Liu wrote:
> Hi, all,
>
> I have a really big matrix that I want to run k-means on.
> I tried:
>>data <-
> read.big.memory('mydata.csv',type='double',backingfile='mydata.bin',descriptorfile='mydata.desc')
> I'm using doMC to register multicore.
>>library(doMC)
>>registerDoMC(cores=8)
>>ans<-bigkmeans(data,k)
>
> In system monitor, it seems only one thread running R. Is there anything I
> did wrong?
> Thanks in advance for any suggestions.
>
> Best,
> Lishu
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] barplot with more than 1 variable
You could have found the solution in
http://had.co.nz/ggplot2/geom_bar.html yourself since all help pages
for ggplot refer you to the web site. But to speed things up for you,
try this for starters:
TUSE2 <- data.frame(country = rep(c("United
States","Italy","Germany","Netherlands"),each=2),
year3 = factor(1:2),
mw = c( 245.8, 255.9, 248.5, 207.4,263.9, 197.7, 174.2, 189.5),
st = c( 200.5, 218.0, 236.1, 237.3, 220.5, 242.7, 221.0, 206.0),
all= c( 446.3, 473.9, 484.6, 444.7,484.5, 440.4, 395.2, 395.5))
# Note year3 is a factor
library(ggplot2)
forplot <- stack(TUSE2,c('st','mw','all'))
ggplot(TUSE2) +
geom_bar(aes(year3,forplot$values,fill=forplot$ind),position="dodge")+
facet_grid(. ~ country)
Hope this helps
Elai
On Sun, Feb 19, 2012 at 5:31 AM, Francesco Sarracino
wrote:
> Dear R listers,
> I am trying to produce a simple (for a stata user) barplot with 4
> countries on the x axis, each country observed in 2 subsequent years
> and 3 variables.
> Basically, I should have three bars for each year for each country. I
> am attaching the chart I made in Stata, but I am not sure you'll
> manage to see it!
>
> I did the following:
> #here I create the data-set TUSE2. The vectors mw, st and all are the
> three variables I'd like to plot in each of the two years of the
> variable year3 for each country.
>
> country <- c("United States","United
> States","Italy","Italy","Germany","Germany","Netherlands","Netherlands")
> year3 <- c(1, 2, 1, 2, 1, 2,1 , 2)
> mw <- c( 245.8, 255.9, 248.5, 207.4,263.9, 197.7, 174.2, 189.5)
> st <- c( 200.5, 218.0, 236.1, 237.3, 220.5, 242.7, 221.0, 206.0)
> all <- c( 446.3, 473.9, 484.6, 444.7,484.5, 440.4, 395.2, 395.5)
> TUSE2 <- data.frame(country, year3, mw, st, all)
> attach(TUSE2)
>
> #I load the library ggplot (But if you know of alternative techniques,
> I'll be happy to try them out)
> library(ggplot2)
>
> #I start producing the chart:
> p <- ggplot(TUSE2, aes(country, all, group=year3))
> p <- p + layer(
> geom="bar",
> position="dodge",
> stat="identity")
> #so far I manage to get what I want for just one variable: all. Let's
> try to add another variable, let's say: mw
>
> f <- p + layer(
> geom="bar",
> aes(x=country, y = mw),
> position="dodge",
> stat="identity")
> f
> #and here my adventure crashes. I get the following message and I
> can't figure out how to get my chart. The error message is:
> Error in pmin(y, 0) : object 'y' not found
> In addition: Warning message:
> In min(diff(sort(x))) : no non-missing arguments to min; returning Inf
>
> Can you please help me understanding how to make my chart? I am trying
> to switch from Stata to R, but so far life is really hard!
> Best,
> f.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] World Heat/Rainbow Map
Figure 13.10 in "Lattice: Multivariate Data Visualization with R" might give you some ideas http://lmdvr.r-forge.r-project.org/figures/figures.html On Sun, Feb 19, 2012 at 9:05 AM, Lorenzo Isella wrote: > Dear All, > I would like to do the following: make a plot of the world and color a few > selected states. > Some states have an associated scalar (i.e a simple number) and I would like > to paint them using a heat map and show legend for the color map in the plot. > One or two states do not have any number associated to them and are simply > colored/filled according to a distinctive pattern that I choose. > To fix the ideas, imagine the following situation in which > > (1) the US are striped (or filled with any pattern not in the heat map) and > no number is associated to them; > (2) China, France and UK have associated values 54, 32 and 14 and I would > like to "paint" them according to a heat map. > (3) The other states can be left blank, or I can pretend that there is a > "zero" associated to them. > > I made some attempts with map(world), but so far I have been banging my head > against the wall. > Any suggestion is appreciated. > Cheers > > Lorenzo > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] foreach %do% and %dopar%
Marcos,
Untested because you didn't provide a reproducible example but my
guess, the problem is in the local assignment of MCPVMP*. The %do%
worked just because it operates in the same (local) environment for
all threads. I'll try to clarify with a a self contained example of
sourcing a script with "calculations" one some variables (take the
sqrt):
cat('sqrt(somevar)',file='myscript.R')
library(doMC)
registerDoMC()
foreach(i = 1:2) %dopar% {somevar <- c(1,9)[i] ; source('myscript.R')$value }
# local assignment fails
foreach(i = 1:2) %dopar% {somevar <<- c(1,9)[i] ; source('myscript.R')$value }
# assign to the global environment works
Bottom line try to change all the equal signs to "<<-" and source the
script in the loop (I don't think passing it to .options.smp will work
in this context).
Two unrelated minor points
1) is NsimT a vector 1:2000 or of class 'iter' ? just a number like
simx=2000 is not right
2) With intel's hyperthreading you may really have only 4 real cores not 8
Hope this helps
Elai
On Fri, Feb 17, 2012 at 11:41 PM, Marcos Larios wrote:
> Hi everyone,
>
> I'm working on a script trying to use foreach %dopar% but without success,
> so I manage to run the code with foreach %do% and looks like this:
>
> The code is part of a MCMC model for projects valuation, returning the most
> important results (VPN, TIR, EVA, etc.) of the simulation.
>
> foreach (simx = NsimT, .combine=cbind, .inorder=FALSE, .verbose=TRUE) %do% {
> MCPVMPA = MCVAMPA[simx] #The *[simx] variables are vectors containing
> 100,000 simulations of each variable.
> MCPVMPB = MCVAMPB[simx] #Wich then I want to parse to the script below
> MCPVMPC = MCVAMPC[simx] #In order for the model to take the values of
> each variable.
> MCPVMPD = MCVAMPD[simx]
>
> source(paste(pathglobal, "Valuation.R", sep="")) #This script does
> everyting
> #Before you suggest making it a function, let me say I CAN'T.
> #This script is about 3000 lines long, and connects to several other
> scripts (8), each about 300-500 lines long.
> #pathglobal is a string variable pointing to the working directory, ej.
> ~/Documents/Valuation/
>
> MCVRERM[simx] <- sum(CIFERRN) #This is my response variable 1
> MCVRWKM[simx] <- WKMTWKZ[12] #This is my response variable 2
> MCVRFNM[simx] <- sum(FNMSDE) #This is my response variable 3
> #The response variables come from another script which is call within
> Valuation.R
>
> }
>
> As I said, the above code works flawlessly with %do%, but it takes about 1
> hour just to run 2,000 simulations, so I want to make use of all the 8
> cores of my Intel i7 machine by using %dopar%, so I changed the code to
> something like this:
>
> foreach (simx = NsimT, .combine=cbind, .inorder=FALSE, .verbose=TRUE,
> .options.smp= source(paste(pathglobal, "Valuation.R", sep=""))) %dopar% {
>
> MCPVMPA = MCVAMPA[simx]
> MCPVMPB = MCVAMPB[simx]
> MCPVMPC = MCVAMPC[simx]
> MCPVMPD = MCVAMPD[simx]
>
> MCVRERM[simx] <- sum(CIFERRN)
> MCVRWKM[simx] <- WKMTWKZ[12]
> MCVRFNM[simx] <- sum(FNMSDE)
> }
>
> The changes makes the script fail, saying somtehing like: abscent value
> where TRUE/FALSE is necessary, the variables that are supposed to pass the
> information to the Valuation.R script are empty, and error comes because it
> does something like this at some point:
>
> VarX = A/B #But because the variables values are not getting into the
> Valuation.R script then
> #A = MCPVMPA = 0 and B = MCMKMPA = 0
>
> So VarX = NaN
>
> Then it does something like:
>
> if (VarZ>VarX) VarY = 0 else VarY = VarX-VarZ #This line generates the error
>
> I have all the libraries installed:
> library(foreach)
> library(doMC)
> registerDoMC()
>
> It's not the installation because all the example code with %dopar% I've
> tested work fine.
> So can you help me please to make my code work?
> Any ideas?
>
>
> --
> L.E. Marcos Larios Sata Rosa
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Neighbour List to Matrix
I think is good to know that list contain more than 60 rows with
around 14000 nodes (participants).
?read.table may be unreliable for large matrices and with 14/600
you'll end up with many NA's. You might do better with
nbrs<- scan('nbrs.txt',skip=1,what=list('integer','integer',double(0)))
names(nbrs) <- c('c1','c2','c3')
xtabs(c3~c1+c2,nbrs,sparse=T)
# returns
4 x 9 sparse Matrix of class "dgCMatrix"
1 2 3 4 5 6 7 8 9
1 . 0.065 0.044 . . . . . .
3 0.071 . . 0.016 0.011 . . . .
4 . . 0.004 . . . 0.004 . 0.004
5 . . . 0.010 . 0.011 . 0.009 0.004
Cheers
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 - error message
I think my problem is that I can't >> incorporate the 'lake' variable in a fixed-effect interaction because it is >> only has one binary observation. But I don't know what to do to be able to >> fit this model. Any help would be greatly appreciated! >> -Sean > > In principle you should be able to fit this model, but the error message > is telling you that there are numeric problems -- it may just be that > your data are a little too sparse in some direction. Yes. Consider the collinearity between Age and Year, i.e. for a given cohort (mos or all captured by "fishID" ?) they are essentially the same variable with different units. So I would suspect the problem is you are over fitting those. A few suggestions: > > * try centering Age, or re-introducing the intercept, to see if you > can get something to work. > * You _might_ try the development version of lme4 (lme4Eigen, on r-forge) > * plot your data to see if you see anything odd about the data > * perhaps try making Year a fixed effect -- 4 levels is fairly small > for a random effect > * Ask further questions on the r-sig-mixed-models mailing list. > > Ben Bolker > Add one more to those: make sure your random effects are indeed crossed. If they are nested (without knowing anything about your data, just given their names that's a possibility), you could try nlme::lme. Elai > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] creating series of vectors
# All days in years 2006 to 2009 by month in 48 (12x4) files.
days <- seq(as.Date("2006/1/1"), as.Date("2009/12/31"),by="day") # one
long vector
out <- paste(rep(format(days,'%d%m%y'),each=2),c('aaa','bbb'),sep='_')
# reformat to style
month <- factor(rep(format(days,'%B%y'),each=2)) # group by month.year
for(i in levels(month))
cat(out[month==i],'\n',file=paste(i,'txt',sep='.')) # write external
files
Cheers
On Thu, Feb 16, 2012 at 9:32 AM, Nino Pierantonio wrote:
> Dear All,
>
> I am pretty new to R and thus my question may sound silly.
>
> Is there a way to automatically generate a series of separate vectors (so
> not arranged in a matrix), without typing and changing every time the
> values, and store them as separate *xlsx file, where the "*" is replaced by
> the name of the vector itself?
>
> What i would like to create is a total of 12 vectors, corresponding to the
> 12 months (January to December), say for the year 2006; thus the name of a
> resulting single vector should be something like "January2006", and the
> final file that will be stored in my WD should have the same name
> ("January2009.xlsx").
>
> The number of the elements of each vector must correspond to the length in
> days of the single months (considering a non-leap-year, 356 days) multiplied
> by 2 (e.g. "January2006" will have 31*2=62 elements, "February2006" will
> have 28*2=56 elements, and so on).
>
> Finally, the elements of the vectors should be named as:
> "010106_aaa","010106_bbb","020106_aaa","020106_bbb", ... ,
> "310106_aaa","310106_bbb".
>
> To sum up, at the end of the process i would like to obtain 12 vectors as it
> follows:
>
> Jauary2006("010106_aaa","010106_bbb","020106_aaa","020106_bbb", ... ,
> "310106_aaa","310106_bbb")
> .
> .
> .
> .
> .
> December2006("010106_aaa","010106_bbb","020106_aaa","020106_bbb", ... ,
> "310106_aaa","310106_bbb")
>
> Any help would be particularly welcome and appreciated.
> Cheers,
>
> NP
>
> * Italiano - rilevata
> * Inglese
> * Italiano
> * Francese
> * Spagnolo
> * Tedesco
>
> * Inglese
> * Italiano
> * Francese
> * Spagnolo
> * Tedesco
>
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Graphing lines of different lengths
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. On Thu, Feb 16, 2012 at 12:42 PM, Jordan Patricia Sinclair wrote: > Hello all. > I need to graph multiple lines of different lengths on the same graph. When > I try to add lines I get an error due to different lengths. The only thing I > could find when reading up on it was that 'if the data are inputted > separately they must be of the same length'. > Could someone please let me know how to input the data together? > Thanks for any help, > Jordan > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Defining a viewport scale in {Grid}
Read the Details section in ?viewport carefully. You are treating
xscale/yscale as if they are xlim/ylim in base graphics. They are not.
It may take some trial and error on your part to figure out how
exactly this works, in general you are setting the size and location
of each polygon relative to the parent viewport (in this case the
plotting region), so you need to make the sizes smaller. Here are a
couple of examples:
pushViewport(viewport(width=unit(.5,'npc'),height=unit(.5,'npc'),just=c('left','bottom')))
# ... draw your 5 polygons
dev.off()
pushViewport(viewport(width=unit(.5,'npc'),height=unit(.5,'npc'),just='bottom',yscale=c(-1,1)))
# ... draw your 5 polygons
Hope this helps
On Thu, Feb 16, 2012 at 6:47 AM, geotheory wrote:
> Am just feeling my way into the grid library, and cannot figure out how to
> define the plot limits. 3/5 of the example polygons below plot in the
> default 0-1 range viewport. But when I try to redefine the viewport the
> polygons plot in the same places. I also get the same result without
> employing push/pop. (As you can see from the scale I'm trying to introduce,
> I want to plot map polygons.) Grateful for an idea where I'm going wrong.
>
> pushViewport(viewport(xscale=c(-180,180),yscale=c(-90,90)))
>
> grid.polygon(x=c(0,5,5,0)/100, y=c(0,0,5,5)/100, draw = T)
> grid.polygon(x=c(100,95,95,100)/100, y=c(100,100,95,95)/100, draw = T)
> grid.polygon(x=c(50,70,70,50)/100, y=c(50,50,70,70)/100, draw = T)
> grid.polygon(x=-c(10,10,40,40)/100, y=c(10,40,40,10)/100, draw = T)
> #NEGATIVE X
> grid.polygon(x=c(50,52,52,50)/100, y=-c(50,50,52,52)/100, draw = T)
> #NEGATIVE Y
>
> popViewport()
>
> Thanks in advance.
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Defining-a-viewport-scale-in-Grid-tp4393974p4393974.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] separating plot output from debug output
If you don't dev.off(), all plots will be sent to the open graphical
device. That usually doesn't impact behavior of other output types:
pdf(file='fooout.pdf')
hist(x <- rnorm(100))
y <- sin(x)
print(str(y))
cat(y,file='fooout.txt')
plot(x,y)
dev.off()
Hope this helps
On Wed, Feb 15, 2012 at 3:43 PM, Tom Roche wrote:
>
> I'm attempting to refactor an R script that does a lot of plotting,
> among other things. Ideally I'd like to do something like
>
> setup # does pdf(...)
> for each part of input {
> plot(process(part))
> }
> cleanup # does dev.off()
>
> but have problems:
>
> 1 I'm plotting to PDF, so everytime I dev.off() creates a new file,
> and I want everything in one file (as does my boss :-)
> 2 I'm doing the work on a cluster, where I very much do not have root,
> and which has a fairly minimal set of installed packages, so I can't
> just call something external like 'pdftk' to merge the PDFs as I go.
>
> 3 As part of the processing, I printf status and debug messages, which
> I don't want in my PDF(s).
>
> The solutions I can imagine are
>
> 1 Append to a single PDF, but I understand this is not feasible, no?
>
> 2 Create a buncha PDFs with code above, download them to my laptop,
> merge them to a single PDF, upload it. Feasible but annoying and
> kludgey.
>
> 3 Separate processing from plotting, e.g.,
>
> setup # but not pdf(...)
> for each part of input {
> write(process1(part), intermediate)
> }
> pdf(...)
> for each part of intermediate {
> plot(process2(part))
> }
> cleanup # does dev.off()
>
> Again, feasible but kludgey.
>
> 4 No status and debug messages. I hope to be that good someday :-)
>
> Am I missing something? Are there clean solutions to this problem?
>
> TIA, Tom Roche
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Filling out a data frame row by row.... slow!
First, in R there is no need to declare the dimensions of your objects
before they are populated so couldn't you reduce some run time by not
going through the double data.frame step ?
> df<- data.frame()
> df
data frame with 0 columns and 0 rows
> for(i in 1:100) for(j in 1:3) df[i,j]<- runif(1)
> str(df)
'data.frame': 100 obs. of 3 variables:
...
Second, about populating an environment ?assign might work better for you
> e<- new.env()
> system.time(for(i in 1:1) e$a[i]<- rnorm(1,i))
user system elapsed
0.970.000.96
> rm(e)
> e<- new.env()
> system.time(for(i in 1:1) assign('a',rnorm(1,i),env=e))
user system elapsed
0.170.000.17
Third, how are you reading in the file? and what does that mean "not
knowing in advance..." ? Bill's suggestion to not populate the
data.frame line by line is probably the "real" solution to your
problem, as otherwise it's a little like kicking a turtle to make it
go faster...try to find a rabbit instead.
Posting a minimal example of your file format would have really
helped. Often using ?scan to read the whole (or big chunks of the)
file into R, followed by a customized formatting function that
utilizes ?grep and ?strsplit to reconstruct the data you want in
columns, solves the NEED to populate a data frame line by line.
Hope this helps
Elai
> One complication is I don't know the names of the columns I'm assigning to
> before I read them off the file. And crazily, if I change this:
> data$x[i] <- i + 0.1
>
> where data is an environment and x a primitive vector, to use a computed
> name instead:
>
> data[[colname]][i] <- i + 0.1
>
> Then I get back to way-superlinear performance. Eventually I found I could
> work around it like:
>
> eval(substitute(var[ix] <- data,
> list(var=as.name(colname), ix=i, data = i+0.1)),
> envir = data)
>
> but... as workarounds go that seems to be on the crazy nuts end of the
> scale. Why does [[]] impose a performance penalty?
>
> Peter
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Plotting function image
Inline
On Wed, Feb 15, 2012 at 3:04 AM, uday wrote:
> Hi ,
>
> Thanks for reply
>
> My latitude and longitude contains 9-10 observations per file
> when I run coords <- expand.grid(lat=1:5,long=1:5) then my computer
>
You don't have to run this part. As your original post did not provide
any details on data, coords is only a toy example.
> *** caught segfault ***
> address 0x6951c20, cause 'memory not mapped'
>
> Traceback:
> 1: .C("spline_eval", z$method, nu = as.integer(n), x =
> as.double(xout), y = double(n), z$n, z$x, z$y, z$b, z$c, z$d, PACKAGE =
> "stats")
> 2: spline(gam.data$x[, col.data], gam.smooths.all$fit[, m], xout =
> gam.results.global[m, , "x.values"], ties = mean)
> 3: eval.with.vis(expr, envir, enclos)
> 4: eval.with.vis(ei, envir)
> 5: source(file.path(getwd(), "Skripte", "r",
> "GAM_hourly", "1_calcs_GAM_all_sites_hourly.R"), echo =
> TRUE, max.deparse.length = 2e+05)
>
> Possible actions:
> 1: abort (with core dump, if enabled)
> 2: normal R exit
> 3: exit R without saving workspace
> 4: exit R saving workspace
>
I doubt expand.grid was the cause of this crash, and I see no
reference to it in traceback.
> so is there any other way to deal with this ?
> My main objective to plot this observations on global map
> the sample data is
> lat <- -11.3082 -11.6041 -11.9002 -12.1961 -12.1461 -12.7881 -12.6657
> -12.8467 -13.7233 -13.6271
> lon<- 135.6423 135.5799 135.5184 135.4558 133.5313 135.3321 134.6688
> 133.9839 132.0651 131.5528
> gas <- 1688.91 1679.24 1677.77 1635.60 1652.77 1663.43 1642.16 1671.84
> 1674.65 1665.54
>
Your lat,long do not seem to be forming a uniform grid. ?image would
probably spit an error about "increasing x and y" or something similar
(no access to R on this machine).
You could try
library(lattice)
levelplot(gas~lat+long,data=yourdata)# data argument is optional
-only if lat,long,gas are columns of a data frame
You may also want to look at the maps package for projections of lat/long
Hope this helps
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Plotting-function-image-tp4387796p4389959.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Plotting function image
In the absence of data coords <- expand.grid(lat=1:5,long=1:5) coords$z <- rnorm(25) Coords<- unstack(coords,z~long) image(as.matrix(Coords)) On Tue, Feb 14, 2012 at 10:36 AM, uday wrote: > I have some data set which has latitude, longitude and Z values. > I would like to plot them on global map. > I am thinking to use image > image(x, y, z, zlim, xlim, ylim, col = heat.colors(12), > add = FALSE, xaxs = "i", yaxs = "i", xlab, ylab, > breaks, oldstyle = FALSE, useRaster = FALSE, ...) > > here I have x = longitude > y= latitude > z = z > > but z suppose to be matrix but my z values are single array , could > somebody please tell me how to convert these single array of z to matrix ? > > > > -- > View this message in context: > http://r.789695.n4.nabble.com/Plotting-function-image-tp4387796p4387796.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] best subset selection on random effects model
I don't know of any package that will do it (or if violating the
marginality principle by having non-nested models even makes sense)
but you could always build your own search through all possible
models. Problem is for a large number of fixed effects this can take a
good bit of time...
Here is a piece of code I had from a few of years back. There may be
more elegant solutions that will reduce run time such as using
update.lme. Basically my code finds all possible combinations of model
terms with ?combn and than fits them.
require(nlme)
fm <- lme(distance ~age*Sex, random = ~ 1|Subject,data=Orthodont,method='ML')
Terms <- terms(fm)
todrop <- 1:length(attr(Terms,'term.labels'))
subs <- unlist(sapply(todrop,function(p)
combn(todrop,p,simplify=F)),recursive =F)
fm.subList <- lapply(subs[-length(subs)],function(s,...){
newf<- formula(drop.terms(terms(fm),s,keep.response = TRUE))
update(fm,newf)
})
names(fm.subList) <- sapply(fm.subList, function(x) paste('fm',attr(
terms(x),'term.labels'),sep='.'))
sort(sapply(fm.subList,BIC)) # 1st is best based on BIC
Cheers
On Tue, Feb 14, 2012 at 11:25 AM, Tao Zhang wrote:
> The models are not nested. I would like to consider all the possible
> subsets.
> I hope to output a table, where each row of the table indicates a best
> subset of the fixed effects for a particular model size.
>
> Thank you,
> Tao
>
>
> 2012/2/13 ilai
>>
>> The question is where do your models come from? Passing nested models
>> to ?anova.lme in nlme package or lme4 results in a likelihood ratio
>> test. Are you looking for something else/more ?
>>
>>
>> On Sun, Feb 12, 2012 at 8:02 PM, Tao Zhang wrote:
>> > Hi,
>> > I know leaps() computes the best subset selection for linear model,
>> > and the bestglm() computes the best subset selection for generalized
>> > linear
>> > model.
>> > Is there any package for best subset selection on random effects model,
>> > or
>> > mixed effects model?
>> >
>> > Thank you so much.
>> >
>> > Tao
>> >
>> > [[alternative HTML version deleted]]
>> >
>> > __
>> > R-help@r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide
>> > http://www.R-project.org/posting-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to suppress the empty plots in xyplot (lattice)
> read ?xyplot
> It takes a skip argument:
> ‘skip’: logical vector (default ‘FALSE’), replicated to be as
> long as the number of panels (spanning all pages). For
> elements that are ‘TRUE’, the corresponding panel
> position is skipped; i.e., nothing is plotted in that
> position. The panel that was supposed to be drawn there
> is now drawn in the next available panel position, and
> the positions of all the subsequent panels are bumped up
> accordingly. This is often useful for arranging plots in
> an informative manner.
>
> Or a 'drop.unused.levels' argument, or a 'subset'.
>
> Any one of these would work .
>
> Cheers
>
> On Tue, Feb 14, 2012 at 10:50 AM, Jun Shen wrote:
>> Thanks, Jeff,
>>
>> It did work in one way if I use
>>
>> xyplot(Y~X|as.factor(ID*PERIOD),data=...)
>>
>> But I would like to do something like
>>
>> xyplot(Y~X|as.factor(paste("ID=",ID)*paste("PERIOD=",PERIOD)),data=...)
>>
>> Then, it didn't work
>>
>> The error message:
>> Error in paste("ID=", ID) * paste("PERIOD=", PERIOD) :
>> non-numeric argument to binary operator
>>
>> On Tue, Feb 14, 2012 at 11:08 AM, Jeff Newmiller
>> wrote:
>>
>>> Set up a single (factor) variable that identifies the combinations that
>>> exist, and plot using that variable.
>>> ---
>>> Jeff Newmiller The . . Go Live...
>>> DCN: Basics: ##.#. ##.#. Live
>>> Go...
>>> Live: OO#.. Dead: OO#.. Playing
>>> Research Engineer (Solar/Batteries O.O#. #.O#. with
>>> /Software/Embedded Controllers) .OO#. .OO#. rocks...1k
>>> ---
>>> Sent from my phone. Please excuse my brevity.
>>>
>>> Jun Shen wrote:
>>>
>>> >Dear all,
>>> >
>>> >In a plot command like
>>> >
>>> >xyplot(Y~X|ID*PERIOD,data=...)
>>> >
>>> >xyplot will generate all the possible ID*PERIOD combinations. But not
>>> >all
>>> >of them have data in there. So I have a lot of empty plots. How can I
>>> >suppress those empty plots and ask xyplot only to generate plots
>>> >actually
>>> >with data. Thanks.
>>> >
>>> >Jun
>>> >
>>> > [[alternative HTML version deleted]]
>>> >
>>> >__
>>> >R-help@r-project.org mailing list
>>> >https://stat.ethz.ch/mailman/listinfo/r-help
>>> >PLEASE do read the posting guide
>>> >http://www.R-project.org/posting-guide.html
>>> >and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] different way for a for loop for several columns?
Inline
On Tue, Feb 14, 2012 at 3:16 AM, Nerak T wrote:
> Dear Ilai,
>
>
>
> Thanks for your answer. I'm indeed kind of a beginner in R, starting
> to discover the endless possibilities in R. My goal for the moment is indeed
> to get rid of the use of loops and to see through the apply family (which
> seems to be an endless maze for the moment). (For some reason, the apply
> function doesn’t seems to be logical for my brains which prefer to think in
> a loop way)
>
?apply ?lapply ?tapply, etc. are just wrappers for building more
efficient loops. If you "think in loops" (which you shouldn't) you are
also thinking in "apply". The reason it may seem like an endless maze
is because you use different wrappers for looping over different
object classes and indices, but at the end, the call to apply() is
similar to calling for().
e.g. consider a matrix with dimensions n x p. To sum rows you could
for(i in 1:n) sum(matrix[i,])
But better apply(matrix, 1 , sum) # where the 1 denotes the 1st
dimension (rows)
Same thing for sum columns
for(j in 1:p) sum(matrix[,j])
But better apply(matrix, 2 , sum) # where the 2 denotes the 2nd
dimension (columns)
The same "stuff" that goes in the loop can go to apply with small
syntax changes. e.g.
ll<- list(1:5,1:10,letters)
out<- list()
for(L in 1:3){
# ...
# ... a bunch of complicated functions/calculations
# ...
out[[ L ]] <- length( ll[[ L ]] )
}
out
Can be replaced with
lapply( ll, function(L) {
# ...
# ... a bunch of complicated function/calculations
# ...
length(L)
} )
This time use lapply since you are looping over L elements of the list ll.
>
>
> Your answer is really helpful.
> Something I found really interesting is your general comment. You say that I
> don’t need to declare variables? The reason I started to do this is because
> if I don’t, I get a message that the object is not found.
???
xx<- c(0,10,100) # declare xx
print(xx)
rnorm(3,xx) # use it
rm(xx) # remove xx
rnorm(3, xx<- c(0,10,100)) define and use it "at the same time"
print(xx)
> If I create a data
> frame before the calculation with the right amount of rows, it seems to
> work. But if there is a way not to have to make them before, would be great?
Only in loops. What happens is you need to create a "storage" for the
result of your loop, since the objects created in the loop are
overwritten at each step:
for(i in 1:10) cat( i, '\n' )
i # i = 10 everything before was overwritten
> Because most of the time, I don’t need that column that I created (but
> couldn't create an empty data frame with right dimensions to solve the
> problem)…
See my lapply example for creating "empty" storage.
>
>> Last, in R you want to avoid loops as much as possible especially for
>> large data sets. Operations are performed on objects so 1:4 + 1 is
>> equivalent to for(i in 1:4) i+1
> This part I don’t understand… where do you put that “ 1:4 + 1 ” ?
>
You don't put it anywhere, it was in answer to your comment: " but
it’s not that I have created a function that has to be applied on a
whole column, calculations are done for the different rows…"
So, no! calculations are done on the object (which has some dimension
or is a list), only in rare cases do you need to loop over each
element (or dimension) of the object itself.
>
> Many thanks, I’m trying to learn as much as possible to be able to use R
> more efficient so I really appreciate your help.
Pleasure. Good luck !
>
>
>
> Kind regards,
>
> Nerak
>
>
>
>
>
>
>
>
>
>> Date: Mon, 13 Feb 2012 23:43:29 -0700
>> Subject: Re: [R] different way for a for loop for several columns?
>> From: ke...@math.montana.edu
>> To: nera...@hotmail.com
>
>>
>> Nerak,
>> Your example could have been done without a loop at all (at least this
>> calculation), or as you already know by calling one of the apply
>> family functions which are more efficient (but are still "loops"):
>>
>> test<- data.frame(
>>
>> Date=c(1980,1980,1980,1980,1981,1981,1981,1981,1982,1982,1982,1982,1983,1983,1983,1983),
>> C = c(0,0,0,0,5,2,0,0,0,15,12,10,6,0,0,0),
>> B = c(0,0,0,0,9,6,2,0,0,24,20,16,2,0,0,0),
>> F = c(0,0,0,0,6,5,1,0,0,18,16,12,10,5,1,0)
>> )
>> test.2 <- test[,-1] > 1
>> aggregate(test.2, list(test$Date), sum)
>>
>> # See ?aggregate for more details. it also has a time series method
>> which may be useful for you.
>>
>> A general comment. if you are or will be using R a bit more, it may
>> benefit you to study the manuals or find a good basic tutorial. You
>> seem to be applying the conventions of some other progr
Re: [R] If (x > 0)
On Tue, Feb 14, 2012 at 6:41 AM, Schmidt, Michael
wrote:
> So the function must come BEFORE the call to the function...I see.
Yes. May be different than what you're used to but in R think of
functions as just another set of "objects". Therefore they must be
declared in the environment from which they are called beforehand,
just like the arguments.
Cheers
Elai.
> Thanks for that info.
> Take care
> Mike
>
> -Original Message-
> From: ila...@gmail.com [mailto:ila...@gmail.com] On Behalf Of ilai
> Sent: Monday, February 13, 2012 6:19 PM
> To: Schmidt, Michael
> Subject: Re: [R] If (x > 0)
>
> After placing res <- conditional1(x) ; cat("result= ",res,"\n") AFTER
> defining the function conditional1 (or R wouldn't know what "conditional1"
> is), your script worked flawlessly for me.
> From the terminal in xubuntu10.04:
> $ R --no-save --slave --args 1 < test1.R [1] 1 result= 1 $ R --no-save
> --slave --args -1e-6 < test1.R [1] -1e-06 result= -1 $ R --no-save --slave
> --args 0 < test1.R [1] 0 result= 0
>
> So maybe you had an old (bad) version of conditional1 which was used by test1
> ?
>
> Cheers
>
>
>
> On Mon, Feb 13, 2012 at 2:59 PM, Schmidt, Michael
> wrote:
>> Hi,
>> I am new to R. I was trying to get a very simple program to run. Take
>> one number from the command line. If the number < 0 return -1. If
>> number > 0 return 1 and if the number == 0 return 0. The code is in a
>> file called test1.R
>>
>>
>> The code:
>>
>> #useage: R --no-save --args 5 < test1.R args = (commandArgs(TRUE)) x =
>> as.numeric(args[1])
>> print(x)
>>
>> res <- conditional1(x)
>> cat("result= ",res,"\n")
>>
>> conditional1 <- function(x){
>> result <- 0
>> if (x > 0) {
>> result <- 1
>> } else if (x < 0) {
>> result <- -1
>> }
>> return(result)
>> }
>>
>> result <- 1
>> } else if (x < 0) {
>> result <- -1
>> }
>> return(result)
>> }
>>
>>
>> The output:
>>>R --no-save --slave --args 1 < test1.R
>> [1] 1
>> result= 1
>>>R --no-save --slave --args -1 < test1.R
>> [1] -1
>> result= -1
>>>] R --no-save --slave --args 0 < test1.R
>> [1] 0
>> result= -1
>>
>>
>> The problem:
>> For arguments 1 and -1 it works as intended. For 0 (zero) it does not. If
>> the 0 value is passed into the function named "conditional1" I would expect
>> both if-statements to evaluate to false and 0 being return. From what I can
>> tell (0 < 0) evaluates to true since -1 is returned. Hm...
>> What is going on? What am I doing wrong? Why is this happening? I am baffled!
>> Any help would be appreciated.
>> Thanks
>> Mike
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] passing an extra argument to an S3 generic
Hi Michael,
Try the attached. The only change to your script is in the first line
where I explicitly tell hatvalues to use methods (the "infmlm" class
stays). I also commented out all your TESTME at the end.
source('mlminfl-testHELP.R')
Now this should have worked for you too. Let me know. Sorry about that
"just remove the class". Had somewhat of a brain glitch when writing
the E-mail and wasn't clear.
Cheers
On Tue, Feb 14, 2012 at 8:05 AM, Michael Friendly wrote:
> On 2/11/2012 12:00 PM, ilai wrote:
>>
>> You are setting a new class ("inflmlm") at the end of mlm.influence.
>> Remove that second to last line and enjoy your new S3 method.
>
>
> Thanks for the suggestion, but it doesn't help -- I still get the same
> behavior whether mlm.influence returns a classed object or not.
> As well, I am defining print.inflmlm() and as.data.frame.inflmlm() methods
> for these objects, so I do need mlm.influence to return a classed object.
>
> My hatvalues.mlm is designed to be similar in structure to
> stats::hatvalues.lm where the S3 generic is defined.
>
>
> hatvalues.mlm <- function(model, m=1, infl, ...)
> {
> if (missing(infl)) {
> infl <- mlm.influence(model, m=m, do.coef=FALSE);
> }
> hat <- infl$H
> m <- infl$m
> names(hat) <- if(m==1) infl$subsets else apply(infl$subsets,1, paste,
> collapse=',')
> hat
> }
>
>> hatvalues
> function (model, ...)
> UseMethod("hatvalues")
>
>
>> hatvalues.lm
>
> function (model, infl = lm.influence(model, do.coef = FALSE),
> ...)
> {
> hat <- infl$hat
> names(hat) <- names(infl$wt.res)
> hat
> }
>
>
>
> The idea is that the infl= argument specifies a call to the computational
> function, mlm.influence() in my case, just as lm.influence() does in the
> stats package.
>
> The logic of UseMethod is that it should dispatch on the class of the
> *first* argument to the function, which in my test case is c("mlm", "lm")
>
>
>> Rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ n+s+ns+na+ss, data=Rohwer2)
>> class(Rohwer.mod)
> [1] "mlm" "lm"
>
>> trace(hatvalues)
>> hatvalues(Rohwer.mod, m=2)
> trace: hatvalues(Rohwer.mod, m = 2)
>
> Error in UseMethod("hatvalues") :
> no applicable method for 'hatvalues' applied to an object of class
> "c('double', 'numeric')"
>> hatvalues(Rohwer.mod)
> trace: hatvalues(Rohwer.mod)
>
> 1 2 3 4 5 6 7
> 8
> 0.16700926 0.21845327 0.14173469 0.07314341 0.56821462 0.15432157 0.04530969
> 0.17661104
> ...
>
>
> I'm still stumped on why with the extra argument m=2, R sees this
> as an object of class c('double', 'numeric'). As well, I can't see
> any way to debug this.
>
>
>
>> I'm not sure, but I think it is just the new class "inflmlm" applied
>> to inf in the formals of hatvalues.mlm confused the dispatch
>> mechanism. You would think the error message will call the offending
>> class not "numeric" double" but that's above my pay grade...
>>
>> You could probably put back the inflmlm class assignment with an
>> explicit call to UseMethod in hatvalues.mlm ?
>>
>> Cheers
>>
>> On Fri, Feb 10, 2012 at 2:35 PM, Michael Friendly
>> wrote:
>>>
>>> On 2/10/2012 4:09 PM, Henrik Bengtsson wrote:
>>>>
>>>>
>>>> So people may prefer to do the following:
>>>>
>>>> hatvalues.mlm<- function(model, m=1, infl, ...)
>>>> {
>>>> if (missing(infl)) {
>>>> infl<- mlm.influence(model, m=m, do.coef=FALSE);
>>>> }
>>>>
>>>> hat<- infl$H
>>>> m<- infl$m
>>>> names(hat)<- if(m==1) infl$subsets else apply(infl$subsets,1,
>>>> paste, collapse=',')
>>>> hat
>>>> }
>>>
>>>
>>> Thanks; I tried exactly that, but I still can't pass m=2 to the mlm
>>> method
>>> through the generic
>>>
>>>> hatvalues(Rohwer.mod)
>>>
>>> 1 2 3 4 5 6
>>> 7
>>> 8
>>> 0.16700926 0.21845327 0.14173469 0.07314341 0.56821462 0.15432157
>>> 0.04530969
>>> 0.17661104
>>> 9 10 11 12 13 14
>>> 15
>>> 1
Re: [R] Writing R-scripts
It seems all you are doing in the if statements is defining functions.
You need to actually "apply" them to some arguments, then you can pass
results.
i.e.
f<- function(x,type,...){
a<- function(...){ 2* x }
b<- function(...) { x^2 }
if(type==1){ ret<- a(x) }
if(type==2){ ret<- b(x) }
ret
}
You may also want to see ?switch or ?local
Cheers
On Mon, Feb 13, 2012 at 2:59 PM, Cem Girit wrote:
> Hello,
>
> This is my first attempt to write a script in R. The program below
> is intended to do some parametric tests on group data. There are subroutines
> for each type of test. The call to the "parametric.tests", routine sets the
> argument "testtype" for the test to be used. How can I transfer the
> calculated values (in "result" below) in each routine to the calling
> parametric.tests routine?
>
> Cem
>
>
>
>
> ## testtype : 1: Duncan
> ## 2: Dunnett
> ## resp: response variable, must be numeric and vector
> ## group: group id for resp, numeric or character
> ## alpha: CL 0.05 or 0.01
> ## vehicle: Control group name for Dunnett
>
> parametric.tests<-function(testtype, resp, group, vehicle, alpha)
> {
> if (testtype==1){
>
> ## resp: response variable, must be numeric and vector
> ## group: group id for resp, numeric or character
> ## alpha: CL 0.05 or 0.01
>
> duncan.test <- function (resp, group, alpha) {
>
> .
>
> result <- data.frame(label=label, estimate=Estimate, alpha=alpha,
> lower=Lower, upper=Upper, p.value=pval, significance=sig)
> return(result)
> }
> }
>
> else if (testtype==2){
>
>
> dunnett.test <- function(resp, group, vehicle, alpha)
> {
>
> .
>
> result <- data.frame(label=label, estimate=Estimate, alpha=alpha,
> lower=Lower, upper=Upper, p.value=pval, significance=sig)
> return(result)
> }
> }
> }
>
> Cem
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] only 0s may be mixed with negative subscripts
Like this ? sum every batch of n rows? N <- 2 t(kronecker(diag(nrow(sens2)/N),rep(1,N))) %*% as.matrix(sens2) Time Y X rownumber [1,] 2657130134 0.3236854 0.1080725 3 [2,] 2657130134 0.4621816 0.1198302 7 [3,] 2657130134 0.4341804 0.106958311 Check rownumber column: 1+2 = 3, 3+4 = 7, 5+6=11 ... Or sum every 3 rows: N <- 3 t(kronecker(diag(nrow(sens2)/N),rep(1,N))) %*% as.matrix(sens2 ) Time Y X rownumber [1,] 3985695201 0.5646581 0.1693338 6 [2,] 3985695202 0.6553893 0.165527215 Hope this helps (more...). On Mon, Feb 13, 2012 at 4:04 PM, Hasan Diwan wrote: > On 13 February 2012 14:46, ilai wrote: >> The function you posted runs without error (on these 6 lines), but >> does not return anything that looks remotely like a sum, or cumsum of >> anything. Can you clarify what you are trying to do? I assume by "sum >> of every other row" you don't mean summing Time, X and Y for rows >> 1,3,5,..., ? > > I'm trying to get a piecewise sum of every n rows in a given column in > this data set. >> For the sum of sens2[c(1,3,5,...),] for every column (assuming no NA's >> in the data) you could > > I do format checking well before getting to this stage in the analysis. >> >> (1:nrow(sens2) %% 2) %*% as.matrix(sens2) > > That does not do what I want... Again, what it should return is a list > of the sum of the current and preceding row. > -- > Sent from my mobile device > Envoyait de mon portable > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] only 0s may be mixed with negative subscripts
The function you posted runs without error (on these 6 lines), but
does not return anything that looks remotely like a sum, or cumsum of
anything. Can you clarify what you are trying to do? I assume by "sum
of every other row" you don't mean summing Time, X and Y for rows
1,3,5,..., ?
For the sum of sens2[c(1,3,5,...),] for every column (assuming no NA's
in the data) you could
(1:nrow(sens2) %% 2) %*% as.matrix(sens2)
Time Y X rownumber
[1,] 3985695201 0.56826 0.1369527 9
Hope this helps
On Mon, Feb 13, 2012 at 11:56 AM, Hasan Diwan wrote:
> I'd like to get the sum of every other row in a data.frame. When I
> actually set about doing this, I get the error in the subject line of
> this message. A sample of my data is below, followed by the function
> call that should give me the results I want:
>
>> dput(head(sens2))
> structure(list(Time = c(1328565067, 1328565067.05, 1328565067.1,
> 1328565067.15, 1328565067.2, 1328565067.25), Y = c(0.0963890795246276,
> 0.227296347215609, 0.240972698811569, 0.221208948983498, 0.230898231782485,
> 0.203282153087549), X = c(0.0245045248243853, 0.0835679411703579,
> 0.0612613120609633, 0.058568910563872, 0.0511868450318788, 0.0557714205674231
> ), rownumber = 1:6), .Names = c("Time", "Y", "X", "rownumber"
> ), row.names = c(NA, 6L), class = "data.frame")
>> speedX <- sapply(sens2[sens2$rownumber %% 2 == 0,], function(row) {
>> cumsum(c(sens2[row+1,3], sens2[row,3]))}, simplify=TRUE)
> Error in xj[i] : only 0's may be mixed with negative subscripts
> Help?
> --
> Sent from my mobile device
> Envoyait de mon portable
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] best subset selection on random effects model
The question is where do your models come from? Passing nested models to ?anova.lme in nlme package or lme4 results in a likelihood ratio test. Are you looking for something else/more ? On Sun, Feb 12, 2012 at 8:02 PM, Tao Zhang wrote: > Hi, > I know leaps() computes the best subset selection for linear model, > and the bestglm() computes the best subset selection for generalized linear > model. > Is there any package for best subset selection on random effects model, or > mixed effects model? > > Thank you so much. > > Tao > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] object not found - Can not figure out why I get this error: Error in NROW(yCoordinatesOfLines) : object 'low' not found
Ah, scoping rules...
Consider:
f <- function(x,...) plot(x,xlim=c(low,high),...)
f(1:10,low=2,high=9) # "Error ... object 'low' not found "
But:
f <- function(x,low,high,...) plot(x,xlim=c(low,high),...)
f(1:10,2,9,col=2) # beautiful red points [low,high]
Sorry I can't be more specific but never heard of quantmod.
Hope this helps anyway.
On Sat, Feb 11, 2012 at 7:43 AM, Samo Pahor wrote:
> Hi,
>
> I have been using R for over a year now. I am a very happy user. Thank you
> for making this happen.
>
> This is my first question to this list.
>
> I trying to add some functions to quantmod that would enable me to draw
> arbitrary lines and text and make sure they are redrawn. I have created
> following function:
>
> require(quantmod)
>
> # Add horizontal line to graph produced by quantmod::chart_Series()
> add_HorizontalLine<-function(yCoordinatesOfLines, on=1, ...) {
> lenv <- new.env()
> lenv$add_horizontalline <- function(x, yCoordinatesOfLines, ...) {
> xdata <- x$Env$xdata
> xsubset <- x$Env$xsubset
>
> x0coords <- rep(1, NROW(yCoordinatesOfLines))
> x1coords <- rep(NROW(xdata[xsubset]), NROW(yCoordinatesOfLines))
>
> if ((NROW(x0coords) > 0) & (NROW(x1coords) > 0)) {
> segments(x0coords,
> yCoordinatesOfLines,
> x1coords,
> yCoordinatesOfLines, ...)
> # abline(h=yCoordinatesOfLines, ...)
> }
> }
> mapply(function(name, value) {assign(name,value,envir=lenv)},
> names(list(yCoordinatesOfLines=yCoordinatesOfLines,...)),
> list(yCoordinatesOfLines=yCoordinatesOfLines,...))
> exp <- parse(text=gsub("list","add_horizontalline",
> as.expression(substitute(list(x=current.chob(),
>
> yCoordinatesOfLines=yCoordinatesOfLines, ..., srcfile=NULL)
> plot_object <- current.chob()
> lenv$xdata <- plot_object$Env$xdata
> # plot_object$set_frame(sign(on)*abs(on)+1L)
> plot_object$set_frame(2*on)
> plot_object$add(exp,env=c(lenv, plot_object$Env),expr=TRUE)
> plot_object
> }
>
> # Short test function that uses add_HorizontalLine
> test<-function(series, low=20, high=80) {
> chart_Series(SPX, subset="2012")
> add_TA(RSI(Cl(SPX)))
> plot(add_HorizontalLine(c(low, high), on=2, col=c('green', 'red'),
> lwd=2))
> }
>
> # Actual test
> SPX <- getSymbols("^GSPC", from="2000-01-01", auto.assign=FALSE)
> dev.new()
> test(SPX)
>
> This gives me the following error:
>> test(SPX)
> Error in NROW(yCoordinatesOfLines) : object 'low' not found
>
> What am I doing wrong here? Any hints highly appreciated.
>
> The funniest thing is that this was working and somehow broke it...
>
> Best,
> Samo
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Longitudinal Factor Analysis
?factanal There is also package sem (structural equations model) by John Fox. I'm sure there are more (maybe more fitting your situation) but these two came to mind first... Cheers On Sun, Feb 12, 2012 at 6:51 AM, Gregory Gilbert wrote: > I have a data set in the format below. I would like to perform a factor > analysis on BM1-BM20 as they are 20 biomechanical measurements of the trial > (hitting a baseball). However, my observations are not independent and, I > assume, I have to account for this. I have consulted the R literature via > RSeek and Google (and I have consulted some very knowledgeable colleagues) > and cannot find any method of doing this in R. Is someone aware of a method > of reducing the 20 biomechanical measurements? > > Thank you for your help. > > Greg > > ID PitchType Location SwingNumber BM1...BM20 Hit > Velocity > 1 Fastball 1 1 xx ... xx 1 > 62 > 1 Fastball 1 2 xx ... xx 1 > 68 > 1 Fastball 1 3 xx ... xx 0 > 1 Fastball 1 4 xx ... xx 0 > 1 Fastball 1 5 xx ... xx 1 > 71 > 1 Fastball 2 1 xx ... xx 0 > 1 Fastball 2 2 xx ... xx 1 > 61 > 1 Fastball 2 3 xx ... xx 1 > 64 > 1 Fastball 2 4 xx ... xx 1 > 59 > 1 Fastball 2 5 xx ... xx 0 > 1 Fastball 2 6 xx ... xx 0 > 1 Fastball 2 7 xx ... xx 0 > 1 Fastball 2 8 xx ... xx 0 > 1 Fastball 3 1 xx ... xx 1 > 78 > 1 Fastball 3 2 xx ... xx 0 > 1 Fastball 3 3 xx ... xx 1 > 74 > 1 Fastball 4 1 xx ... xx 0 > 1 Fastball 4 2 xx ... xx 1 > 69 > 1 Fastball 4 3 xx ... xx 1 > 72 > 1 Fastball 4 4 xx ... xx 1 > 65 > 1 Fastball 4 5 xx ... xx 0 > 1 Fastball 4 6 xx ... xx 0 > 1 Fastball 4 7 xx ... xx 0 > 1 Fastball 4 8 xx ... xx 1 > 74 > 1 Fastball 4 9 xx ... xx 1 > 78 > 1 Fastball 4 10 xx ... xx 1 > 76 > 1 Fastball 5 1 xx ... xx 0 > 1 Fastball 5 2 xx ... xx 1 > 53 > 1 Fastball 5 3 xx ... xx 0 > 1 Fastball 5 4 xx ... xx 1 > 51 > 1 Fastball 5 5 xx ... xx 0 > 1 Fastball 5 6 xx ... xx 1 > 55 > . > . > . > 40 Fastball 5 1 xx ... xx 0 > 40 Fastball 5 2 xx ... xx 0 > 40 Fastball 5 3 xx ... xx 1 > 60 > 40 Fastball 5 4 xx ... xx 0 > 40 Fastball 5 5 xx ... xx 0 > > > > > > -- > Gregory E. Gilbert, EdD, MSPH, PStat® > ASA Accredited Professional Statistician™ > > [[alternative HTML version deleted]] > > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] multiple histograms from a dataframe
On Sat, Feb 11, 2012 at 9:10 PM, David Winsemius wrote: > > On Feb 11, 2012, at 6:25 PM, Adel ESSAFI wrote: > >> >> >> Le 11 février 2012 02:33, David Winsemius a écrit >> : >> >> On Feb 10, 2012, at 7:05 PM, Adel ESSAFI wrote: >> >> >> what can I do to draw to figures together using lattice? > > > You need to describe what you mean by "together". It is possible that the > goup parameter is what you want but that's just a guess. It's also possible > that the formular operator "+" will give you what you desire. Perhaps: > > xyplot( law[,67] + law[,66] ~ law[,3]| > interaction(law[,1],law[,2]),type='l') > > >> >> it draws multiple histograms but by selecting distinct values of law[,1] >> The deal is to make the same thing but for a couple of columns > > > That doesn't make any sense to me. But then I do apologize for the English > language. It's horribly complex and syntactically a mess. > English language aside, do you want histograms (count frequency) or scatterplots (y~x) ? as David suggested the group argument or '+' combined with allow.multiple is useful, but not trivial for histograms. A simple solution is to consider alternatives such density curves: densityplot(~as.numeric(volcano[,c(2:5)])|cut(volcano[,1],4), groups=gl(4,nrow(volcano)), key=simpleKey(as.character(2:5),points=F,lines=T,columns=4)) In the future, a minimal working example and clear statement of your goal may reduce the need for multiple posts. Regards > >> >> Thanks in advance for help >> >> Adel >> >> >> -- >> >> David Winsemius, MD >> West Hartford, CT >> >> -- >> PhD candidate in Computer Science >> Address >> 3 avenue lamine, cité ezzahra, Sousse 4000 >> Tunisia >> tel: +216 97 246 706 (+33640302046 jusqu'au 15/6) >> fax: +216 71 391 166 > > > David Winsemius, MD > West Hartford, CT > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Lattice 3d coordinate transformation
Thank you Deepayan, your answer put me on the path to SOLVED !!!
Actually passing projected corners to panel.rect was the first thing I
tried, but couldn't get it to work. However, panel.3dpolygon in
latticeExtra did the trick.
I'm posting it here for completion.
require(lattice) ; require(latticeExtra)
set.seed(1113)
d <- data.frame(x=runif(30),y=runif(30),g=gl(2,15))
d$z <- with(d,rnorm(30,3*asin(x^2)-5*y^as.integer(g),.1))
d$z <- d$z+min(d$z)^2
surf <- by(d,d$g,function(D){
fit <- lm(z~poly(x,2)*poly(y,2),data=D)
outer(seq(0,1,l=10),seq(0,1,l=10),function(x,y,...)
predict(fit,data.frame(x=x,y=y)))
})
panel.3d.surf <- function(x, y, z, rot.mat, distance, zlim.scaled, ...){
zz <- surf[[packet.number()]] ; n <- nrow(zz)
lp <- level.colors(zz, at = do.breaks(range(zz), 20), col.regions
= heat.colors(20))
s <- seq(-.5,.5,l=n) ; cntrds <- expand.grid(s,s) ; index <- 0
apply(cntrds,1,function(i){
index <<- index+1
xx <- i[1]+c(-.5,-.5,.5,.5)/(n-1) ; yy <- i[2]+c(-.5,.5,.5,-.5)/(n-1)
panel.3dpolygon(xx,yy, zlim.scaled[1], rot.mat, distance,
border=lp[index], col=lp[index],...)
})
panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled, ...)
}
cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud = panel.3d.surf,
zoom = 1,screen=list(z= 21,y=0,x=-60),aspect = c(1,1), panel.aspect = 1,
scales=list(z=list(arrows=F,tck=0),x=list(distance=.75)),
par.box=list(lwd=NA),lwd=3)
## Beautiful !
On Sat, Feb 11, 2012 at 6:00 AM, Deepayan Sarkar
wrote:
> On Fri, Feb 10, 2012 at 12:43 AM, ilai wrote:
>> Hello List!
>> I asked this before (with no solution), but maybe this time... I'm
>> trying to project a surface to the XY under a 3d cloud using lattice.
>> I can project contour lines following the code for fig 13.7 in
>> Deepayan Sarkar's "Lattice, Multivariate Data Visualization with R",
>> but it fails when I try to "color them in" using panel.levelplot.
>> ?utilities.3d says there may be some bugs, and I think
>> ltransform3dto3d() is not precise (where did I hear that?), but is
>> this really the source of my problem? Is there a (simple?) workaround,
>> maybe using 3d.wire but projecting it to XY? How? Please, any insight
>> may be useful.
>
> I don't think this will be that simple. panel.levelplot() essentially
> draws a bunch of colored rectangles. For a "3D" projection, each of
> these will become (four-sided) polygons. You need to compute the
> coordinates of those polygons, figure out their fill colors (possibly
> using ?level.colors) and then draw them.
>
> -Deepayan
>
>
>> Thanks in advance,
>> Elai.
>>
>> A working example:
>>
>> ## data "d" and predicted "surf":
>> set.seed(1113)
>> d <- data.frame(x=runif(30),y=runif(30),g=gl(2,15))
>> d$z <- with(d,rnorm(30,3*asin(x^2)-5*y^as.integer(g),.1))
>> d$z <- d$z+min(d$z)^2
>> surf <- by(d,d$g,function(D){
>> fit <- lm(z~poly(x,2)*poly(y,2),data=D)
>> outer(seq(0,1,l=10),seq(0,1,l=10),function(x,y,...)
>> predict(fit,data.frame(x=x,y=y)))
>> })
>> ##
>> # This works to get contours:
>> require(lattice)
>> cloud(z~x+y|g,data=d,layout=c(2,1), type='h', lwd=3, par.box=list(lty=0),
>> scales=list(z=list(arrows=F,tck=0)),
>> panel.3d.cloud = function(x, y, z,rot.mat, distance,
>> zlim.scaled, nlevels=20,...){
>> add.line <- trellis.par.get("add.line")
>> clines <- contourLines(surf[[packet.number()]],nlevels = nlevels)
>> for (ll in clines) {
>> m <- ltransform3dto3d(rbind(ll$x-.5, ll$y-.5,
>> zlim.scaled[1]), rot.mat,
>> distance)
>> panel.lines(m[1,], m[2,], col = add.line$col, lty = add.line$lty,
>> lwd = add.line$lwd)
>> }
>> panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled =
>> zlim.scaled, ...)
>> }
>> )
>> # But using levelplot:
>> panel.3d.levels <- function(x, y, z,rot.mat, distance, zlim.scaled,...)
>> {
>> zz <- surf[[packet.number()]]
>> n <- nrow(zz)
>> s <- seq(-.5,.5,l=n)
>> m <- ltransform3dto3d(rbind(rep(s,n),rep(s,each=n),zlim.scaled[1]),
>> rot.mat, distance)
>> panel.levelplot(m[1,],m[2,],zz,1:n^2,col.regions=heat.colors(20))
>> panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled,
>> ...)
>> }
>> cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud =
&g
Re: [R] passing an extra argument to an S3 generic
You are setting a new class ("inflmlm") at the end of mlm.influence.
Remove that second to last line and enjoy your new S3 method.
I'm not sure, but I think it is just the new class "inflmlm" applied
to inf in the formals of hatvalues.mlm confused the dispatch
mechanism. You would think the error message will call the offending
class not "numeric" double" but that's above my pay grade...
You could probably put back the inflmlm class assignment with an
explicit call to UseMethod in hatvalues.mlm ?
Cheers
On Fri, Feb 10, 2012 at 2:35 PM, Michael Friendly wrote:
> On 2/10/2012 4:09 PM, Henrik Bengtsson wrote:
>>
>> So people may prefer to do the following:
>>
>> hatvalues.mlm<- function(model, m=1, infl, ...)
>> {
>> if (missing(infl)) {
>> infl<- mlm.influence(model, m=m, do.coef=FALSE);
>> }
>>
>> hat<- infl$H
>> m<- infl$m
>> names(hat)<- if(m==1) infl$subsets else apply(infl$subsets,1,
>> paste, collapse=',')
>> hat
>> }
>
> Thanks; I tried exactly that, but I still can't pass m=2 to the mlm method
> through the generic
>
>> hatvalues(Rohwer.mod)
> 1 2 3 4 5 6 7
> 8
> 0.16700926 0.21845327 0.14173469 0.07314341 0.56821462 0.15432157 0.04530969
> 0.17661104
> 9 10 11 12 13 14 15
> 16
> 0.05131298 0.45161152 0.14542776 0.17050399 0.10374592 0.12649927 0.33246744
> 0.33183461
> 17 18 19 20 21 22 23
> 24
> 0.17320579 0.26353864 0.29835817 0.07880597 0.14023750 0.19380286 0.04455330
> 0.20641708
> 25 26 27 28 29 30 31
> 32
> 0.15712604 0.15333879 0.36726467 0.11189754 0.30426999 0.08655434 0.08921878
> 0.07320950
>
>> hatvalues(Rohwer.mod, m=2)
> Error in UseMethod("hatvalues") :
> no applicable method for 'hatvalues' applied to an object of class
> "c('double', 'numeric')"
>
> ## This works:
>> hatvalues.mlm(Rohwer.mod, m=2)
> ... output snipped
>
>> hatvalues
>
> function (model, ...)
> UseMethod("hatvalues")
>
>
>
>>
>
> -Michael
>
>
> --
> Michael Friendly Email: friendly AT yorku DOT ca
> Professor, Psychology Dept.
> York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
> 4700 Keele Street Web: http://www.datavis.ca
> Toronto, ONT M3J 1P3 CANADA
>
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Formatting Y axis.
yaxt='n' in ?par and ?axis are your friends.
# A plot on log scale labeled with original:
plot(x,log(y),yaxt='n')
axis(2,at=pretty(log(y)),labels=round(exp(pretty(log(y)
Works for qqnorm and boxplots, as well as other top level fun.
By the way this is a FAQ.
On Fri, Feb 10, 2012 at 9:43 AM, sock.o wrote:
> I've looked around and I just can't find anything that will work for my
> needs. This is a bit of a 2 part question but pertaining to the same topic
> so bare with me.
>
> The first is with my qq plot. On the Y axis of my qq plot it'll have my
> sample quantities but because my data is log-normal it'll show numbers
> between 0 - 5 (depending on the data). I'd like to know how to get it,
> instead of displaying 0-5 to display actual values so that way when reading
> it I know where my potential problems are. Would just save a lot of time and
> be an easier more understandable read. Here's sample code to plug in and
> see.
>
> x <- rlnorm(400, 2.3, .8)
> y <- rlnorm(400, 3, 1.6)
>
> plot.new()
> q1 <- qqnorm(log1p(y))
> q2 <- qqnorm(log1p(x))
> points( q1, col = "red")
> points( q2, col = "blue")
> abline( qqline(log1p(y), col = "red") )
> abline( qqline(log1p(x), col = "blue") )
> legend( "bottomright", inset = 0.02, title = "Loc_ID",
> c("Sayincom",
> "Sayout"), fill = c( "red", "blue"), horiz = TRUE)
>
> So what I want is for the y axis. or even secondary y axis if it's lined up
> properly to read the actual values rather than the logged values.
>
> The second part of this is in regards to my boxplots. In order to fit all of
> the data in a readable manner I've been making log = "y". Which is fine, but
> when trying to estimate the difference via notches it's difficult to
> estimate how far it really is based on the axis being logged. You can use
> the same x as above, or y, whichever floats your boat.
>
> boxplot ( x, notch = TRUE, log = "y", boxwex = .50)
>
> Versus
>
> boxplot ( x, notch = TRUE, boxwex = .50)
>
> What I need from this is for the ability to zoom in on a particular
> location. Rather than be forced to look at things in a log format on my y
> axis to just be able to say, ok I have my mean at 25 and my standard
> deviation is 10 so I want to concentrate this boxplot between 5 and 45 to
> view 2 standard deviations. (That whole mean thing and standard deviations
> are arbitrary numbers but the concept remains)
>
> Any help would be much appreciated. If there are particular packages
> required for any suggestions please include their name in your response.
>
> If you need actual imaging of what the actual graphs I'm working with look
> like as opposed to the randomly generated values that I'm supplying then
> just say so and I'll post the actual graphs I'm trying to work with.
>
> Sincerely,
>
> Sock
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Formatting-Y-axis-tp4376843p4376843.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Help needed please
Your script is rather inefficient with spurious cbind calls. Any
particular reason not to use
?ar directly ?
Call:
ar.yw.default(x = simtimeseries, order.max = 4)
Coefficients:
1234
1.9440 -1.9529 0.8450 -0.2154
Order selected 4 sigma^2 estimated as 15.29
To repeat the sim, you could use a for() loop but ?sapply is better:
out<- sapply(1:100,function(...){
simtimeseries <- arima.sim(n=1024,list(order=c(4,0,0),
ar=c(2.7607, -3.8106, 2.6535,
-0.9258),sd=sqrt(1)))
aryule <- ar.yw(simtimeseries,order.max=4)
c( c(aryule$ar,NA)[1:4] , aryule$var.pred )
}
)
rowMeans(out[1:4,]) # mean phi(1),...,4 see ?rowMeans for dealing with NA's
mean(out[5,])# mean sig^2
Cheers
On Fri, Feb 10, 2012 at 6:42 AM, Jaymin Shah wrote:
> I have coded a time series from simulated data:
>
> simtimeseries <- arima.sim(n=1024,list(order=c(4,0,0),ar=c(2.7607, -3.8106,
> 2.6535, -0.9258),sd=sqrt(1)))
> #show roots are outside unit circle
> plot.ts(simtimeseries, xlab="", ylab="", main="Time Series of Simulated Data")
>
> # Yule
>
>
> q1 <- cbind(simtimeseries[1:1024])
> q2 <- t(q1)%*%q1
> s0 <- q2/1204
> r1 <- cbind(simtimeseries[1:1023])
> r2 <- cbind(simtimeseries[2:1024])
> r3 <- t(r1)%*%r2
> s1 <- r3/1204
> t1 <- cbind(simtimeseries[1:1022])
> t2 <- cbind(simtimeseries[3:1024])
> t3 <- t(t1)%*%t2
> s2 <- t3/1204
> u1 <- cbind(simtimeseries[1:1021])
> u2 <- cbind(simtimeseries[4:1024])
> u3 <- t(u1)%*%u2
> s3 <- u3/1204
> v1 <- cbind(simtimeseries[1:1020])
> v2 <- cbind(simtimeseries[5:1024])
> v3 <- t(v1)%*%v2
> s4 <- v3/1204
>
> i0 <- c(s0,s1,s2,s3)
> i1 <- c(s1,s0,s1,s2)
> i2 <- c(s2,s1,s0,s1)
> i3 <- c(s3,s2,s1,s0)
>
> gamma <- cbind(i0,i1,i2,i3)
> eta <-c(s1,s2,s3,s4)
> inversegamma <- solve(gamma)
> phihat <- inversegamma%*%eta
> phihat
>
> Phihat <- cbind(phihat)
> s <- c(s1,s2,s3,s4)
> S <- cbind(s)
> sigmasquaredyule <- s0 - (t(Phihat)%*%S)
> sigmasquaredyule
>
>
>
> I did a yule walker estimate on the simulated data and wanted to work out phi
> hat which is a vector of 4 values and sigmasquaredyule which is one value.
> However, I want to run the simulated data 100 times i.e. in a for loop and
> then take the averages of the phi hat values and sigmasquaredyule value.
>
> How would i repeat this simulated time series lots of times (e.g. a 100
> times) and store the average value of phi hat and sigmasquaredyule.
>
> Thank you
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] making multiple lines using qqplot
Melissa,
par(new=T) works as many times as you use it. You don't provide data,
but (assuming it is not NULL) more likely your n=500 qqplot was just
obscuring the points of the n=50 plot.
Reverse the order (i.e. qqplot 500 first, 50, 5 last) and see if all
three are there (as there are more 500 you should still see green on
the extremes).
Second, you say you want "lines" but used pch=20. replace with
type='l' to get lines (may also help with your main problem). If you
want to stick with points and your device supports it, you can
consider semi-transparent colors.
Cheers
On Thu, Feb 9, 2012 at 9:00 PM, Melrose2012
wrote:
> Hi Everyone,
>
> I want to make 3 lines on the same graph (not as subplots, all in the same
> window, one on top of each other) and I want them to be quantile-quantile
> plots (qqplot). Essentially, I am looking for the equivalent of Matlab's
> "hold on" command to use with qqplot. I know I can use 'points' or 'lines',
> but these do not give me a qqplot (only appear to work as scatter plots). I
> found the syntax 'par(new=TRUE)' but that only seems to work for two lines,
> not for three.
>
> My script currently looks like:
>
> qqplot(nq.n5,tq.n5,col="red",xlab="Normal Distribution Quantiles",ylab="t
> Distribution Quantiles",main="Quantile-Quantile Plot of Normal vs
> t-Distribution for Various Sample Sizes",pch=20)
> par(new=TRUE)
> qqplot(nq.n50,tq.n50,col="blue",xlab="",ylab="",pch=20).
> par(new=TRUE)
> qqplot(nq.n500,tq.n500,col="green",xlab="",ylab="",pch=20)
> legend("topleft",c("n=5","n=50","n=500"),fill=c("red","blue","green"))
>
> I realize that this only plots the first and the third qqplot because by
> doing par(new=TRUE) again, it gets rid of the middle one. I don't know how
> to get around this and get all 3 lines on the same plot.
>
> Can anyone please help me with this syntax?
>
> Thank you very much for your time and advice!
>
> Cheers,
> Melissa
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/making-multiple-lines-using-qqplot-tp4375273p4375273.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to change x-axis ?
plot(1:10, xaxt='n') # Don't plot the x-axis axis(1,at=c(2,5,10)) # Construct your own See ?par ?axis On Thu, Feb 9, 2012 at 2:48 PM, summer wrote: > Hi, > I want to define the x-axis in plot as 2,4,6.100 instead of the default > one. How to do that? Many thanks. > > > -- > View this message in context: > http://r.789695.n4.nabble.com/how-to-change-x-axis-tp4374503p4374503.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Apply pmax to dataframe with different args based on dataframe factor
Your attempt was just overly complicated. All you needed was
threshold <- c( .2 , .4 , .5 )[ df$track ]
df$value <- pmax(threshold, df$value)
df # desired outcome
Cheers
On Thu, Feb 9, 2012 at 3:56 PM, Idris Raja wrote:
> # I have a dataframe in the following form:
>
> track <- c(rep('A', 3), rep('B', 4), rep('C', 4))
> value <- c(0.15, 0.25, 0.35, 0.05, 0.99, 0.32, 0.13, 0.80, 0.75, 0.60, 0.44)
> df <- data.frame(track=factor(track), value=value)
>
> #> print(df)
> #track value
> #1 A 0.15
> #2 A 0.25
> #3 A 0.35
> #4 B 0.05
> #5 B 0.99
> #6 B 0.32
> #7 B 0.13
> #8 C 0.80
> #9 C 0.75
> #10 C 0.60
> #11 C 0.44
>
>
> # If any of the values are below a threshold value, I want to replace it
> with the
> # threshold value. The twist is that there is a different threshold value
> for
> # every track.
>
> # I tried something like this, but it's not working
>
> threshold <- list()
> threshold['A'] <- 0.2
> threshold['B'] <- 0.4
> threshold['C'] <- 0.5
>
>
> for (track in levels(df$track)){
> df[df$track==track,]$outcome <- pmax(df[df$track==track,]$outcome,
> threshold[track])
> }
> # Warning messages:
> # 1: In is.na(mmm) : is.na() applied to non-(list or vector) of type 'NULL'
> # 2: In is.na(mmm) : is.na() applied to non-(list or vector) of type 'NULL'
> # 3: In is.na(mmm) : is.na() applied to non-(list or vector) of type 'NULL'
>
>
> #**
> # Desired Results:
>
> #> print(df)
> #track value
> #1 A 0.20 # value changed
> #2 A 0.25
> #3 A 0.35
> #4 B 0.40 # value changed
> #5 B 0.99
> #6 B 0.40 # value changed
> #7 B 0.40 # value changed
> #8 C 0.80
> #9 C 0.75
> #10 C 0.60
> #11 C 0.50 # value changed
>
>
> # Any ideas? Thanks for reading.
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Tukey HSD
Also you want " for(i in 2:ncol(mydf)) { ..."
Your current setup of 1:ncol(mydf)-1 picks columns 0,1,2,...,ncol-1
(missing the last), which is not what you want...
Setting i<-2 in the first line was overridden by the call to loop.
Cheers
On Thu, Feb 9, 2012 at 2:33 PM, peter dalgaard wrote:
>
> On Feb 9, 2012, at 20:28 , jlpugh wrote:
>
>> Hey Folks:
>>
>> New to R and learning as I go, but really, I know just enough to get myself
>> into trouble. I've waded through everything up till now, and don't see
>> anything in the search that is directly helpful for the POS that I've
>> created.
>>
>> The GOAL: All I want in the world is a program that performs 1-way ANOVA's
>> on every column in a data set (taking the first column as the definition of
>> the groups) and then spits out ONLY those results that were significant (p
>> <= whatever I want), with their respective plots and TukeyHSD results.
>> Sounds simple, right?
>
> Quite likely, this could all be done more elegantly (lapply springs to mind),
> and there could be something in there that won't work at all, but to point
> out the glaring obvious: If the grouping factor is in column 1, then
> shouldn't you have var2~var1 in the aov() call?
>
>> var1 <- mydf[,1]
>> var2 <- mydf[,i]
>> fm1 <- aov(var1 ~ var2)
>> tky <- TukeyHSD(fm1)
>
>
> and, by the way, you _really_ do not want to go via capture.output to get at
> the p value. It is accessible from the summary(fm1) directly. For aov
> objects, it's like:
>
> summary(fm1)[[1]]$"Pr(>F)"[1]
>
> (for "lm" objects, try anova(fm)$"Pr(>F)"[1])
>
>>
>> Data:
>>
>> http://r.789695.n4.nabble.com/file/n4374072/test_text.txt test_text.txt
>>
>> (see upload)
>>
>> CODE: ( which occurs after TEST <- read.table("test_text.txt") )
>>
>> i <- 2;
>> sink (file = "test_output.txt", append = FALSE)
>> mydf <- data.frame(TEST)
>> for (j in 1:ncol(mydf)-1) {
>> var1 <- mydf[,1]
>> var2 <- mydf[,i]
>> fm1 <- aov(var1 ~ var2)
>> tky <- TukeyHSD(fm1)
>> otpt <- capture.output(summary(fm1))
>> i <- i+1;
>> lines <- as.vector(unlist(strsplit(otpt[2]," ")),mode="list") # gets
>> the
>> p-value
>> if (grepl("[1234567890]",lines[14],perl = TRUE)) { #make sure that
>> slot for
>> p-value has a number
>> number <- as.numeric(lines[14]) # make it numeric for logic
>> test
>> if (number <= 0.1) {
>> cat(otpt, sep = "\n\n")
>> cat(tky, sep = "\n\n")
>> quartz(boxplot(mydf[,i] ~ mydf[,1]))
>> }
>> }
>> }
>> sink()
>>
>> The ERROR:
>>
>> Error in TukeyHSD.aov(fm1) : no factors in the fitted model
>> In addition: Warning message:
>> In replications(paste("~", xx), data = mf) : non-factors ignored: var2
>>
>> It works if I leave out " tky <- TukeyHSD(fm1) " and the subsequent cat.
>> I've tried doing the Tukey on the mydf[,1].. itself, changing their classes,
>> etc. Perhaps the whole approach is flawed.
>>
>> THANKS!
>>
>>
>> --
>> View this message in context:
>> http://r.789695.n4.nabble.com/Tukey-HSD-tp4374072p4374072.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd@cbs.dk Priv: pda...@gmail.com
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] passing an extra argument to an S3 generic
You do not provide mlm.influence() so your code can't be reproduced.
Or did you mean to put lm.influence() in the formals to your hatvalues.mlm ?
If yes, then 1) you have a typo 2) lm.influence doesn't allow you to
pass on arguments, maybe try influence.lm instead.
Elai
On Thu, Feb 9, 2012 at 1:42 PM, Michael Friendly wrote:
> I'm trying to write some functions extending influence measures to
> multivariate linear models and also
> allow subsets of size m>=1 to be considered for deletion diagnostics. I'd
> like these to work roughly parallel
> to those functions for the univariate lm where only single case deletion
> (m=1) diagnostics are considered.
>
> Corresponding to stats::hatvalues.lm, the S3 method for class "lm" objects,
>
>> hatvalues <-function (model, ...)
> UseMethod("hatvalues")
>
>> hatvalues.lm <-
> function (model, infl = lm.influence(model, do.coef = FALSE), ...)
> {
> hat <- infl$hat
> names(hat) <- names(infl$wt.res)
> hat
> }
>
> I have, for class "mlm" objects
>
> hatvalues.mlm <- function(model, m=1, infl=mlm.influence(model, m=m, do.coef
> = FALSE), ...)
> {
> hat <- infl$H
> m <- infl$m
> names(hat) <- if(m==1) infl$subsets else apply(infl$subsets,1, paste,
> collapse=',')
> hat
> }
>
> where mlm.influence() does the calculations, but also allows the m= argument
> to specify subset size.
> Yet when I test this I can't seem to pass the m= argument directly, so that
> it gets stuffed in to the infl=
> call to mlm.influence.
>
> # fit an mlm
> library(heplots)
> Rohwer2 <- subset(Rohwer, subset=group==2)
> rownames(Rohwer2)<- 1:nrow(Rohwer2)
> Rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ n+s+ns+na+ss, data=Rohwer2)
>
>> class(Rohwer.mod)
> [1] "mlm" "lm"
>
>
> ## this doesn't work, as I would like it to, calling the hatvalues.mlm
> method, but passing m=2:
>> hatvalues(Rohwer.mod, m=2)
> Error in UseMethod("hatvalues") :
> no applicable method for 'hatvalues' applied to an object of class
> "c('double', 'numeric')"
>
> I don't understand why this doesn't just call hatvalues.mlm, since
> Rohwer.mod is of class "mlm".
>
> # These work -- calling hatvalues.mlm explicitly, or passing the infl=
> argument with the call to
> # mlm.influence
>> hatvalues.mlm(Rohwer.mod, m=2)
>> hatvalues(Rohwer.mod, infl=mlm.influence(Rohwer.mod,m=2))
>
> Can someone help me understand what is wrong and how to make the .mlm method
> allow m= to be passed
> directly to the infl= computation?
>
> thx,
> -Michael
>
> --
> Michael Friendly Email: friendly AT yorku DOT ca
> Professor, Psychology Dept.
> York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
> 4700 Keele Street Web: http://www.datavis.ca
> Toronto, ONT M3J 1P3 CANADA
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Lattice 3d coordinate transformation
Oops. Obviously I mean
"> A working example:
> require(lattice)
> ..."
On Thu, Feb 9, 2012 at 12:13 PM, ilai wrote:
> Hello List!
> I asked this before (with no solution), but maybe this time... I'm
> trying to project a surface to the XY under a 3d cloud using lattice.
> I can project contour lines following the code for fig 13.7 in
> Deepayan Sarkar's "Lattice, Multivariate Data Visualization with R",
> but it fails when I try to "color them in" using panel.levelplot.
> ?utilities.3d says there may be some bugs, and I think
> ltransform3dto3d() is not precise (where did I hear that?), but is
> this really the source of my problem? Is there a (simple?) workaround,
> maybe using 3d.wire but projecting it to XY? How? Please, any insight
> may be useful.
> Thanks in advance,
> Elai.
>
> A working example:
>
> ## data "d" and predicted "surf":
> set.seed(1113)
> d <- data.frame(x=runif(30),y=runif(30),g=gl(2,15))
> d$z <- with(d,rnorm(30,3*asin(x^2)-5*y^as.integer(g),.1))
> d$z <- d$z+min(d$z)^2
> surf <- by(d,d$g,function(D){
> fit <- lm(z~poly(x,2)*poly(y,2),data=D)
> outer(seq(0,1,l=10),seq(0,1,l=10),function(x,y,...)
> predict(fit,data.frame(x=x,y=y)))
> })
> ##
> # This works to get contours:
> require(lattice)
> cloud(z~x+y|g,data=d,layout=c(2,1), type='h', lwd=3, par.box=list(lty=0),
> scales=list(z=list(arrows=F,tck=0)),
> panel.3d.cloud = function(x, y, z,rot.mat, distance,
> zlim.scaled, nlevels=20,...){
> add.line <- trellis.par.get("add.line")
> clines <- contourLines(surf[[packet.number()]],nlevels = nlevels)
> for (ll in clines) {
> m <- ltransform3dto3d(rbind(ll$x-.5, ll$y-.5,
> zlim.scaled[1]), rot.mat,
> distance)
> panel.lines(m[1,], m[2,], col = add.line$col, lty = add.line$lty,
> lwd = add.line$lwd)
> }
> panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled =
> zlim.scaled, ...)
> }
> )
> # But using levelplot:
> panel.3d.levels <- function(x, y, z,rot.mat, distance, zlim.scaled,...)
> {
> zz <- surf[[packet.number()]]
> n <- nrow(zz)
> s <- seq(-.5,.5,l=n)
> m <- ltransform3dto3d(rbind(rep(s,n),rep(s,each=n),zlim.scaled[1]),
> rot.mat, distance)
> panel.levelplot(m[1,],m[2,],zz,1:n^2,col.regions=heat.colors(20))
> panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled, ...)
> }
> cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud =
> panel.3d.levels,
> scales=list(z=list(arrows=F,tck=0)),par.box=list(lty=0),lwd=3)
> # I also tried to "fill" between contours but can't figure out what to
> do with the edges and how to incorporate the x,y limits to 1st and nth
> levels.
> panel.3d.contour <- function(x, y, z,rot.mat, distance,xlim,ylim,
> zlim.scaled,nlevels=20,...)
> {
> add.line <- trellis.par.get("add.line")
> zz <- surf[[packet.number()]]
> clines <- contourLines(zz,nlevels = nlevels)
> colreg <- heat.colors(max(unlist(lapply(clines,function(ll) ll$level
> for (i in 2:length(clines)) {
> ll <- clines[[i]]
> ll0 <- clines[[i-1]]
> m <- ltransform3dto3d(rbind(ll$x-.5, ll$y-.5, zlim.scaled[1]),
> rot.mat, distance)
> m0 <- ltransform3dto3d(rbind(ll0$x-.5, ll0$y-.5,
> zlim.scaled[1]), rot.mat, distance)
> xvec <- c(m0[1,],m[1,ncol(m):1])
> yvec <- c(m0[2,],m[2,ncol(m):1])
> panel.polygon(xvec,yvec,col=colreg[ll$level],border='transparent')
> panel.lines(m[1,], m[2,], col = add.line$col, lty = add.line$lty,
> lwd = add.line$lwd)
> }
> panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled, ...)
> }
> cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud =
> panel.3d.contour,
> scales=list(z=list(arrows=F,tck=0)),par.box=list(lty=0),lwd=3)
>
> #
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Lattice 3d coordinate transformation
Hello List!
I asked this before (with no solution), but maybe this time... I'm
trying to project a surface to the XY under a 3d cloud using lattice.
I can project contour lines following the code for fig 13.7 in
Deepayan Sarkar's "Lattice, Multivariate Data Visualization with R",
but it fails when I try to "color them in" using panel.levelplot.
?utilities.3d says there may be some bugs, and I think
ltransform3dto3d() is not precise (where did I hear that?), but is
this really the source of my problem? Is there a (simple?) workaround,
maybe using 3d.wire but projecting it to XY? How? Please, any insight
may be useful.
Thanks in advance,
Elai.
A working example:
## data "d" and predicted "surf":
set.seed(1113)
d <- data.frame(x=runif(30),y=runif(30),g=gl(2,15))
d$z <- with(d,rnorm(30,3*asin(x^2)-5*y^as.integer(g),.1))
d$z <- d$z+min(d$z)^2
surf <- by(d,d$g,function(D){
fit <- lm(z~poly(x,2)*poly(y,2),data=D)
outer(seq(0,1,l=10),seq(0,1,l=10),function(x,y,...)
predict(fit,data.frame(x=x,y=y)))
})
##
# This works to get contours:
require(lattice)
cloud(z~x+y|g,data=d,layout=c(2,1), type='h', lwd=3, par.box=list(lty=0),
scales=list(z=list(arrows=F,tck=0)),
panel.3d.cloud = function(x, y, z,rot.mat, distance,
zlim.scaled, nlevels=20,...){
add.line <- trellis.par.get("add.line")
clines <- contourLines(surf[[packet.number()]],nlevels = nlevels)
for (ll in clines) {
m <- ltransform3dto3d(rbind(ll$x-.5, ll$y-.5,
zlim.scaled[1]), rot.mat,
distance)
panel.lines(m[1,], m[2,], col = add.line$col, lty = add.line$lty,
lwd = add.line$lwd)
}
panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled =
zlim.scaled, ...)
}
)
# But using levelplot:
panel.3d.levels <- function(x, y, z,rot.mat, distance, zlim.scaled,...)
{
zz <- surf[[packet.number()]]
n <- nrow(zz)
s <- seq(-.5,.5,l=n)
m <- ltransform3dto3d(rbind(rep(s,n),rep(s,each=n),zlim.scaled[1]),
rot.mat, distance)
panel.levelplot(m[1,],m[2,],zz,1:n^2,col.regions=heat.colors(20))
panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled, ...)
}
cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud = panel.3d.levels,
scales=list(z=list(arrows=F,tck=0)),par.box=list(lty=0),lwd=3)
# I also tried to "fill" between contours but can't figure out what to
do with the edges and how to incorporate the x,y limits to 1st and nth
levels.
panel.3d.contour <- function(x, y, z,rot.mat, distance,xlim,ylim,
zlim.scaled,nlevels=20,...)
{
add.line <- trellis.par.get("add.line")
zz <- surf[[packet.number()]]
clines <- contourLines(zz,nlevels = nlevels)
colreg <- heat.colors(max(unlist(lapply(clines,function(ll) ll$level
for (i in 2:length(clines)) {
ll <- clines[[i]]
ll0 <- clines[[i-1]]
m <- ltransform3dto3d(rbind(ll$x-.5, ll$y-.5, zlim.scaled[1]),
rot.mat, distance)
m0 <- ltransform3dto3d(rbind(ll0$x-.5, ll0$y-.5,
zlim.scaled[1]), rot.mat, distance)
xvec <- c(m0[1,],m[1,ncol(m):1])
yvec <- c(m0[2,],m[2,ncol(m):1])
panel.polygon(xvec,yvec,col=colreg[ll$level],border='transparent')
panel.lines(m[1,], m[2,], col = add.line$col, lty = add.line$lty,
lwd = add.line$lwd)
}
panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled, ...)
}
cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud = panel.3d.contour,
scales=list(z=list(arrows=F,tck=0)),par.box=list(lty=0),lwd=3)
#
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Row-wise kronecker product with Matrix package
Maybe one of these will improve:
>help.search('kronecker')
...
spam::kronecker Kronecker Products on Sparse Matrices
spam::spam.classClass "spam"
base::kronecker Kronecker Products on Arrays
Matrix::kronecker-methods
Methods for Function 'kronecker()' in Package
'Matrix'
I doubt it because they will calculate the "full" Kronecker prod and
it will be up to you to index the rows, but you never know...
system.time(A[, rep(seq(ncol(A)), each = ncol(B))] *
B[,rep(seq(ncol(B)),ncol(A))])
user system elapsed
0.016 0.000 0.019
> system.time(kronecker(A,B)[c(1,4),])
user system elapsed
0.008 0.000 0.008
> system.time(spam::kronecker(A,B)[c(1,4),])
user system elapsed
0.008 0.000 0.009
Cheers
On Thu, Feb 9, 2012 at 9:38 AM, Ally wrote:
>
> I'm trying to calculate the row-wise kronecker product A \Box B of two
> sparse matrices A and B, and am struggling to find a quick way to do this
> that takes advantage of sparseness. I thought a good idea would be to use
> "rep" to construct 2 matrices of the same dimension of the end product, and
> multiply these two together:
>
> library(Matrix)
> A<-Matrix(c(1,0,0,0,0,1,2,0), 2, 4)
> B<-Matrix(c(2,5,0,0,0,1,0,0,0,0), 2, 5)
>
> A[, rep(seq(ncol(A)), each = ncol(B))] * B[, rep(seq(ncol(B)),ncol(A))]
>
> This works, but for much larger problems is slow (compared to keeping A and
> B dense). I was wondering why this happens, and whether there might be a
> way around it?
>
> Thanks in advance for any advice,
>
> Alastair
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Row-wise-kronecker-product-with-Matrix-package-tp4373437p4373437.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] basic debugging
?debug will satisfy your curiosity regarding "debug mode" - strictly speaking it is not a "mode", just another function. On Wed, Feb 8, 2012 at 7:30 PM, R. Michael Weylandt wrote: > What is your question? More interestingly, what is "debug mode" in R? > > I'd suggest you look at traceback() -- there's also the powerful, but > possibly advanced, options(error=recover) > > Kerninghan said "The most effective debugging tool is still careful > thought, coupled with judiciously placed print statements." It's good > advice, particularly for debugging something in a loop. > > Michael > > On Wed, Feb 8, 2012 at 5:25 PM, ikuzar wrote: >> Hi, >> >> I have to debug my program. When I execute my function (in debug mode), I >> have got an error but I do not know which line is concerned. I do not want >> to do an infinite "Browse[2]>n with each line in my function... >> >> THanks for your help, >> >> ikuzar >> >> -- >> View this message in context: >> http://r.789695.n4.nabble.com/basic-debugging-tp4371113p4371113.html >> Sent from the R help mailing list archive at Nabble.com. >> >> __ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Setting up infile for R CMD BATCH
I'm going to declare this SOLVED. Yes, if you don't want a separate
script for batch, you will need to modify the original script so it
either readline or skips it. Here is an example:
# Save in file myTest.R
# Add this local function to the beginning of your original "program"
and replace all the readline prompts with myreadline.
# The new function takes on two arguments:
# "what": the original object (in your example it was "type"<-...)
# "prompt": The same string from readline
# All it is doing is searching for "Answers.R", sourcing if available
or prompting if not.
myreadline <- function(what,prompt){
ans <- length(grep('Answers.R',list.files(),fixed=T))>0 # add a
warning for multiple files
if(ans) source('Answers.R')
else{
ret <- as.integer(readline(prompt))
assign(what,ret,1)
}
}
# here is an interactive program to square only negative values
print(x <- rnorm(1))
myreadline('type','x>0 ? 1:T,2:F \n')
print(x^type)
### END myTest.R
Running myTest interactivly:
> source('myTest.R')
[1] -0.3712215
x>0 ? 1:T,2:F
2
[1] 0.1378054 # -.37^2
> source('myTest.R')
[1] 0.3160747
x>0 ? 1:T,2:F
1
[1] 0.3160747 # .316^1
# Create a list of answers
> dump('type',file='Answers.R')
# run again this time with answers available
> source('myTest.R')
[1] -1.088665 # skips prompt
[1] -1.088665 # -1.088^1 (type in Answer.R ==1)
# Now you can also run as batch
$ R CMD BATCH myTest.R out.R
$ cat out.R
# ...
> print(x <- rnorm(1))
[1] -1.248487
> myreadline('type','x>0 ? 1:T,2:F \n')
> print(x^type)
[1] -1.248487
That's it. Only thing is to keep the file names in the working
directory straight.
Enjoy
On Wed, Feb 8, 2012 at 12:39 PM, Gang Chen wrote:
> Sorry Elai for the confusions.
>
> Let me try to reframe my predicament. The main program "myTest.R" has
> been written in interactive mode with many readline() lines embedded.
> Suppose a user has already run the program once before in interactive
> mode with all the answers saved in a text file called answer.R. Now
> s/he does not want to go through the interactive again because it's a
> tedious process, and would like to run the program but in batch mode
> with answer.R. And that's why I tried the following which didn't pan
> out:
>
> R CMD BATCH answer.R output.Rout
>
> Of couse I could rewrite a separate program specifically for batch
> mode as you suggested previously with eval(), for example. However, is
> there an approach to keeping the original program so that the user
> could run both interactive and batch mode? That probably requires
> modifying the readline() part, but how?
>
> Thanks,
> Gang
>
>
> On Wed, Feb 8, 2012 at 2:08 PM, ilai wrote:
>> Gang,
>> Maybe someone here has a different take on things. I'm afraid I have
>> no more insights on this unless you explain exactly what you are
>> trying to achieve, or more importantly why? That may help understand
>> what the problem really is.
>>
>> Do you want to save an interactive session for future runs? then
>> ?save.image and all your "answers" are in the environment. In this
>> case consider putting an "if(!exists('type') | length(type)<1 |
>> is.na(type))" before "type<- readline(...)" in your script so type
>> wouldn't be overwritten in subsequent runs.
>>
>> If your goal is to batch evaluate multiple answer files from users
>> (why else would you ask questions with readline?), then you should
>> have enough to go on with my answer and the examples in ?eval.
>>
>> Elai
>>
>>
>> On Wed, Feb 8, 2012 at 9:04 AM, Gang Chen wrote:
>>> Hi Elai,
>>>
>>> Thanks a lot for the suggestions! I really appreciate it...
>>>
>>> Your suggestion of using eval() and creating those answers in a list
>>> would work, but there is no alternative to readline() with which I
>>> could read the input in batch mode? I'm asking this because I'd like
>>> to have the program work in both interactive and batch mode.
>>>
>>> Thanks again,
>>> Gang
>>>
>>>
>>> On Wed, Feb 8, 2012 at 12:50 AM, ilai wrote:
>>>> Ahh,
>>>> I think I'm getting it now. Well, readlines() is not going to work for
>>>> you. The help file ?readline clearly states "In non-interactive use
>>>> the result is as if the response was RETURN and the value is ‘""’."
>>>> The implication
Re: [R] Setting up infile for R CMD BATCH
Gang,
Maybe someone here has a different take on things. I'm afraid I have
no more insights on this unless you explain exactly what you are
trying to achieve, or more importantly why? That may help understand
what the problem really is.
Do you want to save an interactive session for future runs? then
?save.image and all your "answers" are in the environment. In this
case consider putting an "if(!exists('type') | length(type)<1 |
is.na(type))" before "type<- readline(...)" in your script so type
wouldn't be overwritten in subsequent runs.
If your goal is to batch evaluate multiple answer files from users
(why else would you ask questions with readline?), then you should
have enough to go on with my answer and the examples in ?eval.
Elai
On Wed, Feb 8, 2012 at 9:04 AM, Gang Chen wrote:
> Hi Elai,
>
> Thanks a lot for the suggestions! I really appreciate it...
>
> Your suggestion of using eval() and creating those answers in a list
> would work, but there is no alternative to readline() with which I
> could read the input in batch mode? I'm asking this because I'd like
> to have the program work in both interactive and batch mode.
>
> Thanks again,
> Gang
>
>
> On Wed, Feb 8, 2012 at 12:50 AM, ilai wrote:
>> Ahh,
>> I think I'm getting it now. Well, readlines() is not going to work for
>> you. The help file ?readline clearly states "In non-interactive use
>> the result is as if the response was RETURN and the value is ‘""’."
>> The implication is you cannot use it to "insert" different answers as
>> if you were really there.
>> How about using eval() instead? You will need to make the answers a
>> named list (or just assigned objects).
>>
>> test <- expression({
>> if(a>2) print('+')
>> else print('I got more')
>> b <- b+3 # reassign b in the environment
>> print(b)
>> print(c)
>> d^2
>> })
>> dump('test',file='myTest.R') ; rm(test)
>>
>> # make the answers.R file:
>>
>> a=5 ; b=2 ; c=2 ; d=3
>> source("myTest.R")
>> eval(test)
>>
>> # Now, from the terminal R CMD BATCH answers.R out.R
>> # And here is my $ cat out.R
>> ... flushed
>>> a=5 ; b=2 ; c=2 ; d=3
>>> source("myTest.R")
>>> eval(test)
>> [1] "+"
>> [1] 5
>> [1] 2
>> [1] 9
>>>
>>> proc.time()
>> user system elapsed
>> 0.640 0.048 0.720
>>
>> Would this work?
>> Elai
>>
>>
>>
>>
>> On Tue, Feb 7, 2012 at 4:05 PM, Gang Chen wrote:
>>> Suppose I create an R program called myTest.R with only one line like
>>> the following:
>>>
>>> type <- as.integer(readline("input type (1: type1; 2: type2)? "))
>>>
>>> Then I'd like to run myTest.R in batch mode by constructing an input
>>> file called answers.R with the following:
>>>
>>> source("myTest.R")
>>> 1
>>>
>>> When I ran the following at the terminal:
>>>
>>> R CMD BATCH answer.R output.Rout
>>>
>>> it failed to pick up the answer '1' from the 2nd line in answers.R as
>>> shown inside output.Rout:
>>>
>>>> source("myTest.R")
>>> input type (0: quit; 1: type1; 2: type2)?
>>>> 1
>>> [1] 1
>>>
>>> What am I missing here?
>>>
>>> Thanks in advance,
>>> Gang
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Setting up infile for R CMD BATCH
On Tue, Feb 7, 2012 at 10:50 PM, ilai wrote:
> Ahh,
> I think I'm getting it now. Well, readlines() is not going to work for
> you. The help file ?readline clearly states "In non-interactive use
> the result is as if the response was RETURN and the value is ‘""’."
> The implication is you cannot use it to "insert" different answers as
> if you were really there.
> How about using eval() instead? You will need to make the answers a
> named list (or just assigned objects).
>
> test <- expression({
> if(a>2) print('+')
> else print('I got more')
> b <- b+3 # reassign b in the environment
> print(b)
> print(c)
> d^2
> })
> dump('test',file='myTest.R') ; rm(test)
>
> # make the answers.R file:
>
> a=5 ; b=2 ; c=2 ; d=3
> source("myTest.R")
> eval(test)
>
> # Now, from the terminal R CMD BATCH answers.R out.R
> # And here is my $ cat out.R
> ... flushed
>> a=5 ; b=2 ; c=2 ; d=3
>> source("myTest.R")
>> eval(test)
> [1] "+"
> [1] 5
> [1] 2
> [1] 9
>>
>> proc.time()
> user system elapsed
> 0.640 0.048 0.720
>
> Would this work?
> Elai
>
>
>
>
> On Tue, Feb 7, 2012 at 4:05 PM, Gang Chen wrote:
>> Suppose I create an R program called myTest.R with only one line like
>> the following:
>>
>> type <- as.integer(readline("input type (1: type1; 2: type2)? "))
>>
>> Then I'd like to run myTest.R in batch mode by constructing an input
>> file called answers.R with the following:
>>
>> source("myTest.R")
>> 1
>>
>> When I ran the following at the terminal:
>>
>> R CMD BATCH answer.R output.Rout
>>
>> it failed to pick up the answer '1' from the 2nd line in answers.R as
>> shown inside output.Rout:
>>
>>> source("myTest.R")
>> input type (0: quit; 1: type1; 2: type2)?
>>> 1
>> [1] 1
>>
>> What am I missing here?
>>
>> Thanks in advance,
>> Gang
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Setting up infile for R CMD BATCH
You're not missing anything.
In your output.Rout: the ">1" right after the source('test') is the
"1" inputed from answers.R. the "[1] 1" is the result of test. Remove
the second line from answers.R and see what happens (hint: script ends
after the readline prompt).
Just out of curiosity, why will you use a script that requires user
input (readlines) in batch mode ?
Cheers
On Tue, Feb 7, 2012 at 4:05 PM, Gang Chen wrote:
> Suppose I create an R program called myTest.R with only one line like
> the following:
>
> type <- as.integer(readline("input type (1: type1; 2: type2)? "))
>
> Then I'd like to run myTest.R in batch mode by constructing an input
> file called answers.R with the following:
>
> source("myTest.R")
> 1
>
> When I ran the following at the terminal:
>
> R CMD BATCH answer.R output.Rout
>
> it failed to pick up the answer '1' from the 2nd line in answers.R as
> shown inside output.Rout:
>
>> source("myTest.R")
> input type (0: quit; 1: type1; 2: type2)?
>> 1
> [1] 1
>
> What am I missing here?
>
> Thanks in advance,
> Gang
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] binomial vs quasibinomial
Not really an R question, now is it ? more like pure stats. I'm guessing you didn't get an answer because this list can't tell you how to analyze your data (or in your case, approve an incorrect analysis). Regarding the part of your question that is R related, I think you may be confused on what the dispersion parameter is. In summary.glm it is reported in the line above the null deviance: "(Dispersion parameter for quasibinomial family taken to be ...)". You said "the dispersion factor is the same for the quasibinomial...". I doubt that for any real data it gets estimated at 1 (which was assumed for the binomial family model). Also, since you say your p-values are higher on the second model, that leads me to believe it is indeed taken to be greater than 1. Check the two models again. Regards On Tue, Feb 7, 2012 at 4:11 AM, Jhope wrote: > After looking at 48 glm binomial models I decided to try the quasibinomial > with the top model 25 (lowest AIC). To try to account for overdispersion > (residual deviance 2679.7/68 d.f.) After doing so the dispersion factor is > the same for the quasibinomial and less sectors of the beach were > significant by p-value. While the p-values in the binomial were more > significant for each section of the beach. -- telling me more about the > beach. > > Is this ok? Can I just look at the binomial glm model 25 and look at its > p-values for beach sections and forget about the quasibinomial model 25? > > J > > > Call: > glm(formula = cbind(Shells, TotalEggs - Shells) ~ Sector:Veg:Aeventexhumed, > family = quasibinomial, data = data.to.analyze) > > > -- > View this message in context: > http://r.789695.n4.nabble.com/binomial-vs-quasibinomial-tp4364371p4364371.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] I bet apply has a solution
duplicated(Data..) On Mon, Feb 6, 2012 at 11:34 AM, LCOG1 wrote: > Hi all > For the data below, I would like to return a logical value indicating > differences in the data. > > #Create data > Data..<-data.frame(a=rep(1,10),b=c(rep(1,9),2),c=c(rep(1,8),2,2)) > > a b c > 1 1 1 1 > 2 1 1 1 > 3 1 1 1 > 4 1 1 1 > 5 1 1 1 > 6 1 1 1 > 7 1 1 1 > 8 1 1 1 > 9 1 1 2 > 10 1 2 2 > > > So what I want is to return logical value telling me if all the values are > the same. So the result would be a b c DidChange > 1 1 1 1 FALSE > 2 1 1 1 FALSE > 3 1 1 1 FALSE > 4 1 1 1 FALSE > 5 1 1 1 FALSE > 6 1 1 1 FALSE > 7 1 1 1 FALSE > 8 1 1 1 FALSE > 9 1 1 2 TRUE > 10 1 2 2 TRUE > > I bet apply could handle this elegantly but that family of functions is > still not 100% intuitive to me. Thoughts. Thanks everyone > > Cheers, > Josh > > > -- > View this message in context: > http://r.789695.n4.nabble.com/I-bet-apply-has-a-solution-tp4362294p4362294.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multi-page PDF using dev.copy2pdf(filename, onefile=TRUE)?
Doug,
dev.copy2pdf closes the connection after it's "done", so onefile is
meaningless. To look at each plot before copy to a single pdf, you
could open a pdf(...) but revert between it and your graphic device:
graphics.off()
plot(1:7, 1:7)
x11c<- dev.cur() # your current graphics device
pdf(file="test.pdf")
dev.set(which=x11c) # back from pdf
dev.copy() # copy
dev.set(which=x11c)# back to graphic
plot(1:5, 1:5,col=2,pch=2)
dev.copy()
dev.set(which=x11c)
# ... etc. however many more plots
# don't forget to close the pdf device at the end:
dev.off(which=x11c+1)
# Both plots are in 'test.pdf'
Depending on your application you might be able to simplify things
with dev.next/dev.prev, or wrap this sequence into a little helper
function to be used in a loop.
Enjoy,
Elai
On Mon, Feb 6, 2012 at 6:44 AM, Doug Hill wrote:
> Hi all. I want to generate a sequence of n plots and save them into a single
> PDF file, one plot per page. From the R docs and other sources I gather the
> basic way to do this is save plot 1 into a file then append the 2:n plots to
> the same file.
> This code shows my basic approach, but for some reason only the last plot is
> saved into the pdf. I've tried different variations (e.g. using onefile only
> in the second call, or only in the first), to no avail. The comments show
> what I see if I step through the code one line at a time:
> scratch<-function() {
> graphics.off()
> plot(1:7, 1:7) # Opens a graphics window and displays a 7-point plot in it,
> as expected
> dev.copy2pdf(file="test.pdf", onefile=TRUE) # I see the 7-point plot in Adobe
> reader, as expected
> plot(1:5, 1:5) # Overwrites in the graphics window the 7-point plot with a
> 5-point, as expected
> dev.copy2pdf(file="test.pdf", onefile=TRUE) # Overwrites test.pdf so that it
> contains only the 7-point plot
> }
> A couple things:
> (1) The reason I don't just use something like pdf(filename) plot(...)
> plot(...) dev.off() is that I also want to see the plots before they're saved
> (I pause after each plot() command). But according to the docs for
> dev.copy2pdf(), that function accepts the same args as pdf() does, including
> onefile.
> (2) I wrap my code in a function to be able to use it in the StatEt debugger
> in Eclipse.
> If you know what I'm doing wrong, or know of a different/better way, advise
> away! Thanks, Doug
> [[alternative HTML version deleted]]
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] how do I exort a list of numbers into csv file?
cat will wrap Josh's "overkill" approach in one line:
mylist<- list(1:3,3:9)
lapply(mylist , cat , sep=',' , fill=T , append=T , file='foo.csv')
Cheers
On Sun, Feb 5, 2012 at 7:42 PM, Joshua Wiley wrote:
> Hi Michael,
>
> Part of me imagines this is overkill, but this should be one option:
>
> ## your data
> mylist <- list(1:3, 3:6)
> ## open a writeable connection to a file
> con <- file("test.csv", "w")
> ## first collapse each element of the list to be a comma separated
> string, then write each
> ## element of new character vector to con using writeLines
> writeLines(sapply(mylist, paste, collapse = ", "), con = con)
> ## close the connection
> close(con)
>
> see ?writeLines for details on different ways to indicate end of the line.
>
> Hope this helps,
>
> Josh
>
> On Sun, Feb 5, 2012 at 6:01 PM, Michael wrote:
>> Hi all,
>>
>> I have a list of vector of numbers - the reason I used list of vector was
>> that I each list have different numbers of numbers which I don't know
>> before run-time.
>>
>> mylist[[1]] = c(1, 2, 3)
>> mylist[[2]] = c (3, 4, 5, 6)
>> ...
>> ...
>> etc.
>>
>> Could you please tell me if there is a way to dump all these at once into
>> csv file such that each row correspond to a vector or a cell of the list as
>> shown above?
>>
>> Thanks a lot!
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
>
> --
> Joshua Wiley
> Ph.D. Student, Health Psychology
> Programmer Analyst II, Statistical Consulting Group
> University of California, Los Angeles
> https://joshuawiley.com/
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Lattice: correct use of ltransform3dto3d to plot a surface under a cloud ?
Hello list!
I am trying to project the fitted surface to a 3d plot of the data,
similar to figures 13.7 or 6.5 in Deepayan Sarkar's "Lattice,
Multivariate Data Visualization with R", but replace the contour/map
lines with "levelplot". Problem is I can't get the color regions to
line up after the coordinate transformation. Is there a simple
solution my geometry challenged brain missed? It's been driving me
crazy for 2 days now so any help will be greatly appreciated!
I use lattice for all my other figures and would like to stay
consistent, so solutions of the form "package rgl" don't work. Thank
you all in advance.
Here is a minimal (still long) working example of what I mean, and
what I found out so far:
## make data and predicted surf
set.seed(1718)
d <- data.frame(x=runif(60),y=runif(60),g=gl(2,30))
d$z <- with(d,rnorm(60,2*x^as.numeric(g)-y^3))
d$z <- d$z+abs(min(d$z)) # so 'h' goes to the X-Y plane
surf <- by(d,d$g,function(D){
fit <- lm(z~poly(x,2)+poly(y,2),data=D)
outer(seq(0,1,l=10),seq(0,1,l=10),function(x,y,...)
predict(fit,data.frame(x=x,y=y)))
})
###
require(lattice)
# Modified code for plot 13.7 [changed: build clines from surf, -.5
for xy coords (why? don't know, works :), 3dscatter not wire]
panel.3d.contour <- function(x, y, z,rot.mat, distance,
zlim.scaled,nlevels=20,...)
{
add.line <- trellis.par.get("add.line")
clines <- contourLines(surf[[packet.number()]],nlevels = nlevels)
for (ll in clines) {
m <- ltransform3dto3d(rbind(ll$x-.5, ll$y-.5, zlim.scaled[1]),
rot.mat, distance)
panel.lines(m[1,], m[2,], col = add.line$col, lty = add.line$lty,
lwd = add.line$lwd)
}
panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled, ...)
}
cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud = panel.3d.contour,
zoom = 1,screen=list(z= 21,y=0,x=-60),aspect = c(1,1), panel.aspect = 1,
scales=list(z=list(arrows=F,tck=0)),par.box=list(lty=0),lwd=3)
# This works. But for my data the contours are messy, so I am trying
to use levelplot:
panel.3d.levels <- function(x, y, z,rot.mat, distance, zlim.scaled,...)
{
zz <- surf[[packet.number()]]
n <- nrow(zz)
s <- seq(-.5,.5,l=n)
m <- ltransform3dto3d(rbind(rep(s,n),rep(s,each=n),zlim.scaled[1]),
rot.mat, distance)
panel.levelplot(m[1,],m[2,],zz,1:n^2,col.regions=heat.colors(20))
panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled, ...)
}
cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud = panel.3d.levels,
zoom = 1,screen=list(z= 21,y=0,x=-60),aspect = c(1,1), panel.aspect = 1,
scales=list(z=list(arrows=F,tck=0)),par.box=list(lty=0),lwd=3)
# Unexpected...
# I can use panel.points for centroids and color them in "manually"
but that leaves white space or overlap:
cloud(z~x+y|g,data=d,layout=c(2,1), type='h', par.box=list(lty=0),
lwd=3, scales=list(z=list(arrows=F,tck=0)),
panel.3d.cloud = function(x, y, z,rot.mat, distance, zlim.scaled,...){
zz <- surf[[packet.number()]]
n <- nrow(zz)
s <- seq(-.5,.5,l=n)
m <- ltransform3dto3d(rbind(rep(s,n),rep(s,each=n),zlim.scaled[1]),
rot.mat, distance)
lp <- level.colors(zz, at = do.breaks(range(zz), 20),
col.regions = heat.colors(20))
panel.points(m[1,],m[2,],col=lp,pch=18,cex=2.8)
panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled =
zlim.scaled, ...)
})
#So I try to "make my own" using the lp for panel.rect, but I get the
same behavior as points for the x0,x1,y0,y1 :
panel.3d.levels <- function(x, y, z,rot.mat, distance, zlim.scaled,...)
{
zz <- surf[[packet.number()]]
n <- nrow(zz)
s <- seq(-.5,.5,l=n)
lp <- level.colors(zz, at = do.breaks(range(zz), 20),
col.regions = heat.colors(20))
cntrds <- expand.grid(s,s)
apply(cntrds,1,function(i){
xx <- i[1]+c(-.5,.5)/(n-1) ; yy <- i[2]+c(-.5,.5)/(n-1)
m <- ltransform3dto3d(rbind(xx,yy,zlim.scaled[1]), rot.mat, distance)
panel.rect(m[1,1],m[2,1],m[1,2],m[2,2])
})
panel.3dscatter(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled, ...)
}
cloud(z~x+y|g,data=d,layout=c(2,1), type='h', panel.3d.cloud = panel.3d.levels,
zoom = 1,screen=list(z= 21,y=0,x=-60),aspect = c(1,1), panel.aspect = 1,
scales=list(z=list(arrows=F,tck=0)),par.box=list(lty=0),lwd=3)
# This is as close as I got, but how to get each diagonal of
rectangles "shifted" to cover the space? I thought ltransform3dto3d
will take care of it when I transform every line in the loop. But it
didn't.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lean text label below barplot table
par(xpd=T)
bp<- barplot(matrix(1:50,5,10),names=F)
text(bp,-2,c('these labels','are tooo','looong'),pos=2,srt=45)
Xiaogang Yang wrote:
>
> Hi, everyone:
> I am plotting an graph with bar plot, but the label after every bar is too
> long, I wanna if I can draw the label lean to an angle
> thanks
>
> --
> Xiaogang Yang
> Sensorweb Research Laboratory
> http://sensorweb.vancouver.wsu.edu/
> Washington State University Vancouver
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
--
View this message in context:
http://www.nabble.com/lean-text-label-below-barplot-table-tp25336287p25426596.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
