Re: [R] Graphical presentation of logistic regression

[Anon.]

Beale, Colin wrote:
> Hi,
> 
> I wonder if anyone has written any code to implement the suggestions of
> Smart et al (2004) in the Bulletin of the Ecological Society of America
> for a new way of graphically presenting the results of logistic
> regression (see
> www.esapubs.org/bulletin/backissues/085-3/bulletinjuly2004_2column.htm#t
> ools1 for the full text)? I couldn't find anything relating to this sort
> of graphical representation of logistic models in the archives, but
> maybe someone has solved it already? In short, Smart et al suggest that
> a logistic regression be presented as a combination of the two
> histograms for successes and failures (with one presented upside down at
> the top of the figure, the other the right way up at the bottom)
> overlaid by the probability function (ie logistic curve). It's somewhat
> hard to describe, but is nicely illustrated in the full text version
> above. I think it is a sensible way of presenting these results and am
> keen to do so - at the moment I can only do this by generating the two
> histograms and the logistic curve separately (using hist() and lines()),
> then copying and pasting the graphs out of R and inverting one in a
> graphics package, before overlying the others. I'm sure this could be
> done within R and would be a handy plotting function to develop. Has
> anyone done so, or can anyone give me any pointers to doing this? I
> really nead to know how to invert a histogram and how to overlay this
> with another histogram "the right way up".
> 
> Any thoughts would be welcome.
> 
My reaction was that I had seen some R code in a Bulletin of the ESA 
that someone sent me.  A quick search revealed this:

which has the code.

Bob

-- 
Bob O'Hara

Dept. of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland

Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org

__
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] The Perils of PowerPoint

[Anon.]

Mike Waters wrote:

> And thus to that 'New Age' Management Role, that of the Professional
>PowePoint Ranger. He (invariably he) who culls the fruits of the labours of
>others to present in ever more slick PowerPoint compendia, whilst never
>sullying their hands with 'real' work.
>
>  
>
In academia they're known as "professors".

Bob

-- 
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland

Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org

__
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] Converting characters to numbers in data frames

[Anon.]

I'm sure I'm missing something obvious here, but I can't find the 
solution (including in the FAQ etc.).

I have a vector of names of variables like this: NRes.x.y. where x and y 
are numbers.  I want to extract these numbers as numbers to use 
elsewhere.  I can extract the numbers as a list of characters using 
strsplit(), and convert that to a data frame, e.g.:

NAMES=c("NRes.1.2.", "NRes.1.3.", "NRes.1.4.", "NRes.1.5.", "NRes.1.6.")
NUMBERS=strsplit(gsub("NRes.","", NAMES, perl =T), '.', fixed = TRUE)
NUMBERS.df=t(data.frame(NUMBERS))

But I now want to convert the characters to be numeric.  Using 
as.numeric(NUMBERS.df) converts them, but to a vector.  How can I 
convert and keep as a data frame?  I could use this:

matrix(as.numeric(NUMBERS.df), ncol=dim(NUMBERS.df)[2])

but I seem to be jumping through far too many hoops: there must be an 
easier way.  An suggestions?

Bob

-- 
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland

Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org

__
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] help: how to change the size of a window after it has been created

[Anon.]

Marc Schwartz wrote:
> On Thu, 2005-07-14 at 14:38 +0200, wu sz wrote:
> 
>>Hello,
>>
>>I wish to plot some figures in a window in turn, but the size of these
>>figures is different, so how can I change the size of the window by
>>resetting the parameters before each plotting?
>>
>>Thank you,
>>Shengzhe
> 
> 
> Other than dragging a plot window with a mouse, I do not think that
> there is a way to change the size of an open display device via code
> (though somebody will no doubt correct me if I am wrong).
> 

Perhaps some lateral thinking is called for.  Rather than change the 
size of the window, how about changing the size of the plotting region 
within the window.  For example:

split.screen(matrix(c(0,0.5, 0,1,  0.5,1,0,1), ncol=4, byrow=T))
plot(1:5)
close.screen(all=T)

