[R] P values in non linear regression and singular gradients using nls
Carrasco-Torrecilla, Roman R
roman.carrasco06 at imperial.ac.uk
Wed Aug 6 00:44:04 CEST 2008
Dear all,
We are trying to fit a non linear model to dispersal data.
It seems that sometimes when the fit of the model of the data is not
very good we start getting singular gradient errors. However if we
modify slightly the data this won't occurr. We have also tried changing
the initial parameter values and the algorithm for fitting in nls but
didn't help.
So we ended up programming a software to do the fitting for us using a
searching algorithm in the "map" of the parameters and minimising the
sum of squares (we get the same fits as R in the models that could fit
the data well). We still don't know why R could not fit those models.
We now wonder how can we calculate the standard errors and p-values of
the non linear regression parameters. We have been trying to find
information about it and all I find has to do with linear regression.
Could you tell me please how are they calculated in nls in R?
And also what are we doing wrong with nls that was not able to fit the
model?
Many thanks in advance,
Roman Carrasco.
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of r-help-request at r-project.org
Sent: 04 August 2008 11:00
To: r-help at r-project.org
Subject: R-help Digest, Vol 66, Issue 4
Send R-help mailing list submissions to
r-help at r-project.org
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-help
or, via email, send a message with subject or body 'help' to
r-help-request at r-project.org
You can reach the person managing the list at
r-help-owner at r-project.org
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-help digest..."
Today's Topics:
1. Re: R command history -- can it be like Matlab's? (Peter Dalgaard)
2. Re: output components of GAM (Gavin Simpson)
3. Re: Gaps in time series. (Gabor Grothendieck)
4. Re: drop1() seems to give unexpected results compare to
anova() (Thomas P C Chu)
5. Re: How do I get an inset graph (i.e. graph within a graph)?
(Mark Difford)
6. Re: How do I get an inset graph (i.e. graph within a graph)?
(Gabor Grothendieck)
7. Re: losing row.names in matrix operations (Gabor Grothendieck)
8. neighbors of an element of a matrix (rostam shahname)
9. Will This Computer Run 64 bit Ubuntu/R? (Tom La Bone)
10. Installing R on Apache Web Server (Ajay ohri)
11. access object inside a function (Jonas)
12. Re: access object inside a function (Duncan Murdoch)
13. Re: running strucchange? (Burgs)
14. gmaps, projection, inset (John P. Burkett)
15. Re: access object inside a function (Patrick Burns)
16. Re: access object inside a function (Jonas)
17. (sans objet) (pir2.jv)
18. Re: drop1() seems to give unexpected results compare to
anova() (Prof Brian Ripley)
19. Re: (sans objet) (Dirk Eddelbuettel)
20. Re: Bubble plots (Franz Mueter)
21. Re: Will This Computer Run 64 bit Ubuntu/R? (Prof Brian Ripley)
22. Re: Will This Computer Run 64 bit Ubuntu/R? (Dirk Eddelbuettel)
23. Re: access object inside a function (Duncan Murdoch)
24. missing F statistic in anova.gam (Michael A. Milligan)
25. Re: Will This Computer Run 64 bit Ubuntu/R? (Prof Brian Ripley)
26. Convert date to decimal days (cls59)
27. Re: convert for loop into apply() (Charles C. Berry)
28. Re: Convert date to decimal days (Prof Brian Ripley)
29. Re: Convert date to decimal days (cls59)
30. Determining model parameters (rkevinburton at charter.net)
31. Re: Determining model parameters (Ben Bolker)
32. Changing values (Andrew Ramsey)
33. Re: Changing values (jim holtman)
34. Re: Bubble plots (Michael Bibo)
35. Re: Determining model parameters (Ben Bolker)
36. Re: Changing values (J Dougherty)
37. Sweave and ggplot2 (Sorn.Norng at dpi.vic.gov.au)
38. Decomposing tests of interaction terms in mixed-effects
models (Andrew Robinson)
39. about the 95%CI around the median... (Fernando Marmolejo Ramos)
40. Re: about the 95%CI around the median... (Simon Blomberg)
41. Re: graph (MORNEAU Fran?ois)
42. Major difference in multivariate analyses SPSS and R (Draga, R.)
43. Re: Sweave and ggplot2 (ONKELINX, Thierry)
44. Re: Decomposing tests of interaction terms in mixed-effects
models (Peter Dalgaard)
45. Re: Major difference in multivariate analyses SPSS and R
(Peter Dalgaard)
46. Re: about the 95%CI around the median... (Gavin Simpson)
47. Re: Unexpected nls behaviour: Solved (Keith Jewell)
48. Re: How to format the output file just the way I want ?
(Pierre8rou)
49. Re: source a script file straight from a subversion
repository (Martin Maechler)
----------------------------------------------------------------------
Message: 1
Date: Sun, 03 Aug 2008 12:26:48 +0200
From: Peter Dalgaard <p.dalgaard at biostat.ku.dk>
Subject: Re: [R] R command history -- can it be like Matlab's?
To: Gad Abraham <gabraham at csse.unimelb.edu.au>
Cc: r-help at r-project.org, Prof Brian Ripley <ripley at stats.ox.ac.uk>
Message-ID: <489587E8.70108 at biostat.ku.dk>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>> Not quite true that you can't type anything. What happens (for me) is
>> that you are still in reverse-i-search, so you can get this effect
>> from "^R l ^C d".
>>
>> >
>> (reverse-i-search)`l': ls()
>> (reverse-i-search)`l': ls()
>> >
>> (reverse-i-search)`l': ls()
>> (reverse-i-search)`ld': levels(ftpain3) <-
>> list(none="none",intermediate=c("mild","medium"),severe="severe")
>>
>> It snaps out of it if you press ^C twice.
>>
>
> If try typing after a ctrl-C in reverse search mode , nothing
> displays, and the console beep sounds. It doesn't matter how many
> times I press ctrl-C, the same happens.
I think this depends on what is in your history. At least typing 's'
should do something if you had an 'ls()' previously.
> The only way out of it is to press either enter (which shows the
> previously highlighted command ls and then runs it) or to use the up
> key, and then ctrl-C.
Now that is actually true for me too. I must have hit "up" after the 2nd
^C.
>
> I'm using Ubuntu Hardy.
>
And I Fedora 9, but I'm pretty sure this is generic to systems that use
libreadline.
--
O__ ---- Peter Dalgaard ?ster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
------------------------------
Message: 2
Date: Sun, 03 Aug 2008 12:45:50 +0100
From: Gavin Simpson <gavin.simpson at ucl.ac.uk>
Subject: Re: [R] output components of GAM
To: andreasgwinter at yahoo.com
Cc: r-help at r-project.org
Message-ID: <1217763950.2650.6.camel at graptoleberis.geog.ucl.ac.uk>
Content-Type: text/plain
is this gam::gam or mgcv::gam? If the latter, and if I understand what
you want, if gam.out contains your GAM model fitted using gam() from
mgcv, then:
out.term <- predict(gam.out, type = "terms")
gives you the contributions of each covariate on the response for the
fitted model. read :?predict.gam for details of what type = "terms"
actually returns.
Note that you should be wary of accessing the components of model fits
directly. Instead, get the various bits using accessor functions,
coefs() for $coefficients, fitted() for the fitted values etc. Often the
components don't hold exactly what you think but which are processed
correctly when accessed by the appropriate accessor function.
HTH
G
On Sat, 2008-08-02 at 18:47 -0700, Andreas Winter wrote:
> I would like to request help with the following:
> I am trying to use a Generalized Additive Model (gam) to examine the
> density distribution of fish as a function of latitude and longitude
> as continuous variables, and year as a categorical variable. The model
> is written as:
>
> gam.out <- gam(Density ~ s(Lat) + s(Lon) + as.factor(Year))
>
> The fitted model prediction of the link function is gam.out
> $linear.predictors. Presumably, gam.out$linear.predictors must be
> derived from some combination of the original predictor variables
> (Lat, Lon, Year), their corresponding coefficients and the intercept
> (gam.out$coefficients), and the smooth outputs gam.out$smooth and/or
> gam.out$sp.
>
> By comparison, for a glm model:
>
> glm.out <- glm(Density ~ Lat + Lon + as.factor(Year))
>
> this is simply:
>
> glm.out$linear.predictors = glm.out$coefficients(intercept) + glm.out
> $coefficients (year) + glm.out$coefficients(lat) x Lat + glm.out
> $coefficients(lon) x Lon
>
> My problem is that I cannot figure out how to get the equivalent from
> the gam model. I would like to know how to decompose gam.out
> $linear.predictors into its components so that I can evaluate the
> effects of the different predictor variables separately.
>
> I would appreciate any comments that can help me with this.
>
> Thank you,
>
> Andreas Winter
> Blacksburg, VA
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
------------------------------
Message: 3
Date: Sun, 3 Aug 2008 09:02:24 -0400
From: "Gabor Grothendieck" <ggrothendieck at gmail.com>
Subject: Re: [R] Gaps in time series.
To: rkevinburton at charter.net
Cc: R-help at r-project.org
Message-ID:
<971536df0808030602x6abdb431gaccd0ae3a3251b66 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
If you are using zoo for this you can do:
library(zoo)
set.seed(1)
x <- zoo(rnorm(7), as.numeric(1:7))
y <- zoo(rnorm(4), c(1, 2, 6, 7))
xy <- merge(x,y)
xy$y - xy$x
# or if you just want the difference at times existing in both
y - x
On Sat, Aug 2, 2008 at 5:16 PM, <rkevinburton at charter.net> wrote:
> I like the fact that in subtracting two time series objects that there
is some effort to align the series. So if I have a time series of that
begins at 1 and one that begins at 2 a subtraction operation makes sure
that the proper values are subtracted. But I am unclear as to the best
way to build a time series with "holes". say that I have data for "day"
1,2,6,7 in one time series and 1,2,3,4,5,6,7 in another. How do I
construct the time series and indicate the missing data for 3,4,5 as in
the first series? Eventually I would like to find the difference between
the two series. If the difference is series 2 minues series one, then
for the indexes of 3,4,5 in the resultant difference it would just be
the series 2 value. The resultant series would have a length of 7.
>
> Thank you.
>
> Kevin
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
------------------------------
Message: 4
Date: Sun, 03 Aug 2008 09:53:59 -0400
From: Thomas P C Chu <tchu11 at netscape.net>
Subject: Re: [R] drop1() seems to give unexpected results compare to
anova()
To: r-help at r-project.org
Message-ID: <8CAC38783A1F557-DCC-1DB6 at WEBMAIL-MA04.sysops.aol.com>
Content-Type: text/plain; charset="us-ascii"; format=flowed
Thanks to Mr Dalgaard for his advice and everyone else who has
contributed. Inclusion of an error term at the end of sim.set$y = ...
line did cure my problems with drop1() and step().
I suppose it is my own inexperience in carrying out simulations caused
such gaffe.
Thomas
________________________________________________________________________
AOL Email goes Mobile! You can now read your AOL Emails whilst on the
move. Sign up for a free AOL Email account with unlimited storage today.
------------------------------
Message: 5
Date: Sun, 3 Aug 2008 06:55:41 -0700 (PDT)
From: Mark Difford <mark_difford at yahoo.co.uk>
Subject: Re: [R] How do I get an inset graph (i.e. graph within a
graph)?
To: r-help at r-project.org
Message-ID: <18798809.post at talk.nabble.com>
Content-Type: text/plain; charset=us-ascii
Hi Arthur,
This can be done quite easily using the appropriate arguments listed
under
?par; and there are other approaches. Ready-made functions exist in
several
packages. I tend to use ?add.scatter from package ade4. It's a short
function, so it's easy to customize it, but it works well straight out
of
the box.
HTH, Mark.
Arthur Roberts wrote:
>
> Hi, all,
>
> How do I get an inset graph (i.e. graph within a graph)? Your input
> is greatly appreciated.
>
> Best wishes,
> Art
> University of Washington
> Department of Medicinal Chemistry
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
--
View this message in context:
http://www.nabble.com/How-do-I-get-an-inset-graph-%28i.e.-graph-within-a
-graph%29--tp18796563p18798809.html
Sent from the R help mailing list archive at Nabble.com.
------------------------------
Message: 6
Date: Sun, 3 Aug 2008 10:03:54 -0400
From: "Gabor Grothendieck" <ggrothendieck at gmail.com>
Subject: Re: [R] How do I get an inset graph (i.e. graph within a
graph)?
To: "Arthur Roberts" <aroberts99163 at yahoo.com>
Cc: r-help at stat.math.ethz.ch
Message-ID:
<971536df0808030703u4c74d306ub84cc6a741a670b0 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
On Sun, Aug 3, 2008 at 3:51 AM, Arthur Roberts <aroberts99163 at yahoo.com>
wrote:
> How do I get an inset graph (i.e. graph within a graph)? Your input
is
> greatly appreciated.
See ?subplot in the TeachingDemos package.
------------------------------
Message: 7
Date: Sun, 3 Aug 2008 10:07:24 -0400
From: "Gabor Grothendieck" <ggrothendieck at gmail.com>
Subject: Re: [R] losing row.names in matrix operations
To: rcoder <mpdotbook at gmail.com>
Cc: r-help at r-project.org
Message-ID:
<971536df0808030707y3b83a302l4f56bf653971441d at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Constructing a zoo object is well covered in ?zoo
and in the three vignettes:
vignette("zoo")
vignette("zoo-quickref")
vignette("zoo-faq")
On Sat, Aug 2, 2008 at 10:27 AM, rcoder <mpdotbook at gmail.com> wrote:
>
> Thanks for your reply Gabor.
>
> As I already have a dates column in my dataframe (in column
row.names), is
> it possible to preserve this whilst still making a data set suitable
for
> rollapply()?
>
> Thanks,
>
> rcoder
>
>
>
> Gabor Grothendieck wrote:
>>
>> If its regular you can convert it to ts or zoo.
>> If its irregular convert it to zoo. There is no
>> reason to expect rollapply to work with objects
>> of other classes. Read ?ts and ?zoo. In
>> ts note the start and frequency arguments.
>>
>>
>> On Sat, Aug 2, 2008 at 7:50 AM, rcoder <mpdotbook at gmail.com> wrote:
>>>
>>> Hi everyone,
>>>
>>> I have a data frame, with the following format:
>>>
>>> MatDate->
>>> row.names ID1 ID2 ID3
>>> 1 date1
>>> 2 date1
>>> 3 date3
>>> etc
>>>
>>> but I cannot perform a rollapply() statement on the matrix without
>>> converting the matrix into a time series.
>>> i.e. MatTs<-ts(MatDate)
>>
>> Use the start and frequency arguments. See ?ts
>>
>>>
>>> Only then will my rollapply statement work:
>>> MatMin<-rollapply(MatTs, 2,by=2, min, na.rm=F)
>>>
>>> If I apply the rollapply() statement to the dataframe, I get the
>>> following
>>> error: Error: could not find function "rollapply"
>>>
>>> The problem is that when I convert the data.frame matrix into a time
>>> series
>>> matrix, I lose the dates in the row.names column. I just want to
know if
>>> anyone could suggest a way to get around this problem, i.e. keep the
>>> row.names column in place, and use the rollapply() statement as
above.
>>>
>>> Thanks,
>>>
>>> rcoder
>>> --
>>> View this message in context:
>>>
http://www.nabble.com/losing-row.names-in-matrix-operations-tp18788509p1
8788509.html
>>> Sent from the R help mailing list archive at Nabble.com.
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>
> --
> View this message in context:
http://www.nabble.com/losing-row.names-in-matrix-operations-tp18788509p1
8789688.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
------------------------------
Message: 8
Date: Sun, 3 Aug 2008 12:03:58 -0300
From: "rostam shahname" <rostamepython at gmail.com>
Subject: [R] neighbors of an element of a matrix
To: r-help at r-project.org
Message-ID:
<de90f9680808030803j61884d1br3f98972ce76fd31d at mail.gmail.com>
Content-Type: text/plain
Hi R users,
I wonder if there is any function which would render the value of the
neighbors of the given element [i,j] of a matrix.
Thanks,
Rostam
[[alternative HTML version deleted]]
------------------------------
Message: 9
Date: Sun, 3 Aug 2008 08:04:04 -0700 (PDT)
From: Tom La Bone <booboo at gforcecable.com>
Subject: [R] Will This Computer Run 64 bit Ubuntu/R?
To: r-help at r-project.org
Message-ID: <18799376.post at talk.nabble.com>
Content-Type: text/plain; charset=us-ascii
After doing some reading about 64-bit systems and software I am still
somewhat uncertain about some things. I have a Dell Dimension XPS 400
with a
dual core Intel Pentium D 940 (3.2 GHz) and 4 Gb of memory. I currently
dual
boot the system with XP Professional and Ubuntu 8.04 (32 bit). If I
simply
install Ubuntu 8.04 (64 bit) on this system will it run the 64 bit linux
version of R? I would appreciate any advice you might be willing to
offer on
this.
Tom
--
View this message in context:
http://www.nabble.com/Will-This-Computer-Run-64-bit-Ubuntu-R--tp18799376
p18799376.html
Sent from the R help mailing list archive at Nabble.com.
------------------------------
Message: 10
Date: Sun, 3 Aug 2008 20:45:19 +0530
From: "Ajay ohri" <ohri2007 at gmail.com>
Subject: [R] Installing R on Apache Web Server
To: "r-help at r-project.org" <r-help at r-project.org>
Message-ID:
<4bc14e460808030815y2b6767f0j7d1b2c383e7e6c85 at mail.gmail.com>
Content-Type: text/plain
Hi List,
I have run R from Windows Terminal server and it works fine. I now need
to
install it on apache Web server / Amazon EC 2 server to test R (with a
GUI)
for computing on a cloud (using a framework I wrote about).Intent is to
benchmark R 's superior algorithms on very big files > 1-2 gb (against
mainstream s software)
Also is there any data compression /encryption technique or package in
R.
Intent is reduce bandwidth transfer time and costs.
Regards,
Ajay
decisionstats .com
[[alternative HTML version deleted]]
------------------------------
Message: 11
Date: Sun, 3 Aug 2008 17:29:48 +0200
From: Jonas <jonas73x at gmail.com>
Subject: [R] access object inside a function
To: r-help at r-project.org
Message-ID:
<1875de320808030829u70eadd9ao6b2a3908f78c2930 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
hi,
my apologies if this has been covered numerous times before and it's
only my lack of search and/or R skills that is stopping me from
finding the solution.
i'm trying to access an object defined and created inside a function
(the object is not returned by the function), but i can't seem to get
it to work, i keep getting an "object not found" error message. i
thought the solution where to be found in changing the
frame/environment. but i can't seem to understand how do do it
correctly. what i want to do is something like:
my.function() <- {
function.from.package(x,y)
plot (object.inside.package.function)
}
i'm using R 2.7.1 on linux
sincerely,
jonas
------------------------------
Message: 12
Date: Sun, 03 Aug 2008 11:34:32 -0400
From: Duncan Murdoch <murdoch at stats.uwo.ca>
Subject: Re: [R] access object inside a function
To: Jonas <jonas73x at gmail.com>
Cc: r-help at r-project.org
Message-ID: <4895D008.60405 at stats.uwo.ca>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
On 03/08/2008 11:29 AM, Jonas wrote:
> hi,
>
> my apologies if this has been covered numerous times before and it's
> only my lack of search and/or R skills that is stopping me from
> finding the solution.
> i'm trying to access an object defined and created inside a function
> (the object is not returned by the function), but i can't seem to get
> it to work, i keep getting an "object not found" error message. i
> thought the solution where to be found in changing the
> frame/environment. but i can't seem to understand how do do it
> correctly. what i want to do is something like:
>
> my.function() <- {
> function.from.package(x,y)
> plot (object.inside.package.function)
> }
That pseudo-code doesn't really make sense: you didn't save the result
of function.from.package() anywhere. If you change it to
my.function() <- {
z <- function.from.package(x,y)
plot (z)
}
then it's perfectly standard. So the question is: where are you
expecting R to look to find object.inside.package.function?
Duncan Murdoch
>
> i'm using R 2.7.1 on linux
>
> sincerely,
> jonas
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
------------------------------
Message: 13
Date: Sun, 3 Aug 2008 08:46:22 -0700 (PDT)
From: Burgs <sjh408 at hotmail.com>
Subject: Re: [R] running strucchange?
To: r-help at r-project.org
Message-ID: <18799748.post at talk.nabble.com>
Content-Type: text/plain; charset=us-ascii
Thanks to all for your replies. It is appreciated.
Prof Brian Ripley wrote:
>
> On Sat, 2 Aug 2008, Burgs wrote:
>
>>
>> Yes, I'm using Windows. I'm getting an error when I use the source()
>> command, it seems not to like the path to the .R file I'm putting in,
>> although it is the correct path...thanks for your response anyway.
>
> See rw-FAQ Q2.16.
>
> If you want to send code to the R Console to run it, you just paste it
in.
> Provided it is a complete line (with newline at the end) it will be
run.
> Otherwise as Don says, you need to hit Return/Enter to run it.
>
>>
>>
>> Don MacQueen wrote:
>>>
>>> Most of the time I use the source() command.
>>>
>>> source('path_to_file')
>>>
>>> Maybe if you press the Return key after pasting into the console?
>>> I assume you're on Windows, which I don't use, so basically I have
no
>>> idea...
>>>
>>> -Don
>>>
>>> At 2:45 PM -0700 8/2/08, Burgs wrote:
>>>> Thank you for your post. I think my question is simply how to
"run"
>>>> code
>>>> contained in an .R text file? From the manual I see that it should
be
>>>> highlighted/selected and then run using a Ctrl character however
that
>> simply
>>>> seems to paste the code in the RConsole without actually doing
>>>> anything...
>>>>
>>>>
>>>>
>>>> Don MacQueen wrote:
>>>>>
>>>>> I would suggest imitating an example from the help page for the
>>>>> function you're trying to use.
>>>>>
>>>>> Since you haven't shown any example of what you tried, including
>>>>> error messages, there's really no way anyone can help, beyond
>>>>> generalities like mine. (Please review the posting guide, whose
URL
>>>>> is at the bottom of every email to r-help.)
>>>>>
>>>>> -Don
>>>>>
>>>>> At 2:11 PM -0700 8/2/08, Burgs wrote:
>>>>>> Greetings,
>>>>>>
>>>>>> I'm complety new to "R" and have a question. I've read
through a
>>>>>> couple
>>>>>> of manuals but I'm having a problem with getting something run
>>>>>> properly.
>>>>>> I'd like to attempt to use the "strucchange" package with some
sample
>> data
>>>>>> however I'm having trouble understanding the proper syntax of the
>> commands
>>>>>> from which to do so. I basically want to take some data I have
and
>>>>>> run
>> it
>>>>>> through strucchange so as to be able to observe the plots. Any
>>>> suggestions
>>>>>> greatly appreciated.
>>>>>> --
>>>>>> View this message in context: http:// www.
>>>>>> nabble.com/running-strucchange--tp18791489p18791489.html
>>>>>> Sent from the R help mailing list archive at Nabble.com.
>>>>>>
>>>>>> ______________________________________________
>>>>>> R-help at r-project.org mailing list
>>>>>> https:// stat.ethz.ch/mailman/listinfo/r-help
>>>>>> PLEASE do read the posting guide http:// www.
>>>> R-project.org/posting-guide.html
>>>>>> and provide commented, minimal, self-contained, reproducible
code.
>>>>>
>>>>>
>>>>> --
>>>>> ---------------------------------
>>>>> Don MacQueen
>>>>> Lawrence Livermore National Laboratory
>>>>> Livermore, CA, USA
>>>>> 925-423-1062
>>>>> macq at llnl.gov
>>>>>
>>>>> ______________________________________________
>>>>> R-help at r-project.org mailing list
>>>>> https:// stat.ethz.ch/mailman/listinfo/r-help
>>>>> PLEASE do read the posting guide
>>>>> http:// www. R-project.org/posting-guide.html
>>>>> and provide commented, minimal, self-contained, reproducible
code.
>>>>>
>>>>>
>>>>
>>>> --
>>>> View this message in context: http:// www.
>>>> nabble.com/running-strucchange--tp18791489p18793670.html
>>>> Sent from the R help mailing list archive at Nabble.com.
>>>>
>>>> ______________________________________________
>>>> R-help at r-project.org mailing list
>>>> https:// stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide http:// www.
>> R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>> --
>>> ---------------------------------
>>> Don MacQueen
>>> Lawrence Livermore National Laboratory
>>> Livermore, CA, USA
>>> 925-423-1062
>>> macq at llnl.gov
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>
>> --
>> View this message in context:
>> http://www.nabble.com/running-strucchange--tp18791489p18795602.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
> --
> Brian D. Ripley, ripley at stats.ox.ac.uk
> Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel: +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UK Fax: +44 1865 272595
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
--
View this message in context:
http://www.nabble.com/running-strucchange--tp18791489p18799748.html
Sent from the R help mailing list archive at Nabble.com.
------------------------------
Message: 14
Date: Sun, 03 Aug 2008 11:59:21 -0400
From: "John P. Burkett" <burkett at uri.edu>
Subject: [R] gmaps, projection, inset
To: R-help at r-project.org
Message-ID: <4895D5D9.8030100 at uri.edu>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Running R version 2.6.1 under Gentoo Linux, I'm trying to produce a
thematic map of the USA using the gmaps package. The result thus far
has two problems from my point of view.
First, the projection (Miller cylindrical?) elongates southern states
and flattens northern ones unattractively. I'd prefer Albers conic
projection or something similar.
Second, the inset for Alaska is placed south of the contiguous states.
I'd prefer to put it north of them.
I would be very grateful to anyone who could suggest (a) how to change
the projection and/or (b) relocate the inset.
The most relevant section of my R code is pasted below. If needed, I can
supply the entire program and data.
-John
*****************Code starts here*************************
library(gmaps)
grid.newpage()
grid.frame(name="map")
statenames <- list("Alabama", "Alaska", "Arizona", "Arkansas",
"California", "Colorado", "Connecticut",
"Florida", "Georgia", "Hawaii", "Idaho",
"Illinois", "Indiana", "Iowa", "Kansas",
"Kentucky", "Louisiana", "Maine", "Maryland",
"Massachusetts", "Michigan", "Minnesota", "Mississippi",
"Missouri", "Montana", "Nebraska", "Nevada",
"New Hampshire", "New Jersey", "New Mexico", "New York",
"North Carolina", "North Dakota", "Ohio", "Oklahoma",
"Oregon", "Pennsylvania", "Rhode Island", "South
Carolina",
"South Dakota", "Tennessee", "Texas", "Utah",
"Vermont", "Virginia", "Washington", "West Virginia",
"Wisconsin", "Wyoming", "Delaware") # Delaware is last
because it has a missing value and should not be shaded on map.
grid.pack("map",USALevelPlot(states=statenames,levels=spc,col.fun=greens
,normalize=T),height=unit(1,'null'))
*****************Code ends here***************************
--
John P. Burkett
Department of Environmental and Natural Resource Economics
and Department of Economics
University of Rhode Island
Kingston, RI 02881-0808
USA
phone (401) 874-9195
------------------------------
Message: 15
Date: Sun, 03 Aug 2008 17:05:56 +0100
From: Patrick Burns <pburns at pburns.seanet.com>
Subject: Re: [R] access object inside a function
To: Duncan Murdoch <murdoch at stats.uwo.ca>
Cc: r-help at r-project.org
Message-ID: <4895D764.7010504 at pburns.seanet.com>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Duncan Murdoch wrote:
> On 03/08/2008 11:29 AM, Jonas wrote:
>> hi,
>>
>> my apologies if this has been covered numerous times before and it's
>> only my lack of search and/or R skills that is stopping me from
>> finding the solution.
>> i'm trying to access an object defined and created inside a function
>> (the object is not returned by the function), but i can't seem to get
>> it to work, i keep getting an "object not found" error message. i
>> thought the solution where to be found in changing the
>> frame/environment. but i can't seem to understand how do do it
>> correctly. what i want to do is something like:
>>
>> my.function() <- {
>> function.from.package(x,y)
>> plot (object.inside.package.function)
>> }
>
> That pseudo-code doesn't really make sense: you didn't save the result
> of function.from.package() anywhere. If you change it to
>
> my.function() <- {
> z <- function.from.package(x,y)
> plot (z)
> }
>
> then it's perfectly standard. So the question is: where are you
> expecting R to look to find object.inside.package.function?
I suspect the OP doesn't care where R looks as
long as it finds it. The fact that everything about
'function.from.package' except 'z' has been sucked
down a black hole by the time 'plot' is on the scene
is going to make that search impossible.
One solution is to create 'my.function.from.package'
and change the return value to include
'object.inside.package.function'.
Patrick Burns
patrick at burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of S Poetry and "A Guide for the Unwilling S User")
>
> Duncan Murdoch
>
>>
>> i'm using R 2.7.1 on linux
>>
>> sincerely,
>> jonas
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
------------------------------
Message: 16
Date: Sun, 3 Aug 2008 18:13:40 +0200
From: Jonas <jonas73x at gmail.com>
Subject: Re: [R] access object inside a function
To: "Patrick Burns" <pburns at pburns.seanet.com>
Cc: r-help at r-project.org, Duncan Murdoch <murdoch at stats.uwo.ca>
Message-ID:
<1875de320808030913n372546bfq2c7e7a17800c59bc at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Duncan, Patrick,
Sorry for not being very clear (and having non-working code in my
example). What I'm actually trying to achieve is creating a wrapper
around the NelsonSiegel function in the fBonds package in order to
(among other things) create my own plot function. However, all of the
objects i would like to access inside of my.function are not returned
by the NelsonSiegel function. The are only local inside the
NelsonSiegel function. Is there a way to inside my wrapper function
access these objects?
my.function() <- {
z <- NelsonSiegel(Yield, Maturity, doplot=FALSE)
#here i would like to access objects not returned by NelsonSiegel
}
//Jonas
On Sun, Aug 3, 2008 at 6:05 PM, Patrick Burns <pburns at pburns.seanet.com>
wrote:
> Duncan Murdoch wrote:
>>
>> On 03/08/2008 11:29 AM, Jonas wrote:
>>>
>>> hi,
>>>
>>> my apologies if this has been covered numerous times before and it's
>>> only my lack of search and/or R skills that is stopping me from
>>> finding the solution.
>>> i'm trying to access an object defined and created inside a function
>>> (the object is not returned by the function), but i can't seem to
get
>>> it to work, i keep getting an "object not found" error message. i
>>> thought the solution where to be found in changing the
>>> frame/environment. but i can't seem to understand how do do it
>>> correctly. what i want to do is something like:
>>>
>>> my.function() <- {
>>> function.from.package(x,y)
>>> plot (object.inside.package.function)
>>> }
>>
>> That pseudo-code doesn't really make sense: you didn't save the
result of
>> function.from.package() anywhere. If you change it to
>>
>> my.function() <- {
>> z <- function.from.package(x,y)
>> plot (z)
>> }
>>
>> then it's perfectly standard. So the question is: where are you
>> expecting R to look to find object.inside.package.function?
>
> I suspect the OP doesn't care where R looks as
> long as it finds it. The fact that everything about
> 'function.from.package' except 'z' has been sucked
> down a black hole by the time 'plot' is on the scene
> is going to make that search impossible.
>
> One solution is to create 'my.function.from.package'
> and change the return value to include
> 'object.inside.package.function'.
>
> Patrick Burns
> patrick at burns-stat.com
> +44 (0)20 8525 0696
> http://www.burns-stat.com
> (home of S Poetry and "A Guide for the Unwilling S User")
>>
>> Duncan Murdoch
>>
>>>
>>> i'm using R 2.7.1 on linux
>>>
>>> sincerely,
>>> jonas
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>
------------------------------
Message: 17
Date: Sun, 3 Aug 2008 18:19:54 +0200
From: "pir2.jv" <pir2.jv at wanadoo.fr>
Subject: [R] (sans objet)
To: r-help at r-project.org
Message-ID: <E0320F3B-A3C8-4CE6-AD52-B949F1BF83CF at wanadoo.fr>
Content-Type: text/plain; charset=ISO-8859-1; delsp=yes; format=flowed
- Peut-on poser des questions en fran?ais?
- Tr?s simple. Je fais un programme "truc.r" (sous emacs)
Ce programme r?alise des graphiques.
Je voudrais les enregistrer au fur et ? mesure -- sinon, quand la
session est finie, ces graphiques sont perdus--
J'ai essay?:
pdf("truc.pdf")
Mais le fichier "truc.pdf" cr?? est vide.
Comment fait-on?
Merci et salutations
Jacques Vernin
------------------------------
Message: 18
Date: Sun, 3 Aug 2008 17:23:47 +0100 (BST)
From: Prof Brian Ripley <ripley at stats.ox.ac.uk>
Subject: Re: [R] drop1() seems to give unexpected results compare to
anova()
To: Thomas P C Chu <tchu11 at netscape.net>
Cc: r-help at r-project.org
Message-ID:
<alpine.LFD.1.10.0808031720390.16254 at gannet.stats.ox.ac.uk>
Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
On Sun, 3 Aug 2008, Thomas P C Chu wrote:
> Thanks to Mr Dalgaard for his advice and everyone else who has
contributed.
> Inclusion of an error term at the end of sim.set$y = ... line did cure
my
> problems with drop1() and step().
>
> I suppose it is my own inexperience in carrying out simulations caused
such
> gaffe.
R 2.8.0 will warn if asked to compute F tests that are essentially 0/0
(where 'essentially' is relative to the total sum of squares). It does
so
in your examples, so in future it will be easier to notice this.
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
------------------------------
Message: 19
Date: Sun, 3 Aug 2008 16:24:45 +0000
From: Dirk Eddelbuettel <edd at debian.org>
Subject: Re: [R] (sans objet)
To: "pir2.jv" <pir2.jv at wanadoo.fr>
Cc: r-help at r-project.org
Message-ID: <20080803162444.GA16797 at master.debian.org>
Content-Type: text/plain; charset=iso-8859-1
On Sun, Aug 03, 2008 at 06:19:54PM +0200, pir2.jv wrote:
>
> - Peut-on poser des questions en fran?ais?
Pas vraiment. R-help est une liste organisee en anglais.
> - Tr?s simple. Je fais un programme "truc.r" (sous emacs)
> Ce programme r?alise des graphiques.
> Je voudrais les enregistrer au fur et ? mesure -- sinon, quand la
> session est finie, ces graphiques sont perdus--
> J'ai essay?:
> pdf("truc.pdf")
> Mais le fichier "truc.pdf" cr?? est vide.
>
> Comment fait-on?
C'est une F.A.Q. -- il te faut un 'dev.off()' apres le (ou les)
commandes plot(), lines(), ... donc:
pdf("truc.pdf")
plot(1:10, type='b', main='Ma graphique')
dev.off()
> Merci et salutations
De rien, mais reviens 'en anglais' s'il-te-plait.
Dirk
--
Three out of two people have difficulties with fractions.
------------------------------
Message: 20
Date: Sun, 3 Aug 2008 08:30:47 -0800
From: "Franz Mueter" <fmueter at alaska.net>
Subject: Re: [R] Bubble plots
To: "'Cody Hamilton'" <Cody_Hamilton at Edwards.com>,
<r-help at r-project.org>
Message-ID: <000f01c8f586$496b4a10$dc41de30$@net>
Content-Type: text/plain; charset="us-ascii"
I believe
?symbols
is what you are looking for. Of course, you need to convert categories
to
numeric values for plotting:
For example, adding some actual data to your data frame:
> D <- cbind(D, P = runif(15))
> symbols(as.numeric(D$time), as.numeric(D$y), circles=D$P, inches=0.2,
ann=F, xaxt="n", yaxt="n")
Use axis() to get proper axis labels:
> axis(1, 1:3, levels(D$time))
> axis(2, 1:4, levels(D$y))
Franz
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of Cody Hamilton
Sent: Friday, August 01, 2008 7:07 PM
To: r-help at r-project.org
Subject: [R] Bubble plots
Is there a way to create a 'bubble plot' in R?
For example, if we define the following data frame containing the level
of y
observed for 5 patients at three time points:
time<-c(rep('time 1',5),rep('time 2',5),rep('time 3',5))
y<-c('a','b','c','d','a','b','c','a','d','a','a','a','b','c','d')
D<-data.frame(cbind(y,time))
I would like to display the percentage of subjects in each level of y at
each time point as a bubble whose size is proportional to the percentage
of
subjects in the given level of y at the given time point. Thus, in the
case
of the data frame above the plot would have the levels of y
('a','b','c','d') on the y-axis and the levels of time ('time 1','time
2',
time 3') on the x-axis with four bubbles above each time point (e.g. the
size of the bubble in the bottom left corner of the plot would be
proportional to the percentage of patients with y='a' at time='time 1').
I am running R 2.7.1 under windows.
Regards,
-Cody
________________________________
This message contains information which may be confident...{{dropped:8}}
------------------------------
Message: 21
Date: Sun, 3 Aug 2008 17:38:01 +0100 (BST)
From: Prof Brian Ripley <ripley at stats.ox.ac.uk>
Subject: Re: [R] Will This Computer Run 64 bit Ubuntu/R?
To: Tom La Bone <booboo at gforcecable.com>
Cc: r-help at r-project.org
Message-ID:
<alpine.LFD.1.10.0808031725360.16254 at gannet.stats.ox.ac.uk>
Content-Type: TEXT/PLAIN; format=flowed; charset=US-ASCII
The Pemtium D 940 supports EMT64 according to
http://www.pcstats.com/articleview.cfm?articleID=1988 (I googled it).
That is the key: EMT64 is Intel's name for x86_64 and so your box should
run x86_64 Linux which I believe Debian/Ubuntu refer to as 'amd64'
(which
understandably Intel does not favour and GNU configure does not use on
Linux).
R compiles readily on x86_64 Linux (it is our main development platform)
so in so far as there is 'the 64 bit linux version of R' (not really),
it
will run. More precisely you should be able to compile R from the
sources
and also install an amd64 Ubuntu package for R.
On Sun, 3 Aug 2008, Tom La Bone wrote:
>
> After doing some reading about 64-bit systems and software I am still
> somewhat uncertain about some things. I have a Dell Dimension XPS 400
with a
> dual core Intel Pentium D 940 (3.2 GHz) and 4 Gb of memory. I
currently dual
> boot the system with XP Professional and Ubuntu 8.04 (32 bit). If I
simply
> install Ubuntu 8.04 (64 bit) on this system will it run the 64 bit
linux
> version of R? I would appreciate any advice you might be willing to
offer on
> this.
>
> Tom
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
------------------------------
Message: 22
Date: Sun, 3 Aug 2008 16:43:31 +0000
From: Dirk Eddelbuettel <edd at debian.org>
Subject: Re: [R] Will This Computer Run 64 bit Ubuntu/R?
To: Prof Brian Ripley <ripley at stats.ox.ac.uk>
Cc: Tom La Bone <booboo at gforcecable.com>, r-help at r-project.org
Message-ID: <20080803164331.GA2156 at master.debian.org>
Content-Type: text/plain; charset=us-ascii
On Sun, Aug 03, 2008 at 05:38:01PM +0100, Prof Brian Ripley wrote:
> The Pemtium D 940 supports EMT64 according to
> http://www.pcstats.com/articleview.cfm?articleID=1988 (I googled it).
>
> That is the key: EMT64 is Intel's name for x86_64 and so your box
should
> run x86_64 Linux which I believe Debian/Ubuntu refer to as 'amd64'
(which
> understandably Intel does not favour and GNU configure does not use on
> Linux).
>
> R compiles readily on x86_64 Linux (it is our main development
platform)
> so in so far as there is 'the 64 bit linux version of R' (not really),
it
> will run. More precisely you should be able to compile R from the
sources
> and also install an amd64 Ubuntu package for R.
I'd rather just grab current prebuild binaries from any CRAN mirror, eg
http://cran.r-project.org/bin/linux/ubuntu/
as Vincent and Michael do a really excellent job rebuilding my Debian
packages for Ubuntu x86 and amd64, typically within a day. There are a
few other goodies there too, see the README in the aforementioned
directory.
Dirk
--
Three out of two people have difficulties with fractions.
------------------------------
Message: 23
Date: Sun, 03 Aug 2008 12:43:47 -0400
From: Duncan Murdoch <murdoch at stats.uwo.ca>
Subject: Re: [R] access object inside a function
To: Jonas <jonas73x at gmail.com>
Cc: r-help at r-project.org, Patrick Burns <pburns at pburns.seanet.com>
Message-ID: <4895E043.3070909 at stats.uwo.ca>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
On 03/08/2008 12:13 PM, Jonas wrote:
> Duncan, Patrick,
>
> Sorry for not being very clear (and having non-working code in my
> example). What I'm actually trying to achieve is creating a wrapper
> around the NelsonSiegel function in the fBonds package in order to
> (among other things) create my own plot function. However, all of the
> objects i would like to access inside of my.function are not returned
> by the NelsonSiegel function. The are only local inside the
> NelsonSiegel function. Is there a way to inside my wrapper function
> access these objects?
>
> my.function() <- {
> z <- NelsonSiegel(Yield, Maturity, doplot=FALSE)
> #here i would like to access objects not returned by NelsonSiegel
> }
No, as Patrick said, once the function returns all the locals (other
than z) disappear. You need to get the source to the function and
modify it.
Duncan Murdoch
>
>
> //Jonas
>
> On Sun, Aug 3, 2008 at 6:05 PM, Patrick Burns
<pburns at pburns.seanet.com> wrote:
>> Duncan Murdoch wrote:
>>> On 03/08/2008 11:29 AM, Jonas wrote:
>>>> hi,
>>>>
>>>> my apologies if this has been covered numerous times before and
it's
>>>> only my lack of search and/or R skills that is stopping me from
>>>> finding the solution.
>>>> i'm trying to access an object defined and created inside a
function
>>>> (the object is not returned by the function), but i can't seem to
get
>>>> it to work, i keep getting an "object not found" error message. i
>>>> thought the solution where to be found in changing the
>>>> frame/environment. but i can't seem to understand how do do it
>>>> correctly. what i want to do is something like:
>>>>
>>>> my.function() <- {
>>>> function.from.package(x,y)
>>>> plot (object.inside.package.function)
>>>> }
>>> That pseudo-code doesn't really make sense: you didn't save the
result of
>>> function.from.package() anywhere. If you change it to
>>>
>>> my.function() <- {
>>> z <- function.from.package(x,y)
>>> plot (z)
>>> }
>>>
>>> then it's perfectly standard. So the question is: where are you
>>> expecting R to look to find object.inside.package.function?
>> I suspect the OP doesn't care where R looks as
>> long as it finds it. The fact that everything about
>> 'function.from.package' except 'z' has been sucked
>> down a black hole by the time 'plot' is on the scene
>> is going to make that search impossible.
>>
>> One solution is to create 'my.function.from.package'
>> and change the return value to include
>> 'object.inside.package.function'.
>>
>> Patrick Burns
>> patrick at burns-stat.com
>> +44 (0)20 8525 0696
>> http://www.burns-stat.com
>> (home of S Poetry and "A Guide for the Unwilling S User")
>>> Duncan Murdoch
>>>
>>>> i'm using R 2.7.1 on linux
>>>>
>>>> sincerely,
>>>> jonas
>>>>
>>>> ______________________________________________
>>>> R-help at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide
>>>> http://www.R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
------------------------------
Message: 24
Date: Sun, 3 Aug 2008 10:25:27 -0700 (PDT)
From: "Michael A. Milligan" <pakdke at yahoo.com>
Subject: [R] missing F statistic in anova.gam
To: r-help at r-project.org
Message-ID: <418588.4927.qm at web50305.mail.re2.yahoo.com>
Content-Type: text/plain; charset=utf-8
Hello,
I have encountered results which I am not sure how to interpret when
using anova.gam to compare 2 different models. For certain tests the
results do not include an F- or associated p-statistic. This happens
when comparing certain models and not others, and I do not discern a
patten explaining when the test works and when it does not.
Here is some output for some of my tests (y#, x1, and x2 are each
1-variable vectors, while X is a matrix of several variables).
These results compare a model additively separable in x1 and x2 with a
model in which they are not assumed additively separable:
Model 1: y1 ~ s(x1) + s(x2) + X
Model 2: y1 ~ X + s(x1, x2)
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 3815.5111 29860.6
2 3810.3577 29898.8 5.1534 -38.2
No F statistic is computed, though the statistic is computed when other
dependent variables are used.
Here are some results for a similar analysis with a different dependent
variable:
Model 1: y2 ~ x1 + x2 + X
Model 2: y2 ~ s(x1) + s(x2) + X
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 3822.000 33970
2 3819.535 33921 2.465 49 2.2257 0.09578 .
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> anova(EqQuad,EqGamAS,test="F")
Analysis of Deviance Table
Model 1: y2 ~ x1 + x1sq + x2 + x2sq + x1x2 + X
Model 2: y2 ~ s(x1) + s(x2) + X
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 3819.00000 33955
2 3819.53502 33921 -0.53502 34
Model 1: y2 ~ s(x1) + s(x2) + X
Model 2: y2 ~ X + s(x1, x2)
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 3819.5350 33921
2 3821.8323 33967 -2.2973 -46 2.2391 0.09863 .
An F statistic is reported for comparing the linear model with the
additively separable semiparametric model, and for comparing the
additively separable model with the non-additvely separable model, but
not when comparing the partially quadratic model (x#sq means x#^2) with
the additively separable semiparametric model.
I'm happy to provide more information about my dataset or my estimation,
but I don't know what might be helpful, as I really don't understand at
all the cause of this problem. My dataset is not small (about 3800
observations). I will say that x2 has many observations of value 0.
I appreciate any light anyone can shed on this issue. Thank you very
much.
Michael Milligan
Ph.D. Candidate
University of New Mexico
------------------------------
Message: 25
Date: Sun, 3 Aug 2008 18:29:10 +0100 (BST)
From: Prof Brian Ripley <ripley at stats.ox.ac.uk>
Subject: Re: [R] Will This Computer Run 64 bit Ubuntu/R?
To: Dirk Eddelbuettel <edd at debian.org>
Cc: Tom La Bone <booboo at gforcecable.com>, r-help at r-project.org
Message-ID:
<alpine.LFD.1.10.0808031828210.17908 at gannet.stats.ox.ac.uk>
Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
On Sun, 3 Aug 2008, Dirk Eddelbuettel wrote:
> On Sun, Aug 03, 2008 at 05:38:01PM +0100, Prof Brian Ripley wrote:
>> The Pemtium D 940 supports EMT64 according to
>> http://www.pcstats.com/articleview.cfm?articleID=1988 (I googled it).
>>
>> That is the key: EMT64 is Intel's name for x86_64 and so your box
should
>> run x86_64 Linux which I believe Debian/Ubuntu refer to as 'amd64'
(which
>> understandably Intel does not favour and GNU configure does not use
on
>> Linux).
>>
>> R compiles readily on x86_64 Linux (it is our main development
platform)
>> so in so far as there is 'the 64 bit linux version of R' (not
really), it
>> will run. More precisely you should be able to compile R from the
sources
>> and also install an amd64 Ubuntu package for R.
>
> I'd rather just grab current prebuild binaries from any CRAN mirror,
eg
>
> http://cran.r-project.org/bin/linux/ubuntu/
>
> as Vincent and Michael do a really excellent job rebuilding my Debian
> packages for Ubuntu x86 and amd64, typically within a day. There are a
> few other goodies there too, see the README in the aforementioned
> directory.
That is what I meant by 'an amd64 Ubuntu package for R'.
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
------------------------------
Message: 26
Date: Sun, 3 Aug 2008 11:10:01 -0700 (PDT)
From: cls59 <sharpsteen at mac.com>
Subject: [R] Convert date to decimal days
To: r-help at r-project.org
Message-ID: <18800462.post at talk.nabble.com>
Content-Type: text/plain; charset=us-ascii
Hello all,
I have a quick question about formatting date strings..
I am currently debugging some Matlab code someone else wrote and since
it is
so bad that I have to go through it line by line I figured that I would
just
rewrite the thing in R.
The code produces plots of wave spectra with decimal days since the
Epoch as
the x axis and wave period as the Y axis. I am able to convert the date
stamp in the input file to a POSIXct object with the following call:
date = as.POSIXct(header[[1]][2],tz='GMT',format='%Y%m%d_%k')
date
"2008-07-11 03:00:00 GMT"
So every thing looks groovy, using as.double(date) gives me the number
of
seconds since the epoch which reproduces the correct date when fed to
date
-u -r on the command line.
However, I can't seem to find a clean way of converting the date to
decimal
days since the epoch. Using
as.double(date)/86400
Gives:
14071.12
Which is a catastrophic loss of significance, the number should be
14071.125
I tried playing around with the format command to see if I could get R
to
format the date differently, but no luck.
Any suggestions?
-Charlie
--
View this message in context:
http://www.nabble.com/Convert-date-to-decimal-days-tp18800462p18800462.h
tml
Sent from the R help mailing list archive at Nabble.com.
------------------------------
Message: 27
Date: Sun, 3 Aug 2008 11:23:43 -0700
From: "Charles C. Berry" <cberry at tajo.ucsd.edu>
Subject: Re: [R] convert for loop into apply()
To: Bill.Venables at csiro.au
Cc: R-help at r-project.org, anhnttran at ucla.edu
Message-ID: <Pine.LNX.4.64.0808031115280.4659 at tajo.ucsd.edu>
Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
If length( levels( a1$cat ) ) is not very large, and the
structure suggested in the toy example holds in the actual case (probes
in
a2 do not 'overlap' and are in order according to a2$st and a2$en, then
this might do:
> unsplit(lapply(levels(a1$cat), function(x) {
+ tmp1 <- subset( a1, cat==x )
+ tmp2 <- subset( a2, cat==x)
+ findInterval( tmp1$en, tmp2$st )-findInterval( tmp1$st, tmp2$en+1 )
+ }), a1$cat)
[1] 1 1 2 2 0 1
>
HTH,
Chuck
On Sun, 3 Aug 2008, Bill.Venables at csiro.au wrote:
> Here is a way to speed up your toy example:
>
> _________________
> a1 <- data.frame(id = 1:6,
> cat = paste('cat', rep(1:3, c(2,3,1))),
> st = c(1, 7, 30, 40, 59, 91),
> en = c(5, 25, 39, 55, 70, 120))
>
> a2 <- data.frame(id = paste('probe', 1:8),
> cat = paste('cat', rep(1:3, c(2,3,3))),
> st = c(1, 9, 20, 38, 53, 70, 80, 95),
> en = c(6, 15, 36, 43, 58, 75, 85, 98))
>
> c1 <- outer(a2$st , a1$en , "<=")
> c2 <- outer(a2$en , a1$st , ">=")
> c3 <- outer(a2$cat, a1$cat, "==")
>
> a1$coverage <- colSums(c1*c2*c3)
> __________________
>
> This won't work in one step if a1 has 30000 rows and a2 has 200000,
> unless you memory size is approximately infinite, so you will need a
> loop. Suppose you can handle 1000 probes at a time. You might be
able
> to get away with something like this:
> __________________
>
[[elided Yahoo spam]]
>
> checkCoverage <- function(a1, a2)
> colSums(outer(a2$st , a1$en , "<=") *
> outer(a2$en , a1$st , ">=") *
> outer(a2$cat, a1$cat, "=="))
>
> coverage <- numeric(N <- nrow(a2))
> m2 <- 0
> while(m2 < N) {
> m1 <- m2 + 1
> m2 <- min(m2 + chunk, N)
> coverage[m1:m2] <- checkCoverage(a1, a2[m1:m2, ])
> }
> a1$coverage <- coverage
> __________________
>
> (Warning: untested code.)
>
> Failing that, go to C code and be done with it.
>
> BTW, why did you think apply() was going to be useful here?
>
>
> Bill Venables
> CSIRO Laboratories
> PO Box 120, Cleveland, 4163
> AUSTRALIA
> Office Phone (email preferred): +61 7 3826 7251
> Fax (if absolutely necessary): +61 7 3826 7304
> Mobile: +61 4 8819 4402
> Home Phone: +61 7 3286 7700
> mailto:Bill.Venables at csiro.au
> http://www.cmis.csiro.au/bill.venables/
>
> -----Original Message-----
> From: r-help-bounces at r-project.org
[mailto:r-help-bounces at r-project.org]
> On Behalf Of Anh Tran
> Sent: Saturday, 2 August 2008 4:04 PM
> To: rlist
> Subject: [R] convert for loop into apply()
>
> Hi all,I know this topic has came up multiple times, but I've never
> fully
> understand the apply() function.
>
> Anyway, I'm here asking for your help again to convert this loop to
> apply().
>
> I have 2 data frames with the following information: a1 is the
fragment
> that
> is need to be covered, a2 is the probes that cover the specific
> fragment.
>
> I need to count the number of probes cover every given fragment (they
> need
> to have the same cat ID to be on the same fragment)
>
> a1<-data.frame(id=c(1:6), cat=c('cat 1','cat 1','cat 2','cat 2','cat
> 2','cat
> 3'), st=c(1,7,30,40,59,91), en=c(5,25,39,55,70,120));
> a2<-data.frame(id=paste('probe',c(1:8)), cat=c('cat 1','cat 1','cat
> 2','cat
> 2','cat 2','cat 3','cat 3','cat 3'), st=c(1,9,20,38,53,70,80,95),
> en=c(6,15,36,43,58,75,85,98));
> a1$coverage<-NULL;
>
> I came up with this for loop (basically, if a probe starts before the
> fragment end, and end after a fragment start, it cover that fragment)
>
> for (i in 1:length(a1$id))
> {
>
a1$coverage[i]<-length(a2[a2$st<=a1$en[i]&a2$en>=a1$st[i]&a2$cat==a1$cat
> [i],]$id);
> }
>
>> a1$coverage
> [1] 1 1 2 2 0 1
>
>
> This loop runs awefully slow when I have 200,000 probes and 30,000
> fragments. Is there anyway I can speed this up with apply()?
>
> This is the time for my for loop to scan through the first 20 record
of
> my
> dataset:
> user system elapsed
> 2.264 0.501 2.770
>
> I think there is room for improvement here. Any idea?
>
> Thanks
> --
> Regards,
> Anh Tran
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
Charles C. Berry (858) 534-2098
Dept of Family/Preventive
Medicine
E mailto:cberry at tajo.ucsd.edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego
92093-0901
------------------------------
Message: 28
Date: Sun, 3 Aug 2008 20:53:52 +0100 (BST)
From: Prof Brian Ripley <ripley at stats.ox.ac.uk>
Subject: Re: [R] Convert date to decimal days
To: cls59 <sharpsteen at mac.com>
Cc: r-help at r-project.org
Message-ID:
<alpine.LFD.1.10.0808032050520.21020 at gannet.stats.ox.ac.uk>
Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
I think you will find that the 'catastrophic loss of significance' is in
the printing at the default level. Try print(x, digits=15).
On Sun, 3 Aug 2008, cls59 wrote:
>
> Hello all,
>
> I have a quick question about formatting date strings..
Datetime strings, methinks. If you do just want days, use as.Date not
as.POSIXct.
> I am currently debugging some Matlab code someone else wrote and since
it is
> so bad that I have to go through it line by line I figured that I
would just
> rewrite the thing in R.
>
> The code produces plots of wave spectra with decimal days since the
Epoch as
> the x axis and wave period as the Y axis. I am able to convert the
date
> stamp in the input file to a POSIXct object with the following call:
>
> date = as.POSIXct(header[[1]][2],tz='GMT',format='%Y%m%d_%k')
> date
> "2008-07-11 03:00:00 GMT"
>
> So every thing looks groovy, using as.double(date) gives me the number
of
> seconds since the epoch which reproduces the correct date when fed to
date
> -u -r on the command line.
>
> However, I can't seem to find a clean way of converting the date to
decimal
> days since the epoch. Using
>
> as.double(date)/86400
>
> Gives:
>
> 14071.12
>
> Which is a catastrophic loss of significance, the number should be
14071.125
>
> I tried playing around with the format command to see if I could get R
to
> format the date differently, but no luck.
>
> Any suggestions?
>
> -Charlie
> --
> View this message in context:
http://www.nabble.com/Convert-date-to-decimal-days-tp18800462p18800462.h
tml
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
------------------------------
Message: 29
Date: Sun, 3 Aug 2008 13:19:06 -0700 (PDT)
From: cls59 <sharpsteen at mac.com>
Subject: Re: [R] Convert date to decimal days
To: r-help at r-project.org
Message-ID: <18802117.post at talk.nabble.com>
Content-Type: text/plain; charset=us-ascii
You're right, the problem was with the print function not displaying
enough
precision.
Thanks!
--
View this message in context:
http://www.nabble.com/Convert-date-to-decimal-days-tp18800462p18802117.h
tml
Sent from the R help mailing list archive at Nabble.com.
------------------------------
Message: 30
Date: Sun, 3 Aug 2008 13:55:07 -0700
From: <rkevinburton at charter.net>
Subject: [R] Determining model parameters
To: r-help at r-project.org
Message-ID: <20080803165507.GXQ3D.82014.root at mp17>
Content-Type: text/plain; charset=utf-8
This may be a begining question. If so, please bear with me.
If I have some data that based on the historgram and other plots it
"looks" like a beta distribution. Is there a function or functions
within R to help me determine the model parameters for such a
distirbution? Similarily for other "common" distirbutions,
Poisson(lambda), Chi-Square(degrees of freedom, chi-square statistic),
etc. For a normal distribution it is relatively straight forward as mean
and variance can be computed failrly readily.
So far in my reading such parameters are more or less given but to do
any kind of inference using a prior model it seems key that these
parameters be at least estimated accurately.
Thank you.
Kevin
------------------------------
Message: 31
Date: Sun, 3 Aug 2008 21:06:31 +0000 (UTC)
From: Ben Bolker <bolker at ufl.edu>
Subject: Re: [R] Determining model parameters
To: r-help at stat.math.ethz.ch
Message-ID: <loom.20080803T210608-375 at post.gmane.org>
Content-Type: text/plain; charset=us-ascii
<rkevinburton <at> charter.net> writes:
> If I have some data that based on the historgram and other plots it
"looks"
like a beta distribution. Is there
> a function or functions within R to help me determine the model
parameters for
such a distirbution?
library(MASS)
?fitdistr
------------------------------
Message: 32
Date: Sun, 3 Aug 2008 16:54:43 -0400
From: Andrew Ramsey <aramsey at schoolph.umass.edu>
Subject: [R] Changing values
To: r-help at r-project.org
Message-ID: <58C18421-686A-4D5E-A952-29D547E07434 at schoolph.umass.edu>
Content-Type: text/plain; charset=US-ASCII; delsp=yes; format=flowed
Hello--
I am a relatively new user to R and I cannot find the information I
need. Please help.
I have a very large data set with values including letters, numbers,
and symbols (sometimes within the same vector value [ie X9-].
I've imported the data using read.fwp and it arrives in list format.
I'd like to change the letters and symbols to numbers (ie X9- ->
00911) in every entry.
How would you recommend I try to do so?
Thank you!
------------------------------
Message: 33
Date: Sun, 3 Aug 2008 17:55:17 -0400
From: "jim holtman" <jholtman at gmail.com>
Subject: Re: [R] Changing values
To: "Andrew Ramsey" <aramsey at schoolph.umass.edu>
Cc: r-help at r-project.org
Message-ID:
<644e1f320808031455s3b08c840v3338774f06aed855 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Exactly how is the translation of "X9-" -> 00911 done? Are there
unique mappings of character sequences to numbers? How many different
ones might there be? Why do you have leading zeros on the result?
If they are changed to numeric, then the default printing results in
'911'. Can you provide an idea of what you data looks like - e.g.,
str(yourData). Is this to be done on all, or selected, columns? You
need to provide a bit more background so that we can understand the
problem you are trying to solve.
On Sun, Aug 3, 2008 at 4:54 PM, Andrew Ramsey
<aramsey at schoolph.umass.edu> wrote:
> Hello--
>
> I am a relatively new user to R and I cannot find the information I
need.
> Please help.
>
> I have a very large data set with values including letters, numbers,
and
> symbols (sometimes within the same vector value [ie X9-].
>
> I've imported the data using read.fwp and it arrives in list format.
>
> I'd like to change the letters and symbols to numbers (ie X9- ->
00911) in
> every entry.
>
> How would you recommend I try to do so?
>
> Thank you!
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Jim Holtman
Cincinnati, OH
+1 513 646 9390
What is the problem that you are trying to solve?
------------------------------
Message: 34
Date: Sun, 3 Aug 2008 23:13:46 +0000 (UTC)
From: Michael Bibo <michael_bibo at health.qld.gov.au>
Subject: Re: [R] Bubble plots
To: r-help at stat.math.ethz.ch
Message-ID: <loom.20080803T225600-446 at post.gmane.org>
Content-Type: text/plain; charset=us-ascii
Cody Hamilton <Cody_Hamilton <at> Edwards.com> writes:
>
> Is there a way to create a 'bubble plot' in R?
>
> For example, if we define the following data frame containing the
level of y
observed for 5 patients at three
> time points:
>
> time<-c(rep('time 1',5),rep('time 2',5),rep('time 3',5))
> y<-c('a','b','c','d','a','b','c','a','d','a','a','a','b','c','d')
> D<-data.frame(cbind(y,time))
>
> I would like to display the percentage of subjects in each level of y
at
each time point as a bubble whose size
> is proportional to the percentage of subjects in the given level of y
at the
given time point. Thus, in the
> case of the data frame above the plot would have the levels of y
('a','b','c','d') on the y-axis and the
> levels of time ('time 1','time 2', time 3') on the x-axis with four
bubbles
above each time point (e.g. the
> size of the bubble in the bottom left corner of the plot would be
proportional to the percentage of patients
> with y='a' at time='time 1').
It sounds like function balloonplot from package gplots might be just
what you
are looking for.
There is, however, a bug in the labelling of the plot that I have been
meaning
to contact the maintainer about. If you change lines 364 and 376
respectively
from "labels=ynames," to "labels=ylab, and "labels=xnames," to
"labels=xlab,",
it will work for this specific purpose. I have included the full code
for the
amended balloonplot function to the end of the email. You can just copy
and
paste into your R console.
Then, all that is required (given your dataframe D) is:
attach(D)
balloonplot(prop.table(table(time,y))*100)
This gives percentages. You can play with prop.table syntax to get
percentages of row or column totals, and play with titles, margin totals
etc
from there.
The amended balloonplot function (this is not an all-purpose fix - it
just
works for this specific purpose):
# $Id: balloonplot.R 908 2006-03-02 21:43:24Z warnes $
balloonplot <- function(x,...)
UseMethod("balloonplot",x)
balloonplot.table <- function(x, xlab, ylab, zlab, show.zeros = FALSE,
show.margins = TRUE, ... )
{
obj <- x
tmp <- as.data.frame(x)
x <- tmp[,1]
y <- tmp[,2]
z <- tmp[,3]
tableflag <- TRUE
if(missing(xlab)) xlab <- names(dimnames(obj))[1]
if(missing(ylab)) ylab <- names(dimnames(obj))[2]
if(missing(zlab)) zlab <- "Freq"
balloonplot.default(x, y, z, xlab=xlab, ylab=ylab, zlab=zlab,
show.zeros = show.zeros,
show.margins = show.margins, ...)
}
balloonplot.default <- function(x,y,z,
xlab,
ylab,
zlab=deparse(substitute(z)),
dotsize=2/max(strwidth(19),strheight(19)),
dotchar=19,
dotcolor="skyblue",
main,
label=TRUE,
label.digits=2,
scale.method=c("volume","diameter"),
colsrt=par("srt"),
rowsrt=par("srt"),
colmar=1,
rowmar=2,
show.zeros=FALSE,
show.margins=TRUE,
cum.margins=TRUE,
sorted=TRUE,
label.lines=TRUE,
fun=function(x)sum(x,na.rm=T),
hide.duplicates=TRUE,
... )
{
if(is.null(names(x)))
{
xnames <- as.character(substitute(x))
if(length(xnames)>1) xnames <- xnames[-1]
}
else
xnames <- names(x)
if(is.null(names(y)))
{
ynames <- as.character(substitute(y))
if(length(ynames)>1) ynames <- ynames[-1]
}
else
ynames <- names(y)
if(missing(xlab))
xlab <- paste( xnames, collapse=", " )
if(missing(ylab))
ylab <- paste( ynames, collapse=", " )
####
## Handle arguments
####
scale.method <- match.arg(scale.method)
if( any(z < 0 ) )
warning("z value(s) below zero detected.",
" No balloons will be displayed for these cells.")
if(missing(main))
{
if(scale.method=="volume")
main <- paste("Balloon Plot for ",
paste(xnames, collapse=", "),
" by ",
paste(ynames, collapse=", "),
".\nArea is proportional to ", zlab, ".", sep='')
else
main <- paste("Balloon Plot for ",
paste(ynames, collapse=", "),
" by ",
paste(ynames, collapse=", "),
".\nDiameter is proportional to ", zlab, ".",
sep='')
}
if(length(dotcolor)<length(z))
dotcolor <- rep(dotcolor, length=length(z))
####
## Make sure x and y are lists
####
if(is.list(x))
{
xlabs <- x
x$sep=":"
x <- do.call(paste, x)
}
else
xlabs <- list(x)
if(is.list(y))
{
ylabs <- y
y$sep=":"
y <- do.call(paste, y)
ylab <- paste( names(y) )
}
else
ylabs <- list(y)
####
## sort everything into a useful order
####
if(sorted)
{
ord.x <- do.call(order, xlabs)
ord.y <- do.call(order, ylabs)
}
else
ord.x <- ord.y <- 1:length(x)
forceOrder <- function(X, sord, lord)
factor(X[sord], levels=unique(X[lord]))
x <- forceOrder(x, ord.y, ord.y)
y <- forceOrder(y, ord.y, ord.y)
z <- as.numeric(z[ord.y])
dotcolor <- dotcolor[ord.y]
xlabs <- unique(data.frame(lapply(xlabs, forceOrder,
sord=ord.y, lord=ord.y)))
ylabs <- unique(data.frame(lapply(ylabs, forceOrder,
sord=ord.y, lord=ord.y)))
####
## Function to scale circles to fill the containing box
####
scale <- function(X, min=0, max=16, scale.method)
{
if(scale.method=="volume")
{
X[X<0] <- 0
X <- sqrt(X)
}
X <- min + (X/max(X, na.rm=TRUE) * (max - min) )
cin.x <- par("cin")[1]
cin.y <- par("cin")[2]
if(cin.x < cin.y) X <- X * cin.x/cin.y
X
}
nlabels.y <- length(ylabs)
nlabels.x <- length(xlabs)
####
## Combine duplicate entries
####
# Do twice, once for data, once for colors
tab1 <- split( data.frame(z,dotcolor,x,y), f=list(x,y) )
ztab <- do.call(rbind,
lapply(
tab1,
FUN=function(X) cbind(z=fun(X[,1]),X[1,-1])
)
)
####
## Do the plotting
###
oldpar <- par("xpd","mar")
on.exit( par(oldpar) )
#par(xpd=NA, mar=c(1,1,5,1)+0.1) # clip drawing to device region
if(!show.margins)
{
xlim=c(-0.5,nlevels(x)+nlabels.y*rowmar-0.25) # extra space on
either
# end of plot for
labels
ylim=c(0.50,nlevels(y)+nlabels.x*colmar+1) # and so dots don't
cross
# into margins,
}
else
{
xlim=c(-0.5,nlevels(x)+nlabels.y*rowmar+1) # extra space on
either
# end of plot for
labels
ylim=c(0,nlevels(y)+nlabels.x*colmar+1) # and so dots don't cross
# into margins,
}
plot(x=nlabels.y*rowmar+0.25 + as.numeric(ztab$x) - 1,
y=nlevels(y) - as.numeric(ztab$y) + 1,
cex=scale(ztab$z, max=dotsize, scale.method=scale.method),
pch=dotchar, # plot character
col=as.character(ztab$dotcolor), # dot color
xlab="",
ylab="",
xaxt="n", # no x axis lables
yaxt="n", # no y axis lables
bty="n", # no box around the plot
xaxs = "i",
yaxs = "i",
xlim=xlim,
ylim=ylim,
...
)
ny <- nlevels(ztab$y)
nx <- nlevels(ztab$x)
sumz <- sum(ztab$z, na.rm=TRUE)
colsumz <- sapply(split( ztab$z, ztab$y), sum, na.rm=TRUE) # works
rowsumz <- sapply(split( ztab$z, ztab$x), sum, na.rm=TRUE) # broken
if(show.margins)
{
## column totals
text(
x=(1:nx) + nlabels.y*rowmar + 0.25 -1,
y=0.25,
labels=format(c(sumz, rowsumz), digits=label.digits)[-1],
font=1,
cex=par("cex")*0.75,
adj=c(0.5,0.0)
)
## row totals
rowlabs <- format(c(sumz, colsumz), digits=label.digits)[-1]
width <- max(strwidth(rowlabs),na.rm=TRUE)
text(
x=nx + nlabels.y*rowmar-0.25+width,
y= (ny:1),
labels=rowlabs,
font=1,
cex=par("cex")*0.75,
adj=c(1.0,0.5)
)
## overall total
text(
x=nx + nlabels.y*rowmar-0.25+width,
y=0.25,
labels=sumz,
font=1,
cex=par("cex")*0.75,
adj=c(1.0,0.0)
)
}
if(cum.margins)
{
## Row Sums at left
cx <- c(0, cumsum(rowsumz) / sumz)
rect(xleft = nlabels.y*rowmar - 1 - 0.25 + 1:nx,
xright = nlabels.y*rowmar - 1 + 0.75 + 1:nx,
ybottom = ny+0.75+cx[1:nx]*colmar*nlabels.x,
ytop = ny+0.75+cx[2:(nx+1)]*colmar*nlabels.x,
col = "lightgray",
border = NA)
## Col Sums at top
cy <- c(0, cumsum(colsumz) / sumz)
rect(xleft = -0.5 +rowmar*cy[ny:1]*nlabels.y,
xright = -0.5 +rowmar*cy[(ny+1):2]*nlabels.y,
ybottom = 1:ny-0.5,
ytop = 1:ny+0.5,
col = "lightgray",
border = NA)
tx <- paste(levels(x),"\n[",rowsumz,"]")
ty <- paste(levels(y),"\n[",colsumz,"]")
}
###
## Horizontal borders between cells
###
segments(
x0=nlabels.y*rowmar-0.25,
x1=nx+nlabels.y*rowmar-0.25,
y0=(0:ny)+0.5,
y1=(0:ny)+0.5
)
###
## Vertical borders between cells
###
segments(
x0=(0:nx)+nlabels.y*rowmar-0.25,
x1=(0:nx)+nlabels.y*rowmar-0.25,
y0= 0.5,
y1=ny+0.5,
)
if(hide.duplicates)
undupe <- function(X)
{
# convert duplicates into blanks
X <- as.character(X)
c(X[1], ifelse(X[-1] == X[-length(X)], "", X[-1]))
}
else
undupe <- function(X) X
###
## Column labels
###
for(i in 1:nlabels.x)
{
y <- ny + 0.75 + (nlabels.x - i + .5)*colmar
text(
x= (1:nx) + nlabels.y*rowmar + 0.25 - 1,
y= y,
labels=undupe(xlabs[,i]),
srt=colsrt,
font=1
)
}
###
## Row labels
###
for(i in 1:length(ylabs))
{
text(
y=ny:1,
x= (i-0.5)*rowmar-0.5,
labels=undupe(ylabs[,i]),
srt=rowsrt,
font=1
)
}
####
## Column headers for row labels
####
text(
x=((1:length(ylabs))-0.5)*rowmar-0.5,
y=ny+0.5,
labels=ylab,
srt=colsrt,
font=2,
adj=c(0.5,0.0)
)
####
## Row headers for column labels
####
text(
x= nlabels.y*rowmar - 0.25 - strwidth(','),
y= ny + 0.75 + ((nlabels.x:1) - 1 + .5)*colmar,
labels=xlab,
srt=colsrt,
font=2,
adj=c(1,0.5)
)
###
## add borders to row and column headers
###
if(label.lines)
{
segments( # left: vertical lines
x0=(0:nlabels.y)*rowmar-0.5,
x1=(0:nlabels.y)*rowmar-0.5,
y0=0.5,
y1=ny+0.5
)
segments(
x0=nlabels.y*rowmar-0.25, # top: horizontal lines
x1=nlabels.y*rowmar + nx - 0.25,
y0=(0:nlabels.x)*colmar +ny+0.75,
y1=(0:nlabels.x)*colmar +ny+0.75
)
}
####
## annotate cells with actual values
####
if(label){
if(show.zeros)
indiv <- 1:length(ztab$y)
else
indiv <- which(ztab$z != 0)
text(x=as.numeric(ztab$x[indiv])+ nlabels.y*rowmar - 0.75, #
as.numeric give numeric values
y=ny - as.numeric(ztab$y[indiv]) + 1,
labels=format(ztab$z[indiv], digits=label.digits), #
label value
col="black", # text color
font=2,
adj=c(0.5,0.5)
)
}
# put a nice title
title(main=main)
}
------------------------------
Message: 35
Date: Sun, 03 Aug 2008 20:39:06 -0400
From: Ben Bolker <bolker at ufl.edu>
Subject: Re: [R] Determining model parameters
To: rkevinburton at charter.net
Cc: r-help at r-project.org
Message-ID: <48964FAA.4020203 at ufl.edu>
Content-Type: text/plain; charset=UTF-8; format=flowed
~ The beta distribution only applies to data that are bounded
between 0 and 1 (in some cases strictly bounded, i.e. 0<x<1 not
just 0<=x<=1).
~ As a start, try something like
d2 = (Diff-min(Diff)+0.001)/(max(Diff)-min(Diff)+0.002)
fitdistr(d2,densfun="beta",start=list(shape1=3,shape2=2))
~ [cc'ing back to R-help for posterity ...]
~ Ben Bolker
rkevinburton at charter.net wrote:
| OK. I see in the documentation that "start" is required for certain
distributions. Now I get:
|
| fitdistr(Diff, densfun="beta", list(shape1=3,shape2=2))
| Error in optim(x = c(0, 516, 968, 703, 960, 32, 956, 1417, -92, 471,
531, :
| initial value in 'vmmin' is not finite
|
| Any idea what vmmin is? And what the cause of the error is?
|
| Thanks again.
|
| Kevin
| ---- Ben Bolker <bolker at ufl.edu> wrote:
|> <rkevinburton <at> charter.net> writes:
|>
|>> If I have some data that based on the historgram and other plots it
"looks"
|> like a beta distribution. Is there
|>> a function or functions within R to help me determine the model
parameters for
|> such a distirbution?
|>
|>
|> library(MASS)
|> ?fitdistr
|>
------------------------------
Message: 36
Date: Sun, 3 Aug 2008 20:00:15 -0700
From: J Dougherty <jwd at surewest.net>
Subject: Re: [R] Changing values
To: r-help at r-project.org
Message-ID: <200808032000.15726.jwd at surewest.net>
Content-Type: text/plain; charset="iso-8859-1"
On Sunday 03 August 2008 01:54:43 pm Andrew Ramsey wrote:
> Hello--
>
> I am a relatively new user to R and I cannot find the information I
> need. Please help.
>
> I have a very large data set with values including letters, numbers,
> and symbols (sometimes within the same vector value [ie X9-].
>
> I've imported the data using read.fwp and it arrives in list format.
>
> I'd like to change the letters and symbols to numbers (ie X9- ->
> 00911) in every entry.
>
> How would you recommend I try to do so?
>
> Thank you!
>
Andrew,
Your information is a little sparse, which limits how helpful list
members can
be. What is the purpose to recoding the values? By "read.fwp" did you
mean "read.fwf" as in fixed-wdith format? Looking at your "before"
and "after" example, the leading zeros in the "after" example value
would
seem to mean the value is STILL an alphanumeric string rather than a
"number"
per se. It might order more neatly in a printout than common ascii, but
it
is really no different and should still be a list. As a rule, I would
recommend handling tasks related to data coding in a dbms such as
PostgreSQL
or Access. It is afterall what they are for.
John
------------------------------
Message: 37
Date: Mon, 4 Aug 2008 14:54:45 +1000
From: Sorn.Norng at dpi.vic.gov.au
Subject: [R] Sweave and ggplot2
To: r-help at r-project.org
Message-ID:
<OFE0099DDC.32610562-ONCA25749B.001A7FF8-CA25749B.001AFC3E at nre.vic.gov.a
u>
Content-Type: text/plain
Hi all,
I've been trying to run Sweave with R code embedded - using the ggplot2
package and in particular the qplot command. There appears to be a
problem
in Sweave not picking up that qplot is a function. Has anybody else
tried
to use qplot in Sweave and have you been successful? Any help would be
very much appreciated.
Kind Regards,
Sorn
Notice:
This email and any attachments may contain information
t...{{dropped:14}}
------------------------------
Message: 38
Date: Mon, 4 Aug 2008 16:32:14 +1000
From: Andrew Robinson <A.Robinson at ms.unimelb.edu.au>
Subject: [R] Decomposing tests of interaction terms in mixed-effects
models
To: R-Help Discussion <r-help at stat.math.ethz.ch>
Message-ID: <20080804063214.GN3754 at ms.unimelb.edu.au>
Content-Type: text/plain; charset=us-ascii
Dear R colleagues,
a friend and I are trying to develop a modest workflow for the problem
of decomposing tests of higher-order terms into interpretable sets of
tests of lower order terms with conditioning.
For example, if the interaction between A (3 levels) and C (2 levels)
is significant, it may be of interest to ask whether or not A is
significant at level 1 of C and level 2 of C.
The textbook response seems to be to subset the data and perform the
tests on the two subsets, but using the RSS and DF from the original
fit.
We're trying to answer the question using new explanatory variables.
This approach (seems to) work just fine in aov, but not lme.
For example,
########################################################################
##
# Build the example dataset
set.seed(101)
Block <- gl(6, 6, 36, lab = paste("B", 1:6, sep = ""))
A <- gl(3, 2, 36, lab = paste("A", 1:3, sep = ""))
C <- gl(2, 1, 36, lab = paste("C", 1:2, sep = ""))
example <- data.frame(Block, A, C)
example$Y <- rnorm(36)*2 + as.numeric(example$A)*as.numeric(example$C)^2
+
3 * rep(rnorm(6), each=6)
with(example, interaction.plot(A, C, Y))
# lme
require(nlme)
anova(lme(Y~A*C, random = ~1|Block, data = example))
summary(aov(Y ~ Error(Block) + A*C, data = example))
# Significant 2-way interaction. Now we would like to test the effect
# of A at C1 and the effect of A at C2. Construct two new variables
# that decompose the interaction, so an LRT is possible.
example$AC <- example$AC1 <- example$AC2 <- interaction(example$A,
example$C)
example$AC1[example$C == "C1"] <- "A1.C1" # A is constant at C1
example$AC2[example$C == "C2"] <- "A1.C2" # A is constant at C2
# lme fails (as does lmer)
lme(Y ~ AC1 + AC, random = ~ 1 | Block, data = example)
# but aov seems just fine.
summary(aov(Y ~ Error(Block) + AC1 + AC, data = example))
## AC was not significant, so A is not significant at C1
summary(aov(Y ~ Error(Block) + AC2 + AC, data = example))
## AC was significant, so A is significant at C2
## Some experimentation suggests that lme doesn't like the 'partial
## confounding' approach that we are using, rather than the variables
## that we have constructed.
lme(Y ~ AC, random = ~ 1 | Block, data = example)
lme(Y ~ A + AC, random = ~ 1 | Block, data = example)
########################################################################
##
Are we doing anything obviously wrong? Is there another approach to
the goal that we are trying to achieve?
Many thanks,
Andrew
--
Andrew Robinson
Department of Mathematics and Statistics Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
------------------------------
Message: 39
Date: Mon, 4 Aug 2008 16:19:28 +0930
From: Fernando Marmolejo Ramos
<fernando.marmolejoramos at adelaide.edu.au>
Subject: [R] about the 95%CI around the median...
To: r-help at r-project.org
Message-ID: <1217832568.4896a678c5c5e at webmail.adelaide.edu.au>
Content-Type: text/plain; charset=ISO-8859-1
Dear people
I've learnt that by using the "boxplot.stats" command in the "grDevices"
library
I can get the 5-number summaries of a boxplot, plus other important
information,
like the confidence interval around the median.
I'm interested in knowing the actual formula to used in that package to
calculate that confidence interval.
Can someone help me with this?
Cheers,
Fernando
------------------------------
Message: 40
Date: Mon, 04 Aug 2008 17:00:09 +1000
From: Simon Blomberg <s.blomberg1 at uq.edu.au>
Subject: Re: [R] about the 95%CI around the median...
To: Fernando Marmolejo Ramos <fernando.marmolejoramos at adelaide.edu.au>
Cc: r-help <r-help at r-project.org>
Message-ID: <1217833209.17713.40.camel at sib-sblomber01d.sib.uq.edu.au>
Content-Type: text/plain
See ?fivenum in the stats package. If you just type
stats::fivenum
you will get the code. The crucial calculations are in the last few
lines.
Simon.
On Mon, 2008-08-04 at 16:19 +0930, Fernando Marmolejo Ramos wrote:
> Dear people
>
> I've learnt that by using the "boxplot.stats" command in the
"grDevices" library
> I can get the 5-number summaries of a boxplot, plus other important
information,
> like the confidence interval around the median.
>
> I'm interested in knowing the actual formula to used in that package
to
> calculate that confidence interval.
>
> Can someone help me with this?
>
> Cheers,
>
> Fernando
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Simon Blomberg, BSc (Hons), PhD, MAppStat.
Lecturer and Consultant Statistician
Faculty of Biological and Chemical Sciences
The University of Queensland
St. Lucia Queensland 4072
Australia
Room 320 Goddard Building (8)
T: +61 7 3365 2506
http://www.uq.edu.au/~uqsblomb
email: S.Blomberg1_at_uq.edu.au
Policies:
1. I will NOT analyse your data for you.
2. Your deadline is your problem.
The combination of some data and an aching desire for
an answer does not ensure that a reasonable answer can
be extracted from a given body of data. - John Tukey.
------------------------------
Message: 41
Date: Mon, 4 Aug 2008 09:26:40 +0200
From: MORNEAU Fran?ois <francois.morneau at ifn.fr>
Subject: Re: [R] graph
To: "elyakhlifi mustapha" <elyakhlifi_mustapha at yahoo.fr>
Cc: r-help at r-project.org
Message-ID: <5E3D22A4869BB94AA1138AB97660D8B6013CF672 at POPULUS.ifn.fr>
Content-Type: text/plain; charset="iso-8859-1"
Hello,
Is ?segments what you are looking for ?
Fran?ois
-----Message d'origine-----
De?: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
De la part de elyakhlifi mustapha
Envoy??: vendredi 1 ao?t 2008 17:15
??: r-help at r-project.org
Objet?: [R] graph
Hello.
?
I don't know how to do to ouput segments between points like the chart
attached.
Can you help me please?
Thanks.
________________________________________________________________________
____
.yahoo.fr
------------------------------
Message: 42
Date: Mon, 4 Aug 2008 09:56:02 +0200
From: "Draga, R." <R.Draga at umcutrecht.nl>
Subject: [R] Major difference in multivariate analyses SPSS and R
To: <r-help at r-project.org>
Message-ID:
<C222D04679FF794EA98B823041802F1CE27C20 at EXV4.ds.umcutrecht.nl>
Content-Type: text/plain
Dear colleagues,
I know SPSS can not compute linear mixed models. I used 'R' before for
computing multivariate analyses. But, I never encountered such a major
difference in outcome between SPSS and 'R':
In SPSS the Pearson correlation between variable 1 and variable 2 is 31%
p<0.001.
In SPSS binary logistic regression gives us an Odds Ratio (OR)=4.9 (95%
CI 2.7-9.0), p<0.001, n=338.
OR lower upper
gender 1,120 0,565 2,221
age 0,985 0,956 1,015
variable 2 4,937 2,698 9,032
In R multilevel logistic regression using statistical package 'lmer'
gives us an Odds Ratio=10.2 (95% CI 6.3-14), p=0.24, n=338, groups:
group 1, 98; group 2 84.
OR lower upper
gender 2,295 -2,840 7,430
age 0,003 -70,047 70,054
variable 2 10,176 6,295 14,056
The crosstabs gives us:
variable A
Var B 0 1
0 156 108
1 17 57
Would somebody know how it is possible that in SPSS we get p<0.001 and
in R we get p=0.24?
And, in 'R' the 95% CI of the Odds Ratio is 6.2-14.1. Why is the
p-value=0.24?
Thanks,
Ronald
[[alternative HTML version deleted]]
------------------------------
Message: 43
Date: Mon, 4 Aug 2008 10:02:34 +0200
From: "ONKELINX, Thierry" <Thierry.ONKELINX at inbo.be>
Subject: Re: [R] Sweave and ggplot2
To: <Sorn.Norng at dpi.vic.gov.au>, <r-help at r-project.org>
Message-ID:
<2E9C414912813E4EB981326983E0A1040527D64B at inboexch.inbo.be>
Content-Type: text/plain; charset="us-ascii"
Dear Sorn,
It's hard to guess what your problem is, as you don't provide any sample
code. My guess is that the graphics are empty. Did you use
print(qplot(...)) or just qplot(). The latter won't work. You need
print(qplot(...))
HTH,
Thierry
------------------------------------------------------------------------
----
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey
-----Oorspronkelijk bericht-----
Van: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
Namens Sorn.Norng at dpi.vic.gov.au
Verzonden: maandag 4 augustus 2008 6:55
Aan: r-help at r-project.org
Onderwerp: [R] Sweave and ggplot2
Hi all,
I've been trying to run Sweave with R code embedded - using the ggplot2
package and in particular the qplot command. There appears to be a
problem
in Sweave not picking up that qplot is a function. Has anybody else
tried
to use qplot in Sweave and have you been successful? Any help would be
very much appreciated.
Kind Regards,
Sorn
Notice:
This email and any attachments may contain information\
...{{dropped:10}}
------------------------------
Message: 44
Date: Mon, 04 Aug 2008 10:17:38 +0200
From: Peter Dalgaard <p.dalgaard at biostat.ku.dk>
Subject: Re: [R] Decomposing tests of interaction terms in
mixed-effects models
To: Andrew Robinson <A.Robinson at ms.unimelb.edu.au>
Cc: R-Help Discussion <r-help at stat.math.ethz.ch>
Message-ID: <4896BB22.70302 at biostat.ku.dk>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Andrew Robinson wrote:
> Dear R colleagues,
>
> a friend and I are trying to develop a modest workflow for the problem
> of decomposing tests of higher-order terms into interpretable sets of
> tests of lower order terms with conditioning.
>
> For example, if the interaction between A (3 levels) and C (2 levels)
> is significant, it may be of interest to ask whether or not A is
> significant at level 1 of C and level 2 of C.
>
> The textbook response seems to be to subset the data and perform the
> tests on the two subsets, but using the RSS and DF from the original
> fit.
>
> We're trying to answer the question using new explanatory variables.
> This approach (seems to) work just fine in aov, but not lme.
>
> For example,
>
>
########################################################################
##
>
> # Build the example dataset
>
> set.seed(101)
>
> Block <- gl(6, 6, 36, lab = paste("B", 1:6, sep = ""))
> A <- gl(3, 2, 36, lab = paste("A", 1:3, sep = ""))
> C <- gl(2, 1, 36, lab = paste("C", 1:2, sep = ""))
> example <- data.frame(Block, A, C)
> example$Y <- rnorm(36)*2 +
as.numeric(example$A)*as.numeric(example$C)^2 +
> 3 * rep(rnorm(6), each=6)
>
> with(example, interaction.plot(A, C, Y))
>
> # lme
>
> require(nlme)
> anova(lme(Y~A*C, random = ~1|Block, data = example))
>
> summary(aov(Y ~ Error(Block) + A*C, data = example))
>
> # Significant 2-way interaction. Now we would like to test the effect
> # of A at C1 and the effect of A at C2. Construct two new variables
> # that decompose the interaction, so an LRT is possible.
>
> example$AC <- example$AC1 <- example$AC2 <- interaction(example$A,
example$C)
>
> example$AC1[example$C == "C1"] <- "A1.C1" # A is constant at C1
> example$AC2[example$C == "C2"] <- "A1.C2" # A is constant at C2
>
> # lme fails (as does lmer)
>
> lme(Y ~ AC1 + AC, random = ~ 1 | Block, data = example)
>
> # but aov seems just fine.
>
> summary(aov(Y ~ Error(Block) + AC1 + AC, data = example))
>
> ## AC was not significant, so A is not significant at C1
>
> summary(aov(Y ~ Error(Block) + AC2 + AC, data = example))
>
> ## AC was significant, so A is significant at C2
>
> ## Some experimentation suggests that lme doesn't like the 'partial
> ## confounding' approach that we are using, rather than the variables
> ## that we have constructed.
>
> lme(Y ~ AC, random = ~ 1 | Block, data = example)
> lme(Y ~ A + AC, random = ~ 1 | Block, data = example)
>
>
########################################################################
##
>
> Are we doing anything obviously wrong? Is there another approach to
> the goal that we are trying to achieve?
>
This looks like a generic issue with lme/lmer, in that they are not
happy with rank deficiencies in the design matrix.
Here's a "dirty" trick which you are not really supposed to know about
because it is hidden inside the "stats" namespace:
M <- model.matrix(~AC1+AC, data=example)
M2 <- stats:::Thin.col(M)
summary(lme(Y ~ M2 - 1, random = ~ 1 | Block, data = example)
(Thin.col(), Thin.row(), and Rank() are support functions for
anova.mlm() et al., but maybe we should document them and put them out
in the open.)
--
O__ ---- Peter Dalgaard ?ster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
------------------------------
Message: 45
Date: Mon, 04 Aug 2008 10:44:57 +0200
From: Peter Dalgaard <p.dalgaard at biostat.ku.dk>
Subject: Re: [R] Major difference in multivariate analyses SPSS and R
To: "Draga, R." <R.Draga at umcutrecht.nl>
Cc: r-help at r-project.org
Message-ID: <4896C189.7080501 at biostat.ku.dk>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Draga, R. wrote:
> .....
>
>
> Would somebody know how it is possible that in SPSS we get p<0.001 and
> in R we get p=0.24?
>
> And, in 'R' the 95% CI of the Odds Ratio is 6.2-14.1. Why is the
> p-value=0.24?
>
>
>
>
You asked this before and you are still not providing sufficient
information.
Harold already pointed out that excluding an important variance
component in SPSS could easily explain a difference in p values.
However, the discrepancy between CI and p value suggests that you could
be misreading the output from lmer(), and/or misspecifying the model.
Now, you are simply not telling us
- what you did
- what the data are like (esp. which variables vary within and between
levels)
- what the output was
- which version of lmer() was used
Also you are confusing us with "variable 1", "variable 2", "variable A",
and "Var B"
--
O__ ---- Peter Dalgaard ?ster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
------------------------------
Message: 46
Date: Mon, 04 Aug 2008 09:57:25 +0100
From: Gavin Simpson <gavin.simpson at ucl.ac.uk>
Subject: Re: [R] about the 95%CI around the median...
To: Simon Blomberg <s.blomberg1 at uq.edu.au>
Cc: r-help <r-help at r-project.org>
Message-ID: <1217840246.2655.10.camel at graptoleberis.geog.ucl.ac.uk>
Content-Type: text/plain
On Mon, 2008-08-04 at 17:00 +1000, Simon Blomberg wrote:
> See ?fivenum in the stats package. If you just type
>
> stats::fivenum
>
> you will get the code. The crucial calculations are in the last few
> lines.
That will only give the code to calculate the five number summary, but
Fernando wants to know how the confidence interval is calculated in
boxplot.stats.
To see the code just type boxplot.stats followed by return at the
command line in R.
The relevant line is:
conf <- if (do.conf)
stats[3] + c(-1.58, 1.58) * iqr/sqrt(n)
Which is working on the median (stats[3]). Details of this computation
are in ?boxplot.stats which should have been Fernando's first port of
call. Two works are cited in support of the calculation with full
references in the References section of that help page.
HTH
G
>
> Simon.
>
> On Mon, 2008-08-04 at 16:19 +0930, Fernando Marmolejo Ramos wrote:
> > Dear people
> >
> > I've learnt that by using the "boxplot.stats" command in the
"grDevices" library
> > I can get the 5-number summaries of a boxplot, plus other important
information,
> > like the confidence interval around the median.
> >
> > I'm interested in knowing the actual formula to used in that package
to
> > calculate that confidence interval.
> >
> > Can someone help me with this?
> >
> > Cheers,
> >
> > Fernando
> >
> > ______________________________________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
------------------------------
Message: 47
Date: Mon, 4 Aug 2008 09:58:14 +0100
From: "Keith Jewell" <k.jewell at campden.co.uk>
Subject: Re: [R] Unexpected nls behaviour: Solved
To: r-help at stat.math.ethz.ch
Message-ID: <g76gbc$h01$1 at ger.gmane.org>
Hi Everyone,
I'd omitted the non-optional 'parameters' argument to selfStart. Making
this
change to SSbatch gives the same (successful) result from the two calls
to
nls.
SSbatch<-selfStart(
model=function(Batch, Coeffs)
{
Coeffs[Batch]
}
,initial=function(mCall, data, LHS)
{
# Estimate coefficients as mean of each batch
xy <- sortedXyData(mCall[["Batch"]], LHS, data)
Batch <- data[[as.character(mCall[["Batch"]])]]
# check Batch is successive integers starting at 1
if ((min(xy$x) !=1) | (any(diff(xy$x)!=1))) stop(
"Batch is not a successive integers sequence")
Lval <- list(xy$y)
names(Lval) <- mCall["Coeffs"]
Lval
}
, parameters = c("Coeffs")
)
Sorry for wasting anyones time.
Keith Jewell
-------------------------------------------------
"Keith Jewell" <k.jewell at campden.co.uk> wrote in message news:...
> Hi everyone,
>
> I thought that for a selfStart function, these two should be exactly
> equivalent
>> nls(Aform, DF)
>> nls(Aform, DF, start=getInitial(Aform, DF))
> but in this example that is not the case in R (although it is in
S-plus
> V6.2)
> ------------------------------
> SSbatch<-selfStart(
> model=function(Batch, Coeffs)
> {
> Coeffs[Batch]
>
> }
> ,initial=function(mCall, data, LHS)
> {
> # Estimate coefficients as mean of each batch
> xy <- sortedXyData(mCall[["Batch"]], LHS, data)
> Batch <- data[[as.character(mCall[["Batch"]])]]
> # check Batch is successive integers starting at 1
> if ((min(xy$x) !=1) | (any(diff(xy$x)!=1))) stop(
> "Batch is not a successive integers sequence")
> Lval <- list(xy$y)
> names(Lval) <- mCall["Coeffs"]
> Lval
> }
> )
> DF <- data.frame(A=c(0.9, 1.1, 1.9, 2.0, 2.1, 2.9, 3.0),
> Batch=c(1,1,2,2,2,3,3))
> Aform <- formula(A~SSbatch(Batch,cA))
> nls(Aform, DF, start=getInitial(Aform, DF))
> nls(Aform, DF)
> ------------------------------------
> Don't ask why I'd want such a silly selfStart, that's a long story.
> I guess wherever I would have used nls(Aform, DF)
> I could use nls(Aform, DF, start=getInitial(Aform, DF))
> but that seems clumsy.
>
> Can anyone point out my mistake? Or is this a limitation of nls in R
(I
> hesitate to use the b*g word).
>
> Thanks in advance,
>
> Keith Jewell
> ----------------------------------
>
> I don't think it's relevant but, for completeness:
>
>> sessionInfo()
>
> version 2.7.0 (2008-04-22)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
> Kingdom.1252;LC_MONETARY=English_United
> Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] stats graphics grDevices datasets tcltk utils
methods
> base
>
> other attached packages:
> [1] xlsReadWrite_1.3.2 svSocket_0.9-5 svIO_0.9-5
R2HTML_1.58
> svMisc_0.9-5 svIDE_0.9-5
>
> loaded via a namespace (and not attached):
> [1] tools_2.7.0 VGAM_0.7-7
> --------------------------------
>
>
------------------------------
Message: 48
Date: Mon, 4 Aug 2008 02:40:44 -0700 (PDT)
From: Pierre8rou <pierre8r-nabble at yahoo.fr>
Subject: Re: [R] How to format the output file just the way I want ?
To: r-help at r-project.org
Message-ID: <18808276.post at talk.nabble.com>
Content-Type: text/plain; charset=us-ascii
I think I have found the answer myself.
If you have something better, please write it.
Pierre8r
My R code :
-----------
library(quantmod)
library(xts)
Lines <-
"2008.07.01,02:00,1.5761,1.5766,1.5760,1.5763,65
2008.07.01,02:15,1.5762,1.5765,1.5757,1.5761,95
2008.07.01,02:30,1.5762,1.5765,1.5758,1.5759,58
2008.07.01,02:45,1.5758,1.5758,1.5745,1.5746,91
2008.07.01,03:00,1.5745,1.5751,1.5744,1.5744,117
2008.07.01,03:15,1.5742,1.5742,1.5727,1.5729,100
2008.07.01,03:30,1.5730,1.5736,1.5730,1.5736,61
2008.07.01,03:45,1.5735,1.5740,1.5735,1.5739,55"
quotes <- read.csv(textConnection(Lines), header=FALSE)
x <- as.xts(quotes[,-(1:2)],
as.POSIXct(paste(quotes[,1],quotes[,2]),format='%Y.%m.%d %H:%M'))
colnames(x) <- c('Open','High','Low','Close','Volume')
x$Open <- sprintf("%5.04f", x$Open)
x$High <- sprintf("%5.04f", x$High)
x$Low <- sprintf("%5.04f", x$Low)
x$Close <- sprintf("%5.04f", x$Close)
x$Volume <- sprintf("%.0f", x$Volume)
write.table(x,file = "K:\\OutputFile.csv", quote=FALSE,
col.names=FALSE,
row.names=format(index(x),"%Y.%m.%d,%H:%M"), sep=",")
--
View this message in context:
http://www.nabble.com/How-to-format-the-output-file-just-the-way-I-want-
--tp18790926p18808276.html
Sent from the R help mailing list archive at Nabble.com.
------------------------------
Message: 49
Date: Mon, 4 Aug 2008 11:50:53 +0200
From: Martin Maechler <maechler at stat.math.ethz.ch>
Subject: Re: [R] source a script file straight from a subversion
repository
To: "Steven McKinney" <smckinney at bccrc.ca>
Cc: marc_schwartz at comcast.net, r-help at r-project.org,
murdoch at stats.uwo.ca
Message-ID: <18582.53501.790713.250345 at stat.math.ethz.ch>
Content-Type: text/plain; charset=us-ascii
>>>>> "SM" == Steven McKinney <smckinney at bccrc.ca>
>>>>> on Fri, 1 Aug 2008 18:38:54 -0700 writes:
SM> Thanks to Duncan Murdoch and Marc Schwartz for their
SM> excellent help.
SM> As we don't yet have the apache http: interface to svn
SM> running yet, 'svn export' is the access mechanism.
yes.
Note that svn.r-project.org for the R sources (and some
recommended packages' ones) *only* runs the apache secure http
interface 'https:/.........' for security reasons.
It has the big advantage that there exist no (!) user logins on
that machine at all, and hence no compromised user ssh
information can be misused to break in.
[[elided Yahoo spam]]
Martin Maechler, ETH Zurich
------------------------------
_______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
End of R-help Digest, Vol 66, Issue 4
More information about the R-help
mailing list