Re: [R] Direction of panel plots in trellis graphics
On 13 Jul 2007, at 17:11, [EMAIL PROTECTED] wrote: > > Another high level option is to change the rule determining how > packets are chosen for a given panel in the layout. > > print(tmp.tr3, > packet.panel = function(layout, row, column, ...) { > layout <- layout[c(2, 1, 3)] > packet.panel.default(layout = layout, > row = column, > column = row, ...) > }) > > This effectively transposes the layout, which (along with > as.table=TRUE) is what you want. Thank you Deepayan, that does exactly what I want. Also thanks for the other suggestions by Gabor and Richard. Unfortunately these don't quite work in my case, as I am printing to multiple pages (so Richard's suggestion of transposing becomes tricky), and my 2 conditioning variables are on the rhs and lhs of the formula (i.e. y1 + y2 = x | v) (so I can't easily create a new combination factor as in Gabor's situation). Cheers Yan __ 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] Direction of panel plots in trellis graphics
Hi, Using library(lattice), is there any way to tell xyplot to plot panels top to bottom, then left to right (i.e. panels are appended vertically, then horizontally). as.table changes the plot direction from left-to-right then top-to-bottom, to right-to-left then bottom- to-top, but that's not quite what I want to do. Thanks Yan __ 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] Suggestions for help & weighted.mean
On 11 Aug 2006, at 12:49, Duncan Murdoch wrote: > It makes sense in this case, but in the case where there is more > than one infinite weight, the result has to be NaN. > ... it would be a lot more complicated if it were to handle this > very special case. Yes - I see that it may not be worth slowing down the code to cope with this one particular case. I suppose it really comes down to a question of completeness versus speed. > On the other hand, if you know that in your situation there is at > most one infinite weight, then you could use > > if (any(inf <- is.infinite(w))) x[inf] > else weighted.mean(x, w) Thanks. I think I do something like that already, but your code is cleaner than mine! >> p.s. a while ago I suggested using '??xxx' as a shortcut for >> help.search("xxx"), much like '?xxx' is a shortcut for help >> ("xxx"). I was just wondering if anyone had any more thoughts on >> the matter? > Suggestions like this (and probably the one above) belong more in > the R-devel list than here. OK. Thanks. > I think your ?? suggestion is reasonable; why don't you write up > the necessary patch to implement it, and see if it's feasible? > Include that in your post to R-devel, and it will be easier for > others to see the pitfalls (if there are any). Great. I'll do that when I have time (and if I can work out how the codebase works), then try posting it to R-devel. Cheers Yan __ 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] Suggestions for help & weighted.mean
Hi, just a quick question: Should weighted.mean be able to cope with the specific case where one weight is Inf? I came across this when trying to implement a simple weighted moving average algorithm for surface smoothing: these algorithms often result in a single infinite weight when predicting the surface value at known data points. e.g. > weighted.mean(c(77,88,99), c(Inf, 1, 2)) #should this return 77? [1] NaN Cheers Yan p.s. a while ago I suggested using '??xxx' as a shortcut for help.search("xxx"), much like '?xxx' is a shortcut for help("xxx"). I was just wondering if anyone had any more thoughts on the matter? __ 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] extractAIC using surf.ls
On 11 Aug 2006, at 08:17, Prof Brian Ripley wrote: > Roger, thank you for looking into this. Yes, and thanks to both of you for the corrections. > However, the posting guide asked the poster to contact the > maintainer. Had > (s)he done so, I would have pointed out that spatial 7.28-2 (the > current > version for R-devel) has this corrected (in a slightly simpler way). Thanks. There have been a few times when what I thought was a bug was due to a misunderstanding on my part. It seems better to me to check on the R-help list that it really is a bug before bothering the maintainers. On that note, I see that Prof. Ripley is the author of the loess package. When trying to adjust some of the control parameters for a loess fit, I tried the following (wrong) incantation. loess(dist ~ speed, cars, control = list(statistics = "exact")) Although it is wrong (should be control = loess.control(statistics = "exact")), entering this as a command actually crashes R on the 2 systems I have tried (OS X, linux, both R-2.3.1), with the error below. I'm not sure if you consider this a bug, as the command I typed was invalid. --- *** caught segfault *** address (nil), cause 'memory not mapped' Traceback: 1: .C(R_loess_raw, as.double(y), as.double(x), as.double (weights), as.double(robust), as.integer(D), as.integer(N), as.double(span), as.integer(degree), as.integer(nonparametric), as.integer(order.drop.sqr), as.integer(sum.drop.sqr), as.double (span * cell), as.character(surf.stat), fitted.values = double (N), parameter = integer(7), a = integer(max.kd), xi = double (max.kd), vert = double(2 * D), vval = double((D + 1) * max.kd), diagonal = double(N), trL = double(1), delta1 = double (1), delta2 = double(1), as.integer(surf.stat == "interpolate/ exact")) 2: simpleLoess(y, x, w, span, degree, parametric, drop.square, normalize, control$statistics, control$surface, control$cell, iterations, control$trace.hat) 3: loess(dist ~ speed, cars, control = list(statistics = "exact")) Possible actions: 1: abort (with core dump) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace __ 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] extractAIC using surf.ls
Although the 'spatial' documentation doesn't mention that extractAIC works, it does seem to give an output. I may have misunderstood, but shouldn't the following give at least the same d.f.? > library(spatial) > data(topo, package="MASS") > extractAIC(surf.ls(2, topo)) [1] 46. 437.5059 > extractAIC(lm(z ~ x+I(x^2)+y+I(y^2)+x:y, topo)) [1] 6. 357.5059 # and if the AIC values differ, shouldn't they do so by an additive constant? > (extractAIC(surf.ls(2, topo))-extractAIC(lm(z ~ x+I(x^2)+y+I(y^2) +x:y, topo)))[2] [1] 80 > (extractAIC(surf.ls(1, topo))-extractAIC(lm(z ~ x+y, topo)))[2] [1] 92 # Using R 2.3.1 (OS X), spatial version 7.2-27.1 Thanks Yan __ 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] Slight fault in error messages
Just a quick point which may be easy to correct. Whilst typing the wrong thing into R 2.2.1, I noticed the following error messages, which seem to have some stray quotation marks and commas in the list of available families. Perhaps they have been corrected in the latest version (sorry, I don't want to upgrade yet, but it should be easy to check)? > glm(1 ~ 2, family=quasibinomial(link=foo)) Error in quasibinomial(link = foo) : ‘foo’ link not available for quasibinomial family, available links are "logit", ", ""probit" and "cloglog" > glm(1 ~ 2, family=binomial(link=foo)) Error in binomial(link = foo) : link "foo" not available for binomial family, available links are "logit", ""probit", "cloglog", "cauchit" and "log" I hope this is helpful, Yan __ 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] multiple plots with par mfg
On 23 May 2006, at 21:47, Romain Francois wrote: > Hi, > > An other possibility might be to use two devices and use dev.set to > go from one to another : Thanks. Actually I did try that, but there are quite a lot of points to plot, and the switching between plots slowed the whole simulation down a bit. Thanks too to Henrik Bengtsson for the code snippet. I've managed to get it working now, using par(mfg). Yan __ 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] multiple plots with par mfg
On 23 May 2006, at 15:57, Greg Snow wrote: > The best thing to do is to create the first plot, add everything to > the > first plot that you need to, then go on to the 2nd plot, etc. Yes, I realise that. The problem is that the data are being simulated on the fly, and I wish to display multiple plots which are updated as the simulation progresses. So I do need to return to each plot on every generation of the simulation. > If you > really need to go back to the first plot to add things after plotting > the 2nd plot then here are a couple of ideas: > > Look at the examples for the cnvrt.coords function in the > TeachingDemos > package (my quick test showed they work with layout as well as > par(mfrow=...)). > > The other option is when you use par(mfg) to go back to a previous > plot > you also need to reset the usr coordinates, for example: Aha. I didn't realise that the usr coordinates could be stored and reset using par. > Hope this helps, I think that's exactly what I need. Thank you very much. Yan __ 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] multiple plots with par mfg
On 23 May 2006, at 12:48, Prof Brian Ripley wrote: > On Tue, 23 May 2006, Yan Wong wrote: >> if not, can anyone suggest a way of appending to 2 >> separate plots on the fly. > > No, it is user error. par(mfg=) specifies where the next figure > will the drawn, and points() does not draw a figure but adds to > one. As the help page says: > > 'mfg' A numerical vector of the form 'c(i, j)' where 'i' and 'j' > indicate which figure in an array of figures is to be drawn > next (if setting) or is being drawn (if enquiring). > OK. I didn't appreciate the distinction between drawing and adding. > You need to use screen() or layout() to switch back to an existing > plot. Thanks, but the help page for screen says "The behavior associated with returning to a screen to add to an existing plot is unpredictable and may result in problems that are not readily visible." I assume this to mean that I shouldn't do it using screen(). I can't find any description of how to add points to several different plots generated after a layout() call. Is there a way? Cheers Yan __ 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] multiple plots with par mfg
Hi, I'm trying to add points to 2 plots on the fly using par(mfg=vector) so switch between them. However, the appropriate scales aren't switched when changing from one plot to another, e.g. par(mfcol=c(2,1)) plot(1,1, col="blue")# blue plot plot(1.2,1.2, col="red") # red plot points(1.1,1.1) # appears to bottom left of red point par(mfg=c(1,1)) # switch plots points(1.1,1.1) # should appear at top of blue point, but appears as on red plot # > version _ platform powerpc-apple-darwin7.9.0 arch powerpc os darwin7.9.0 system powerpc, darwin7.9.0 status Patched major2 minor2.1 year 2006 month03 day 02 svn rev 37488 language R Is this a bug? if not, can anyone suggest a way of appending to 2 separate plots on the fly. Thanks Yan Wong Leeds __ 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] Drop1 and weights
On 2 Mar 2006, at 17:56, Prof Brian Ripley wrote: > On Thu, 2 Mar 2006, Yan Wong wrote: > >> >> Indeed, I realised this, but assumed that the problem could be >> understood even without my dataset. I'll test my data on the patched >> version when it becomes available, and re-post then. I've just tried R-patched, and it now works as expected. Thanks for the quick fix. As you point out, the glm and lm versions differ in AIC reported, but only by an additive constant, and so the tests are not affected. Yan > weighted.lm <- lm(trSex ~ (river+length +depth)^2-length:depth, dno2, weights=males+females) > d1.lm <- drop1(weighted.lm, test="F"); > weighted.glm <- glm(trSex ~ (river+length +depth)^2-length:depth, dno2, weights=males+females, family="gaussian") > d1.glm <- drop1(weighted.glm, test="F"); > > #differ by a constant > d1.lm$AIC - d1.glm$AIC [1] 628.1381 628.1381 628.1381 > > Anova(weighted.lm) Anova Table (Type II tests) Response: trSex Sum Sq Df F value Pr(>F) river 1.471 1 3.3775 0.06843 . length1.002 1 2.2999 0.13187 depth 2.974 1 6.8295 0.01005 * river:length 0.075 1 0.1733 0.67790 river:depth 2.020 1 4.6397 0.03313 * Residuals55.303 127 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > d1.lm Single term deletions Model: trSex ~ (river + length + depth)^2 - length:depth Df Sum of Sq RSS AIC F value Pr(F) 55.303 -104.710 river:length 1 0.075 55.379 -106.529 0.1733 0.67790 river:depth 1 2.020 57.324 -101.938 4.6397 0.03313 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > d1.glm Single term deletions Model: trSex ~ (river + length + depth)^2 - length:depth Df Deviance AIC F value Pr(F) 55.30 -732.85 river:length 155.38 -734.67 0.1733 0.67790 river:depth 157.32 -730.08 4.6397 0.03313 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 __ 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] Drop1 and weights
On 2 Mar 2006, at 13:25, Prof Brian Ripley wrote: > It looks like at some point weighted.residuals() was changed, and > broke this. The models are fitted correctly, but AIC is wrong. > However, note that it does work correctly for a glm() fit which > gives you a workaround. Thank you. > You example is not reproducible, so I was unable to test the > workaround nor the correction which will short;y be in R-patched. Indeed, I realised this, but assumed that the problem could be understood even without my dataset. I'll test my data on the patched version when it becomes available, and re-post then. Thanks again Yan __ 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] Drop1 and weights
Hi, If I used drop1 in a weighted lm fit, it seems to ignore the weights in the AIC calculation of the dropped terms, see the example below. Can this be right? Yan library(car) > unweighted.model <- lm(trSex ~ (river+length +depth)^2- length:depth, dno2) > Anova(unweighted.model) Anova Table (Type II tests) Response: trSex Sum Sq Df F value Pr(>F) river0.1544 1 3.3151 0.07100 . length 0.1087 1 2.3334 0.12911 depth0.1917 1 4.1145 0.04461 * river:length 0.0049 1 0.1049 0.74652 river:depth 0.3008 1 6.4567 0.01226 * Residuals5.9168 127 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > drop1(unweighted.model) Single term deletions Model: trSex ~ (river + length + depth)^2 - length:depth Df Sum of Sq RSS AIC 5.92 -401.97 river:length 1 0.0048895.92 -403.86 river:depth 1 0.306.22 -397.37 ### Both drop1 & Anova suggest dropping river:length ### Compare with the following: > weighted.model <- lm(trSex ~ (river+length +depth)^2-length:depth, dno2, weights=males+females) > Anova(weighted.model) Anova Table (Type II tests) Response: trSex Sum Sq Df F value Pr(>F) river 1.471 1 3.3775 0.06843 . length1.002 1 2.2999 0.13187 depth 2.974 1 6.8295 0.01005 * river:length 0.075 1 0.1733 0.67790 river:depth 2.020 1 4.6397 0.03313 * Residuals55.303 127 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > drop1(weighted.model) Single term deletions Model: trSex ~ (river + length + depth)^2 - length:depth Df Sum of Sq RSS AIC 55.30 -104.71 river:length 1-49.335.97 -402.79 river:depth 1-49.046.26 -396.39 __ 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] iPlots and mouse clicks
Hi, I'm investigating the iPlots package as a teaching aid, and I'd like to be able to detect mouse clicks (with x-y position) in the plotting window. I can't see how to do this - the only event loop that is illustrated in the examples uses iset.sel.changed(). Is there anyone out there who has made use of the iPlots event loop functionality and might be able to give me some tips? Thanks, and apologies if this is not the correct list, Yan Wong Leeds __ 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] Possible bug in lmer nested analysis with factors
On 18 Sep 2005, at 16:27, Douglas Bates wrote: > I have a couple of other comments. You can write the nested grouping > factors as Sundar did without having to explicitly evaluate the > interaction term and drop unused levels. The lmer function, like most > modeling functions in R, uses the optional argument drop.unused.levels > = TRUE when creating the model frame. In other words, the use of "b:c" in a model formula, where both b and c are factors, results in an internal call to evalq(b:c)[,drop=T] (or equivalent), which is treated as a factor in a temporary model data frame? I know little of the internals to R - that is new to me, but does make sense for factors. Thus I could use |a:b and |a:b:c as random terms in lme or lmer, even though 'a' is a fixed, unnested factor. Notation like this in the model formula does indeed aid clarity. By the way, I noticed that in your mlmRev vignette you recommended this as good practice for lea:school (page 3), but then omitted to do it on page 4. > John Maindonald has already suggested the use of > > (1|b/c) => (1|b:c) + (1|b) > > as "syntactic sugar" for the lmer formula and it is a reasonable > request. This is, indeed, the behaviour I was expecting. > Unfortunately, implementing this is not high on my priority > list at present. (We just made a massive change in the sparse matrix > implementation in the Matrix package and shaking the bugs out of that > will take a while.) All your efforts in these areas are, I'm sure, much appreciated. I'm certainly very interested in learning to use lmer, and welcome all the improvements that are being made. > In any case the general recommendation for nested grouping factors is > first to ensure that they are stored as factors and then to create the > sequence of interaction terms. As a brief aside, I know that some people assume that lme treats random effects as factors even if they are of a numeric type. It might be worth doing a check in lmer (and even lme) that random effects are factors, producing a warning if not. Again, this is a non-vital suggestion, and I don't wish to take up any more of your time! Thanks Yan __ 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] Possible bug in lmer nested analysis with factors
On 18 Sep 2005, at 16:04, Douglas Bates wrote: > You are correct that good documentation of the capabilities of lmer > does not currently exist. lmer is still under active development and > documentation is spread in several places. The vignette in the mlmRev > package explores some of the capabilities of lmer. Also see the > examples in that package. Yes. Thanks for this, and indeed for the development of the package. I'm currently trying to do GLMMs (binary response), so I thought that I should learn mixed modelling using a library with these capabilities. > You are correct that the denominator degrees of freedom associated > with terms in the fixed effects is different between lme and lmer. > ... > Some arguments on > degrees of freedom can be made for nested grouping factors but the > question of testing fixed effects terms for models with partially > crossed grouping factors is difficult. Would it not be possible to recognise when the model is fully nested, and make this a special case? I was imagining using lmer as a replacement for lme, so finding that they differ in this way came as some surprise. When learning to use a new, relatively undescribed routine, I usually try to see if I can reproduce known results. This is where I was coming unstuck when trying to reproduce lme results using lmer. I suspect that many people (I know of one other in my group) will use lmer as a drop-in replacement for lme specifically for its GLMM capabilities rather than for its partial nesting. I realise, however, that this might not be your priority. > This area could be a very fruitful research area for people with > strong mathematical and implementation skills. That's not me, I'm afraid. I am only just working through Chapter 1 of your (excellent) "mixed effects models in S" book. > There are already some facilities for lmer models such as mcmcsamp and > simulate which can be used for evaluating the posterior distribution > of a single coefficient or for a parametric bootstrap of the reference > distribution of a quantity like the likelihood ratio statistic for a > hypothesis test. This, again, is beyond me at the moment. But I do hope that someone else can respond to the call, especially for "textbook" as well as more complex examples of lmer usage. Best wishes Yan Wong __ 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] Possible bug in lmer nested analysis with factors
On 16 Sep 2005, at 17:21, Sundar Dorai-Raj wrote: > My guess is he wants this: > > c1 <- factor(c) > d1 <- factor(d) > m <- lmer(a ~ b + (1|c1:d1)+(1|c1)) > > which assumes d1 is nested within c1. > > Take a look at Section 3 in the "MlmSoftRev" vignette: > > library(mlmRev) > vignette("MlmSoftRev") Ah, that vignette is extremely useful: it deserves to be more widely known. I mainly intended this reply to be a thank you to yourself and Harold. Unfortunately, there is (as always), one last thing that is still puzzling me: the df for fixed factors seems to vary between what I currently understand to be equivalent calls to lme and lmer, e.g: --- a<-rnorm(36); b<-factor(rep(1:3,each=12)) c<-factor(rep(1:2,each=6,3)) d<-factor(rep(1:3,each=2,6)) c <- evalq(b:c)[,drop=T] #make unique factor levels d <- evalq(c:d)[,drop=T] #make unique factor levels summary(lme(a ~ b, random=~1|c/d)) # output includes estimates for fixed effects such as #Value Std.Error DFt-value p-value # (Intercept) 0.06908901 0.3318330 18 0.2082041 0.8374 # b2 0.13279084 0.4692828 3 0.2829655 0.7956 # b3 -0.26146698 0.4692828 3 -0.5571630 0.6163 # I understand the above lme model to be equivalent to summary(lmer(a ~ b + (1|c) +(1|c:d)) #but this gives fixed effects estimates with differing DF: # Estimate Std. Error DF t value Pr(>|t|) # (Intercept) 0.069089 0.331724 33 0.2083 0.8363 # b2 0.132791 0.469128 33 0.2831 0.7789 # b3 -0.261467 0.469128 33 -0.5573 0.5811 Again, many thanks for your help: even more so if you or anyone else can answer this last little niggle of mine. Yan __ 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] Possible bug in lmer nested analysis with factors
On 16 Sep 2005, at 17:12, Doran, Harold wrote: > I think you might have confused lme code with lmer code. Why do you > have > c/d in the random portion? Apologies. I obviously have done something of the sort. I assumed that the 'random' assignment in lme could just be incorporated into an lmer call by placing it in brackets and removing the ~, so that lme(a ~ b, random= ~ 1|c/d) would be equivalent to lmer(a ~ b + (1|c/d)) Is there a good guide somewhere to lmer calling conventions? I obviously don't understand them. As you can see, I would like to nest d within c, (and actually, c is nested in b too). Perhaps it would be better with some explanation of the Crawley data. There are 3 fixed drug treatments ('b') given to 2 rats (6 rats in all: 'c'), and 3 samples ('d') are taken from each of the rat's livers, with some response variable recorded for each sample ('a': here just allocated a Normal distribution for testing purposes). I.e. c and d are random effects, with d %in% c and c %in% b. He analyses it via aov(a ~ b+c+d+Error(a/b/c)) I'm wondering what the lme and lmer equivalents are. I've been told that a reasonable form of analysis using lme is a<-rnorm(36);b<-rep(1:3,each=12);d<-rep(1:3,each=2,6) c <- rep(1:6, each=6) #use unique labels for each rat ## I got this wrong in my previous example model1 <- lme(a ~ b, random= ~ 1|c/d) Which gives what looks to be a reasonable output (but I'm new to all this mixed modelling stuff). How would I code the same thing using lmer? > I think what you want is > >> lmer(a ~ b + (1 | c)+(1|d)) >> > > Which gives the following using your data I'm not sure this is what I wanted to do. But thanks for the all the help. Yan __ 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] Possible bug in lmer nested analysis with factors
On 16 Sep 2005, at 16:57, Yan Wong wrote: > Hello, > > Is this a bug in the lmer routine? > > ... > Error in lmer(a ~ b + (1 | c/d)) : entry 0 in matrix[0,0] has row > 2147483647 and column 2147483647 > In addition: Warning message: > / not meaningful for factors in: Ops.factor(c, d) Sorry, I forgot to specify: R 2.1.1, Mac OS X 10.4.2, lme4 version 0.98-1, Matrix version 0.98-7. Yan __ 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] Possible bug in lmer nested analysis with factors
Hello, Is this a bug in the lmer routine? > library(lme4) > ### test case based on rats data from Crawley > a<-rnorm(36);b<-rep(1:3,each=12);c<-rep(1:2,each=6,3);d<-rep (1:3,each=2,6) > > ### mixed model works when c & d are numeric, lmer assumes they are factors > m <- lmer(a ~ b + (1|c/d)) > > ### but bails out when they are actually specified as factors > c<-factor(c); d<-factor(d) > m <- lmer(a ~ b + (1| c / d)) Error in lmer(a ~ b + (1 | c/d)) : entry 0 in matrix[0,0] has row 2147483647 and column 2147483647 In addition: Warning message: / not meaningful for factors in: Ops.factor(c, d) Cheers Yan __ 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] Simple suggestion for improvement
On 3 Mar 2005, at 17:17, Adaikalavan Ramasamy wrote: How will you deal with multiple word searches such as help.search("eps dev") One way to implement would be ??"eps dev" but this looks awkward to me. That's what you have to do with the normal help function sometimes anyway, e.g. ?"+" ?"base-defunct" etc. __ 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] Simple suggestion for improvement
On 3 Mar 2005, at 10:08, Duncan Murdoch wrote: That's not a bad suggestion, but it might not be trivial to implement. Right now the "?" is an operator that is parsed like other operators such as "+": it becomes a function call . To have "??" mean something special would mean changes to the parser, or a special case to the .helpForCall function that the "?" function calls. OK, I can see that. Adding another operator might be seen as too great a change to consider "trivial". But the second way (changing the function) seems a little "hacky" to me. Anyway, it is just a suggestion that would save me (and others) some typing time. Thanks for replying to my original post, Yan __ 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] Simple suggestion for improvement
Hello, Being relatively new to R, I often find myself searching for functions using help.search("term"). Why not have the command ??term invoke it in the same way as ?topic invokes index.search("topic")? Using a double question mark to invoke a wider search for a term seems relatively intuitive to me, and presumably would be trivial to implement. Cheers Yan Wong Leeds University __ 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