split.screen(matrix(c(0,0.5, 0,0.5,  0.5,1,0.5,1), ncol=4, byrow=T))
plot(1:5)
close.screen(all=T)

Bob

-- 
Bob O'Hara

Dept. of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland

Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org

__
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] Very Slow Gower Similarity Function

[Anon.]

Jari Oksanen wrote:
On 18 Apr 2005, at 19:10, Tyler Smith wrote:
Hello,
I am a relatively new user of R. I have written a basic function to 
calculate
the Gower similarity function. I was motivated to do so partly as an 
excercise
in learning R, and partly because the existing option (vegdist in the 
vegan
package) does not accept missing values.

Speed is the reason to use C instead of R. It should be easy, almost 
trivial, to modify the vegdist.c  so that it handles missing values. I 
guess this handling means ignoring the value pair if one of the values 
is missing -- which is not so gentle to the metric properties so dear 
to Gower. Package vegan is designed for ecological community data 
which generally do not have missing values (except in environmental 
data), but contributions are welcome.

The only reason you never see ecological community data with missing 
values is because the ecologists remove those species/sites from their 
Excel sheets before they give it to you to sort out their mess.  This is 
actually one of the few things they know how to do in Excel - I'm 
dreading the day when a paper appears in JAE saying that you can use 
Excel to produce P-values.

To be slightly more serious, as an exercise the OP could consider 
writing a wrapper function in R that removes the missing data and then 
calls vegdist to calculate his Gower similarity index.

Bob
--
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org
__
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] Different graph type can coexisti??

[Anon.]

Landini Massimiliano wrote:
Please consider a data frame
A   B   C   D
1   4   5   0
2   3   2   75
3   4   1   84
4   5   1   90
5   3   0   100
Is there a way to plot column B and C as barplot *and* D as line on the same
graph??
 

It's not totally clear to me how you want to plot B and C, but the 
following code, and ?barplot should get you what you want:

Bars=matrix(cbind(B,C), nrow=2, byrow=T)
PLOT=barplot(Bars, ylim=c(0,100))
lines(PLOT,D, col=2)
axis(1)
The only tricky thing is dealing with the way R scales the x-axis, but 
fortunately the function barplot() returns the mid-points: hence this 
construction.

Bob
--
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org
__
[EMAIL PROTECTED] 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] How about a mascot for R?

[Anon.]

David Scott wrote:
On Fri, 3 Dec 2004, Patrick Connolly wrote:
On Thu, 02-Dec-2004 at 02:08PM -0600, Erin Hodgess wrote:
|> How about an R-madillo?
I'm for R-gnu which already has a song (as pointed out by Murray
Jorgensen some time back).  I'm sure an arrangement could be worked
out with the Emacs people.
The Flanders and Swann song, the exact title of which I don't recall,
was a response to someone voicing the opinion that the said animal was
an Helk.  Rather witty in my view and has a spirit in keeping with
this list.
Knowledge of and a liking for Flanders and Swan does rather date one, 
but I must confess my partiality for their humour. For the words of 
the Gnu Song:

http://www.members.optusnet.com.au/pennywyatt/Interests/FlandersSwann/DropOfaHat/At%20the%20Drop%20of%20a%20Hat04.html 


Indeed, they also had a song about an R-madillo:

And no, this doesn't date you - they had split up before I was even 
born.  I think it just shows that you have the good taste to be partial 
to humour of timeless quality.

Whilst I'm wasting bandwidth, I found the original platypus and kangaroo 
(Kangaroo?  Another F & S composition!) reference:

 (see the bottom)
The paper should be available on JSTOR, for those of you with access to it.
Bob
--
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org
__
[EMAIL PROTECTED] 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] Running sum

[Anon.]

Prof Brian Ripley wrote:
?cumsum
On Fri, 19 Nov 2004, Sean Davis wrote:
I have vector X of length N that I want to have a running sum for 
(called Y). I just need max(Y).  I do this with a "for" loop like so:

   Y <- vector(length=N)
   Y[1] <- X[1]
   for (i in 2:N) {
 Y[i] <- Y[i-1]+X[i]
   }
   return(max(Y))
