Re: [R] x,y,z table to matrix with x as rows and y as columns
> I am sure I am missing something obvious but I cannot find the > function I am looking for. I have a data frame with three columns: X, > Y and Z, with X and Y being grid coordinates and Z the value > associated with these coordinates. I want to transform this data > frame in a matrix of Z values, on the grid defined by X and Y (and, > as a plus, fill the X.Y combinations which do no exist in the > original data frame with NAs in the resulting matrix). I could do > this manually but I guess the appropriate function should be > somewhere around. I just can't find it. Immediately after my last post I realized there was a much better solution > mydat <- expand.grid(x=1:5, y=1:5) > mydat <- data.frame(mydat, z=rnorm(25)) > mydat$z[sample(1:25,4)] <- NA > data2mat <- function(x, y, z) + { + out <- matrix(unlist(split(z, interaction(x,y))), ncol=length(unique(y))) + dimnames(out) <- list(unique(x), unique(y)) + out + } > with(mydat, data2mat(x, y, z)) Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] x,y,z table to matrix with x as rows and y as columns
jiho gmail.com> writes: > > Hello all, > > I am sure I am missing something obvious but I cannot find the > function I am looking for. I have a data frame with three columns: X, > Y and Z, with X and Y being grid coordinates and Z the value > associated with these coordinates. I want to transform this data > frame in a matrix of Z values, on the grid defined by X and Y (and, > as a plus, fill the X.Y combinations which do no exist in the > original data frame with NAs in the resulting matrix). I could do > this manually but I guess the appropriate function should be > somewhere around. I just can't find it. I have used xtabs in the past, but xtabs returns 0 instead of NA, which makes great for cross-tabulation, the "offending" line is x[is.na(x)] <- 0. So you would need to either modify the xtabs function or trust that z is never 0 and replace all 0's with NA. > mydat <- expand.grid(x=1:5, y=1:5) > mydat <- data.frame(mydat, z=rnorm(25)) > mydat$z[sample(1:25,4)] <- NA > mytab <- xtabs(z~x+y, mydat) Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] problems with character objects and calls to list()
> Really I'd like the call to list() to behave as though the text had > been entered directly so that you get > > > list(1:2, 3:4, 5:6) > [[1]] > [1] 1 2 > > [[2]] > [1] 3 4 > > [[3]] > [1] 5 6 > > eval(parse(text=paste("list(",to.convert,")",sep=""))) [[1]] [1] 1 2 [[2]] [1] 3 4 [[3]] [1] 5 6 [[4]] [1] 7 8 [[5]] [1] 9 10 [[6]] [1] 11 12 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Adding points to a wireframe with conditioning variable
Deepayan Sarkar gmail.com> writes: > > On 7/5/07, Mark Lyman gmail.com> wrote: > > I would like to add points to a wireframe but with a conditioning variable. I > > found a solution for this without a conditioning variable here, > > http://finzi.psych.upenn.edu/R/Rhelp02a/archive/65321.html. Does anyone know > > how to plot a wireframe conditioned on a variable and add the points > > conditioned on the same variable, similar to the solution at the link above? > > Depends on what form you have your points in. Inside a panel function, > packet.number() will give you the packet number in sequential order > (column major order if you think of the conditioning variables as > defining a multiway cross-tabulation), and which.packet() will give > you a vector with the current level of each conditioning variable. You > can use these to extract the appropriate subset of points. > > -Deepayan Thank you Deepayan. I modified the panel function you wrote in the post referenced above according to your suggestion. Here is the panel function, I came up with panel.3dwire.points <- function(x, y, z, xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, pts, ...) { panel.3dwire(x=x, y=y, z=z, xlim=xlim, ylim=ylim, zlim=zlim, xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled, zlim.scaled=zlim.scaled, ...) pts.sub <- subset(pts, g==levels(g)[which.packet()]) xx <- xlim.scaled[1] + diff(xlim.scaled)*(pts.sub$x - xlim[1])/diff (xlim) yy <- ylim.scaled[1] + diff(ylim.scaled)*(pts.sub$y - ylim[1])/diff (ylim) zz <- zlim.scaled[1] + diff(zlim.scaled)*(pts.sub$z - zlim[1])/diff (zlim) panel.3dscatter(x=xx, y=yy, z=zz, xlim=xlim, ylim=ylim, zlim=zlim, xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled, zlim.scaled=zlim.scaled, ...) } where pts is a dataframe with the three dimensional points (x,y,z) and the conditioning variable g. Mark __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Adding points to a wireframe with conditioning variable
I would like to add points to a wireframe but with a conditioning variable. I found a solution for this without a conditioning variable here, http://finzi.psych.upenn.edu/R/Rhelp02a/archive/65321.html. Does anyone know how to plot a wireframe conditioned on a variable and add the points conditioned on the same variable, similar to the solution at the link above? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] barchart (lattice) with text labels
Deepayan Sarkar wrote: > On 2/24/07, Mark and Heather Lyman <[EMAIL PROTECTED]> wrote: >> I would like to place the value for each bar in barchart (lattice) at >> the top of each bar. Something like the following code produces. >> >> library(lattice) >> >> mypanelfunc <- function(x, y, ...) >> { >> panel.barchart(x, y, ...) >> panel.text(x, y, labels=as.character(round(x,2)), ...) >> } >> >> myprepanelfunc <- function(x, y, ...) list(xlim=c(0, max(x)+.1)) >> >> mydata <- expand.grid(a=factor(1:5), b=factor(1:3), c=factor(1:2)) >> mydata$x <- runif(nrow(mydata)) >> >> barchart(a~x|b, mydata, groups=c, panel=mypanelfunc, >> prepanel=myprepanelfunc, adj=c(-0.1,0.5)) >> >> However, I cannot figure out how to shift the values to correspond with >> their respective grouped bar. > > You should look at panel.barchart and try to reproduce the > calculations done there. > > Deepayan > This is the panel function that I ended up using. As Deepayan suggested, I borrowed heavily from panel.barchart. mypanelfunc <- function(x, y, groups, box.ratio, ...) { panel.barchart(x, y, groups=groups, box.ratio=box.ratio, ...) origin <- current.panel.limits()$xlim[1] groupSub <- function(groups, subscripts, ...) groups[subscripts] groups <- as.numeric(groupSub(groups, ...)) vals <- sort(unique(groups)) nvals <- length(vals) height <- box.ratio/(1+ nvals * box.ratio) for (i in unique(y)) { ok <- y == levels(y)[i] nok <- sum(ok, na.rm = TRUE) panel.text(x = x[ok], y = (i + height * (groups[ok] - (nvals + 1)/2)), labels=as.character(signif(x[ok],2)), ...) } } Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] sapply again return value
Antje yahoo.de> writes: > Hello, > > I used an sapply to get some data back (s <- sapply(...) ). The output > of s would then deliver something like this: > > B06_lamp.csv C06_lamp.csv D06_lamp.csv > [1,] NULL NULL Numeric,512 > [2,] NULL NULL Numeric,512 > [3,] NULL NULL 2 > > mode(s) > [1] "list" > > dim(s) > [1] 3 3 > > > > Now, I'd like to remove the columns which contain NULL (it's alway the > whole column). > How can I do this??? > > Antje As long as it is always the whole column, you can just test the first row, like this: > s<-matrix(list(NULL,NULL,NULL,NULL,NULL,NULL,numeric(512),numeric (512),2),ncol=3) > s[,!apply(s,2,sapply,is.null)[1,],drop=FALSE] [,1] [1,] Numeric,512 [2,] Numeric,512 [3,] 2 If you would like to make sure that the whole column is NULL, you could do something like the following: > s[,!apply(apply(s,2,sapply,is.null),2,sum)==nrow(s),drop=FALSE] [,1] [1,] Numeric,512 [2,] Numeric,512 [3,] 2 I use the "drop=FALSE" here for display purposes, but you may want to leave it off depending on desired format. Mark __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] plotting derived values not raw
James Root gmail.com> writes: > > I am trying to plot the mean and standard error of three separate > conditions. For various reasons, I do not have access to the raw data from > which the mean and error were derived and would like to make error bar plots > utilizing only the actual mean and standard error values. Is there a way to > do this in R? Thanks for any help in advance. > james I believe that xYplot in the Hmisc package will do what you want Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using variable names in for loops - Generating plot s semi-automatically from a series of variables Partly solved
Anthony Staines gmail.com> writes: For each > disease, I want to pull out the column of data containing the word > 'Male' and plot this against age, and then add a line to the plot for > the corresponding column containing 'Female'. > -- > attach(data) > > Diseases <- c("Cardiovascular disease","Road Traffic Injury", ... > ,"All causes") > Male <- names(data)[grep("Male",names(data))] > Female <- names(data)[grep("Female",names(data))] > #Disease contains disease labels in the correct order, and Male and > Female now hold the (correct) variables. > > for (i in seq(1,length(Diseases))) > { > jpeg(paste(Diseases[i],".jpg")) #This works fine! > plot(Male[i]~Age) #This does not > lines(Female[i]~Age) > dev.off() > } > detach(data) [SNIP] There are a few ways you can do this. Using plot.formula like you do here you can use this: > mydata <- as.data.frame(matrix(runif(300),30,10)) > attach(mydata) > plot(formula(paste(mynames[1],mynames[2],sep='~'))) A way to do this without a formula would be: > plot(eval(parse(text=mynames[1])),eval(parse(text=mynames[2]))) > detach(mydata) Take a look at the help files for eval and parse. I still do not have a firm grasp on how to use them and other related functions, like substitute, but what I have been able figure out has been very useful. Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] LM Model
Simon P. Kempf web.de> writes: [SNIP] > > But I want to skip the lm function and specify my own regression equation > RENT= 15 -0.15*AGE1 and then use the predict.lm function. However, in order > to use the predict.lm function I need an object of class lm. Is there any > way to do so? Or maybe somebody has another solution? > > Thanks in advance, > > Simon Here is one way. Take a look at help(offset), help(lm), and help(lm.predict). > xx <- runif(30) > yy <- rnorm(30) > mydata<-data.frame(xx,yy) > lm(yy~offset(15*rep(1,30))+offset(-0.15*xx)-1,mydata) Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem with fitted.value function
I don't know how to use the fitted > values function with a given function and given input-variables but yet > unknown result-values. Take a look at the predict function, ?predict. Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] applying cor.test in two different but matched dataframes
Gabriele Stocco units.it> writes: > > Dear R help, > I have two dataframes with the same number of rows and columns. Each > column is for a patient, each row for a variable. The columns and rows > are matched in the two dataframes, so that each patient (column) and > each variable (row) has the same position in the two dataframes. The > values are different since each dataframe reports data from a > different way to measure the same variables in the sampe patients. > > How can I make a cor.test, applying the function row by row in the two > different dataframes, possibly without merging them? > > Thanks, > > Gabriele Stocco > University of Trieste Try: # Create Data Frames With Rows as Variables x<-as.data.frame(t(as.data.frame(matrix(rnorm(100),10 y<-as.data.frame(t(as.data.frame(matrix(rnorm(100),10 # Convert to Lists So that Each Element of List is a variable t.x<-as.list(as.data.frame(t(x))) t.y<-as.list(as.data.frame(t(y))) # Use mapply to apply function to the two lists; see ?mapply mapply(cor.test,t.x,t.y) Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] update index in
Kim Milferstedt uiuc.edu> writes: > > Hello, > > I have a time series of data as a data.frame. Occasionally there is > one or more days missing (e.g. data available for days 2, 3, 4, 8, 9, > 10 --> missing days between 4 and 8). The experimental time > information can be found in the 2nd column of "data". I would like to > have a continuous time line with one time point per day. Therefore I > try to insert lines for the missing days that contain zeros for the > data categories just to fill the columns. > > > Thanks already, > > Kim > I believe this will also do what you want: > days<-c(1:10)[-5:-7] > xx<-rnorm(7) > data<-data.frame(xx,days) > new.data<-merge(data,data.frame(days=1:10),all.y=TRUE) It usually is not a good idea to use zeroes as placeholders for missing values. Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Plotmath expression
Christos Hatzis nuverabio.com> writes: > > Hello, > > I've been trying to plot a subscript in a text formula using plotmath but I > haven't been able to do so. > > In my example below I would like the text label to show > X[min] = 10.1 +/- 5.5 > > Here is the code: > > ll <- c(x=10.1, sde=5.5) > plot(1:10) > text(x=9, y=2, pos=2, expression(paste(X[min], "=", paste(ll, > collapse="+/-" How about the following? text(x=9, y=2, pos=2,substitute(X[min]==x%+-%sde,as.list(ll))) Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] write data to pdf
Franco Mendolia gmx.de> writes: > > Hello! > > Is there a possibility in R to save data in pdf-format? > I do not want to save a plot but some lines of simple text. > > Regards, > > Franco Mendolia You could also use pdf() and textplot() in the gplots package Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How can I delete components in a column ?
Yen Ngo yahoo.se> writes: > > Hi all R-helpers, > > i am a new R-user and have problem with deleting some components in a column. I have a dataset like > > Name Idx >empty 2 >empty 3 > anone2 > bnone3 > d none 2 > ad cfh 4 > bf cdt 5 >empty 2 >empty 2 > gf cdh 4 > d none 5 > > and want to eliminate all components that have id=none and empty . The remaining data should be Take a look at the subset function. Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Incompatability of lme4 and Matrix packages
I just installed R 2.4 on a windows machine and installed the lme4 package version 0.9975-1. Unfortunately this version requires the 0.9975-2 version of the Matrix package which I cannot seem to find. It seems the newest version I can find is 0.9975-1. Is there somewhere I could find this newer version. Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Help on apply() function
Anusha hotmail.com> writes: > > Respected Sir/Madam, > > I have a problem with apply function. I have to two matrices of dimension of > one column but n rows. I have to check whether one matrix is greater than > other by going thru each row (ie) using if condition to check one matrix > with another matrix. > > I like to use apply() function to this approach. That is apply function > between two matrices. I searched for examples online but I couldn't find > any. > > I don't know how to loop thru the matrices. > > Please help with this problem. > > Any help is appreciated. > > My mailid is monkponu yahoo.com. > > Thanks, > Anusha. > You can use the functions all.equal with the function isTRUE (see ?all.equal) to check if two objects are nearly equal (within a certain tolerance). Or you can use identical (see ?identical) to check if they are exactly the same. See the examples in the help for identical. Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Type II and III sum of square in Anova (R, car package)
> 1. First of all, more general question. Standard anova() function for lm() > or aov() models in R implements Type I sum of squares (sequential), which > is not well suited for unbalanced ANOVA. Therefore it is better to use > Anova() function from car package, which was programmed by John Fox to use > Type II and Type III sum of squares. Did I get the point? > > 2. Now more specific question. Type II sum of squares is not well suited > for unbalanced ANOVA designs too (as stated in STATISTICA help), therefore > the general rule of thumb is to use Anova() function using Type II SS > only for balanced ANOVA and Anova() function using Type III SS for > unbalanced ANOVA? Is this correct interpretation? > > 3. I have found a post from John Fox in which he wrote that Type III SS > could be misleading in case someone use some contrasts. What is this about? > Could you please advice, when it is appropriate to use Type II and when > Type III SS? I do not use contrasts for comparisons, just general ANOVA > with subsequent Tukey post-hoc comparisons. There are many threads on this list that discuss this issue. Not being a great statistician myself, I would suggest you read through some of these as a start. As I understand, the best philosophy with regards to types of sums of squares is to use the type that tests the hypothesis you want. They were developed as a convenience to test many of the hypotheses a person might want "automatically," and put it into a nice, neat little table. However, with an interactive system like R it is usually even easier to test a full model vs. a reduced model. That is if I want to test the significance of an interaction, I would use anova(lm.fit2, lm.fit1) where lm.fit2 contains the interaction and lm.fit2 does not. The anova function will return the appropriate F-test. The danger with worrying about what type of sums of squares to use is that often we do not think about what hypotheses we are testing and if those hypotheses make sense in our situation. Mark Lyman __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to iteratively extract elements out of a list
lapply(m,function(x)x[x>2]) [[1]] [1] 3 4 [[2]] [1] 4 5 [[3]] [1] 4 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Specify Z matrix with lmer function
Is there a way to specify a Z matrix using the lmer function, where the model is written as y = X*Beta + Z*u + e? I am trying to reproduce smoothing methods illustrated in the paper "Smoothing with Mixed Model Software" my Long Ngo and M.P. Wand. published in the /Journal of Statistical Software/ in 2004 using the lme4 and Matrix packages. The code and data sets used can be found at http://www.jstatsoft.org/v09/i01/. There original code did not work for me without slight modifications here is the code that I used with my modifications noted. x <- fossil$age y <- 10*fossil$strontium.ratio knots <- seq(94,121,length=25) n <- length(x) X <- cbind(rep(1,n),x) Z <- outer(x,knots,"-") Z <- Z*(Z>0) # I had to create the groupedData object with one group to fit the model I wanted grp <- rep(1,n) grp.dat<-groupedData(y~Z|grp) fit <- lme(y~-1+X,random=pdIdent(~-1+Z),data=grp.dat) I would like to know how I could fit this same model using the lmer function. Specifically can I specify a Z matrix in the same way as I do above in lme? Thanks, Mark __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Post Hoc Groupings
Looking at the errors your code produces, it looks like you need to make Dock and Slip factors. dock_2004_data$Dockf<-factor(dock_2004_data$Dock) dock_2004_data$Slipf<-factor(dock_2004_data$Slip) rich.aov <- aov(X.open ~ Dockf*Slipf, data=dock_2004_data) TukeyHSD(rich.aov, c("Dockf", "Slipf")) Jarrett Byrnes wrote: >Indeed, the following works as well >On Oct 26, 2005, at 5:23 PM, P Ehlers wrote: > > > >>fm1 <- aov(breaks ~ wool*tension, data = warpbreaks) >>TukeyHSD(fm1, c("wool","tension", "wool:tension")) >> >> > >However, when working with my own dataset, I get the following errors. > I have some inkling this may be due to a slightly unbalanced sample >size, but am uncertain of this. > > > rich.aov <- aov(X.open ~ Dock*Slip, data=dock_2004_data) > > TukeyHSD(rich.aov, c("Dock", "Slip")) > >Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep' >In addition: Warning messages: >1: non-factors ignored: Slip in: replications(paste("~", xx), data = mf) >2: non-factors ignored: Dock, Slip in: replications(paste("~", xx), >data = mf) > >I am pleased to know that these errors are not quite the > >__ >R-help@stat.math.ethz.ch mailing list >https://stat.ethz.ch/mailman/listinfo/r-help >PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > > > > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Bug in lmer?
I am using version 0.98-7 of the Matrix package. I used the RGui "Install Packages..." menu option to get the lme4 package from CRAN and this version of the Matrix was automatically downloaded as well. Martin Maechler wrote: >>>>>>"Mark" == Mark Lyman <[EMAIL PROTECTED]> >>>>>>on Thu, 29 Sep 2005 14:44:38 -0600 writes: >>>>>> >>>>>> > >Mark> I am relatively new to R so I am not confident enough in what I am > doing >Mark> to be certain this is a bug. > > >Mark> I am running R 2.1.1 on a Windows XP >Mark> machine and the lme4 package version 0.98-1. > >lme4 nowadays is heavily based on "Matrix" which version are you >using there? > >Mark> The following code fits the model I want using the >Mark> nlme package version 3.1-60. > > < .. > {see a script at the end} > >Mark> The problem is that when I try fitting the model using >Mark> the lmer function with the following code: > >Mark> lmer(adg~trt+(1|loc)+(1|block:loc)+(1|loc:trt),mltloc) > >Mark> I get this message from Windows and R closes. > >Mark> >> R for Windows GUI front-end has encountered a problem and needs > to >Mark> >> close. We are sorry for the inconvenience. > >That definitely means there is a bug. >The question is *where* the bug is: "lme4", "Matrix", "R", "Windows". > >One first thing coming to mind is a mismatch of "lme4" and "Matri > >Mark> This same code works on a Macintosh. So it doesn't >Mark> seem that I have made an error in my code. > >correct; I can also fit the model nice and quickly on Linux, >and summary() confirms the same fit {with the "usual problem" of >different estimates for the degrees of freedoms 'df'}. > >So currently the bug only shows on the Windows platform. >Could it be that you have a mismatching package "Matrix" version >there, but not on the Mac? > >Mark> Also if anyone of the random effect terms is removed there is >Mark> no problem. Is this something that is being looked at? > >not yet, AFAIK. > >Mark> Or I have I made a mistake somewhere? I have included >Mark> the data that I am using below. > >I'm putting the data and an R script up for FTP, >so that you or others can run this ``from anywhere'' via > > source("ftp://stat.ethz.ch/U/maechler/R/mltloc-ex.R";, echo = TRUE) > >Maybe this helps diagnosis, >Martin Maechler > > > > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Bug in lmer?
I am relatively new to R so I am not confident enough in what I am doing to be certain this is a bug. I am running R 2.1.1 on a Windows XP machine and the lme4 package version 0.98-1. The following code fits the model I want using the nlme package version 3.1-60. mltloc$loc <- factor(mltloc$loc) mltloc$block <- factor(mltloc$block) mltloc$trt <- factor(mltloc$trt) Mltloc <- groupedData(adg~trt|loc,mltloc) plot(Mltloc) plot(Mltloc,outer=~trt) lme(adg~trt, random=pdBlocked(list(pdIdent(~1),pdIdent(~block-1),pdIdent(~trt-1))),Mltloc) The problem is that when I try fitting the model using the lmer function with the following code: lmer(adg~trt+(1|loc)+(1|block:loc)+(1|loc:trt),mltloc) I get this message from Windows and R closes. R for Windows GUI front-end has encountered a problem and needs to close. We are sorry for the inconvenience. This same code works on a Macintosh. So it doesn't seem that I have made an error in my code. Also if anyone of the random effect terms is removed there is no problem. Is this something that is being looked at? Or I have I made a mistake somewhere? I have included the data that I am using below. X obs loc block trt adg fe 1 1 3 A 1 3 3.16454 7.1041 2 2 4 A 1 4 3.12500 6.6847 3 3 6 A 1 2 3.15944 6.8338 4 4 7 A 1 1 3.25000 6.5254 5 5 9 A 2 2 2.71301 8.2505 6 6 10 A 2 1 3.20281 7.5922 7 7 12 A 2 3 3.02423 7.3894 8 8 16 A 2 4 2.87245 7.4604 9 9 19 A 3 1 2.68878 8.2785 10 10 20 A 3 2 2.86862 7.9470 11 11 21 A 3 3 2.89923 7.9739 12 12 22 A 3 4 3.02806 8.4331 13 13 25 B 1 3 2.18131 6.6691 14 14 27 B 1 4 2.51914 5.6281 15 15 28 B 1 2 1.88739 7.0723 16 16 31 B 1 1 2.34685 6.0295 17 17 33 B 2 4 2.45608 5.6195 18 18 35 B 2 1 2.25225 6.3978 19 19 36 B 2 3 2.23649 6.1799 20 20 40 B 2 2 2.47523 5.9985 21 21 41 B 3 2 1.94200 7.2975 22 22 44 B 3 3 2.43243 6.4350 23 23 47 B 3 4 2.30180 6.3339 24 24 48 B 3 1 2.53378 6.1564 25 25 50 C 1 4 2.96014 7.5110 26 26 51 C 1 2 3.23551 7.4762 27 27 54 C 1 3 3.24638 7.2063 28 28 56 C 1 1 3.04710 7.6389 29 29 58 C 2 3 3.26449 7.5466 30 30 59 C 2 2 2.71377 9.0895 31 31 61 C 2 1 3.06522 7.8723 32 32 62 C 2 4 2.71739 8.2318 33 33 66 C 3 4 3.03623 7.9426 34 34 67 C 3 3 3.10507 8.4608 35 35 69 C 3 1 3.16304 8.5549 36 36 70 C 3 2 3.02899 8.5038 37 37 74 D 1 1 2.49164 9.5758 38 38 77 D 1 3 2.51833 9.5121 39 39 79 D 1 2 2.35631 10.3264 40 40 80 D 1 4 2.30331 9.7715 41 41 81 D 2 3 2.72688 9.5628 42 42 83 D 2 2 2.59512 9.9414 43 43 85 D 2 1 2.56516 9.3887 44 44 88 D 2 4 2.91523 8.3158 45 45 89 D 3 3 2.57943 10.4416 46 46 90 D 3 4 2.98159 8.7710 47 47 93 D 3 2 2.35370 11.0148 48 48 94 D 3 1 2.21953 11.2417 49 49 99 E 1 3 2.84158 8.7886 50 50 100 E 1 4 2.65264 8.6946 51 51 102 E 1 2 2.47112 9.7143 52 52 103 E 1 1 2.89769 9.2401 53 53 105 E 2 2 2.57343 9.5353 54 54 106 E 2 1 2.99752 8.7538 55 55 110 E 2 4 2.95380 8.8210 56 56 112 E 2 3 3.08663 8.9427 57 57 114 E 3 1 2.72525 9.4308 58 58 115 E 3 2 2.75825 9.7721 59 59 116 E 3 3 3.08333 8.9010 60 60 117 E 3 4 3.12129 8.4852 61 61 122 F 1 1 3.20600 6.3983 62 62 123 F 1 2 2.89500 6.6569 63 63 125 F 1 4 3.36900 6.0821 64 64 126 F 1 3 3.12000 6.5349 65 65 130 F 2 2 3.19300 6.6729 66 66 131 F 2 1 3.29800 6.5488 67 67 133 F 2 4 3.09700 6.6598 68 68 135 F 2 3 3.38500 6.2998 69 69 139 F 3 3 3.44900 6.2849 70 70 140 F 3 2 3.05000 6.9957 71 71 141 F 3 1 3.43500 6.7302 72 72 143 F 3 4 3.60600 6.3827 73 73 145 G 1 2 2.58669 8.1394 74 74 147 G 1 1 3.17892 7.0972 75 75 148 G 1 4 2.95284 7.3140 76 76 151 G 1 3 3.17924 6.9430 77 77 154 G 2 4 2.62344 7.5150 78 78 155 G 2 3 2.64286 8.0237 79 79 157 G 2 1 3.12760 7.3169 80 80 160 G 2 2 2.54993 8.1957 81 81 163 G 3 4 2.58322 7.9687 82 82 164 G 3 3 2.84813 7.9284 83 83 166 G 3 2 2.69279 8.5303 84 84 167 G 3 1 3.14424 7.3564 85 85 169 H 1 3 3.39974 6.5945 86 86 173 H 1 2 3.12370 6.7530 87 87 175 H 1 4 3.17969 6.4279 88 88 176 H 1 1 3.70052 6.4830 89 89 177 H 2 2 2.95192 7.3809 90 90 179 H 2 3 3.44661 6.7929 91 91 182 H
[R] lmer random effect model matrix question
I have one fixed effect, sor, with two levels. I have eight lots and three wafers from each lot. I have included the data below. I would like to fit a mixed model that estimates a covariance parameter for wafer, which is nested in lot, and two covariance parameters for lot, one for each level of sor. The following command fits the model that I want, except for it estimates the correlation between the two covariance parameters for lot. Is there anyway to make R not estimate this correlation? Thank you. lmer(y~sor+(sor-1|lot)+(1|wafer:lot),wafer) For those familiar with proc mixed the following SAS code fits the model that I want: proc mixed scoring=4; class sor lot wafer site; model y= sor/ddfm=satterth; random lot(sor)/group=sor; random wafer(lot); run; sor lot wafer sitey 11 1 11 2006 21 1 12 1999 31 1 13 2007 41 1 21 1980 51 1 22 1988 61 1 23 1982 71 1 31 2000 81 1 32 1998 91 1 33 2007 10 1 2 11 1991 11 1 2 12 1990 12 1 2 13 1988 13 1 2 21 1987 14 1 2 22 1989 15 1 2 23 1988 16 1 2 31 1985 17 1 2 32 1983 18 1 2 33 1989 19 1 3 11 2000 20 1 3 12 2004 21 1 3 13 2004 22 1 3 21 2001 23 1 3 22 1996 24 1 3 23 2004 25 1 3 31 1999 26 1 3 32 2000 27 1 3 33 2002 28 1 4 11 1997 29 1 4 12 1994 30 1 4 13 1996 31 1 4 21 1996 32 1 4 22 2000 33 1 4 23 2002 34 1 4 31 1987 35 1 4 32 1990 36 1 4 33 1995 37 2 5 11 2013 38 2 5 12 2004 39 2 5 13 2009 40 2 5 21 2023 41 2 5 22 2018 42 2 5 23 2010 43 2 5 31 2020 44 2 5 32 2023 45 2 5 33 2015 46 2 6 11 2032 47 2 6 12 2036 48 2 6 13 2030 49 2 6 21 2018 50 2 6 22 2022 51 2 6 23 2026 52 2 6 31 2009 53 2 6 32 2010 54 2 6 33 2011 55 2 7 11 1984 56 2 7 12 1993 57 2 7 13 1993 58 2 7 21 1992 59 2 7 22 1992 60 2 7 23 1990 61 2 7 31 1996 62 2 7 32 1993 63 2 7 33 1987 64 2 8 11 1996 65 2 8 12 1989 66 2 8 13 1996 67 2 8 21 1997 68 2 8 22 1993 69 2 8 23 1996 70 2 8 31 1990 71 2 8 32 1989 72 2 8 33 1992 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html