Re: [R] Direction of panel plots in trellis graphics

2007-07-13 Thread Yan Wong

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

2007-07-13 Thread Yan Wong
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

2006-08-11 Thread Yan Wong

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

2006-08-11 Thread Yan Wong
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

2006-08-11 Thread Yan Wong

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

2006-08-06 Thread Yan Wong
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

2006-06-13 Thread Yan Wong
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

2006-05-24 Thread Yan Wong

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

2006-05-23 Thread Yan Wong

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

2006-05-23 Thread Yan Wong

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

2006-05-23 Thread Yan Wong
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

2006-03-08 Thread Yan Wong
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

2006-03-02 Thread Yan Wong

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

2006-03-01 Thread Yan Wong
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

2006-02-17 Thread Yan Wong
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

2005-09-20 Thread Yan Wong

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

2005-09-20 Thread Yan Wong

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

2005-09-18 Thread Yan Wong

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

2005-09-16 Thread Yan Wong

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

2005-09-16 Thread Yan Wong

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

2005-09-16 Thread Yan Wong
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

2005-03-04 Thread Yan Wong
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

2005-03-03 Thread Yan Wong
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

2005-03-03 Thread Yan Wong
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