Is there a faster way to do this?

My apologies, I am not being entirely serious, but...
Please do read the posting guide! 
(particularly the third point in the "Responding to other posts" section).

Bob
--
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org
__
[EMAIL PROTECTED] 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] (no subject)

[Anon.]

fang lai wrote:
Dear all,
I have several questions regarding fisher.test() in
R, and I'd highly appreciate any help with it.
I have a group of observations, each having people's
income, and an indicator of whether selected in or out
a program. I want to test the difference between
income of people who are in and out.
Because the distribution is far from normal, I decide
to use the fisher's exact test, using either mean or
rank as statistics.
Question 0 is: Can I do this test using fisher.test()
in R?
 

The answer is "no".  I think you have mis-understood the purpose of 
Fisher's exact test: read (and understand!) the description in the help 
page.

If your data were normall distributed, then you could use a t-test 
(t.test()).  As you are not happy with the normality assumption, you 
could try an equivalent non-paramteric test, such as the Wilcoxon test 
(wilcox.test()), also known as a Mann-Whitney test.  I would recommend 
that you make sure you understand how the test works first.

Bob
--
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org
__
[EMAIL PROTECTED] 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] Odd behaviour with scale()

[Anon.]

Moi!
A student here has been getting a bit irritated with some side effects 
of scale() (OS is Windows XP, the behaviour occurs in R 2.0.0, but not 
1.7.1).  The problem is that she scales a variable in a data frame, then 
does a regression, and tries to get some predictions for some new data.  
However, at this point she gets an error (see the example below).  This 
seems to be because the scaled variable in the new data frame does not 
have the center and scale attributes, but the one in the old data frame 
does.

The work-around is to put the scaled variable intro a new data frame, 
which again won't have the attributes.  But it seems odd to me that 
whether a scale()'d variable has attributes depends on where it's 
placed.  I presume that this is because I'm not understanding something 
about the way R is working, rather than it being a bug.  Would anyone 
care to enlighten me?

> Data1=data.frame(xx=1:10, yy=2.1:12)
> Data1$xx=scale(Data1$xx)
>
> reg1=lm(yy~xx, data=Data1)
> New=data.frame(xx=2:4)
> b=predict(reg1, New, se.fit=T)
Error: variable 'xx' was fitted with nmatrix.1 but numeric was supplied
>
> New=data.frame(xx=scale(2:4, center=5.5, scale=3.02))
> b=predict(reg1, New, se.fit=T)
Error: variable 'xx' was fitted with nmatrix.1 but numeric was supplied
>
Bob
--
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org
__
[EMAIL PROTECTED] 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] [Fwd: OpenBUGS]

[Anon.]

This is slightly off-topic, but there was a discussion about MCMC a 
couple of weeks ago.  OpenBUGS can operate from R, at least in Windows 
(there are some problems in Linux at the moment).

The interface with R is one are that there are plans to work on: it's a 
bit basic at the moment.

Bob
--
Bob O'Hara
Dept. of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org
 Original Message 
Subject: OpenBUGS
Date: Wed, 29 Sep 2004 13:53:30 +0100
From: Andrew Thomas <[EMAIL PROTECTED]>
Reply-To: Andrew Thomas <[EMAIL PROTECTED]>
To: [EMAIL PROTECTED]

We would like to announce that OpenBUGS an open source version of BUGS
is now available and has its own web page at
http://mathstat.helsinki.fi/openbugs/
Eventually the OpenBUGS code base will replace v1.4 of WinBUGS
Enjoy
Andrew Thomas
__
[EMAIL PROTECTED] 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] BUGS and OS X

[Anon.]

Jari Oksanen wrote:
On Wed, 2004-09-15 at 21:29, Tamas K Papp wrote:
 

On Wed, Sep 15, 2004 at 02:21:18PM -0400, Liaw, Andy wrote:
   

That's more of a question for the BUGS developers.  BUGS is not open
 

source,
 

so whatever binary is provided, that's all you can use.  If I'm not
mistaken, WinBUGS is the only version under development.
 

I found something called JAGS, and I am still exploring it.  It
appears to be an open-source BUGS replacement, thought with
limitations.
   

MacOS X is a kind of unix (where the emphasis is on the "kind of"), so
you can get and compile any source code developed for unix -- with some
luck. One alternative is Bassist available at
http://www.cs.helsinki.fi/research/fdk/bassist/. I just tried and found
out that you can compile and install it in MacOS X in the usual way
(./configure && make && sudo make install). That's all I can say about
it. It may not be easiest to use. The current version seems to be a bit
oldish and not quite complete, but somebody claimed that they may start
developing Bassist again. Actually, Bob O'Hara (who usually calls
himself "Anon." in this list) should know more, and hopefully this
message will prompt him to tell us, too.
 

I don't think Andrew has even tried installing WinBUGS it on a Mac 
(first problem: find a Mac).

I haven't heard about any plans to develop BASSIST further, but as it's 
open source, someone else could take it up.  The BASSIST language is 
very BUGS-like, so shouldn't be a big problem to learn.

To answer Ramon about BUGS and R, be patient: I think an announcement 
will come soon, and I'd rather not pre-empt that.

We're drifting away from R, so can I suggest either ending this, 
sticking to R, or shifting it to the WinBUGS list?  That's probably a 
better forum for getting feedback on Bayesian stuff anyway.

Oh, and Jari, are you likely to submit your R: opas ekologielle to CRAN 
at any time?  Hint, hint.

Bob
--
Bob O'Hara
Department of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: www.jnr-eeb.org
__
[EMAIL PROTECTED] 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] getting started on Bayesian analysis

[Anon.]

Jari Oksanen wrote:
On Wed, 2004-09-15 at 03:27, HALL, MARK E wrote:
 I've found 
 
Bayesian Methods: A Social and Behavioral Sciences Approach
by Jeff Gill 

useful as an introduction.  The examples are written in R and S with generalized scripts for doing 
a variety of problems.  (Though I never got change-point analysis to successfully in R.)

Change point analysis? I haven't seen the book, but I read lecture
handouts of one Bayesian course over here in Finland (Antti Penttinen,
Jyväskylä), and translated his example to R during one (rare) warm
summer day in a garden. So do you mean this (binary case):



This was off-topic. So something about business: isn't the (Win)BUGS
author working with a R port?
Yes, he is, and he's got it working in Windows.  But if anyone wants to 
discuss how R does memory management with Andrew, he'll be all ears.

In the mean time, there is the R2WinBUGS package on CRAN for you to enjoy.
Bob
--
Bob O'Hara
Dept. of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org
__
[EMAIL PROTECTED] 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] Population simulation.

[Anon.]

[EMAIL PROTECTED] wrote:
>
> Hello,
>
> can anyone tell me if R has any special function for simulating the 
structure
> of human populations? Something like the genetic algorithm?
> I need to simulate a sample of a population with a specific 
structure. Is
> there something on R that can help me?
>
As others have pointed out, you need to give us some more details. 
There is an artile about simulation of ecological dynamics in last 
December's R News, which may be of some help:


If you're simulating something like a coalescence process, then it might 
be easier to use some specialist software.  In that case, googling with 
a few relevant terms might be the way to go.

Bob
--
Bob O'Hara
Dept. of Mathematics and Statistics
P.O. Box 68 (Gustaf Hällströmin katu 2b)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 51479
Mobile: +358 50 599 0540
Fax:  +358-9-191 51400
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Odd behaviour of step (and stepAIC)?

[Anon.]

I can only assume I'm betraying my ignorance here, but this is not what 
I would expect.

I'm getting the following from a stepwise selection (with both step and 
stepAIC):

> step(lm(sqrt(Grids)~ SE + Edge + NH), scope=~ (Edge + SE + NH)^2)
Start:  AIC= 593.56
 sqrt(Grids) ~ SE + Edge + NH
  Df Sum of SqRSSAIC
 2147.0  593.6
+ Edge:NH  1   3.0 2143.9  595.1
+ SE:NH4  23.2 2123.8  598.4
- NH   1  75.8 .8  601.6
- Edge 1 448.7 2595.7  646.4
- SE   41033.7 3180.6  699.1
My problem is that the SE:Edge term is not added.  Now, I know that I've 
specified the terms in a different order in the start model and the 
scope, but this shouldn't matter, should it?  Aren't these model 
specifications commutative?

Am I missing something important, or is this just an obscure feature?

Bob

--
Bob O'Hara
Dept. of Mathematics and Statistics
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Probit predictions outside (0,1) interval

[Anon.]

Arnab mukherji wrote:
Hi!

I was trying to implement a probit model on a dichotomous outcome variable and found that the predictions were outside the (0,1) interval that one should get. I later tried it with some simulated data with a similar result. 

Here is a toy program I wrote and I cant figure why I should be getting such odd predictions.

Let me be the first to write "read the help file".

There are several scales th at you can predict on in a GLM.  The 
relevant part of the help file is this:

 type: the type of prediction required.  The default is on the scale
  of the linear predictors; the alternative `"response"' is on
  the scale of the response variable.  Thus for a default
  binomial model the default predictions are of log-odds
  (probabilities on logit scale) and `type = "response"' gives
  the predicted probabilities.
Your predictions are on the probit scale.

Bob

--
Bob O'Hara
Dept. of Mathematics and Statistics
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Problem with Integrate

[Anon.]

Peter Dalgaard wrote:
"Anon." <[EMAIL PROTECTED]> writes:


I assume that this is because for much of the range, the integral is
basically zero.
The help page for integrate() says
When integrating over infinite intervals do so explicitly,
rather
than just using a large number as the endpoint.  This increases
the chance of a correct answer - any function whose integral over
an infinite interval is finite must be near zero for most of that
interval.
That is, if you want an integral from 0 to Inf, do that.
Sorry, I forgot to put that in my message.  It gives the same error as
a large value.  I assume it's all a result of the NaN's being returned.


It does seem to help a bit if you modify the integrand to


PLN1
function(lam, Count, mu, sigma2) {
   t1 <- dpois(Count, exp(lam), log=F)
   t2 <- dnorm(lam, mu, sqrt(sigma2))
   ifelse(t1==0|t2==0,0,t1*t2)
}



Mange tak!

And thank you to everyone else for their help too.  The generic solution 
will help me with the next stage of my work.  I can see that it might 
break down, but only under rather bizarre circumstances.

Bob

--
Bob O'Hara
Dept. of Mathematics and Statistics
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Problem with Integrate

[Anon.]

Thomas Lumley wrote:
On Tue, 2 Mar 2004, Anon. wrote:


The background: I'm trying to fit a Poisson-lognormal distrbutuion to
some data.  This is a way of modelling species abundances:
N ~ Pois(lam)
log(lam) ~ N(mu, sigma2)
The number of individuals are Poisson distributed with an abundance
drawn from a log-normal distrbution.
To fit this to data, I need to integrate out lam.  In principle, I can
do it this way:
PLN1 <- function(lam, Count, mu, sigma2) {
  dpois(Count, exp(lam), log=F)*dnorm(LL, mu, sqrt(sigma2))
}
and integrate between -Inf and Inf.  For example, with mu=2, and
sigma2=2.8 (which are roughly right for the data), and Count=73, I get this:
> integrate(PLN1, -10, 10, Count=73, mu=2, sigma2=2.8)
0.001289726 with absolute error < 2.5e-11
> integrate(PLN1, -20, 20, Count=73, mu=2, sigma2=2.8)
0.001289726 with absolute error < 2.5e-11
> integrate(PLN1, -100, 100, Count=73, mu=2, sigma2=2.8)
2.724483e-10 with absolute error < 5.3e-10
> integrate(PLN1, -500, 500, Count=73, mu=2, sigma2=2.8)
1.831093e-73 with absolute error < 3.6e-73
> integrate(PLN1, -1000, 1000, Count=73, mu=2, sigma2=2.8)
Error in integrate(PLN1, -1000, 1000, Count = 73, mu = 2, sigma2 = 2.8):
non-finite function value
In addition: Warning message:
NaNs produced in: dpois(x, lambda, log)
So, the integral gets smaller, and then gives an error.


I assume that this is because for much of the range, the integral is
basically zero.


The help page for integrate() says

 When integrating over infinite intervals do so explicitly, rather
 than just using a large number as the endpoint.  This increases
 the chance of a correct answer - any function whose integral over
 an infinite interval is finite must be near zero for most of that
 interval.
That is, if you want an integral from 0 to Inf, do that.

Sorry, I forgot to put that in my message.  It gives the same error as a 
large value.  I assume it's all a result of the NaN's being returned.

Bob

--
Bob O'Hara
Dept. of Mathematics and Statistics
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Problem with Integrate

[Anon.]

The background: I'm trying to fit a Poisson-lognormal distrbutuion to 
some data.  This is a way of modelling species abundances:
N ~ Pois(lam)
log(lam) ~ N(mu, sigma2)
The number of individuals are Poisson distributed with an abundance 
drawn from a log-normal distrbution.

To fit this to data, I need to integrate out lam.  In principle, I can 
do it this way:

PLN1 <- function(lam, Count, mu, sigma2) {
  dpois(Count, exp(lam), log=F)*dnorm(LL, mu, sqrt(sigma2))
}
and integrate between -Inf and Inf.  For example, with mu=2, and 
sigma2=2.8 (which are roughly right for the data), and Count=73, I get this:

> integrate(PLN1, -10, 10, Count=73, mu=2, sigma2=2.8)
0.001289726 with absolute error < 2.5e-11
> integrate(PLN1, -20, 20, Count=73, mu=2, sigma2=2.8)
0.001289726 with absolute error < 2.5e-11
> integrate(PLN1, -100, 100, Count=73, mu=2, sigma2=2.8)
2.724483e-10 with absolute error < 5.3e-10
> integrate(PLN1, -500, 500, Count=73, mu=2, sigma2=2.8)
1.831093e-73 with absolute error < 3.6e-73
> integrate(PLN1, -1000, 1000, Count=73, mu=2, sigma2=2.8)
Error in integrate(PLN1, -1000, 1000, Count = 73, mu = 2, sigma2 = 2.8):
non-finite function value
In addition: Warning message:
NaNs produced in: dpois(x, lambda, log)
So, the integral gets smaller, and then gives an error.

I then tried entering the formula directly:
PLN2 <- function(LL, Count, mu, sigma2) {
exp(-LL-(log(LL)-mu)^2/(2*sigma2))*LL^(Count-1)/(gamma(Count+1)*sqrt(2*pi*sigma2))
}
> integrate(PLN2, 0, 100, Count=73, mu=2, sigma2=2.8)
0.001287821 with absolute error < 2.6e-10
> integrate(PLN2, 0, 1000, Count=73, mu=2, sigma2=2.8)
0.001289726 with absolute error < 2.9e-08
> integrate(PLN2, 0, 1, Count=73, mu=2, sigma2=2.8)
0.001289726 with absolute error < 9.7e-06
> integrate(PLN2, 0, 19100, Count=73, mu=2, sigma2=2.8)
1.160307e-08 with absolute error < 2.3e-08
> integrate(PLN2, 0, 19200, Count=73, mu=2, sigma2=2.8)
Error in integrate(PLN2, 0, 19200, Count = 73, mu = 2, sigma2 = 2.8) :
non-finite function value
And the same thing happens.

I assume that this is because for much of the range, the integral is 
basically zero.

Can anyone suggest a fix?  Preferably one that will work with Count=320 
and Count=0 (both of which I have in the data).

Bob

--
Bob O'Hara
Dept. of Mathematics and Statistics
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Using split.screen

[Anon.]

I want to draw a figure with several panels of unequal size, so i 
thought I would try using screen().  However, I can't figure out how to 
define the sizes as a matrix.  I've tried this:

split.screen(matrix(c(0,0.5,0,0.5,  0.5,1,0.5,1), byrow=F, ncol=4))

and a couple of variants on it, but get the same error:

Error in par(.split.screens[[cur.screen]]) :
invalid value specified for graphics parameter "fig".
The help usefully says that they are defined in NDC units, but I don't 
know what an NDC unit is, and there isn't any example.  The code in 
kjetil brinchmann halvorsen's message on R-help on  Mar 31 2003 (do a 
search for "NDC units"!) didn't work either, it gives the same message:

> split.screen( matrix( c(0, 0.3, 0.5, 1, 0.3, 0.7, 0.5, 1,
+ + 0.7, 1, 0.5, 1, 0, 0.5, 0, 0.5,
+ + 0.5, 1, 0, 0.5), 5, 4, byrow=TRUE))
Error in par(.split.screens[[cur.screen]]) :
invalid value specified for graphics parameter "fig".
I get the same in R-1.8.1 on Windows, and R-1.5.1 on Linux.
As Kjetil pointed out then, "NDC" is not explained in the help pages, 
and I don't have my copy of MASS with me.

Bob

--
Bob O'Hara
Rolf Nevanlinna Institute
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Offsets in glmmPQL?

[Anon.]

I've got a colleague who's using a GLMM to analyse her data, and I've 
told her that she needs to include an offset.  However, glmmPQL doesn't 
seem to allow one to be included.  Is there anyway of doing this?

Bob

--
Bob O'Hara
Rolf Nevanlinna Institute
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help

[R] Maximisation of likelihood of a discrete parameter

[Anon.]

Moi!

I have a problem where I want to find the ML estimate of a discrete 
parameter.  I just want a function like optim that finds the max/min 
value for a function.  Does anyone know of such a function for R?

Thanks.

Bob

--
Bob O'Hara
Rolf Nevanlinna Institute
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help

Re: [R] Plot of Canonical Correlation Analysis

[Anon.]

[EMAIL PROTECTED] wrote:
Dear all,

I didn't find any graphical solution in the package "mva" to plot the
canonical scores from a CCA (canonical correlation analysis).
Does anybody knows how to plot or has anybody already programmed :
  - the map of the canonical scores,
  - the graph of the canonical weights,
  - the correlation circle i.e. the canonical loadings ?
Thank you for help ...
Have a look at Jari Oksanen's vegan package.  It's on CRAN, but it also 
has its own page: 

Bob

--
Bob O'Hara
Rolf Nevanlinna Institute
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help

[R] Strange axis labels?

[Anon.]

Moi!

I'm trying to add a rather long label to a y-axis, and it's so long that
it won't fit into one line.  However, things get strange when I try and
split it over 2 lines.  The problem seems to be the plus/minus symbol,
which means I have to use expression().  I'm using R1.6.1 on Windoze
2000.

As an example:

thing1 <- expression(paste("log odds of survival probability ( " %+-%
"95% C.I.)", sep=" ")
thing2 <- expression(paste("log odds of survival probability \n( " %+-%
"95% C.I.)"))

plot(c(1,2), c(1,2), ylab=thing1)
par(mar =c(2.1, 6.1, 2.1, 2.1)); plot(c(1,2), c(1,2), ylab=thing2)

The first plot has the label on one line, with only a single space
between the bracket and the +/-.  Tbis is what I want, but I wonder why
paste has a sep=, if it doesn't do anything (I assume I've missed
something here).

The second plot has the label on 2 lines (you need the par() to get it
all in!), but the open bracket, "(", lines up with the start of the
first lines, and there is then a large gap, and the string continues
with the plus/minus below and to the right of the final part of the
first line.  Not what I expected at all.  Can anyone explain what's
going on, and how to remove this large space?

Thanks!

Bob

-- 
Bob O'Hara

Rolf Nevanlinna Institute
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland

tel: +358 9 191 23743  mobile: +358 50 599 0540
fax: +358 9 191 22779email: [EMAIL PROTECTED]

It is being said of a certain poet, that though he tortures the English
language, he has still never yet succeeded in forcing it to reveal his
meaning
- Beachcomber

__
[EMAIL PROTECTED] mailing list
http://www.stat.math.ethz.ch/mailman/listinfo/r-help