Re: [R] Wichmann-Hill Random Number Generator and the Birthday Problem

[Peter Dalgaard]

Prof Brian Ripley wrote:
Your coincidence calculations may be correct for _independent_ draws 
from a discrete distribution on M values, but independence is not 
satisfied.
Yet again, you are trying to do things that any good text on 
simulation would warn you against, and which (in a thread on R-devel) 
you have already been told a good way to do.


A good example are the "good" old linear congruential random generators. 
These will start repeating at the first coincidence, for the pretty 
obvious reason that the next random number is a function only of the 
previous one. So the number of distinct values in N draws is min(N, 
cycle_length). In particular, the number of coincidences is 0 when N is 
less than cycle_length.


--
  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
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

__
R-help@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] nonlinear constrained optimization

[Chua Siang Li]

   Hi.  I need some advises on how to use R to find pi (i is the index) with
   the following objective function and constraint:

   max (sum i)[ f(ai, bi, pi) * g(ci, di, pi) * Di ]

   s.t.  (sum i)[ f(ai, bi, pi) * Di * pi] / (sum i)[ f(ai, bi, pi) * Di ]  <=
   constant

   f and g are diffentiable.
   So, I am thinking of optim with method = "BFGS"? But wonder how to include
   the constraint.
   I also saw constrOptim but it says for linear constraint only.
   So, what will be the suitable function to use then?
   Also, how should I write (sum i) in R?
   Many thanks and your help is much appreciated.
   
   Chua Siang Li
   Consultant - Operations Research
   Acceval Pte Ltd
   Tel: 6297 8740
   Email: [EMAIL PROTECTED]
   Website: www.acceval-intl.com
   This message and any attachments (the "message"...{{dropped:13}}
__
R-help@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.

Re: [R] Open directory within a menu in tcltk

[Peter Dalgaard]

Judith Flores wrote:

Hello,

I am trying to setup a menu to open files and directories. I have the 
following code:




opendir<-function() {
  dirname<<-tclvalue(tkchooseDirectory())
}

openfile<-function() {
  filename<<-tclvalue(tkgetOpenFile())
}

require(tcltk)
t1<-tktoplevel()
topMenu<-tkmenu(t1)
tkconfigure(t1,menu=topMenu)

plotMenu<-tkmenu(topMenu, tearoff=FALSE)

tkadd(plotMenu,"command",label="Open file", command=openfile)
tkadd(plotMenu, "command",label="Select directory",command=opendir)
tkadd(plotMenu,"command",label="Quit", command=function() tkdestroy(t1))
tkadd(topMenu, "cascade", label="Menu",menu=plotMenu) 


tkfocus(t1)

  I can extract a file, but not a directory, I obtain the following error 
message:

Error in function ():
cannot change value of locked binding for 'dirname'.

  I would appreciate very much any help or explanation on this.
  


This is mainly the effect of an unfortunate choice of variable name. 
Beware the semantics of "<<-". You need to ensure that you have control 
of where the assignment goes to.


It is in general unfortunate to make such assignments to the global 
environment, where they can mask and clobber other objects. I'd consider 
a redesign along the lines of


myGui <- function() {
   opendir<-function() {
   dirname<<-tclvalue(tkchooseDirectory())
  }
   dirname <- NA
   tl <- tktoplevel()
  ..
}



As always, thank you,

Judith

__
R-help@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.
  



--
  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
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

__
R-help@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] Writing R Extensions : A new R package for Gini Index decomposition to prupose

[Ndoye Souleymane]

Dear All,

I have developed a programme the anable the decomposition of the Gini index, it 
complets tha valuable work of Achim Zeileis, the author of the ineq package. 
I would like to make it to be part of all R package. How should I proceed.
Must I sent it to the the Core developement team ?
The proogramme is written in R. 

Many thanks for your advice,

Best regards,

Souleymane
_
Retouchez, classez et partagez vos photos gratuitement avec le logiciel Galerie 
de Photos !

__
R-help@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] Exponential smoothing?

[Öhagen Patrik]

Dear List,

I have used all my resources (i.e. "help.search) and I still havn't been able 
to figure out if there is an Exponential Smoothing command in R.

Thank you in advance!

__
R-help@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.

Re: [R] Exponential smoothing?

[Søren Højsgaard]

Try ?HoltWinters and ?filter
Regards
Søren
 


-Oprindelig meddelelse-
Fra: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] På vegne af Öhagen Patrik
Sendt: 19. august 2008 10:00
Til: R hELP
Emne: [R] Exponential smoothing?

Dear List,

I have used all my resources (i.e. "help.search) and I still havn't been able 
to figure out if there is an Exponential Smoothing command in R.

Thank you in advance!

__
R-help@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@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] No object with name: PDE

[Birgitle]

Hello R-User!

When I plot my rpart-object:

 plot(data.rpart); text(data.rpart,cex=0.2, use.n=T)

I get this error message

No object with name:  PDE

Can somebody tell me why and what the meaning of the message is?

Many thanks in advance

B.



-
The art of living is more like wrestling than dancing.
(Marcus Aurelius)
-- 
View this message in context: 
http://www.nabble.com/No-object-with-name%3A--PDE-tp19046286p19046286.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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.

Re: [R] variance components

[Pablo G Goicoechea]

   Dear John:
   Weir, BS 1996 Genetic Data Analysis II . Sinaur Associates, Sunderland, MA,;
   should get you started for methods in population genetics
   (otherreferencecouldbetheArlequin'smanual:
   [1]http://cmpg.unibe.ch/software/arlequin3/)
   However, you are probably interested in variance components for quantitative
   genetics (half-full sibs, mother-progeny correlations, etc.). In such a
   case,you   could   take   a   look   at   Lynch   &   Walsh   book
   [2]http://www.amazon.com/Genetics-Analysis-Quantitative-Traits-Michael/dp/08
   78934812
   Hope this helps
   Pablo
   John Sorkin escribió:

Please forgive the cross posting. I sent this to Medstats and have not received
 any response.

I am looking for a good written explanation of components of variance (a.k.a va
riance components), in particular as they are used in genetic analyses. Does an
yone have a suggestion? I am interested in books, articles and on-line material
. I intent to give the material to some of my students none of whom are statist
icians, they are all physicians or Ph.D.s in life sciences.
Thanks
John

John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)

Confidentiality Statement:
This email message, including any attachments, is for th...{{dropped:6}}

__
[EMAIL PROTECTED] mailing list
[4]https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide [5]http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



   --

   Pablo G Goicoechea

   Bioteknología Saila / Dpto Biotecnología

   NEIKER-Tecnalia

   Apdo 46

   01080 Vitoria-Gasteiz (SPAIN)

   Phone: +34 902 540 546 Fax: +34 902 540 547

   [EMAIL PROTECTED]

References

   1. http://cmpg.unibe.ch/software/arlequin3/
   2. 
http://www.amazon.com/Genetics-Analysis-Quantitative-Traits-Michael/dp/0878934812
   3. mailto:R-help@r-project.org
   4. https://stat.ethz.ch/mailman/listinfo/r-help
   5. http://www.R-project.org/posting-guide.html
   6. mailto:[EMAIL PROTECTED]
__
R-help@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.

Re: [R] Exponential smoothing?

[Hans W. Borchers]

Öhagen Patrik  mpa.se> writes:
> 
> Dear List,
> 
> I have used all my resources (i.e. "help.search) and I still havn't been
> able to figure out if there is an Exponential Smoothing command in R.

A few weeks ago the book "Forecasting with Exponential Smoothing" by 
Hyndman et al. has appeared .

In connection with this book the 'forecasting' bundle is available
from the CRAN package repository, containing functions and data sets
described in the book.

//  Hans Werner Borchers

> Thank you in advance!
>

__
R-help@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] Scaling lattice tiffs

[Dieter Menne]

Dear ListeRs,

Since I normally produce output by Sweave, I never had to care about pixel
resolution. Now someone asked me for 600dpi tiffs of lattice.

I did not find an easy way to simply scale a tiff-image that looks good in
default display. See example below.

Dieter


library(lattice)

p = xyplot(Sepal.Length + Sepal.Width ~ Petal.Length + Petal.Width |
Species, 
   data = iris)

# looks good, but low resolution
trellis.device("tiff",file="a.tif")
print(p)
dev.off()

# do I really have to scale all fonts to get this right?
trellis.device("tiff",file="a1.tif",res=600,width=6000,height=4000)
print(p)
dev.off()

__
R-help@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.

Re: [R] changing plot font for Times new roman

[Prof Brian Ripley]

For X11 the family is "Times" (not "times").  See ?X11Fonts for the 
alternative "serif".


Font values > 5 apply only to the Windows family of devices, and can be 
changed there via the Rdevga file.  But family="serif" works there too 
(although its 'Times New Roman' may not be identical to that on X11, but 
then what you get on X11 depends on the support software installed, in 
particular fontconfig's knowledge of fonts -- on my system it is using 
the URW clone).


On Mon, 18 Aug 2008, milton ruser wrote:


Dear all,

I know that it is a know issue, but I would like to change
the type of font on my plot, and I am not sure the
rigth way. I would like to use Times New Roman font,
but according to the "par()" help, some device allow
we choose an family of fonts. I tryed par(family="times")
without success.

Surfing on R archieve I got an suggestion of
use par(font.lab=6), but when I go to the
par(font.lab) help the highest value there is 5.
How can I have sure that font.lab=6 is the Times
New Roman?

Thanks in advance,

miltinho astronauta
brazil
---

op<-par()

x11(800,500)
par(mfrow=c(1,2))

x<-plot(runif(100),rnorm(100) , main="standard font")

par(font.lab=6)
par(font.axis=6)

x<-plot(runif(100),rnorm(100) , main="font=6")
par<-op

[[alternative HTML version deleted]]

__
R-help@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,  [EMAIL PROTECTED]
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, UKFax:  +44 1865 272595

__
R-help@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.

Re: [R] "nested" getInitial calls; variable scoping problems: Solved??

[Keith Jewell]

Hi All,

I've found a solution to my problems which I think is optimal, but I've been 
wrong before (!!) so comments are welcome.

To avoid scoping problems in the 'inner' getInitial I modify the formula so 
that all variable values are specified explitly in the formula, the inner 
getInitial need not search for any values so there are no scoping problems:
---
ssA <- selfStart(
  model = function(x, Coeff, A)
{
paste(x, Coeff, A)
},
  initial = function(mCall, data, LHS)
{
x <- eval(mCall[["x"]], data, parent.frame())
A <- eval(mCall[["A"]], data, parent.frame())
paste("CoeffA", x, A)
},
  parameters = c("Coeff")
   )
ssB <- selfStart(
  model = function(x, Coeff, A)
{
paste(x, Coeff, A)
},
  initial = function(mCall, data, LHS)
{
A <- eval(mCall[["A"]], data, parent.frame())
Amod <- paste(A, "mod in B")
#---
# Doing a getInitial for (for example)  y' ~ ssA(x, Coeff, Amod)
# to avoid scoping problems, specify values in the formula for LHS and
# all ssA arguments except Coeff; i.e. all except attr(ssA, "pnames")
#--
object <- y ~ ssA(x, Coeff, A)# define a formula object
object[[2]] <- eval(LHS, data, parent.frame())# left hand side 
values
# name RHS arguments by match.call, then assign values  to x and A
object[[3]] <- match.call(get(as.character(object[[3]][[1]])), 
object[[3]])
object[[3]]$x <-  as.character(eval(mCall[["x"]], data, parent.frame())) 
# x
object[[3]]$A <-  Amod  # A
getInitial(object, data) # do the getInitial
 },
  parameters = c("Coeff")
 )
getInitial(y ~ ssA("this", "that", "other"), data= data.frame(x=c("test"), 
y=1))
getInitial(y ~ ssB("this", "that", "other"), data= data.frame(x=c("test"), 
y=1))
getInitial(y ~ ssB(x, "that", "other"), data= data.frame(x=c("test1", 
"test2"), y=1:2))
--

Hope that helps someone,

Keith Jewell

"Keith Jewell" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
> Hi All,
>
> Another nls related problem (for background, I'm migrating a complicated 
> modelling package from S-plus to R).
>
> Below I've reduced this to the minimum necessary to demonstrate my problem 
> (I think); the real situation is more complicated.
>
> Two similar selfStart functions, ssA and ssB.
> The 'initial' function for ssB modifies its arguments a little and then 
> calls getInital for ssA.
> In addition to the "x" and the fitted coefficients ("Coeff"), ssA and ssB 
> have arguments ("A") which are not the same length as the dataframe, so 
> cannot be passed as columns of that dataframe.
>
> The initial function for ssA uses eval(... parent.frame()) to find "A", 
> which is fine when called from .GlobalEnv. But when called from the 
> initial function of ssB it can't find (the modified version of) "A" .
>
> If A is assigned as a parameter of the dataframe then getInitial.formula 
> returns "A" as initial estimates of "Coeff", so that doesn't work.
> I've explored the evaluation frame parent-child structure created by the 
> nested getInitial calls, but it doesn't seem helpful, and I certainly 
> couldn't trust my understanding of it.
> I'm considering making up the matched call and calling the 'initial' 
> function of ssA directly, in the same manner as getInitial.selfStart
>attr(ssA, "initial"))(mCall = mCall, data = data, LHS = LHS)
> or perhaps if I call getInitial thus...
>getInitial(ssA, data, mCall, LHS = NULL, ...)
> it will call getInitial.selfStart without calling getInitial.formula.
> In that case perhaps eval(... parent.frame())  will find the argument?
> Or, at least I could attach the argument to the dataframe as a parameter.
>
> But this all seems a but clumsy, and I feel there must be a better way. 
> Any comments and/or advice will be welcome.
>
> Thanks in advance;
>
> Keith Jewell
> ---
> code showing the problem:
> ssA <- selfStart(
>  model = function(x, Coeff, A)
>{
>paste(x, Coeff, A)
>},
>  initial = function(mCall, data, LHS)
>{
>x <- eval(mCall[["x"]], data, parent.frame())
>A <- eval(mCall[["A"]], data, parent.frame())
>paste("CoeffA", x, A)
>},
>  parameters = c("Coeff")
>   )
> ssB <- selfStart(
>  model = function(x, Coeff, A)
>{
>paste(x, Coeff, A)
>},
>  initial = function(mCall, data, LHS)
>{
>x <- eval(mCall[["x"]], data, parent.frame())
>A <- eval(mCall[["A"]], data, parent.frame())
>Amod <- paste(A, "mod in B")
>getInitial(y ~ ssA(x, Coeff, Amod), data)
>},
>  parameters = c("Coeff")
> )
> getInitial(y ~ ssA("this", "that", "other"), data= 
> data.frame(x=c("test")))
> getInitial(y ~ ssB("this", "that", "other"), data= 
> data.frame(x=c("test")))
> ---
>> sessionInfo()
> R version 2.7.1 Patched (

Re: [R] Fwd: Parsing XML or KML into CSV /Using R for geo coding , OR problem

[Roger Bivand]

Ajay ohri  gmail.com> writes:

> 
> Hi,
> 
> I have a data file in a KML format which is Google Earth's format for
> geographic data. I need to import it into a csv file . How can I do that.
> KML format is just like XML format .example below

The easiest way is most likely to use the rgdal KML driver, now also reading KML
files. Install the rgdal package, and use the readOGR() function. Then coerce 
the
imported object to a data frame, and export to csv with write.csv().

Roger

> 
> Which R module with deal with an OR problem (like transportation problem
> using geo coded data as in the example below)
> 
> Regards,
> Ajay
> 
> www.decisionstats.com
> -

__
R-help@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] gam.check in gam (mgcv)

[Van Wyk, Jaap]

Hallo
I need some help with the output provided by gam.check after a gam fit (using 
the package mgcv).
To give a brief description of my data, I have
claims: a vector of values, which include NA's and one large negative value - 
otherwise all positive (55 values in total that are not NA).
origin: a factor with 10 levels
j : taking the values 1, 2, , 10
I have fitted a gam, with:
 
> library(mgcv)
> modelg <- gam(claims ~ origin + s(log(col)), family = quasipoisson(),
 +  subset=!is.na(claims))
 
(I have obtained a work-around of the quasipoisson family to deal with the 
negative value.)
 
After this I tried the following:
 
> gam.check(modelg)
 
This gives a plot in the usual way (that looks fine, except for the "outlier" 
which was a negative value), but returns the following message:
 
fit method: deviance based outer iter. - newton, exact hessian.
step failed after 1 iteration.
gradient range [9.136047,9.136047] (score 994.1401 & scale 733.0711).
Hessian positive definite, eigenvalue range [42.98903,42.98903].

What does this mean? Is there a problem?
 
Any help is much appreciated.
 
Thanks.
Jacob L van Wyk
Dept of Statistics
University of Johannesburg
P O Box 524
Auckland Park 2006
South Africa
Tel: +27-11-559-3080
Fax: +27-11-559-2832
Cell: +27-82-859-2031

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Converting monthly data to quarterly dataMonday, August 18, 2008 11:38 AM

[Denise Xifara]

Dear Gavin,

This is really great, thank you!  I created some long loops to get rid of
extra months at the beginning and the end of my data but your code is great
for putting it then together quarterly.
thanks again,
Denise

On Mon, 2008-08-18 at 14:31 +0100, Denise Xifara wrote:
> Thank you very much Stephen, but how will aggregate deal with months that
> fall outside annual quarters? eg, one extra month at the end of the
dataset?

[Without your data I'm kind of guessing at the exact format and problem,
but the example below shows one way to deal with quarters that span
years so should be adaptable to your problem. You also have a different
idea of quarters to my ecological one...]

## first some dummy data
set.seed(12345)
dat <- data.frame(price = cumsum(rnorm(450)),
  date = seq(from = as.Date("2008-01-01"),
 length = 450, by = "days"))

## get the months of observations
dat$month <- factor(format(dat$date, format = "%b"), levels = month.abb)

## and then format this for the quarters
dat$quarter <- character(length = NROW(dat))
## I'm sure these 4 steps can be simplified but just
## how escapes me at the moment
dat$quarter[dat$month %in% month.abb[c(12,1,2)]] <- "Winter"
dat$quarter[dat$month %in% month.abb[c(3:5)]] <- "Spring"
dat$quarter[dat$month %in% month.abb[c(6:8)]] <- "Summer"
dat$quarter[dat$month %in% month.abb[c(9:11)]] <- "Autumn"
dat$quarter <- factor(dat$quarter,
  levels = c("Spring","Summer","Autumn","Winter"))

## look at the fruits of our labour
head(dat)

## create period
runs <- rle(as.numeric(dat$quarter))$lengths
dat$period <- factor(rep(seq_along(runs), times = runs))

## aggregate
with(dat, aggregate(price, list(quarter = quarter, period = period),
FUN = mean))

I use the rle() function (run length encoding) calculate the number of
observations where the 'quarter' remains the same:

> rle(as.numeric(dat$quarter))
Run Length Encoding
  lengths: int [1:6] 60 92 92 91 90 25
  values : num [1:6] 4 1 2 3 4 1

The 'values' here are the numeric representation of the quarter factor.
The most interesting for us is the second 4 - this is the winter 2008/9.
I use the lengths to replicate a period number (1,2,...,n) the correct
number of times. Now we have the period correctly calculated, we just
aggregate by quarter and period to give the averages you want.

If you are working on months 1-3 as quarter 1, 4-6 as quarter 2 etc,
then it is much easier, just aggregate by quarter and year:

## copy the data above
dat2 <- dat
## change meaning of quarter
dat2$quarter <- character(length = NROW(dat2))
dat2$quarter[dat2$month %in% month.abb[c(1:3)]] <- "Q1"
dat2$quarter[dat2$month %in% month.abb[c(4:6)]] <- "Q2"
dat2$quarter[dat2$month %in% month.abb[c(7:9)]] <- "Q3"
dat2$quarter[dat2$month %in% month.abb[c(10:12)]] <- "Q4"
dat2$quarter <- factor(dat2$quarter, levels = c("Q1","Q2","Q3","Q4"))
## year variable
dat2$year <- factor(format(dat2$date, format = "%Y"))

## drop the first 40 days to simulate a late starting record
## and aggregate
with(dat2[-(1:40), ], aggregate(price, list(quarter = quarter, year =
year), FUN = mean))

Which gives:
  quarter yearx
1  Q1 2008 13.58644
2  Q2 2008 24.16523
3  Q3 2008 28.56004
4  Q4 2008 32.60900
5  Q1 2009 44.86594

Do these examples help solve your problem?

G
>
> 2008/8/18 stephen sefick <[EMAIL 
> PROTECTED]
>
>
> > ?aggregate
> > may do what you want
> >
> > On Mon, Aug 18, 2008 at 8:19 AM, Denise Xifara
> > <[EMAIL PROTECTED]>
wrote:
> > > Dear R users,
> > >
> > > I have a dataframe where column is has countries, column 2 is dates
> > > (monthly) for each countrly, the next 10 columns are my factors where
I
> > have
> > > measurements for each country and  for each date.  I have attached a
> > sample
> > > of the data in csv format with the data for 3 countries.
> > >
> > > I would like to convert my monthly data into quarterly data, finding
the
> > > mean over 3 month periods for factors a-i, and the sum for factor j.
My
> > > problem is that not all countries have starting date at the beginning
of
> > a
> > > quarter for a particular year, ie some countries start in May or
> > September,
> > > and also some countries have one extra month, some have two extra
months
> > so
> > > there's no way of deleting some rows with a simple command (I want to
get
> > > rid of all extra data that does not fall into the quarters for each
> > > country), since the amount of data to get rid of for each country
varies.
> > >
> > > I tried for example:
> > > i=1
> > > denise<-data[((data$country)==unique(data$country[i]),]
> > > denise[,2]<- as.Date(denise$date, "%Y-%m-%d")
> > > denise2<-denise[order(denise[,2],decreasing=FALSE),]
> > > len<-length(denise[,1])
> > > limit<-floor(len/3)+1
> > > splitter<-rep(1:limit,each=3)
> > > spl.dat<-split(denise2,sp

Re: [R] Jitter in xYplot?

[Jim Lemon]

On Mon, 2008-08-18 at 19:25 -0400, John Poulsen wrote:
> Hello,
> 
> The simple example below plots three species by year and site.  However, 
>   I would like the points to occur in groups horizontally, rather than 
> on top of each other.
> 
> I thought jitter might work, but haven't been successful in getting it 
> to do what I want.
> 
Hi John,
I don't know whether this is what you want, but the cluster.overplot
function (plotrix) turns coincident points into groups of up to nine
points. You may be able to adapt the code for use in xYplot.

Jim

__
R-help@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] atomic value

[Laura Cordero Llana]

Hi all,
 
I am trying to run a code using the function "tapply", and after test it and 
see that the value is atomic, when I run the script it gives an error:
 
Error in sort(unique.default(x), na.last = TRUE) : 'x' must be atomic
for (i in 1:50){
year.name <- paste("year",i,sep="")
year.tmp <- rep(x=periodtime[i], each=length(period[i]))
assign(year.name,year.tmp)
 
grab.p <- get(paste("period",i,sep=""))
daf <- data.frame(year.name, grab.p, value=rnorm(length(year.name)))
repl <- tapply(daf$value, daf$period, mean, na.rm=T)
not <- which(is.na(daf$value))daf$value[not] <- repl[daf$period[not]]
data.frame(xtime, daf$value[not])
}
 
Where is the error if x is atomic?
 
 
All the best,
 
Laura
_


[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Converting monthly data to quarterly data

[Gabor Grothendieck]

Another possibility is to use the yearmon and yearqtr classes in the
zoo package.  Assuming quarters end in Feb, May, Aug and Nov
we have (using dat from the prior post):

> library(zoo)
> z <- zoo(dat$price, dat$date)
> aggregate(z, as.yearqtr(as.yearmon(time(z)) + 1/12), mean)
 2008 Q1  2008 Q2  2008 Q3  2008 Q4  2009 Q1  2009 Q2
 5.19913 21.86151 26.06295 29.68677 41.92547 48.10799

or to use calendar quarters its just:

aggregate(z, as.yearqtr, mean)

In the zoo package see ?zoo, ?aggregate.zoo, ?yearqtr, ?yearmon
and the three zoo vignettes.

On Mon, Aug 18, 2008 at 2:38 PM, Gavin <[EMAIL PROTECTED]> wrote:
> On Mon, 2008-08-18 at 14:31 +0100, Denise Xifara wrote:
>> Thank you very much Stephen, but how will aggregate deal with months that
>> fall outside annual quarters? eg, one extra month at the end of the dataset?
>
> [Without your data I'm kind of guessing at the exact format and problem,
> but the example below shows one way to deal with quarters that span
> years so should be adaptable to your problem. You also have a different
> idea of quarters to my ecological one...]
>
> ## first some dummy data
> set.seed(12345)
> dat <- data.frame(price = cumsum(rnorm(450)),
>  date = seq(from = as.Date("2008-01-01"),
> length = 450, by = "days"))
>
> ## get the months of observations
> dat$month <- factor(format(dat$date, format = "%b"), levels = month.abb)
>
> ## and then format this for the quarters
> dat$quarter <- character(length = NROW(dat))
> ## I'm sure these 4 steps can be simplified but just
> ## how escapes me at the moment
> dat$quarter[dat$month %in% month.abb[c(12,1,2)]] <- "Winter"
> dat$quarter[dat$month %in% month.abb[c(3:5)]] <- "Spring"
> dat$quarter[dat$month %in% month.abb[c(6:8)]] <- "Summer"
> dat$quarter[dat$month %in% month.abb[c(9:11)]] <- "Autumn"
> dat$quarter <- factor(dat$quarter,
>  levels = c("Spring","Summer","Autumn","Winter"))
>
> ## look at the fruits of our labour
> head(dat)
>
> ## create period
> runs <- rle(as.numeric(dat$quarter))$lengths
> dat$period <- factor(rep(seq_along(runs), times = runs))
>
> ## aggregate
> with(dat, aggregate(price, list(quarter = quarter, period = period),
>FUN = mean))
>
> I use the rle() function (run length encoding) calculate the number of
> observations where the 'quarter' remains the same:
>
>> rle(as.numeric(dat$quarter))
> Run Length Encoding
>  lengths: int [1:6] 60 92 92 91 90 25
>  values : num [1:6] 4 1 2 3 4 1
>
> The 'values' here are the numeric representation of the quarter factor.
> The most interesting for us is the second 4 - this is the winter 2008/9.
> I use the lengths to replicate a period number (1,2,...,n) the correct
> number of times. Now we have the period correctly calculated, we just
> aggregate by quarter and period to give the averages you want.
>
> If you are working on months 1-3 as quarter 1, 4-6 as quarter 2 etc,
> then it is much easier, just aggregate by quarter and year:
>
> ## copy the data above
> dat2 <- dat
> ## change meaning of quarter
> dat2$quarter <- character(length = NROW(dat2))
> dat2$quarter[dat2$month %in% month.abb[c(1:3)]] <- "Q1"
> dat2$quarter[dat2$month %in% month.abb[c(4:6)]] <- "Q2"
> dat2$quarter[dat2$month %in% month.abb[c(7:9)]] <- "Q3"
> dat2$quarter[dat2$month %in% month.abb[c(10:12)]] <- "Q4"
> dat2$quarter <- factor(dat2$quarter, levels = c("Q1","Q2","Q3","Q4"))
> ## year variable
> dat2$year <- factor(format(dat2$date, format = "%Y"))
>
> ## drop the first 40 days to simulate a late starting record
> ## and aggregate
> with(dat2[-(1:40), ], aggregate(price, list(quarter = quarter, year =
> year), FUN = mean))
>
> Which gives:
>  quarter yearx
> 1  Q1 2008 13.58644
> 2  Q2 2008 24.16523
> 3  Q3 2008 28.56004
> 4  Q4 2008 32.60900
> 5  Q1 2009 44.86594
>
> Do these examples help solve your problem?
>
> G
>>
>> 2008/8/18 stephen sefick <[EMAIL PROTECTED]>
>>
>> > ?aggregate
>> > may do what you want
>> >
>> > On Mon, Aug 18, 2008 at 8:19 AM, Denise Xifara
>> > <[EMAIL PROTECTED]> wrote:
>> > > Dear R users,
>> > >
>> > > I have a dataframe where column is has countries, column 2 is dates
>> > > (monthly) for each countrly, the next 10 columns are my factors where I
>> > have
>> > > measurements for each country and  for each date.  I have attached a
>> > sample
>> > > of the data in csv format with the data for 3 countries.
>> > >
>> > > I would like to convert my monthly data into quarterly data, finding the
>> > > mean over 3 month periods for factors a-i, and the sum for factor j.  My
>> > > problem is that not all countries have starting date at the beginning of
>> > a
>> > > quarter for a particular year, ie some countries start in May or
>> > September,
>> > > and also some countries have one extra month, some have two extra months
>> > so
>> > > there's no way of deleting some rows with a simple command (I want to get
>> > > rid of all extra data that does 

Re: [R] Converting monthly data to quarterly data

[Denise Xifara]

Wonderful!
Thank you very much Gabor! I'll check it out now.

Kind regards,
Denise

2008/8/19 Gabor Grothendieck <[EMAIL PROTECTED]>

> Another possibility is to use the yearmon and yearqtr classes in the
> zoo package.  Assuming quarters end in Feb, May, Aug and Nov
> we have (using dat from the prior post):
>
> > library(zoo)
> > z <- zoo(dat$price, dat$date)
> > aggregate(z, as.yearqtr(as.yearmon(time(z)) + 1/12), mean)
>  2008 Q1  2008 Q2  2008 Q3  2008 Q4  2009 Q1  2009 Q2
>  5.19913 21.86151 26.06295 29.68677 41.92547 48.10799
>
> or to use calendar quarters its just:
>
> aggregate(z, as.yearqtr, mean)
>
> In the zoo package see ?zoo, ?aggregate.zoo, ?yearqtr, ?yearmon
> and the three zoo vignettes.
>
> On Mon, Aug 18, 2008 at 2:38 PM, Gavin <[EMAIL PROTECTED]> wrote:
> > On Mon, 2008-08-18 at 14:31 +0100, Denise Xifara wrote:
> >> Thank you very much Stephen, but how will aggregate deal with months
> that
> >> fall outside annual quarters? eg, one extra month at the end of the
> dataset?
> >
> > [Without your data I'm kind of guessing at the exact format and problem,
> > but the example below shows one way to deal with quarters that span
> > years so should be adaptable to your problem. You also have a different
> > idea of quarters to my ecological one...]
> >
> > ## first some dummy data
> > set.seed(12345)
> > dat <- data.frame(price = cumsum(rnorm(450)),
> >  date = seq(from = as.Date("2008-01-01"),
> > length = 450, by = "days"))
> >
> > ## get the months of observations
> > dat$month <- factor(format(dat$date, format = "%b"), levels = month.abb)
> >
> > ## and then format this for the quarters
> > dat$quarter <- character(length = NROW(dat))
> > ## I'm sure these 4 steps can be simplified but just
> > ## how escapes me at the moment
> > dat$quarter[dat$month %in% month.abb[c(12,1,2)]] <- "Winter"
> > dat$quarter[dat$month %in% month.abb[c(3:5)]] <- "Spring"
> > dat$quarter[dat$month %in% month.abb[c(6:8)]] <- "Summer"
> > dat$quarter[dat$month %in% month.abb[c(9:11)]] <- "Autumn"
> > dat$quarter <- factor(dat$quarter,
> >  levels = c("Spring","Summer","Autumn","Winter"))
> >
> > ## look at the fruits of our labour
> > head(dat)
> >
> > ## create period
> > runs <- rle(as.numeric(dat$quarter))$lengths
> > dat$period <- factor(rep(seq_along(runs), times = runs))
> >
> > ## aggregate
> > with(dat, aggregate(price, list(quarter = quarter, period = period),
> >FUN = mean))
> >
> > I use the rle() function (run length encoding) calculate the number of
> > observations where the 'quarter' remains the same:
> >
> >> rle(as.numeric(dat$quarter))
> > Run Length Encoding
> >  lengths: int [1:6] 60 92 92 91 90 25
> >  values : num [1:6] 4 1 2 3 4 1
> >
> > The 'values' here are the numeric representation of the quarter factor.
> > The most interesting for us is the second 4 - this is the winter 2008/9.
> > I use the lengths to replicate a period number (1,2,...,n) the correct
> > number of times. Now we have the period correctly calculated, we just
> > aggregate by quarter and period to give the averages you want.
> >
> > If you are working on months 1-3 as quarter 1, 4-6 as quarter 2 etc,
> > then it is much easier, just aggregate by quarter and year:
> >
> > ## copy the data above
> > dat2 <- dat
> > ## change meaning of quarter
> > dat2$quarter <- character(length = NROW(dat2))
> > dat2$quarter[dat2$month %in% month.abb[c(1:3)]] <- "Q1"
> > dat2$quarter[dat2$month %in% month.abb[c(4:6)]] <- "Q2"
> > dat2$quarter[dat2$month %in% month.abb[c(7:9)]] <- "Q3"
> > dat2$quarter[dat2$month %in% month.abb[c(10:12)]] <- "Q4"
> > dat2$quarter <- factor(dat2$quarter, levels = c("Q1","Q2","Q3","Q4"))
> > ## year variable
> > dat2$year <- factor(format(dat2$date, format = "%Y"))
> >
> > ## drop the first 40 days to simulate a late starting record
> > ## and aggregate
> > with(dat2[-(1:40), ], aggregate(price, list(quarter = quarter, year =
> > year), FUN = mean))
> >
> > Which gives:
> >  quarter yearx
> > 1  Q1 2008 13.58644
> > 2  Q2 2008 24.16523
> > 3  Q3 2008 28.56004
> > 4  Q4 2008 32.60900
> > 5  Q1 2009 44.86594
> >
> > Do these examples help solve your problem?
> >
> > G
> >>
> >> 2008/8/18 stephen sefick <[EMAIL PROTECTED]>
> >>
> >> > ?aggregate
> >> > may do what you want
> >> >
> >> > On Mon, Aug 18, 2008 at 8:19 AM, Denise Xifara
> >> > <[EMAIL PROTECTED]> wrote:
> >> > > Dear R users,
> >> > >
> >> > > I have a dataframe where column is has countries, column 2 is dates
> >> > > (monthly) for each countrly, the next 10 columns are my factors
> where I
> >> > have
> >> > > measurements for each country and  for each date.  I have attached a
> >> > sample
> >> > > of the data in csv format with the data for 3 countries.
> >> > >
> >> > > I would like to convert my monthly data into quarterly data, finding
> the
> >> > > mean over 3 month periods for factors a-i, and the sum for factor j.
>  My

[R] how can i get hessian matrix at "constrOptim"

[Kanak Choudhury]

Hi,
i have made a code for optimizing a function using "constrOptim". i need
hessain matrix of the parameters. how could i get hessain matrix when i will
use "constrOptim"? May i get get any help from anyone?

thank you in advance.

Kanak Choudhury.

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Format with n significant figures / digits

[wannymahoots]

In case anyone has the same problem and wants a solution, here is mine
(albeit with some horrible hacks):

"my.signif" <- function(vec, digits=6) {
# get correct number of significant figures
vec = signif(vec, digits)

# now get correct number of trailing zeros by calling sprintf and
specifying
# the correct number of decimal places seperately for each
number
om = floor(log10(abs(vec)))   # order of magnitude
dp = digits-om-1  # no. decimal points
dp[which(dp<0)] = 0   # use zero when dp is negative
for(i in 1:length(vec)) {
if(is.na(vec[i])) { next }
vec[i] = sprintf(paste("%.",dp[i],"f", sep=""), vec[i])
}
return(vec)
}


On Aug 14, 3:23 pm, [EMAIL PROTECTED] wrote:
> I am not sure that format solves the problem (although again I may  
> well be missing something)
>
> # trailing zeros on large numbers
>  >format(vec, digits=4, scientific=F)
> [1] "      0.8000" "    123.4567" "      0.1235" "      7.6543"  
> "7654321."
>
> # misses trailing zeros
>  > format(vec[1], digits=4, scientific=F)
> [1] "0.8"
>
> I guess I could resort to writing a routine that called format with  
> different arguments depending on the magnitude of each number in the  
> vector, but I was hoping for a cleaner solution.  Any thoughts?
>
> On 14 Aug 2008, at 14:57, Prof Brian Ripley wrote:
>
>
>
> > I think you missed the function called 'format'.  R's internal print  
> > routines (used by format) calculate the format passed to (C level)  
> > sprintf based on the input, including 'digits'.  Just make sure you  
> > pass one number at a time to format() if you don't want a common  
> > layout for all the numbers.
>
> > On Thu, 14 Aug 2008, [EMAIL PROTECTED] wrote:
>
> >> Hi everyone,
>
> >> I can't figure out how to format numbers to have a certain number  
> >> of significant figures (as opposed to decimal places) without using  
> >> scientific notation and including trailing zeros if necessary.
>
> >> e.g. I would like to achieve the following:
> >> 0.81  --->  0.8000
> >> 123.4567    --->    123.4
> >> 0.1234567  --->    0.1234
> >> 7.654321    --->    7.654
> >> 7654321     --->    7654000
>
> >> It seems like it should be simple, but I can't seem to do this  
> >> using sprintf. Specifically, I get the following outputs:
>
> >>> vec = c(0.81, 123.4567, 0.1234567, 7.654321, 7654321)
>
> >>> #  this specifies precision after decimal point, rather than  
> >>> significant figures
> >>> sprintf("%.4f", vec)
> >> [1] "0.8000"       "123.4567"     "0.1235"       "7.6543"  
> >> "7654321."
>
> >>> #  uses scientific notation, but significant figures are correct
> >>> sprintf("%.3e", vec)
> >> [1] "8.000e-01" "1.235e+02" "1.235e-01" "7.654e+00" "7.654e+06"
>
> >>> #  correct number of significant figures, and without scientific  
> >>> notation, but has trailing zeros on large numbers
> >>> sprintf("%.4f", sprintf("%.3e", vec))
> >> [1] "0.8000"       "123.5000"     "0.1235"       "7.6540"  
> >> "7654000."
>
> >> Am I missing something obvious and can someone help me?
>
> >> Many thanks
>
> >> __
> >> [EMAIL PROTECTED] mailing list
> >>https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
>
> > --
> > Brian D. Ripley,                 [EMAIL PROTECTED]
> > 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
>
> __
> [EMAIL PROTECTED] mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@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.

Re: [R] R vs Stata on generalized linear mixed models: glmer and xtmelogit

[Douglas Bates]

On Mon, Aug 18, 2008 at 7:55 PM,  <[EMAIL PROTECTED]> wrote:
> Hello,
> I have compared the potentials of R and Stata about GLMM, analysing the 
> dataset 'ohio' in the package 'faraway' (the same dataset is analysed with 
> GEE in the book 'Extending the linear model with R' by Julian Faraway).
> Basically, I've tried the 2 commands 'glmmPQL' and 'glmer' of R and the 
> command 'xtmelogit' of Stata. If I'm not wrong, 'glmer' uses the Laplacian 
> approximation as default, corresponding to adaptive Gauss-Hermite 
> approximation with only 1 point, while 'xtmelogit' uses 7 points. In order to 
> compare them, I tried also to change the corresponding parameters.
>
> This is the code for R:
>
> rm(list=ls())
> library(faraway); library(lme4); library(MASS)
> data <- ohio
> pql  <- glmmPQL(resp~smoke+factor(age), random=~1|id, family=binomial,data)
> summary(pql)$tTable["smoke",1:2]
> lap <- glmer(resp~smoke+factor(age)+(1|id), family=binomial,data)
> attributes(summary(lap))$coefs["smoke",1:2]
> agq7 <- glmer(resp~smoke+factor(age)+(1|id),nAGQ=7,family=binomial,data)
> write.csv(data,file="data.csv")
>
> This is the code for Stata:
>
> clear
> insheet using data.csv
> xi: xtmelogit resp smoke i.age, || id:, covariance(independent) laplace
> xi: xtmelogit resp smoke i.age, || id:, covariance(independent)
>
> Results:
> - Both the point estimate and the standard error for the fixed effect, and 
> the standard deviation for random effect of 'glmmPQL' are lower than 'glmer'
> - 'glmer' and 'xtmelogit' with Laplacian approximation give very similar 
> results. 'xtmelogit' with 7 points gives similar point estimates for fixed 
> effects, but a different (lower) estimate for the standard deviation of the 
> random effect (as expected)
> - glmer doesn't work with the parameter 'nAGO' (number of points) set to 7, 
> returning 'Code not yet written'
>
> My questions:
> 1) Is the difference between 'glmmPQL' and 'glmer' expected? Which is more 
> reliable?

Yes, the difference is to be expected.  The results from glmer should
better approximate the maximum likelihood estimates.  The Laplace
approximation is a direct approximation to the log-likelihood of the
model at the given parameter estimates.  The PQL method is an indirect
method of determining parameter estimates.  Also, glmmPQL is based on
the lme function from the nlme package and the model representation
and optimization methods in the lme4 package are considerably more
advanced than those in the nlme package.

> 2) Is there a way to set the parameter 'nAGO' to 7 or perform the same 
> analysis in another way?

You have to wait until Bin Dai incorporates the AGQ code from his
Google Summer of Code (SoC) project into the release version of the
lme4 package.  Officially the SoC projects ended yesterday so that
should occur soon.

> 3) Is there a tutorial on the use of lme4, especially to handle the results 
> (summary, coef, etc.)?

The documentation is still quite scattered.  I am better at changing
code than I am at updating the documentation to reflect the current
state of the code.

I think the best current documentation by me on the use of lme4 and
some of the background on mixed models are the slides for my
presentations in a recent workshop at the University of Potsdam.  They
are available at http://www.stat.wisc.edu/~bates/PotsdamGLMM/ as files
LMMD.pdf and GLMMD.pdfThe Sweave sources are also available there
as are the scripts extracted from the Sweave sources (in the scripts
directory).



>
> Thank you for your time
>
> Antonio Gasparrini
> Public and Environmental Health Research Unit (PEHRU)
> London School of Hygiene & Tropical Medicine
> Keppel Street, London WC1E 7HT, UK
> Office: 0044 (0)20 79272406 - Mobile: 0044 (0)79 64925523
> http://www.lshtm.ac.uk/people/gasparrini.antonio ( 
> http://www.lshtm.ac.uk/pehru/ )
>
> __
> R-help@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@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] R-code: for those who like a challenge; HELP

[Jan Akko Eleveld]

I'll try to be more clear, but will have to give much more info. Here we go

As part of my Masters in Public Health I writing a dissertation reviewing the 
mortality trend of internally displaced persons in camp settings in Sub Saharan 
Africa. I have 50 surveys with
- crude mortality rate (cmr), CMR lower confidence interval(cmrlci), CMR upper 
confidence interval(cmruci),
- recall period minimum and recall period maximum (recallmin, recallmax)
- time of displacement mimimum and time of displacement maximum (todmin, todmax)
- population data at time of survey and proportion of population present in 
camp over time of displacement (data[9] - data[170])

Because of the uncertainty in these variables I am using the Monte Carlo method 
with e.g. 10,000 iterations and am writing a function that can run through each 
of the 50 studies 10,000 times. The aim is to analyse the trend in IDP 
mortality in relation to time of displacement. I'm aiming at a graph with 
x=time, y=CMR and an attached weight for each coordinate depending on the size 
of the population it represents at randomly choosen time of displacement.

I have the following code that gives me problems and I hope one of you is able 
to help me find what I do wrong.

--

### 1 setting iterations
iterations<-1

### 2 creating results array
results<-array(0,dim=c(3,2500))

### 3 reading in survey dataset
data<-read.csv("DATA2.csv",header=TRUE)

### 4 FUNCTION THAT GENERATES x y AND weight FOR ONE SURVEY
f1<-function(data){

### 5 Reading in important variables from data
cmr<-data[[1]]
cmrlci<-data[[2]]
cmruci<-data[[3]]
recallmin<-data[[4]]
recallmax<-data[[5]]
todmin<-data[[6]]
todmax<-data[[7]]
pop<-data[[8]]

### 6 Generating random CMR or y coordinate (from a triangular distribution)
y<-0
rand<-runif(n=1,min=0,max=1)
prob<-((cmr-cmrlci)^2)/((cmruci-cmrlci)*(cmr-cmrlci))
if (rand<=prob) y<-sqrt(rand*(cmruci-cmrlci)*(cmr-cmrlci)) + cmrlci
if (rand>prob) y<-cmruci-sqrt((1-rand)*(cmruci-cmrlci)*(cmruci-cmr))

### 7 rounding y to 1 sign digit and then multiplying by ten so as to get a 
whole number per 100,000 per day (easier for array)
y<-round(y,digits=1)*10

### 8 Generating random time or x coordinate

### 9 First generating random time within recall period
randomrecall<-round(runif(n=1,min=recallmin,max=recallmax))

### 10 Second generating random time of displacement (this
todrange<-todmax-todmin
proprand<-runif(n=1,min=0,max=1)
randomtod<-0
for(i in 1:todrange){
if (proprand>data[[9+i]] & proprand<=data[[10+i]]) randomtod=i
}

### 11 Thus
x<-randomtod-randomrecall

### 12 Generating population weight
weight<-round(pop*data[[9+randomtod]])/10

### 13 Function returns a vector onesurvey consisting of three values
onesurvey<-c(x,y,weight)
invisible(onesurvey)
}

### 14 Applying f1 to each row of the survey dataset, for the set number of 
iterations
replicate(iterations,{

### 15 Applying f1 to each row of the survey dataset, one iteration
### this should return a list allsurveys with three cols and as many rows as 
there are surveys
allsurveys<-apply(data,1,FUN=f1)

### 16 safer to convert list into a dataframe
allsurveys<-as.data.frame(allsurveys)

### 17 Updating "results" array with random values from each survey, after one 
iteration
### loop that picks up each row of allsurveys and updates the array with values 
from that row
for(i in 1:length(allsurveys[[1]])){
x<-allsurveys[[1]][i]
y<-allsurveys[[2]][i]
weight<-allsurveys[[3]][i]
results[x,y]<-results[x,y]+weight
}
results<<-results
}
)

--
If I go through the code step by step from only 1-13, I run into the following 
problems with R warnings and errors.

> y<-0
> rand<-runif(n=1,min=0,max=1)
> prob<-((cmr-cmrlci)^2)/((cmruci-cmrlci)*(cmr-cmrlci))
> if (rand<=prob) y<-sqrt(rand*(cmruci-cmrlci)*(cmr-cmrlci)) + cmrlci
Warning message:
In if (rand <= prob) y <- sqrt(rand * (cmruci - cmrlci) * (cmr -  :
  the condition has length> 1 and only the first element will be used

> if (rand>prob) y<-cmruci-sqrt((1-rand)*(cmruci-cmrlci)*(cmruci-cmr))
Warning message:
In if (rand> prob) y <- cmruci - sqrt((1 - rand) * (cmruci - cmrlci) *  :
  the condition has length> 1 and only the first element will be used

y<-round(y,digits=1)*10
> y
 [1] 32 16 13 13 44 30  9  7  8  4  5  9 23  6  5  6  8  7 12 15  8 20 39  9  9 
12  9  9  7  4 58 13  8 30  8 12  8 13 30  9  7  7 24  8  9 12 13 19 19 12

? What is the effect of the warning message? What do they mean and how can I 
prevent them? It does not seem to influence the calculation
--
> randomrecall<-round(runif(n=1,min=recallmin,max=recallmax))
> randomrecall
[1] 2

> tod

Re: [R] Writing R Extensions : A new R package for Gini Index decomposition to prupose

[stephen sefick]

See if there is interest.  If there is not make your own package or
see if someone else would like to include it into a package that is
complementary.

Stephen Sefick

On Tue, Aug 19, 2008 at 3:46 AM, Ndoye Souleymane <[EMAIL PROTECTED]> wrote:
>
> Dear All,
>
> I have developed a programme the anable the decomposition of the Gini index, 
> it complets tha valuable work of Achim Zeileis, the author of the ineq 
> package.
> I would like to make it to be part of all R package. How should I proceed.
> Must I sent it to the the Core developement team ?
> The proogramme is written in R.
>
> Many thanks for your advice,
>
> Best regards,
>
> Souleymane
> _
> Retouchez, classez et partagez vos photos gratuitement avec le logiciel 
> Galerie de Photos !
>
> __
> R-help@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.
>



-- 
Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods. We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@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] [R-pkgs] New packages bpca: biplot (2d and 3d) and diagnostic tools

[Jose Claudio Faria]

Dear R Community,

We am pleased to announce the release of the 'bcpa' package.  The
source code and binaries is now on CRAN
(http://cran.r-project.org/web/packages/bpca/index.html).

The bpca package implements biplot (2d and 3d) and diagnostic tools of
the quality of the reduction.
It depends R (>= 2.6.0), scatterplot3d, rgl and MASS packages.

Regards,
-- 
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\
Jose Claudio Faria
Estatistica - prof. Titular
UESC/DCET/Brasil
[EMAIL PROTECTED]
[EMAIL PROTECTED]
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\

___
R-packages mailing list
[EMAIL PROTECTED]
https://stat.ethz.ch/mailman/listinfo/r-packages

__
R-help@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.

Re: [R] ARMA(0,2) & GARCH(1,1) - code & hessian

[Yohan Chalabi]

Hi,

As far as I can tell, your code looks very similar to the example of
the paper "Parameter Estimation of ARMA Models with GARCH/APARCH
Errors" available at the rmetrics website. In this paper you can also
find an example how to calculate the hessian matrix.

What is the dataset and the parameters you are trying to estimate?

regards,
Yohan  

regards,
Yohan

 "DK" == Desislava Kavrakova <[EMAIL PROTECTED]>
 on Mon, 18 Aug 2008 02:03:05 -0700 (PDT)

   DK> Hello R-list-members,
   DK> I'm trying to model ARMA(0,2) & GARCH(1,1) process using the code below, 
but according to my textbook, the estimated parameters are wrong. The 
MA-parameters should be negative. (I've got the same problem using garchFit()). 
Can anyone tell me what I'm doing wrong? And how can I calculate the hessian 
matrix?
   DK> 
   DK> Many thanks,
   DK> Desislava Kavrakova
   DK> 
   DK> 
   DK> Code:
   DK> 
   DK> garch<-function(x){
   DK> x<<-ts(x)
   DK> n<-length(x)
   DK> 
   DK> Mean = mean(x); Var = var(x); S = 1e-6
   DK> param = c(a = Mean, b1 = S, b2 = S, alpha0 = 0.1*Var,alpha = 0.1, beta = 
0.8)
   DK> lowerB = c(a = -10*abs(Mean), b1 = S-1, b2 = S-1, alpha0 = S^2, alpha = 
S, beta = S)
   DK> upperB = c(a = 10*abs(Mean), b1 = 1-S, b2 = 1-S, alpha0 = 100*Var, alpha 
= 1-S, beta = 1-S)
   DK> 
   DK> llh<-function(p){
   DK> a<-p[1]
   DK> b1<-p[2]
   DK> b2<-p[3]
   DK> alpha0<-p[4]
   DK> alpha<-p[5]
   DK> beta<-p[6]
   DK> 
   DK> res<-array(length(x))
   DK> hh<-array(length(x))
   DK> res[1]<-x[1]-a
   DK> res[2]<-x[2]-a-b1*res[1]
   DK> for (i in 3:n){
   DK> res[i]<-x[i]-a-b1*res[i-1]-b2*res[i-2]
   DK> }
   DK> res<-ts(res)
   DK> hh[1]<-alpha0
   DK> for (i in 2:n){
   DK> hh[i]<-alpha0+alpha*res[i-1]^2+beta*(hh[i-1]-alpha0)
   DK> }
   DK> hh<-ts(hh)
   DK> h<-sqrt(abs(hh))
   DK>  
   DK> -sum(log(dnorm(x=res/h)/h))
   DK> }
   DK> 
   DK> fit<-nlminb(param, llh, lower=lowerB, upper=upperB)
   DK> fit$par
   DK> }

-- 
PhD student
Swiss Federal Institute of Technology
Zurich

www.ethz.ch

__
R-help@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.

Re: [R] Wichmann-Hill Random Number Generator and the Birthday Problem

[Shengqiao Li]

In fact, that is what I saw in your RNG help in R 2.7.1:

 '"Wichmann-Hill"' The seed, '.Random.seed[-1] == r[1:3]' is an
  integer vector of length 3, where each 'r[i]' is in '1:(p[i]
  - 1)', where 'p' is the length 3 vector of primes, 'p =
  (30269, 30307, 30323)'. The Wichmann-Hill generator has a
  cycle length of 6.9536e12 (= 'prod(p-1)/4', see _Applied
  Statistics_ (1984) *33*, 123 which corrects the original
  article).

You have to admit that Wichmann-Hill's correction appeared here! If it's correct, then 2.8e13, 
which is claimed by Duncan Murdoch in previous email, is wrong!



Shengqiao

On Tue, 19 Aug 2008, Prof Brian Ripley wrote:


Please don't lie!  You falsely claimed:

In R version 2.7.1, Brian Ripley adopted Wichmann's 1984 correction in RNG 
document.


The period in the help has been unchanged since at least Jan 2000 (since 
Random.Rd has not had any changes to those lines since then).


Your coincidence calculations may be correct for _independent_ draws from a 
discrete distribution on M values, but independence is not satisfied.
Yet again, you are trying to do things that any good text on simulation would 
warn you against, and which (in a thread on R-devel) you have already been 
told a good way to do.


On Sun, 17 Aug 2008, Shengqiao Li wrote:



On Sun, 17 Aug 2008, Duncan Murdoch wrote:


Shengqiao Li wrote:

Dear all,

Recently I am generating large random samples (10M) and any duplicated 
numbers are not desired.
We tried several RNGs in R and found Wichmann-Hill did not produce 
duplications.


The duplication problem is the interesting birthday problem. If there are 
M possible numbers, randomly draw N numbers from them,

the average number of dupilcations D = N(N-1)/2/M.

For Knuth-TAOCP and Knuth-TAOCP-2002, M=2^30, since this modulus is used. 
D = 46566.12 for N=10M samples.


For Marsaglia-Multicarry, Super-Duper and Mersene-Twister, M=2^32. D = 
11641.53 for N = 10M samples.


My testing results (see below) agree with above analysis. But for 
Wichmann-Hill, it wasn't. Wichmann-Hill's cycle is 6.9536e12 (refer to 
RNG help by ?RNG and Whichmann's correction in 1984). Thus M <= 
6.9536e12. D 
= 7.19052 when N=1e7 and D>= 179.763 for N=5e7. 

But in my tests, duplications were surprisingly not observed.

It seems that Wichmann-Hill has a much larger cycle than the one 
documented!


Anybody can solve this puzzle?



As I told you, look at the code. The cycle length of Wichmann-Hill in R 
appears to be 2.8e13, and that also appears to be M in your notation. 
Your birthday problem calculation does not apply here.
You could probably get a good approximation to it by rounding the 
Wichmann-Hill output (e.g. look at round(2^30 * runif()), or maybe more 
severe rounding).


M is larger than what was previously documented, but Brian Ripley has 
corrected the documentation.


Wichmann claimed 2.78X10^13 in his 1982 original paper. He made correction 
in 1984, "The period of the generator was incorrectly given as 2.78 X 
10^13. In fact its period is only a quarter of that value: 6.95X10^12." In 
R version 2.7.1, Brian Ripley adopted Wichmann's 1984 correction in RNG 
document. That is the currently documented cycle is 6.95X10^12 not 
2.78X10^13 nor your approxmation 2.8e13.


So, is the cycle claimed in 1982 correct or the one claimed in 1984?

Even if the larger one 2.78e13 claimed in 1982 is correct,  around 45 
duplications were expected for 50M samples. But I got 0.  My testing method 
is based on the example code in the last third line in RNG help.


Shengqiao Li



Duncan Murdoch

Regards,

Shengqiao Li

Department of Statistics
West Virgina Unversity

==Testing===



RNGkind(kind="Knuth-TAOCP");
RNGkind();


[1] "Knuth-TAOCP" "Inversion"


sum(duplicated(runif(1e7)));


[1] 46216


RNGkind(kind="Knuth-TAOCP-2002");
RNGkind();


[1] "Knuth-TAOCP-2002" "Inversion"


sum(duplicated(runif(1e7)));


[1] 46373



RNGkind(kind="Marsaglia-Multicarry");
RNGkind();


[1] "Marsaglia-Multicarry" "Inversion"


sum(duplicated(runif(1e7)));


[1] 11671


RNGkind(kind="Super-Duper");
RNGkind();


[1] "Super-Duper" "Inversion"


sum(duplicated(runif(1e7)));


[1] 11704


RNGkind(kind="Mersenne-Twister");
RNGkind();


[1] "Mersenne-Twister" "Inversion"


sum(duplicated(runif(1e7)));


[1] 11508


RNGkind(kind="Wichmann-Hill");
RNGkind();


[1] "Wichmann-Hill" "Inversion"


sum(duplicated(runif(1e7)));


[1] 0



gc()


   used (Mb) gc trigger  (Mb)  max used   (Mb)
Ncells  199157  5.4 646480  17.3  10149018  271.1
Vcells 4497151 34.4  108714629 829.5 169585390 1293.9



sum(duplicated(runif(5e7)))


[1] 0




sessionInfo()


R version 2.7.1 (2008-06-23)
i386-pc-mingw32

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United 
States.1252;LC_MONETARY=English_United 
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252


attached base packages:
[1] stats graphics  grDevices utils 

[R] jpeg, dev.off() inside function

[marciarr]

Dear all,
I am trying to plot and save the plots. I have written a small function for
that in which I use the jpeg and dev.off() commands. Strangely, when these
commands are inside the function, they are not working anymore.
For instance:

dataframe<- data.frame (x=c(1:20),y=seq(0,9.5, length=20))

for (i in seq(0,15, length=4)){
setwd ("C:/R")
jpeg(filename = paste("test ",deparse(substitute(series)),
i, ".jpg"), width = 1200,
height = 800,pointsize = 12, quality = 400, bg = "white")
plot(dataframe[1+i:((nrow(z)/4)+i),], type="l")
dev.off()} 

This works perfectly. But if I try to make a function out of it, such as:

X<-function (x, y = NULL, plot = TRUE, na.action = na.fail)
{
   if (NCOL(x) == 2) {
series <- deparse(substitute(x))
}
for (i in seq(0,15, length=4)){
setwd ("F:/R/P to B manuscript/area under the curve")
jpeg(filename = paste("test ",deparse(substitute(series)),
i, ".jpg"), width = 1200,
height = 800,pointsize = 12, quality = 400, bg = "white")
plot(z[1+i:((nrow(z)/4)+i),], type="l")
dev.off()}
return(dataframe)}

and try to run:
X(dataframe)
I get this error message:

Error in jpeg(filename = paste("test ", deparse(substitute(series)), i,  : 
  unable to start device devWindows
In addition: Warning message:
In jpeg(filename = paste("test ", deparse(substitute(series)), i,  :
  Unable to open file 'test  "dataframe" 0 .jpg' for writing

Could anyone help me???
Thank you and have a nice day,
Marcia
-- 
View this message in context: 
http://www.nabble.com/jpeg%2C-dev.off%28%29-inside-function-tp19048365p19048365.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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] Querry

[Basnet, Ram Kumar]

Dear all,
 
I am student and doing cluster analysis using "pvclust" package. I found
this package very useful. After bootstrapping, i found the edges number
88, 16, 6 have very high standard error. I got this information "seplot"
fuction and "identify ()". Now i am interested to see the name list of
genotypes (samples) which fall in those edges number(88,16,6). Can you
help me by providing the appropriate R-command (same like list of
genotypes lies in each cluster number when i use pvpick(x,alpha = 0.98)
function but in this case, i am intereseted to know the name list based
on edges number). I will be very very greatful with you.
Looking forward to getting your prompt response.
Thank you.
 
Sincerely,
Ram Kumar Basnet,
MSc student,
Department of plant breeding and genetics.
Wageningen University.
The Netherlands.

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Wichmann-Hill Random Number Generator and the Birthday Problem

[Shengqiao Li]

Prof Brian Ripley wrote:
Your coincidence calculations may be correct for _independent_ draws from a 
discrete distribution on M values, but independence is not satisfied.
Yet again, you are trying to do things that any good text on simulation 
would warn you against, and which (in a thread on R-devel) you have already 
been told a good way to do.


A good example are the "good" old linear congruential random generators. 
These will start repeating at the first coincidence, for the pretty obvious 
reason that the next random number is a function only of the previous one. So 
the number of distinct values in N draws is min(N, cycle_length). In 
particular, the number of coincidences is 0 when N is less than cycle_length.


So, we can conclude that the programming implementation for Wichmann-Hill 
does not reduce the number of possible values but it does for other 
generators. Moreover, the cycle_lenght for other generators in R are quite 
long and these generators are more random than those old, like the W-H 
genrator. When wrapped to 32-bit or 30-bit integers, ties are easily and 
randomly produced for these generators of long cycles. Since they are 
more random, the Poisson appxoimation for the number of ties is appliable and quite accurate.


This sugguest that W-H is not good as others in the terms of randomness?

The following  paragraph in RNG help is misleading since it does not apply 
Wichman-Hill generator. It should be updated.


"All the supplied uniform generators return 32-bit integer values
 that are converted to doubles, so they take at most 2^32 distinct
 values and long runs will return duplicated values."


Shengqiao Li



--
 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
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

__
R-help@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] how to import from SPSS without shortening variable names

[Timo Stolz]

Hello,

as I import '.sav' files from SPSS, the variable names are shortened to
8 uppercase characters:

"sex_of_therapist" will become "SEX_OF_TH"

Is there a way around this? How can I retrieve the full names?

Greets from Southern Germany,
Timo Stolz

__
R-help@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.

Re: [R] another GARCH problem

[Yohan Chalabi]

 "C" == collonil <[EMAIL PROTECTED]>
 on Mon, 18 Aug 2008 07:35:09 -0700 (PDT)

   C> Hallo,
   C> i want to fit a GARCH model with a extern regressor (without arma
   C> components), so i found the following function in package fGarch. I tryed
   C> out a lot of things but usually I get this Error.
   C> 
   C> > garchFit(formula=y~x, formula.var=~garch(1,1),data=w)  
   C> Error in .garchFit(formula.mean, formula.var, series = x, init.rec, 
delta, 
   C> : 
   C>   Algorithm only supported for mci Recursion
   C> 
   C> I think i use the function the wrong way, but cannot find my mistake. Can
   C> somebody help me?

it looks like one of the depended package of fGarch has a change which
breaks garchFit. I would recommend you to use the development version.

If you have R (>=2.7.0) and everything installed in order to compile a
package, you can use the following command to install the development
version of the fGarch package.

source("http://rmetrics.r-forge.r-project.org/installRmetrics.R";)
installRmetrics("fGarch", repos = "http://r-forge.r-project.org";)

but bear in mind that this is development code. So if you have any
suggestion how to improve it, please let us know it!

hope this helps,
Yohan

-- 
PhD student
Swiss Federal Institute of Technology
Zurich

www.ethz.ch

__
R-help@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.

Re: [R] Exponential smoothing?

[S Ellison]

An r hel search picks up
http://tolstoy.newcastle.edu.au/R/help/01a/1162.html 

from 2001, which points towards the ts package.



>>> Öhagen Patrik <[EMAIL PROTECTED]> 19/08/2008 09:00:12 >>>
Dear List,

I have used all my resources (i.e. "help.search) and I still havn't been able 
to figure out if there is an Exponential Smoothing command in R.

Thank you in advance!

__
R-help@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.

***
This email and any attachments are confidential. Any use...{{dropped:8}}

__
R-help@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] S.O.S.

[Mauro Belleggia]

Hi, My name is Mauro Belleggia, I m doing an ANOSIM (using "vegan" package)
to test statistically whether there is a significant difference between 3
(1,2 and 3) groups of sampling units. My results were:

Call:
anosim(dis = distancias, grouping = Lt)
Dissimilarity: bray

ANOSIM statistic R: 0.1313
  Significance: < 0.001

Based on  1000  permutations

Empirical upper confidence limits of R:
   90%95%  97.5%99%
0.0331 0.0464 0.0569 0.0693

Dissimilarity ranks between and within classes:
0%25%50%75%   100%N
Between  2 950.75 1786.5 2641.5 3402.5 2220
1   13 958.50 1927.0 2660.5 3402.5  595
21 629.50 1318.0 2050.5 3258.0  435
3   10 356.00  840.0 1599.0 3367.0  153

I know that there where dissimilarities among grops but i can t interpret
between which. 1 Vs 2?, 2 Vs. 3? Or 1 Vs. 3?

I will be very very greatful with you.
Looking forward to getting your prompt response.


-- 
[EMAIL PROTECTED]

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] how to import from SPSS without shortening variable names

[Peter Dalgaard]

Timo Stolz wrote:
> Hello,
>
> as I import '.sav' files from SPSS, the variable names are shortened to
> 8 uppercase characters:
>
> "sex_of_therapist" will become "SEX_OF_TH"
>
> Is there a way around this? How can I retrieve the full names?
>
>   
Well as long as it didn't become sex_of_the_rapist...

This is actually changed in the most recent version of the foreign
package, which is contained in R 2.7.2RC:

> X <- read.spss("~/Desktop/Biobank_SPSS.sav", to.data.frame=TRUE)
Warning message:
In read.spss("~/Desktop/Biobank_SPSS.sav", to.data.frame = TRUE) :
  ~/Desktop/Biobank_SPSS.sav: File-indicated character representation
code (1252) looks like a Windows codepage
> names(X)
 [1] "Pt_id"   "Popul"   "Age" "Sex" "SKScode"
 [6] "Sys1""Sys2""Sys_pro" "Dia1""Dia2"
[11] "Dia_pro" "MAP1""MAP2""MAP_pro" "MAP_p"
[16] "MAP_30_55"   "Anaesthesia" "Efedrin" "Efedrin1""Efedrin2"
[21] "Metaoxedrin" "Nr"  "ACE1""ACE2""ACE_ID"
[26] "ACE3""ACE4""ADRA1A"  "ADRA2A"  "ADRB1_1"
[31] "ADRB1_2" "ADRB2_1" "ADRB2_2" "ADRB2_3" "ADRB2_4"
[36] "AGT_1"   "AGT_2"   "AGTR1"   "GNB3""AGTR2"
[41] "GNAS""NOS3_1"  "NOS3_2"  "NOS3_3"  "ADRB2_5"
[46] "ADRB2HA" "ADRB2HB" "ADRB2HP" "ADRB1HA" "ADRB1HB"
[51] "ADRB1HP" "AGTHA"   "AGTHB"   "AGTHP"   "ACEHA"
[56] "ACEHB"   "ACEHP"   "Højde"   "Vægt""BMI"
[61] "BMIgroupA"   "BMIgroupB"   "BMIgroupC"

So try the prerelease or wait till Monday (and miss the opportunity to
have any bugs you find fixed before release).

-- 
   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
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

__
R-help@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.

Re: [R] jpeg, dev.off() inside function

[Henrik Bengtsson]

In your function, split up your code in two lines as:

pathname <- paste("test ", deparse(substitute(series)), i, ".jpg");
jpeg(filename=pathname, width=1200, height=800, pointsize=12,
quality=400, bg="white");

The focus on what 'pathname' is/becomes when you put it inside your
function.  That will help you troubleshoot.  It has most likely
nothing to do with the jpeg() function.

My $.02

/Henrik

PS. 'series' is not set inside the function if there are not two columns. DS.


On Tue, Aug 19, 2008 at 6:09 AM, marciarr <[EMAIL PROTECTED]> wrote:
>
> Dear all,
> I am trying to plot and save the plots. I have written a small function for
> that in which I use the jpeg and dev.off() commands. Strangely, when these
> commands are inside the function, they are not working anymore.
> For instance:
>
> dataframe<- data.frame (x=c(1:20),y=seq(0,9.5, length=20))
>
> for (i in seq(0,15, length=4)){
>setwd ("C:/R")
>jpeg(filename = paste("test ",deparse(substitute(series)),
>i, ".jpg"), width = 1200,
>height = 800,pointsize = 12, quality = 400, bg = "white")
>plot(dataframe[1+i:((nrow(z)/4)+i),], type="l")
>dev.off()}
>
> This works perfectly. But if I try to make a function out of it, such as:
>
> X<-function (x, y = NULL, plot = TRUE, na.action = na.fail)
> {
>   if (NCOL(x) == 2) {
>series <- deparse(substitute(x))
>}
>for (i in seq(0,15, length=4)){
>setwd ("F:/R/P to B manuscript/area under the curve")
>jpeg(filename = paste("test ",deparse(substitute(series)),
>i, ".jpg"), width = 1200,
>height = 800,pointsize = 12, quality = 400, bg = "white")
>plot(z[1+i:((nrow(z)/4)+i),], type="l")
>dev.off()}
>return(dataframe)}
>
> and try to run:
> X(dataframe)
> I get this error message:
>
> Error in jpeg(filename = paste("test ", deparse(substitute(series)), i,  :
>  unable to start device devWindows
> In addition: Warning message:
> In jpeg(filename = paste("test ", deparse(substitute(series)), i,  :
>  Unable to open file 'test  "dataframe" 0 .jpg' for writing
>
> Could anyone help me???
> Thank you and have a nice day,
> Marcia
> --
> View this message in context: 
> http://www.nabble.com/jpeg%2C-dev.off%28%29-inside-function-tp19048365p19048365.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@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@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] spatial probit/logit for prediction

[Freeman, Robert]

Hello all,

I am wondering if there is a way to do a spatial error probit/logit model in R? 
 I can't seem to find it in any of the packages.  I can do it in MATLAB with 
Gibbs sampling, but would like to confirm the results.  Ideally I would like to 
use this model to predict probability of parcel conversion in a future time 
period.  This seems especially difficult in a binary outcome model setting.  
Does anyone have any guidance on this aspect?  

Thanks.
Robb

__
R-help@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.

Re: [R] Restricted Least Squares

[Greg Snow]

Nelson,

It sounds like you should create a factor (or some factors) for your dummy 
variables, then use contr.sum to automatically create the dummy variables for 
you such that the coeficients would sum to 0.  See ?C, ?contrasts,
?contr.sum, and ?dummy.coef for details.

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111



> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Nelson Villoria
> Sent: Saturday, August 16, 2008 9:20 AM
> To: Prof Brian Ripley
> Cc: r-help@r-project.org
> Subject: Re: [R] Restricted Least Squares
>
> Thank you very much Dr. Brian -- I appreciate your response.
> I understand the idea of building the constraints into the
> equation but my problem is that i have to deal with
> potentially hundreds of dummy variables, and I want the sum
> of their estimates to equal zero, so LS can be fitted with an
> intercept and a singular X'X matrix - moreover, in different
> regressions the number of dummies is also potentially
> different so i was hoping to use some restrictions of the
> form R*beta = 0, where R is the matrix of restrictions and
> beta the dummy parameters. I have written my own estimator
> based on Greene and Seeks (REView of Economics and
> STATistics, 1991), but wanted to know if there was a package
> out there making things easier.
>
> Nelson Villoria
>
> On Sat, Aug 16, 2008 at 2:24 AM, Prof Brian Ripley
> <[EMAIL PROTECTED]> wrote:
> > On Sat, 16 Aug 2008, Nelson Villoria wrote:
> >
> >> Dear R experts:
> >>
> >> Is there any package that estimates Restricted Least Squares?
> >>
> >> Specifically, If I want to fit:
> >>
> >> G = b0 + b1(Y) + a1(X1) + a2(X2) + a3(X3) + a4(X4) where
> Y, X1 to X4
> >> are variables and b's and a's parameters to be estimated.
> >>
> >> I want to impose a1 + a2 + a3 + a4 = 0.
> >
> > You don't need a package to do that, just re-parametrize.  It is
> >
> > G ~ Y + I(X1-X4) + I(X2-X4) + I(X3-X4)
> >
> > --
> > Brian D. Ripley,  [EMAIL PROTECTED]
> > 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, UKFax:  +44 1865 272595
> >
>
> __
> R-help@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@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.

Re: [R] how can i get hessian matrix at "constrOptim"

[Roland Rau]

Hi,

Kanak Choudhury wrote:

Hi,
i have made a code for optimizing a function using "constrOptim". i need
hessain matrix of the parameters. how could i get hessain matrix when i will
use "constrOptim"? May i get get any help from anyone?



the function fdHess from package 'nlme' can help you?

library(nlme)
?fdHess

I hope this helps,
Roland

__
R-help@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.

Re: [R] To package or not to package?

[Greg Snow]

Not a full function to do this, but one of the examples for the clipplot 
function in the TeachingDemos package shows one way to do something similar to 
what you describe.

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111



> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Carl Witthoft
> Sent: Saturday, August 16, 2008 11:54 AM
> To: r-help@r-project.org
> Subject: [R] To package or not to package?
>
> I'd like a little advice as to when it's appropriate to
> create an R-package and submit it, as opposed to just
> providing the source to some simple code.
> In this case, I've written a function that draws a line plot
> (with options for points, etc) where the color of the line
> changes at specified values of the y-data (e.g., it's one
> color below -10, another color between -10 and -5, etc).
> It's pretty clean, and has a few error-checks,
> self-correctors, etc., so I would be happy to provide it to
> the community as a whole.
> So, is this worthy of a package, or should I just post the
> function code (well commented)?
> Or is this feature available deep inside some graphing
> package I haven't found yet? :-(
>
> thanks for all advice.
> Carl
>
> __
> R-help@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@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.

Re: [R] jpeg, dev.off() inside function

[Prof Brian Ripley]

On Tue, 19 Aug 2008, Henrik Bengtsson wrote:


In your function, split up your code in two lines as:

pathname <- paste("test ", deparse(substitute(series)), i, ".jpg");
jpeg(filename=pathname, width=1200, height=800, pointsize=12,
quality=400, bg="white");

The focus on what 'pathname' is/becomes when you put it inside your
function.  That will help you troubleshoot.  It has most likely
nothing to do with the jpeg() function.


It isn't:


 Unable to open file 'test  "dataframe" 0 .jpg' for writing

  ^^ ^^ ^
and " is invalid in a Windows filename.  (I also think sep="" was intended 
in paste().  This is because 'series' is 'x' deparsed twice.




My $.02

/Henrik

PS. 'series' is not set inside the function if there are not two columns. DS.


On Tue, Aug 19, 2008 at 6:09 AM, marciarr <[EMAIL PROTECTED]> wrote:


Dear all,
I am trying to plot and save the plots. I have written a small function for
that in which I use the jpeg and dev.off() commands. Strangely, when these
commands are inside the function, they are not working anymore.
For instance:

dataframe<- data.frame (x=c(1:20),y=seq(0,9.5, length=20))

for (i in seq(0,15, length=4)){
   setwd ("C:/R")
   jpeg(filename = paste("test ",deparse(substitute(series)),
   i, ".jpg"), width = 1200,
   height = 800,pointsize = 12, quality = 400, bg = "white")
   plot(dataframe[1+i:((nrow(z)/4)+i),], type="l")
   dev.off()}

This works perfectly. But if I try to make a function out of it, such as:

X<-function (x, y = NULL, plot = TRUE, na.action = na.fail)
{
  if (NCOL(x) == 2) {
   series <- deparse(substitute(x))
   }
   for (i in seq(0,15, length=4)){
   setwd ("F:/R/P to B manuscript/area under the curve")
   jpeg(filename = paste("test ",deparse(substitute(series)),
   i, ".jpg"), width = 1200,
   height = 800,pointsize = 12, quality = 400, bg = "white")
   plot(z[1+i:((nrow(z)/4)+i),], type="l")
   dev.off()}
   return(dataframe)}

and try to run:
X(dataframe)
I get this error message:

Error in jpeg(filename = paste("test ", deparse(substitute(series)), i,  :
 unable to start device devWindows
In addition: Warning message:
In jpeg(filename = paste("test ", deparse(substitute(series)), i,  :
 Unable to open file 'test  "dataframe" 0 .jpg' for writing

Could anyone help me???
Thank you and have a nice day,
Marcia
--
View this message in context: 
http://www.nabble.com/jpeg%2C-dev.off%28%29-inside-function-tp19048365p19048365.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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@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,  [EMAIL PROTECTED]
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, UKFax:  +44 1865 272595

__
R-help@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.

Re: [R] how to import from SPSS without shortening variable names

[Prof Brian Ripley]

On Tue, 19 Aug 2008, Timo Stolz wrote:


Hello,

as I import '.sav' files from SPSS, the variable names are shortened to
8 uppercase characters:

"sex_of_therapist" will become "SEX_OF_TH"

Is there a way around this? How can I retrieve the full names?


Have you updated your 'foreign' package recently?  I think that was 
changed in 0.8-27, and 0.8-29 is current.


[sessionInfo(), which the posting guide asked you to paste in the output 
of, would have told us.]


--
Brian D. Ripley,  [EMAIL PROTECTED]
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, UKFax:  +44 1865 272595

__
R-help@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.

Re: [R] nonlinear constrained optimization

[Paul Smith]

On Tue, Aug 19, 2008 at 8:35 AM, Chua Siang Li
<[EMAIL PROTECTED]> wrote:
>
>   Hi.  I need some advises on how to use R to find pi (i is the index) with
>   the following objective function and constraint:
>
>   max (sum i)[ f(ai, bi, pi) * g(ci, di, pi) * Di ]
>
>   s.t.  (sum i)[ f(ai, bi, pi) * Di * pi] / (sum i)[ f(ai, bi, pi) * Di ]  <=
>   constant
>
>   f and g are diffentiable.
>   So, I am thinking of optim with method = "BFGS"? But wonder how to include
>   the constraint.
>   I also saw constrOptim but it says for linear constraint only.
>   So, what will be the suitable function to use then?
>   Also, how should I write (sum i) in R?
>   Many thanks and your help is much appreciated.

Up to my best knowledge, R cannot deal with optimization problems with
nonlinear constraints, unless one uses the penalty method. Outside R,
Ipopt and Algencan can solve problems like yours, but one needs to
program in AMPL and/or C/Fortran.

Paul

__
R-help@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.

Re: [R] Exponential smoothing?

[Giovanni Petris]

See also package forecast.

Giovanni

> Date: Tue, 19 Aug 2008 15:25:35 +0100
> From: S Ellison <[EMAIL PROTECTED]>
> Sender: [EMAIL PROTECTED]
> Precedence: list
> 
> An r hel search picks up
> http://tolstoy.newcastle.edu.au/R/help/01a/1162.html 
> 
> from 2001, which points towards the ts package.
> 
> 
> 
> >>> Öhagen Patrik <[EMAIL PROTECTED]> 19/08/2008 09:00:12 >>>
> Dear List,
> 
> I have used all my resources (i.e. "help.search) and I still havn't been able 
> to figure out if there is an Exponential Smoothing command in R.
> 
> Thank you in advance!
> 
> __
> R-help@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.
> 
> ***
> This email and any attachments are confidential. Any use...{{dropped:8}}
> 
> __
> R-help@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.
> 
> 
> 

-- 

Giovanni Petris  <[EMAIL PROTECTED]>
Associate Professor
Department of Mathematical Sciences
University of Arkansas - Fayetteville, AR 72701
Ph: (479) 575-6324, 575-8630 (fax)
http://definetti.uark.edu/~gpetris/

__
R-help@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.

Re: [R] A doubt about "lm" and the meaning of its summary

[Charles C. Berry]

Alberto,

Moshe was actually on the right track.

Your model does not satisfy enough of the details in the usual 
Gauss-Markov setup to guarantee unbiasedness of the LS solution.


In particular, the errors must be orthogonal to the design columns for the 
unbiasedness of the LS solution.


But, in your case, they are not.

This is easily seen when n == 2, m==0, c==0.

Take eps <- rnorm(2)

Then, x = c( 0, eps[1] ), y <- c( eps[1], eps[2] ).

The expectation of y is evidently 0, so the observed y is the error 
vector.


With basic linear algebra, you can verify that the vector rbind( 1, x ) 
%*% y has a non-zero expectation and that the regression coefficient for 
'x' has negative expectation.


And just for fun:

The correlation when n == 2 is either +1 or -1. You can easily deduce from 
a simple graphical argument (or equivalent arguments invoking symmetry in 
the bivariate distribution of eps) that -1 correlations are exactly 
3 times as likely as +1 correlations.


HTH,

Chuck


On Mon, 18 Aug 2008, Moshe Olshansky wrote:


Hi Alberto,

Please disregard my previous note - I probably had a black-out!!!


--- On Tue, 19/8/08, Alberto Monteiro <[EMAIL PROTECTED]> wrote:


From: Alberto Monteiro <[EMAIL PROTECTED]>
Subject: [R] A doubt about "lm" and the meaning of its summary
To: r-help@r-project.org
Received: Tuesday, 19 August, 2008, 4:31 AM
I have a conceptual problem (sort of). Maybe there's a
simple solution,
maybe not.

First, let me explain the test case that goes ok.

Let x be a (fixed) n-dimensional vector. I simulate a lot
of linear models,
y = m x + c + error, then I do a lot of regressions. As
expected, the
estimated values of m (let's call them m.bar) are
distributed according to a
Student's t distribution.

More precisely (test case, it takes a few minutes to run):
#
# start with fixed numbers
#
  m <- sample(c(-1, -0.1, 0.1, 1), 1)
  c <- sample(c(-1, -0.1, 0, 0.1, 1), 1)
  sigma <- sample(c(0.1, 0.2, 0.5, 1), 1)
  n <- sample(c(4,5,10,1000), 1)
  x <- 1:n # it's possible to use other x
  NN <- sample(c(3000, 4000, 5000), 1)
  m.bar <- m.std.error <- 0  # these vectors will
hold the simulation output

#
# Now let's repeat NN times:
# simulate y
# regress y ~ x
# store m.bar and its error
#
  for (i in 1:NN) {
y <- m * x + c + sigma * rnorm(n)
reg <- lm(y ~ x)
m.bar[i] <- summary(reg)$coefficients['x',
'Estimate']
m.std.error[i] <-
summary(reg)$coefficients['x', 'Std. Error']
  }
#
# Finally, let's analyse it
# The distribution of (m.bar - m) / m.std.error is
# a Student's t with n - 2 degrees of freedom.
# Is it? Let's test!
#
  print(ks.test((m.bar - m) / m.std.error, rt(NN, n - 2)))

Now, this goes ok, with ks.test returning a number
uniformly distributed in
the 0-1 interval.

Next, I did a slightly different model. This model is a
reversion model,
where the simulated random variable lp goes along the
equation:
lp[t + 1] <- (1 + m) * lp[t] + c + error

I tried to use the same method as above, defining
x = lp[1:n], y = lp[2:(n+1}] - lp[1:n], with equation
y <- m * x + c + error

And now it breaks. Test case (it takes some minutes to
run):

#
# start with fixed numbers
#
  m <- runif(1, -1, 0)  # m must be something in the
(-1,0) interval
  c <- sample(c(-1, -0.1, 0, 0.1, 1), 1)
  sigma <- sample(c(0.1, 0.2, 0.5, 1), 1)
  n <- sample(c(4,5,10,1000), 1)
  NN <- sample(c(3000, 4000, 5000), 1)
  m.bar <- m.std.error <- 0  # these vectors will
hold the simulation output

#
# Now let's repeat NN times:
# simulate y
# regress y ~ x
# store m.bar and its error
#
  for (i in 1:NN) {
lp <- 0
lp[1] <- 0
for (j in 1:n) {
  lp[j+1] = (1 + m) * lp[j] + c + sigma * rnorm(1)
}
x <- lp[1:n]
y <- lp[-1] - x
reg <- lm(y ~ x)
m.bar[i] <- summary(reg)$coefficients['x',
'Estimate']
m.std.error[i] <-
summary(reg)$coefficients['x', 'Std. Error']
  }
#
# Finally, let's analyse it
# The distribution of (m.bar - m) / m.std.error should be
# a Student's t with n - 2 degrees of freedom.
# Is it? Let's test!
#
  print(ks.test((m.bar - m) / m.std.error, rt(NN, n - 2)))

... and now it's not.

What's wrong with this? I suspect that the model y ~ x
does only give
meaningful values when x is deterministic; in this case x
is stochastic. Is
there any correct way to estimate this model giving
meaningful outputs?

Alberto Monteiro

__
R-help@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@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
D

Re: [R] nonlinear constrained optimization

[Ravi Varadhan]

Hi,

The "Tango project" website (from Brazil) has some R functions to implement
Algencan.  I haven't used it, but I would be interested to hear others'
experience if someone has already used it or will be using it.


Ravi.

-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Paul Smith
Sent: Tuesday, August 19, 2008 12:57 PM
To: r-help@r-project.org
Subject: Re: [R] nonlinear constrained optimization

On Tue, Aug 19, 2008 at 8:35 AM, Chua Siang Li
<[EMAIL PROTECTED]> wrote:
>
>   Hi.  I need some advises on how to use R to find pi (i is the index)
with
>   the following objective function and constraint:
>
>   max (sum i)[ f(ai, bi, pi) * g(ci, di, pi) * Di ]
>
>   s.t.  (sum i)[ f(ai, bi, pi) * Di * pi] / (sum i)[ f(ai, bi, pi) * Di ]
<=
>   constant
>
>   f and g are diffentiable.
>   So, I am thinking of optim with method = "BFGS"? But wonder how to
include
>   the constraint.
>   I also saw constrOptim but it says for linear constraint only.
>   So, what will be the suitable function to use then?
>   Also, how should I write (sum i) in R?
>   Many thanks and your help is much appreciated.

Up to my best knowledge, R cannot deal with optimization problems with
nonlinear constraints, unless one uses the penalty method. Outside R, Ipopt
and Algencan can solve problems like yours, but one needs to program in AMPL
and/or C/Fortran.

Paul

__
R-help@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@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.

Re: [R] how can i get hessian matrix at "constrOptim"

[Ravi Varadhan]

Hi,

Hessian for a consttarined problem can be tricky, when the converged
parameter vector lies on the boundary, i.e. when one or more of the
constraints are active at convergence.  This can be handled by
reparametrization when the constraints are linear, but is difficult
otherwise.  

Ravi.

-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Roland Rau
Sent: Tuesday, August 19, 2008 11:51 AM
To: Kanak Choudhury
Cc: r-help@r-project.org
Subject: Re: [R] how can i get hessian matrix at "constrOptim"

Hi,

Kanak Choudhury wrote:
> Hi,
> i have made a code for optimizing a function using "constrOptim". i 
> need hessain matrix of the parameters. how could i get hessain matrix 
> when i will use "constrOptim"? May i get get any help from anyone?
> 

the function fdHess from package 'nlme' can help you?

library(nlme)
?fdHess

I hope this helps,
Roland

__
R-help@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@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] Calling a windows program

[Benjamin E. Alexander-Eitzman]

Greetings,

 

How do I make a call to a windows program from R? For example in
windows, I would run this in a  shortcut:

 

"C:\Program Files\WinBUGS14\WinBUGS14.exe" /PAR "myownscript.odc"

 

Thanks,

 

Ben Alexander-Eitzman


[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Calling a windows program

[Henrique Dallazuanna]

See shell command

On Tue, Aug 19, 2008 at 3:48 PM, Benjamin E. Alexander-Eitzman
<[EMAIL PROTECTED]> wrote:
> Greetings,
>
>
>
> How do I make a call to a windows program from R? For example in
> windows, I would run this in a  shortcut:
>
>
>
> "C:\Program Files\WinBUGS14\WinBUGS14.exe" /PAR "myownscript.odc"
>
>
>
> Thanks,
>
>
>
> Ben Alexander-Eitzman
>
>
>[[alternative HTML version deleted]]
>
> __
> R-help@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.
>



-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O

__
R-help@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.

Re: [R] Calling a windows program

[Uwe Ligges]

See ?shell

Uwe Ligges



Benjamin E. Alexander-Eitzman wrote:

Greetings,

 


How do I make a call to a windows program from R? For example in
windows, I would run this in a  shortcut:

 


"C:\Program Files\WinBUGS14\WinBUGS14.exe" /PAR "myownscript.odc"

 


Thanks,

 


Ben Alexander-Eitzman


[[alternative HTML version deleted]]

__
R-help@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@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.

Re: [R] nonlinear constrained optimization

[Hans W. Borchers]

Paul Smith  gmail.com> writes:
> 
> Up to my best knowledge, R cannot deal with optimization problems with
> nonlinear constraints, unless one uses the penalty method. Outside R,
> Ipopt and Algencan can solve problems like yours, but one needs to
> program in AMPL and/or C/Fortran.
> 
> Paul
> 

Please have a look at the 'Rdonlp2' package at .
It provides non-linear optimization with nonlinear constraints.

For more information on advanced optimization procedures available in R take
a short trip to the Optimization Task View.

// Hans Werner Borchers

__
R-help@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] panel function with zoo

[stephen sefick]

n <-  c("m", "j", "j", "a", "s", "o", "n", "d", "j", "f", "m", "a",
"m", "j", "j", "a","s", "o", "n", "d", "j")
pnl.xaxis <- function(...) {
lines(...)
panel.number <- parent.frame()$panel.number
nser <- parent.frame()$nser
# if bottom panel
if (!length(panel.number) || panel.number == nser) {
rng <- range(time(y.zoo))
 axis(1, at = seq(rng[1], rng[2], 1/12), labels = n, tcl = -0.3)
 }
}
plot(y.zoo[,5:9],  panel=pnl.xaxis, type="b", main="215", ylab=names, xaxt="n")

I know this is not reproducible, but this does exactly what I want
except that it only labels the second column of graphs x-axis.  There
are two columns of graphs three in the, first stacked on top of each
other, and two in the second, stacked on top of each other.  Am I
missing something simple?

thanks

Stephen Sefick

-- 
Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods. We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@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] environment problems

[Matias Salibian-Barrera]

Hello,

I realize this is not a new problem, but I'm hoping somebody can point me in 
the right direction. 

The function silly() below calls gam() (from package gam) with different values 
of the span parameter for lo(), which is used inside the formula of gam. The 
values for span are stored in a vector that resides in the function silly(), 
and that lo(), for some reason, needs to see.

A hack to get around this is provided in silly2(), where the current value for 
span is forcefully assigned to a variable in the "top" environment, that is 
always visible (I believe?) Of course, this is not good programming practice 
(to say the least).

I'd very much appreciate it if somebody could help me find a better way to make 
silly() work.

Thanks a lot in advance.

Matias

-
> version
   _   
platform   i386-redhat-linux-gnu   
arch   i386
os linux-gnu   
system i386, linux-gnu 
status 
major  2   
minor  7.1 
year   2008
month  06  
day23  
svn rev45970   
language   R   
version.string R version 2.7.1 (2008-06-23)
> 
> silly <- function(x, y, alphas = seq(.1, .9, length=9)) {
+ a <- rep(0, 9)
+ for(i in 1:9) {
+ a[i] <- predict(gam(y~lo(x, span = alphas[i]), family=poisson), 
type='response')[23]
+ }
+ return(a)
+ }
> 
> silly2 <- function(x, y, alphas = seq(.1, .9, length=9)) {
+ a <- rep(0, 9)
+ for(i in 1:9) {
+ assign("tt", alphas[i], envir=.GlobalEnv)
+ a[i] <- predict(gam(y~lo(x, span = tt), family=poisson), type='response')[23]
+ }
+ return(a)
+ }
> 
> library(gam)
Loading required package: splines
> set.seed(321)
> x <- runif(50)
> y <- rpois(50, lambda=10)
> silly(x=x,y=y)
Error in lodummy(span = alphas[i]) : object "alphas" not found
> silly2(x=x,y=y)
[1] 11.56997 10.98681 10.73250 10.63455 10.62692 10.61285 10.61779 10.63234
[9] 10.66311
> 
> 


  __
[[elided Yahoo spam]]

__
R-help@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] converting coordinates from utm to longitude / latitude

[Werner Wernersen]

Hi,

is there a function in R to convert data read with read.shape and which is 
originally in UTM coordinates into longitude / latitude coordinates?
I found the convUL() function from the PBSmapping package but I have no idea 
how I could apply that to the read.shape object.

Many thanks,
  Werner


__
Do
sragenden Schutz gegen Massenmails. 
http://mail.yahoo.com

__
R-help@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] conditional with and operators

[Altaweel, Mark R.]

Hi,

I have a problem in which I am parsing data from a list.  I am simply trying to 
return data that has several conditions being true.  Here is my syntax below:

d<-sapply(res,function(.df){(.df$TimesVisited[.df$Tick>912 && .df$Id>0])})   
#res is the list, and I am trying to return a result that has two true 
conditions (that is the variable Tick should be greate than 912 and Id should 
be greater than 0

This returns an array of 10 with integer values of 0. This is the incorrect 
result

However, if I do the same syntax except I remove the && statement (i.e. the 
second conditional), then the result producing something that makes sense, 
which is all values that are Tick and greater than 912.

Can someone let me know how I can setup my data to be parsed so I can have 2 or 
multiple conditionals in my function that is looking at an array.

Thanks in advance.

Mark

__
R-help@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] Histogram binning

[Ben Cox]

I am trying to produce frequencies in defined intervals however I can't seem
to figure out how to get R to bin my data the way I want it to.

 

I have several thousand lengths of fish that I want to be binned as follows:


 

Ex.

 

Length  Bin

 

209 200

219 210

431 430

727 720

 

That is, bins with any length equal to or greater than the lower cutoff, but
smaller than the next higher cutoff should be included.  I have tried using
the right=FALSE and include.lowest in the hist function but it still puts
lengths that is greater than the cutoff into the next bin.

All software seems to bin right, but fish nerds do it a little different.

 

 

 

Thanks!  

Ben Cox

Graduate Research Assistant (M.S.)

Montana Cooperative Fishery Research Unit

Montana State University 

PO Box 173460

Bozeman, MT 59717

(406)994-6643 

 


[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Histogram binning

[jim holtman]

Is this what you want:

cuts <- c(0,209,429,719,1500)
x <- runif(1000,0,1447)
barplot(table(cut(x, breaks=cuts)))



On Tue, Aug 19, 2008 at 4:17 PM, Ben Cox <[EMAIL PROTECTED]> wrote:
> I am trying to produce frequencies in defined intervals however I can't seem
> to figure out how to get R to bin my data the way I want it to.
>
>
>
> I have several thousand lengths of fish that I want to be binned as follows:
>
>
>
>
> Ex.
>
>
>
> Length  Bin
>
>
>
> 209 200
>
> 219 210
>
> 431 430
>
> 727 720
>
>
>
> That is, bins with any length equal to or greater than the lower cutoff, but
> smaller than the next higher cutoff should be included.  I have tried using
> the right=FALSE and include.lowest in the hist function but it still puts
> lengths that is greater than the cutoff into the next bin.
>
> All software seems to bin right, but fish nerds do it a little different.
>
>
>
>
>
>
>
> Thanks!
>
> Ben Cox
>
> Graduate Research Assistant (M.S.)
>
> Montana Cooperative Fishery Research Unit
>
> Montana State University
>
> PO Box 173460
>
> Bozeman, MT 59717
>
> (406)994-6643
>
>
>
>
>[[alternative HTML version deleted]]
>
> __
> R-help@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?

__
R-help@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.

Re: [R] conditional with and operators

[Steven McKinney]

Did you try it with the vector '&' and operator?

d<-sapply(res,function(.df){(.df$TimesVisited[.df$Tick>912 & .df$Id>0])}) 

(The '&&' operator is designed for use in
e.g. if() clauses where you want a scalar logical
answer)

HTH

Steve McKinney

-Original Message-
From: [EMAIL PROTECTED] on behalf of Altaweel, Mark R.
Sent: Tue 8/19/2008 1:10 PM
To: r-help@r-project.org
Subject: [R] conditional with and operators
 
Hi,

I have a problem in which I am parsing data from a list.  I am simply trying to 
return data that has several conditions being true.  Here is my syntax below:

d<-sapply(res,function(.df){(.df$TimesVisited[.df$Tick>912 && .df$Id>0])})   
#res is the list, and I am trying to return a result that has two true 
conditions (that is the variable Tick should be greate than 912 and Id should 
be greater than 0

This returns an array of 10 with integer values of 0. This is the incorrect 
result

However, if I do the same syntax except I remove the && statement (i.e. the 
second conditional), then the result producing something that makes sense, 
which is all values that are Tick and greater than 912.

Can someone let me know how I can setup my data to be parsed so I can have 2 or 
multiple conditionals in my function that is looking at an array.

Thanks in advance.

Mark

__
R-help@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@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.

Re: [R] conditional with and operators

[Christos Hatzis]

Mark,

It is not clear whether you are dealing with a single list or something more
complex, such as a list of lists or a list of data frames.  In the case of a
single list, you don't really need any of the 'apply' functions.  The main
problem in your code is the use of '&&' instead of '&':

> test.1 <- list(id = 1:50, tick = rnorm(50, 900, 50), tvis = sample(1:100,
50))
> test.1
$id
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
25 26 27 28 29
[30] 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

$tick
 [1]  891.5450  997.4262  884.0224  892.8317  903.7194  862.7283  913.8126
895.6072
 [9]  923.2599  823.0804  828.0575  919.2257  879.8593  833.6082  924.6253
889.6797
[17]  880.2397 1022.9587  898.3806  953.0213  911.8132  932.3018  960.7257
905.0871
[25] 1018.8407  991.3910  931.3846  960.0569  793.5923  899.3268  921.0783
929.6885
[33]  940.9550  882.9791  914.9050  897.7252  970.7575  867.9848  901.2766
912.0851
[41]  856.8671  878.3230  906.9869  903.8044  882.1902  972.3030  923.9107
869.6903
[49]  836.6934  912.6101

$tvis
 [1] 67  8 89 36 91 30 63 38 22 44 18 79 48 90 82  6 81 39 69  7 45 86 32 43
99 23 34 54 37
[30]  1 15  5 14 87 77 27 93 94 71 60 83 62 92 72 26 80 58 98 84 78

> test.1$id > 10 & test.1$tick > 910
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
FALSE FALSE
[15]  TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE
TRUE  TRUE
[29] FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE FALSE  TRUE
FALSE FALSE
[43] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE

> test.1$tvis[test.1$id > 10 & test.1$tick > 910]
 [1] 79 82 39  7 45 86 32 99 23 34 54 15  5 14 77 93 60 80 58 78 

If this is not what you have, it would be helpful to provide an example of
the data structure that you are dealing with, as the posting guide suggests.

-Christos

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Altaweel, Mark R.
> Sent: Tuesday, August 19, 2008 4:10 PM
> To: r-help@r-project.org
> Subject: [R] conditional with and operators
> 
> Hi,
> 
> I have a problem in which I am parsing data from a list.  I 
> am simply trying to return data that has several conditions 
> being true.  Here is my syntax below:
> 
> d<-sapply(res,function(.df){(.df$TimesVisited[.df$Tick>912 && 
> .df$Id>0])})   #res is the list, and I am trying to return a 
> result that has two true conditions (that is the variable 
> Tick should be greate than 912 and Id should be greater than 0
> 
> This returns an array of 10 with integer values of 0. This is 
> the incorrect result
> 
> However, if I do the same syntax except I remove the && 
> statement (i.e. the second conditional), then the result 
> producing something that makes sense, which is all values 
> that are Tick and greater than 912.
> 
> Can someone let me know how I can setup my data to be parsed 
> so I can have 2 or multiple conditionals in my function that 
> is looking at an array.
> 
> Thanks in advance.
> 
> Mark
> 
> __
> R-help@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@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.

Re: [R] conditional with and operators

[Altaweel, Mark R.]

Hi,

Ok that worked...im from a Java background so I guess Im used to the && rather 
than &

Thanks!

Mark


-Original Message-
From: Steven McKinney [mailto:[EMAIL PROTECTED]
Sent: Tue 8/19/2008 3:33 PM
To: Altaweel, Mark R.; r-help@r-project.org
Subject: RE: [R] conditional with and operators
 


Did you try it with the vector '&' and operator?

d<-sapply(res,function(.df){(.df$TimesVisited[.df$Tick>912 & .df$Id>0])}) 

(The '&&' operator is designed for use in
e.g. if() clauses where you want a scalar logical
answer)

HTH

Steve McKinney

-Original Message-
From: [EMAIL PROTECTED] on behalf of Altaweel, Mark R.
Sent: Tue 8/19/2008 1:10 PM
To: r-help@r-project.org
Subject: [R] conditional with and operators
 
Hi,

I have a problem in which I am parsing data from a list.  I am simply trying to 
return data that has several conditions being true.  Here is my syntax below:

d<-sapply(res,function(.df){(.df$TimesVisited[.df$Tick>912 && .df$Id>0])})   
#res is the list, and I am trying to return a result that has two true 
conditions (that is the variable Tick should be greate than 912 and Id should 
be greater than 0

This returns an array of 10 with integer values of 0. This is the incorrect 
result

However, if I do the same syntax except I remove the && statement (i.e. the 
second conditional), then the result producing something that makes sense, 
which is all values that are Tick and greater than 912.

Can someone let me know how I can setup my data to be parsed so I can have 2 or 
multiple conditionals in my function that is looking at an array.

Thanks in advance.

Mark

__
R-help@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@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.

Re: [R] Survey Design / Rake questions

[Farley, Robert]

While I'm trying to catch up on the statistical basis of my task, could
someone point me to how I should fix my R error? 


Thanks




> library(survey)
> SurveyData <- read.spss("C:/Data/R/orange_delivery.sav",
use.value.labels=TRUE, max.value.labels=Inf, to.data.frame=TRUE)
>
#===

> temp <- sub(' +$', '', SurveyData$direction_) 
> SurveyData$direction_ <- temp
>
#===

>
SurveyData$NumStn=abs(as.numeric(SurveyData$lineon)-as.numeric(SurveyDat
a$lineoff))
> EBSurvey <- subset(SurveyData, direction_ == "EASTBOUND" )
> XTTable <- xtabs(~direction_ , EBSurvey)
> XTTable
direction_
EASTBOUND 
  345 
> WBSurvey <- subset(SurveyData, direction_ == "WESTBOUND" )
> XTTable <- xtabs(~direction_ , WBSurvey)
> XTTable
direction_
WESTBOUND 
  307 
> #
> EBDesign <- svydesign(id=~sampn, weights=~expwgt, data=EBSurvey)
> #   svytable(~lineon+lineoff, EBDesign)
> OnLabels<- c( "Warner Center", "De Soto", "Pierce College",
"Tampa", "Reseda", "Balboa", "Woodley", "Sepulveda", "Van Nuys",
"Woodman", "Valley College", "Laurel Canyon", "North Hollywood")
> EBOnNewTots <- c(1000,   600, 1200,
500, 1000,  500,   200, 250,   1000,   300,
100,  50,73.65 )
> EBNumStn <- c(673.65, 800, 1000, 1000,  800,  700,  600, 500, 400,
200,  50, 50 )
> ByEBOn <- data.frame(OnLabels,EBOnNewTots)
> ByEBNum <- data.frame(c(1:12),EBNumStn)
> RakedEBSurvey <- rake(EBDesign, list(~ByEBOn, ~ByEBNum),
list(EBOnNewTots, EBNumStn ) )
Error in model.frame.default(margin, data = design$variables) : 
  invalid type (list) for variable 'ByEBOn'
>

###
sessionInfo()
R version 2.7.1 (2008-06-23) 
i386-pc-mingw32 

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] graphics  grDevices utils datasets  stats methods   base


other attached packages:
[1] survey_3.8 fortunes_1.3-5 moonsun_0.1prettyR_1.3-2
foreign_0.8-28
>


Robert Farley
Metro
www.Metro.net 

-Original Message-
From: Farley, Robert 
Sent: Monday, August 18, 2008 16:18
To: 'r-help@r-project.org'
Subject: RE: [R] Survey Design / Rake questions

Thank you for the list of references.  Do you know of any "free"
references available online?  I'll have to find my library card :-)


My motivation is to try to correct for a "time on board" bias we see in
our surveys.  Not surprisingly, riders who are only on board a short
time don't attempt/finish our survey forms.  We're able to weight our
survey to the "bus stop-on by bus run" level.  I want to keep that, and
rake on new (imposed?) marginals, like an estimate of how many minutes
they were on-board derived from their origin-destination.  In practice,
we'll have thousands of observations on hundreds of runs.  As I see it,
my work-plan involves: 

Running rake successfully on test data
Preparing "bus stop-on by run" marginals automatically
Plus any other "pre-existing" marginals to be kept.
Appending "time on bus" estimates
Determining the "time on bus" distribution (second survey?)
Implementing the raking adjustment for a production (large)
dataset


As of yet, I cannot get the first step to work  :-(



I hope there are no "fatal flaws" in this concept





Robert Farley
Metro
www.Metro.net 

-Original Message-
From: Stas Kolenikov [mailto:[EMAIL PROTECTED] 
Sent: Monday, August 18, 2008 10:32
To: Farley, Robert
Cc: r-help@r-project.org
Subject: Re: [R] Survey Design / Rake questions

Your reading, in increasing order of difficulty/mathematical details,
might be Lohr's "Sampling"
(http://www.citeulike.org/user/ctacmo/article/1068825), Korn &
Graubard's "Health Surveys"
(http://www.citeulike.org/user/ctacmo/article/553280), and Sarndal et.
al. Survey Math Bible
(http://www.citeulike.org/user/ctacmo/article/716032). You certainly
should try to get a hold of the primary concepts before collecting
your data (or rather before designing your survey... so it might
already be too late!). Post-stratification is not that huge topic, for
some reason; a review of mathematical details is given by Valliant
(1993) (http://www.citeulike.org/user/ctacmo/article/1036976). On
raking, the paper on top of Google Scholar search by Deville, Sarndal
and Sautory (1993)
(http://www.citeulike.org/user/ctacmo/article/3134001) is certainly
coming from the best people in the field.

I am not aware of general treatment of transportation survey sampling,
although I suspect such references do exist in tr

[R] address (nil), cause 'memory not mapped'

[Juan Manuel Barreneche]

Dear users,
I got this problem and i don't have a clue of what it could be happening...

The context: i'm running a loop in which i extract information from a raster
map (I work with GRASS and R, using spgrass6 package), and rearrange it to
create a matrix. I've tried it with small maps and it works smoothly, but in
large maps I have to resort to saving data in the hard disk regularly in
order to avoid the "cannot allocate vector of size..." problem (builted in
the loop itself).

Even so, the message "address (nil), cause 'memory not mapped'" appeared,
and i really don't have a clue of what it means.

I run R 2.7.0 on openSUSE 10.2

If anyone can be of any help, it'll be greatly appreciated... below I'm
writing the output i got (and the R.Version() output...)

thank you,

Juan


Null cells excluded from cost evaluation.
Source map is: Floating point (double) cell type
189 rows, 178 cols
Creating some temporary files
Reading riosypondas3000
 100%
Initializing output
 100%
Finding cost path
Writing nombre
Double cell type.
Writing...
 100%
Peak cost value: 13.614127
Creating BIL support files...
Header File =
/home/mastermind/GRASSDATA//yacare1/jmb/.tmp/SpongeBob/nombre.hdr
World File =
/home/mastermind/GRASSDATA//yacare1/jmb/.tmp/SpongeBob/nombre.wld
Exporting raster as double values (bytes=8)
Using the current region settings...
north=6681799.617727
south=6115500.667945
east=756812.903880
west=223200.856418
r=189
c=178
 100%


 *** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 1: system("g.proj -j -f", intern = TRUE, ignore.stderr = ignore.stderr)
 2: ifelse(.Platform$OS.type == "windows", projstr <- system("g.proj -j
-f", intern = TRUE), projstr <- system("g.proj -j -f", intern =
TRUE, ignore.stderr = ignore.stderr))
 3: getLocationProj()
 4: nchar(projargs)
 5: CRS(getLocationProj())
 6: readRAST6("nombre")
 7: dist.grass(matrizPaisaje.riosypondas3000, "riosypondas3000", direct
= "/home/mastermind/r/distancias", division = 17)
 8: eval(expr, envir, enclos)
 9: eval(parse(text = sprintf("dist.%s <<- dist.grass(matrizPaisaje.%s,
'%s', direct='%s', division=%s)", i, i, i, direct, division)))
10: dist.grass.loop(c("riosypondas3000", "maxent3000", "alturas3000.2"),
direct = "/home/mastermind/r/distancias", division = 17)
11: eval(expr, envir, enclos)
12: eval(expr, envir = loc.frame)
13: system.time(dist.grass.loop(c("riosypondas3000", "maxent3000",
"alturas3000.2"), direct = "/home/mastermind/r/distancias", division =
17))

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace

> R.Version()
$platform
[1] "i686-pc-linux-gnu"

$arch
[1] "i686"

$os
[1] "linux-gnu"

$system
[1] "i686, linux-gnu"

$status
[1] ""

$major
[1] "2"

$minor
[1] "7.0"

$year
[1] "2008"

$month
[1] "04"

$day
[1] "22"

$`svn rev`
[1] "45424"

$language
[1] "R"

$version.string
[1] "R version 2.7.0 (2008-04-22)"

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] address (nil), cause 'memory not mapped'

[Juan Manuel Barreneche]

just to mention that the loop was already running for about 20 hours...

Juan

On Tue, Aug 19, 2008 at 6:31 PM, Juan Manuel Barreneche
<[EMAIL PROTECTED]>wrote:

> Dear users,
> I got this problem and i don't have a clue of what it could be happening...
>
> The context: i'm running a loop in which i extract information from a
> raster map (I work with GRASS and R, using spgrass6 package), and rearrange
> it to create a matrix. I've tried it with small maps and it works smoothly,
> but in large maps I have to resort to saving data in the hard disk regularly
> in order to avoid the "cannot allocate vector of size..." problem (builted
> in the loop itself).
>
> Even so, the message "address (nil), cause 'memory not mapped'" appeared,
> and i really don't have a clue of what it means.
>
> I run R 2.7.0 on openSUSE 10.2
>
> If anyone can be of any help, it'll be greatly appreciated... below I'm
> writing the output i got (and the R.Version() output...)
>
> thank you,
>
> Juan
>
> 
> Null cells excluded from cost evaluation.
> Source map is: Floating point (double) cell type
> 189 rows, 178 cols
> Creating some temporary files
> Reading riosypondas3000
>  100%
> Initializing output
>  100%
> Finding cost path
> Writing nombre
> Double cell type.
> Writing...
>  100%
> Peak cost value: 13.614127
> Creating BIL support files...
> Header File =
> /home/mastermind/GRASSDATA//yacare1/jmb/.tmp/SpongeBob/nombre.hdr
> World File =
> /home/mastermind/GRASSDATA//yacare1/jmb/.tmp/SpongeBob/nombre.wld
> Exporting raster as double values (bytes=8)
> Using the current region settings...
> north=6681799.617727
> south=6115500.667945
> east=756812.903880
> west=223200.856418
> r=189
> c=178
>  100%
> 
>
>  *** caught segfault ***
> address (nil), cause 'memory not mapped'
>
> Traceback:
>  1: system("g.proj -j -f", intern = TRUE, ignore.stderr = ignore.stderr)
>  2: ifelse(.Platform$OS.type == "windows", projstr <- system("g.proj -j
> -f", intern = TRUE), projstr <- system("g.proj -j -f", intern =
> TRUE, ignore.stderr = ignore.stderr))
>  3: getLocationProj()
>  4: nchar(projargs)
>  5: CRS(getLocationProj())
>  6: readRAST6("nombre")
>  7: dist.grass(matrizPaisaje.riosypondas3000, "riosypondas3000", direct
> = "/home/mastermind/r/distancias", division = 17)
>  8: eval(expr, envir, enclos)
>  9: eval(parse(text = sprintf("dist.%s <<- dist.grass(matrizPaisaje.%s,
> '%s', direct='%s', division=%s)", i, i, i, direct, division)))
> 10: dist.grass.loop(c("riosypondas3000", "maxent3000",
> "alturas3000.2"), direct = "/home/mastermind/r/distancias", division =
> 17)
> 11: eval(expr, envir, enclos)
> 12: eval(expr, envir = loc.frame)
> 13: system.time(dist.grass.loop(c("riosypondas3000", "maxent3000",
> "alturas3000.2"), direct = "/home/mastermind/r/distancias", division =
> 17))
>
> Possible actions:
> 1: abort (with core dump, if enabled)
> 2: normal R exit
> 3: exit R without saving workspace
> 4: exit R saving workspace
>
> > R.Version()
> $platform
> [1] "i686-pc-linux-gnu"
>
> $arch
> [1] "i686"
>
> $os
> [1] "linux-gnu"
>
> $system
> [1] "i686, linux-gnu"
>
> $status
> [1] ""
>
> $major
> [1] "2"
>
> $minor
> [1] "7.0"
>
> $year
> [1] "2008"
>
> $month
> [1] "04"
>
> $day
> [1] "22"
>
> $`svn rev`
> [1] "45424"
>
> $language
> [1] "R"
>
> $version.string
> [1] "R version 2.7.0 (2008-04-22)"
>
>

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Multiple Plotting help (lines don't always connect)

[Carl Witthoft]

Just my view on the lines/interpolate/don't interpolate controversy here:

I was taught once upon a time to plot actual data as points (with or 
without error bars), and plot fitted curves as lines.
Only when the density of data points is so high as to make a 'mess' of 
the chart is it ok to plot a line representing the actual data.


YMMV


Carl

__
R-help@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] matrix question

[Altaweel, Mark R.]

Hi,

I have a vector and a list, with data I would like to multiply together.

So for instance I have a vector s:

[[1]]
[1] 44308

[[2]]
[1] 4371

Also, I have a list d:

[[1]]
 [1] 1201 6170 2036 2927 1625 1391 2074 1453 3172 3027 4691 3719 1185  320 2071 
1027 1046 1186 1403  580 1382  4408  174

[[2]]
 [1] 6521  688 2678 3409 3033 1608 3618 1461 1836 2104  879 1095 2630 1591 2986 
 703 2548  913 1426  753  256  869  106


I want to multiply them together and put the results in a matrix.

This is my syntax:

for(i in 1:length(s))
 for(j in 1:length(d))
   m<-d[[j]][j]/s[[i]]  #m is the matrix of dimensions set to X(e.g. 10X10)

However, it seems I only get one result when i get m, which is the last value 
of d/j, which is m[1]=0.04 in this case.

I am sure Im doing something wrong here, but can't quite find the solution.

Thanks.

Mark  

__
R-help@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.

Re: [R] converting coordinates from utm to longitude / latitude

[Werner Wernersen]

It would be nicer to convert directly the entire shapefile object to long/lat 
coordinates but if that is not possible, I will convert the other points to UTM.
Hence, I am playing around with rgdal.

library(rgdal)
SP <- SpatialPoints(cbind(32.29252, -0.3228500), 
  proj4string=CRS("+proj=longlat"))
spTransform(SP, CRS("+proj=utm +zone=36"))


> spTransform(SP, CRS("+proj=utm +zone=36"))
SpatialPoints:
 coords.x1 coords.x2
[1,]  421274.4 -35687.37
Coordinate Reference System (CRS) arguments: +proj=utm +zone=36 +ellps=WGS84 

This result corresponds with what I get when using convUL() but my map of that 
area in UTM coordinates does not extend to the negative.
An external program converts the point to x=420994   y=9964407 which also seems 
correct with respect to the map. Fore sure, I am using the function wrongly 
somehow. Can anyone give me a hint?

That's very much appreciated!

Thanks,
   Werner



- Ursprüngliche Mail 
Von: Werner Wernersen <[EMAIL PROTECTED]>
An: [EMAIL PROTECTED]
Gesendet: Dienstag, den 19. August 2008, 20:28:29 Uhr
Betreff: converting coordinates from utm to longitude / latitude

Hi,

is there a function in R to convert data read with read.shape and which is 
originally in UTM coordinates into longitude / latitude coordinates?
I found the convUL() function from the PBSmapping package but I have no idea 
how I could apply that to the read.shape object.

Many thanks,
  Werner


__
fügt über einen herausragenden Schutz gegen Massenmails. 
http://mail.yahoo.com


__
Do 
ragenden Schutz gegen Massenmails. 
http://mail.yahoo.com

__
R-help@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.

Re: [R] converting coordinates from utm to longitude / latitude

[Dylan Beaudette]

If you would like to convert the entire shapfile check out GDAL:

http://www.gdal.org/

http://casoilresource.lawr.ucdavis.edu/drupal/node/98

Cheers,

Dylan


On Tuesday 19 August 2008, Werner Wernersen wrote:
> It would be nicer to convert directly the entire shapefile object to
> long/lat coordinates but if that is not possible, I will convert the other
> points to UTM. Hence, I am playing around with rgdal.
>
> library(rgdal)
> SP <- SpatialPoints(cbind(32.29252, -0.3228500),
>   proj4string=CRS("+proj=longlat"))
> spTransform(SP, CRS("+proj=utm +zone=36"))
>
> > spTransform(SP, CRS("+proj=utm +zone=36"))
>
> SpatialPoints:
>  coords.x1 coords.x2
> [1,]  421274.4 -35687.37
> Coordinate Reference System (CRS) arguments: +proj=utm +zone=36
> +ellps=WGS84
>
> This result corresponds with what I get when using convUL() but my map of
> that area in UTM coordinates does not extend to the negative. An external
> program converts the point to x=420994   y=9964407 which also seems correct
> with respect to the map. Fore sure, I am using the function wrongly
> somehow. Can anyone give me a hint?
>
> That's very much appreciated!
>
> Thanks,
>Werner
>
>
>
> - Ursprüngliche Mail 
> Von: Werner Wernersen <[EMAIL PROTECTED]>
> An: [EMAIL PROTECTED]
> Gesendet: Dienstag, den 19. August 2008, 20:28:29 Uhr
> Betreff: converting coordinates from utm to longitude / latitude
>
> Hi,
>
> is there a function in R to convert data read with read.shape and which is
> originally in UTM coordinates into longitude / latitude coordinates? I
> found the convUL() function from the PBSmapping package but I have no idea
> how I could apply that to the read.shape object.
>
> Many thanks,
>   Werner
>
>
> __
> fügt über einen herausragenden Schutz gegen Massenmails.
> http://mail.yahoo.com
>
>
> __
> Do
> ragenden Schutz gegen Massenmails.
> http://mail.yahoo.com
>
> __
> R-help@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.



-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341

__
R-help@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.

Re: [R] matrix question

[Rolf Turner]

On 20/08/2008, at 10:08 AM, Altaweel, Mark R. wrote:


Hi,

I have a vector and a list, with data I would like to multiply  
together.


No you haven't.  You have two ***lists***.  Lists and vectors
are not the same thing.  If you don't distinguish between them
you will confuse everyone, including yourself.


So for instance I have a vector s:

[[1]]
[1] 44308

[[2]]
[1] 4371

Also, I have a list d:

[[1]]
 [1] 1201 6170 2036 2927 1625 1391 2074 1453 3172 3027 4691 3719  
1185  320 2071 1027 1046 1186 1403  580 1382  4408  174


[[2]]
 [1] 6521  688 2678 3409 3033 1608 3618 1461 1836 2104  879 1095  
2630 1591 2986  703 2548  913 1426  753  256  869  106



I want to multiply them together and put the results in a matrix.

This is my syntax:

for(i in 1:length(s))
 for(j in 1:length(d))
   m<-d[[j]][j]/s[[i]]  #m is the matrix of dimensions set to X 
(e.g. 10X10)


However, it seems I only get one result when i get m, which is the  
last value of d/j, which is m[1]=0.04 in this case.


I am sure Im doing something wrong here, but can't quite find the  
solution.


	The object m has been repeatedly assigned the value ``d[[j]][j]/s 
[[i]]'' and

so naturally ends up having the last value that it was assigned.

	If you want m to be a matrix, make it one, and assign values to its  
entries!


m <- matrix(NA,length(s),length(d))

for(i in 1:length(s)) {
for(j in 1:length(d)) {
m[i,j] <-d[[j]][j]/s[[i]]
}
}

It's not at all clear if the foregoing is really what you are  
intending to do, since it's
not at all clear what you are intending to do.  Did you really want  
``d[[j]][j]'' (i.e. 2 j-s)???
Are you aware that the two entries of d are of different lengths (24  
and 23 respectively) and
that you only use the first entry of the first and the second entry  
of the second.  Why are you

using lists rather than vectors and matrices?

You could achieve the same result as the convoluted shaganappi  
foregoing code by doing



u <- matrix(c(1201,688),2,2,byrow=TRUE)
v <- c(44308,4371)
m <- u/v

But as I indicated, that's probably not what you really want to do.   
You need to get

your thoughts organized.

cheers,

Rolf Turner

##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

__
R-help@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.

Re: [R] Writing R Extensions : A new R package for Gini Index decomposition to prupose

[Robert Duval]

there are many useful decompositions of the Gini... which one did you implement?

On Tue, Aug 19, 2008 at 8:42 AM, stephen sefick <[EMAIL PROTECTED]> wrote:
> See if there is interest.  If there is not make your own package or
> see if someone else would like to include it into a package that is
> complementary.
>
> Stephen Sefick
>
> On Tue, Aug 19, 2008 at 3:46 AM, Ndoye Souleymane <[EMAIL PROTECTED]> wrote:
>>
>> Dear All,
>>
>> I have developed a programme the anable the decomposition of the Gini index, 
>> it complets tha valuable work of Achim Zeileis, the author of the ineq 
>> package.
>> I would like to make it to be part of all R package. How should I proceed.
>> Must I sent it to the the Core developement team ?
>> The proogramme is written in R.
>>
>> Many thanks for your advice,
>>
>> Best regards,
>>
>> Souleymane
>> _
>> Retouchez, classez et partagez vos photos gratuitement avec le logiciel 
>> Galerie de Photos !
>>
>> __
>> R-help@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.
>>
>
>
>
> --
> Let's not spend our time and resources thinking about things that are
> so little or so large that all they really do for us is puff us up and
> make us feel like gods. We are mammals, and have not exhausted the
> annoying little problems of being mammals.
>
>-K. Mullis
>
> __
> R-help@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@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] bug in lme4?

[Antonio.Gasparrini]

Dear all,
 
I found a problem with 'lme4'. Basically, once you load the package 'aod' 
(Analysis of Overdispersed Data), the functions 'lmer' and 'glmer' don't work 
anymore:
 
library(lme4)
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
  family = binomial, data = cbpp))
install.packages("aod")
library(aod)
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
  family = binomial, data = cbpp))
 
Taking into account that this package is used to perform similar analyses, this 
could be a problem.
 
All the best
 
Antonio Gasparrini
Public and Environmental Health Research Unit (PEHRU)
London School of Hygiene & Tropical Medicine
Keppel Street, London WC1E 7HT, UK
Office: 0044 (0)20 79272406 - Mobile: 0044 (0)79 64925523
http://www.lshtm.ac.uk/people/gasparrini.antonio ( 
http://www.lshtm.ac.uk/pehru/ )
--- Begin Message ---
Dear all,
 
I found a problem with 'lme4'. Basically, once you load the package 'aod' 
(Analysis of Overdispersed Data), the functions 'lmer' and 'glmer' don't work 
anymore:
 
library(lme4)
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
  family = binomial, data = cbpp))
install.packages("aod")
library(aod)
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
  family = binomial, data = cbpp))

Taking into account that this package is used to perform similar analyses, this 
could be a problem.

All the best

--- End Message ---
__
R-help@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.

Re: [R] converting coordinates from utm to longitude / latitude

[Jim Regetz]

Werner Wernersen wrote:

It would be nicer to convert directly the entire shapefile object to long/lat 
coordinates but if that is not possible, I will convert the other points to UTM.
Hence, I am playing around with rgdal.

library(rgdal)
SP <- SpatialPoints(cbind(32.29252, -0.3228500), 
  proj4string=CRS("+proj=longlat"))

spTransform(SP, CRS("+proj=utm +zone=36"))



spTransform(SP, CRS("+proj=utm +zone=36"))

SpatialPoints:
 coords.x1 coords.x2
[1,]  421274.4 -35687.37
Coordinate Reference System (CRS) arguments: +proj=utm +zone=36 +ellps=WGS84 


This result corresponds with what I get when using convUL() but my map of that 
area in UTM coordinates does not extend to the negative.
An external program converts the point to x=420994   y=9964407 which also seems 
correct with respect to the map. Fore sure, I am using the function wrongly 
somehow. Can anyone give me a hint?


It looks like you are specifying 36S in your external program, and 
(implicitly) 36N in R. Using your SP from above, note the following:


> spTransform(SP, CRS("+proj=utm +zone=36 +north"))
SpatialPoints:
 coords.x1 coords.x2
[1,]  421274.4 -35687.37
Coordinate Reference System (CRS) arguments: +proj=utm +zone=36 
+ellps=WGS84


> spTransform(SP, CRS("+proj=utm +zone=36 +south"))
SpatialPoints:
 coords.x1 coords.x2
[1,]  421274.4   9964313
Coordinate Reference System (CRS) arguments: +proj=utm +zone=36 +south 
+ellps=WGS84


The latter gets in the ballpark of output from your external program. 
I'd speculate that the remaining differences owe to different 
assumptions about e.g. the datum (WGS84 for the R statements given 
above), but I must admit I'm not a geographer.


Hope that helps,
Jim


That's very much appreciated!

Thanks,
   Werner



- Ursprüngliche Mail 
Von: Werner Wernersen <[EMAIL PROTECTED]>
An: [EMAIL PROTECTED]
Gesendet: Dienstag, den 19. August 2008, 20:28:29 Uhr
Betreff: converting coordinates from utm to longitude / latitude

Hi,

is there a function in R to convert data read with read.shape and which is 
originally in UTM coordinates into longitude / latitude coordinates?
I found the convUL() function from the PBSmapping package but I have no idea 
how I could apply that to the read.shape object.

Many thanks,
  Werner


__
fügt über einen herausragenden Schutz gegen Massenmails. 
http://mail.yahoo.com



__
Do 
ragenden Schutz gegen Massenmails. 
http://mail.yahoo.com


__
R-help@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@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] Confidence Interval

[Raphael Saldanha]

Hi!

With the following script, I'm trying to make a demonstration of a
Confidence Interval, but I'm observing some differences on tails.

# Teste de média entre uma amostra e uma população normal
# Autor: Raphael de Freitas Saldanha
# Agosto de 2008

n<- 200# Sample size
xbar <- 100# Sample mean
s<- 2  # Sample SD
nc   <- 0.95   # Confidence level (95% -> 0.95)
rep  <- 1000   # Loops

###

y <- NULL# Vetor com as médias da amostra
for (i in 1:rep){# Loop
x <- rnorm(n,xbar,s) # Gere uma amostra normal n elementos, xbar média e
s desvio-padrão
x <- mean(x) # Calcule a média (exata) dessa amostra
y <- c(y,x)  # Coloque essa média em um registro em y
}

y <- sort(y) # Ordene todas as médias geradas

LI <- y[((1-nc)/2)*rep] # Limite inferior,
LS <- y[rep-(((1-nc)/2)*rep)]   # Limite superior

### PARTE GRÁFICA ###

x <- mean(y)

xvals <- seq(-LI, LI, length.out=5000)
dvals <- dnorm(xvals,mean(y), sd(y))[1:5000]

xbvals <- seq(LS, LS*2, length.out=5000)
dbvals <- dnorm(xbvals,mean(y), sd(y))[1:5000]

ahist <- hist(y, freq=FALSE, col="lightblue", main="Intervalo de confiança")

polygon(c(xvals,LI,LI), c(dvals,dnorm(LI,mean(y), sd(y)),min(dvals)),
col="orange", lty=0)
polygon(c(LS,LS,xbvals), c(min(dbvals),dnorm(LS,mean(y), sd(y)),dbvals),
col="orange", lty=0)
curve(dnorm(x,mean(y), sd(y)),add=TRUE, lty=1, col="red",lwd=2)

### Intervalo de Confiança ###

LI # Limite inferior
LS # Limite superior

-- 
Atenciosamente,

Raphael Saldanha
[EMAIL PROTECTED]

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Confidence Interval

[Prof Brian Ripley]

On Tue, 19 Aug 2008, Raphael Saldanha wrote:


Hi!

With the following script, I'm trying to make a demonstration of a
Confidence Interval, but I'm observing some differences on tails.


You need to tell us what you are trying to do.  You seem to be computing a 
parametric bootstrap interval, but incorrectly (you need to reflect 
percentile intervals).  See Davison & Hinkley (1997) 'The Bootstrap and 
its Application' for more details.


In any case, your simulation is not repeatable (you set no seed), so we 
don't know what you saw and what 'differences' disturbed you.


Your calculation of quantiles is not quite correct: use quantile().
(The indices go from 1 to 1000, not 0 to 1000.)  E.g.


quantile(1:1000, c(0.025, 0.975))

   2.5%   97.5%
 25.975 975.025

and not 25, 975.

As your results show, a histogram is not a good summary of the shape of a 
distribution on 1000 points.  We can do much better with an ecdf or a 
density estimate.




# Teste de m?dia entre uma amostra e uma popula??o normal
# Autor: Raphael de Freitas Saldanha
# Agosto de 2008

n<- 200# Sample size
xbar <- 100# Sample mean
s<- 2  # Sample SD
nc   <- 0.95   # Confidence level (95% -> 0.95)
rep  <- 1000   # Loops

###

y <- NULL# Vetor com as m?dias da amostra
for (i in 1:rep){# Loop
x <- rnorm(n,xbar,s) # Gere uma amostra normal n elementos, xbar m?dia e
s desvio-padr?o
x <- mean(x) # Calcule a m?dia (exata) dessa amostra
y <- c(y,x)  # Coloque essa m?dia em um registro em y
}

y <- sort(y) # Ordene todas as m?dias geradas

LI <- y[((1-nc)/2)*rep] # Limite inferior,
LS <- y[rep-(((1-nc)/2)*rep)]   # Limite superior

### PARTE GR?FICA ###

x <- mean(y)

xvals <- seq(-LI, LI, length.out=5000)
dvals <- dnorm(xvals,mean(y), sd(y))[1:5000]

xbvals <- seq(LS, LS*2, length.out=5000)
dbvals <- dnorm(xbvals,mean(y), sd(y))[1:5000]

ahist <- hist(y, freq=FALSE, col="lightblue", main="Intervalo de confian?a")

polygon(c(xvals,LI,LI), c(dvals,dnorm(LI,mean(y), sd(y)),min(dvals)),
col="orange", lty=0)
polygon(c(LS,LS,xbvals), c(min(dbvals),dnorm(LS,mean(y), sd(y)),dbvals),
col="orange", lty=0)
curve(dnorm(x,mean(y), sd(y)),add=TRUE, lty=1, col="red",lwd=2)

### Intervalo de Confian?a ###

LI # Limite inferior
LS # Limite superior

--
Atenciosamente,

Raphael Saldanha
[EMAIL PROTECTED]

[[alternative HTML version deleted]]




--
Brian D. Ripley,  [EMAIL PROTECTED]
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, UKFax:  +44 1865 272595

__
R-help@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] R CMD INSTALL works but not install.packages ...

[Tae-Hoon Chung]

Hi, All;

I tried to install exon.pmcdf_1.1.tar.gz package using install.packages()
but it didn't with the following messages:

> install.packages('/Users/thchung/Downloads/exon.pmcdf_1.1.tar.gz',
repos=NULL)
Error in gzfile(file, "r") : cannot open the connection
In addition: Warning message:
In gzfile(file, "r") :
  cannot open compressed file 'exon.pmcdf_1.1.tar.gz/DESCRIPTION', probable
reason 'No such file or directory'
> sessionInfo()
R version 2.7.2 beta (2008-08-16 r46368)
i386-apple-darwin9.4.0

locale:
C

attached base packages:
[1] splines   tools stats graphics  grDevices utils datasets
[8] methods   base

other attached packages:
 [1] exonmap_1.6.05   plier_1.10.0 RMySQL_0.6-0
 [4] DBI_0.2-4RColorBrewer_1.0-2   simpleaffy_2.16.0
 [7] gcrma_2.12.1 matchprobes_1.12.0   genefilter_1.20.0
[10] survival_2.34-1  affy_1.18.2  preprocessCore_1.2.1
[13] affyio_1.8.0 Biobase_2.0.1

loaded via a namespace (and not attached):
[1] AnnotationDbi_1.2.2 RSQLite_0.6-9   annotate_1.18.0

Surprisingly, the R CMD INSTALL mechanism worked:

BiobankMacBookPro:Downloads thchung$ R CMD INSTALL exon.pmcdf_1.1.tar.gz
* Installing to library '/Library/Frameworks/R.framework/Resources/library'
* Installing *source* package 'exon.pmcdf' ...
** R
** data
** help
 >>> Building/Updating help pages for package 'exon.pmcdf'
 Formats: text html latex example
  exon.pmcdftexthtmllatex
  geometry  texthtmllatex   example
** building package indices ...
* DONE (exon.pmcdf)

Can anyone tell me what the problem is? Thanks in advance,

Tae-Hoon

-- 
Tae-Hoon Chung

Korea Centers for Disease Control & Prevention (KCDC)
Korea National Institute of Health
Center for Genome Sciences :: Biobank Team

(T) +82-2-380-2252
(F) +82-2-354-1078
(M) +82-10-3511-1012

[[alternative HTML version deleted]]

__
R-help@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] Conversion - lowercase to Uppercase letters

[suman Duvvuru]

I would like to know how to convert a string with characters to all
uppercase or all lowercase? If anyone could let me know if there exists a
function in R for the conversion, that will be very helpful.

Regards,
Suman

[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Conversion - lowercase to Uppercase letters

[Rolf Turner]

On 20/08/2008, at 4:19 PM, suman Duvvuru wrote:


I would like to know how to convert a string with characters to all
uppercase or all lowercase? If anyone could let me know if there  
exists a

function in R for the conversion, that will be very helpful.


?tolower

##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

__
R-help@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.

Re: [R] Conversion - lowercase to Uppercase letters

[Moshe Olshansky]

Use toupper or tolower (see ?toupper, ?tolower)


--- On Wed, 20/8/08, suman Duvvuru <[EMAIL PROTECTED]> wrote:

> From: suman Duvvuru <[EMAIL PROTECTED]>
> Subject: [R] Conversion - lowercase to Uppercase letters
> To: r-help@r-project.org
> Received: Wednesday, 20 August, 2008, 2:19 PM
> I would like to know how to convert a string with characters
> to all
> uppercase or all lowercase? If anyone could let me know if
> there exists a
> function in R for the conversion, that will be very
> helpful.
> 
> Regards,
> Suman
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@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@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] Random sequence of days?

[Lauri Nikkinen]

Dear list,

I tried to find a solution for this problem from the archives but
couldn't find any. I would like sample sequence of ten days from
vector d

d <- seq(as.Date("2007-02-12"), as.Date("2008-08-18"), by="days")

so that the days follow each other (sample(d, 10) is not the
appropriate solution). Any ideas?

Best regards,
Lauri

__
R-help@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.

Re: [R] Random sequence of days?

[Moshe Olshansky]

How about

d[sample(length(d),10)]


--- On Wed, 20/8/08, Lauri Nikkinen <[EMAIL PROTECTED]> wrote:

> From: Lauri Nikkinen <[EMAIL PROTECTED]>
> Subject: [R] Random sequence of days?
> To: [EMAIL PROTECTED]
> Received: Wednesday, 20 August, 2008, 4:04 PM
> Dear list,
> 
> I tried to find a solution for this problem from the
> archives but
> couldn't find any. I would like sample sequence of ten
> days from
> vector d
> 
> d <- seq(as.Date("2007-02-12"),
> as.Date("2008-08-18"), by="days")
> 
> so that the days follow each other (sample(d, 10) is not
> the
> appropriate solution). Any ideas?
> 
> Best regards,
> Lauri
> 
> __
> R-help@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@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.

Re: [R] Random sequence of days?

[Jeff Newmiller]

Lauri Nikkinen wrote:

Dear list,

I tried to find a solution for this problem from the archives but
couldn't find any. I would like sample sequence of ten days from
vector d

d <- seq(as.Date("2007-02-12"), as.Date("2008-08-18"), by="days")

so that the days follow each other (sample(d, 10) is not the
appropriate solution). Any ideas?


"Follow each other" is not clear to me.  Do you want 10 consecutive days?
(Pick a random day between the first and the last less ten, and then
extract the consecutive days.)  Do you want random days, but in
sorted order? (sample ten days and sort the result.)

--
---
Jeff NewmillerThe .   .  Go Live...
DCN:<[EMAIL PROTECTED]>Basics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k

__
R-help@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.

Re: [R] Random sequence of days?

[Lauri Nikkinen]

Thanks for quick reply. I meant that the sequence should be like this:

"2007-06-16" "2007-06-17" "2007-06-18" "2007-06-19" "2007-06-20"
"2007-06-21" "2007-06-22" "2007-06-23"
"2007-06-24" "2007-06-25"

so that the the days are adjacent...

-Lauri

2008/8/20, Moshe Olshansky <[EMAIL PROTECTED]>:
> How about
>
> d[sample(length(d),10)]
>
>
> --- On Wed, 20/8/08, Lauri Nikkinen <[EMAIL PROTECTED]> wrote:
>
> > From: Lauri Nikkinen <[EMAIL PROTECTED]>
> > Subject: [R] Random sequence of days?
> > To: [EMAIL PROTECTED]
> > Received: Wednesday, 20 August, 2008, 4:04 PM
> > Dear list,
> >
> > I tried to find a solution for this problem from the
> > archives but
> > couldn't find any. I would like sample sequence of ten
> > days from
> > vector d
> >
> > d <- seq(as.Date("2007-02-12"),
> > as.Date("2008-08-18"), by="days")
> >
> > so that the days follow each other (sample(d, 10) is not
> > the
> > appropriate solution). Any ideas?
> >
> > Best regards,
> > Lauri
> >
> > __
> > R-help@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@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] cmprsk and a time dependent covariate in the model

[phguardiol]

 Dear R users,





I d like to assess the effect of "treatment" covariate on a disease relapse 
risk with the package cmprsk. 

However, the effect of this covariate on survival is time-dependent
(assessed with cox.zph): no significant effect during the first year of 
follow-up,
then after 1 year a favorable effect is observed on survival (step
function might be the correct way to say that ?). 
For overall survival analysis
I have used a time dependent Cox model which has confirmed this positive effect 
after
1 year.


Now I m moving to disease relapse incidence and a similar time dependency seems 
to be present. 


what I d like to have is that: for
patients without "treatment" the code for "treatment" covariate is
always 0, and for patients who received "treatment" covariate I d like
to have it = 0 during time interval 0 to 1 year, and equal to 1 after 1
year. Correct me if I m wrong in trying to do so.




First, I have run the following script (R2.7.1 under XPpro) according to 
previous advices:



library(cmprsk)

attach(LAMrelapse)

fit1<- crr(rel.t, rel.s, treatment, treatment, function(uft)
cbind(ifelse(uft<=1,1,0),ifelse(uft>1,1,0)), failcode=1,
cencode=0, na.action=na.omit, gtol-06, maxiter)

fit1



where:

rel.t = time to event (in years)

rel.s = status , =1 if disease relapse, =2 if death from non disease
related cause (toxicity of previous chemotherapy), =0 if alive &
not in relapse

treatment = binary covariate (value: 0 or 1) representing the treatment to test 
(different from chemotherapy above, with no known toxicity)

I have not yet added other covariates in the model.



this script gave me the following result:

> fit1 <- crr(relcmp.t, relcmp.s, treatment, treatment, function(uft) 
> cbind(ifelse(uft <= 1, 1, 0), ifelse(uft > 1, 1, 0)), failcode = 1, cencode = 
> 0, 

??? na.action = na.omit, gtol = 1e-006, maxiter = 10)

> fit1

convergence:? TRUE 

coefficients:

[1] -0.6808? 0.7508

standard errors:

[1] 0.2881 0.3644

two-sided p-values:

[1] 0.018 0.039


...That I dont understand at all since it looks like if "treatment" covariate 
had also a significant effect of the first period of time !? 
This is absolutely not the case. 
So I m surely wrong with a part of this script... cov2 and tf are pretty 
obscure for me in the help file of the package. I would really appreciate 
advices regarding these 2 "terms". 



I was thinking that I might changed : cbind(ifelse(uft <= 1, 1, 0), ifelse(uft 
> 1, 1, 0) ? ? ? ? ? ? ? ? ? into:??? cbind(ifelse(uft <= 1, 0, 1), 
ifelse(uft > 1, 1, 0)



But since I only have one covariate (treatment) to test, shouldnt I only write 
the following:

fit1<- crr(rel.t,
rel.s, treatment, treatment, function(uft) ifelse(uft<=1,0,1)), failcode=1,
cencode=0, na.action=na.omit, gtol-06, maxiter)


which gives me :

> fit1


convergence:? TRUE 


coefficients:


[1]? 0.06995 -0.75080


standard errors:


[1] 0.2236 0.3644


two-sided p-values:


[1] 0.750 0.039




which, if I understand things
correctly (I m not sure at all !) confirms that before 1 year, the effect of 
"treatment" covariate is not
significant, but is significant after 1 year of follow up. But there I m again 
not sure of the result I obtain...


any help would be greatly appreciated with cov2 and tf



thanks for? if you have some time for this,







Philippe Guardiola



 




[[alternative HTML version deleted]]

__
R-help@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] Understanding output of summary(glm(...))

[Daren Tan]

Simple example of 5 groups of 4 replicates.
 
>set.seed(5)
 
>tmp <- rnorm(20)
 
>gp <- as.factor(rep(1:5,each=4))
 
>summary(glm(tmp ~ -1 + gp, data=data.frame(tmp, gp)))$coefficients  
>Estimate Std. Error   t value  Pr(>|t|)gp1 -0.1604613084  0.4899868 
>-0.3274809061 0.7478301gp2  0.0002487984  0.4899868  0.0005077655 0.9996016gp3 
> 0.0695463698  0.4899868  0.1419352018 0.8890200gp4 -0.6121682841  0.4899868 
>-1.2493567852 0.2306791gp5 -0.6999545014  0.4899868 -1.4285171713 0.1736348
 
>m <- data.frame(tmp, gp)
>sapply(gp, function(x) sd(m[m[,"gp"]==x,1])) [1] 1.169284 1.169284 1.169284 
>1.169284 1.142974 1.142974 1.142974 1.142974 [9] 0.862423 0.862423 0.862423 
>0.862423 0.535740 0.535740 0.535740 0.535740[17] 1.047538 1.047538 1.047538 
>1.047538
Why doesn't the standard deviation of each group correlates with the Pr e.g., 
gp = 4 has the smallest sd of 0.535740, but its Pr is not the lowest (i.e., 
only 0.23 vs 0.1736 of gp = 5). 
 
Another example with new tmp1
 
>tmp1
 [1]  9.577969  9.310792  9.666767  9.610164 10.181692 10.155899 10.025943 [8]  
9.971243 10.177766  9.265793  9.415818 10.099874 10.238829  9.575591[15]  
9.560879  9.617891  9.617891 10.158160 10.592377 10.068443
 
>summary(glm(tmp1 ~ -1 + age, data=data.frame(as.vector(as.matrix(tmp1)), 
>age)))$coefficients  Estimate Std. Error  t value Pr(>|t|)age1  
>9.541423  0.1611603 59.20456 3.380085e-19age2 10.083694  0.1611603 62.56935 
>1.479781e-19age3  9.739813  0.1611603 60.43557 2.485380e-19age4  9.748297  
>0.1611603 60.48821 2.453251e-19age5 10.109218  0.1611603 62.72773 1.424913e-19
m1 <- data.frame(tmp1, gp)
 
>sapply(age, function(x) sd(m1[m1[,"age"]==x,1])) [1] 0.1580745 0.1580745 
>0.1580745 0.1580745 0.1013207 0.1013207 0.1013207 [8] 0.1013207 0.4658736 
>0.4658736 0.4658736 0.4658736 0.3279128 0.3279128[15] 0.3279128 0.3279128 
>0.3995426 0.3995426 0.3995426 0.3995426
 
Can I conclude from the Pr of summary that tmp1 are of better "quality" than 
tmp, given that its Pr. values are signficantly smaller ? 
 
_


[[alternative HTML version deleted]]

__
R-help@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.

Re: [R] Random sequence of days?

[Lauri Nikkinen]

Thanks Jeff for pointing me the right direction. "Consecutive" was the
right word. This will do it

d.samp <- sample(d[1:(length(d)-10)], 1)
seq(d.samp, d.samp+10, by="days")

Regards,
Lauri

2008/8/20, Lauri Nikkinen <[EMAIL PROTECTED]>:
> Thanks for quick reply. I meant that the sequence should be like this:
>
> "2007-06-16" "2007-06-17" "2007-06-18" "2007-06-19" "2007-06-20"
> "2007-06-21" "2007-06-22" "2007-06-23"
> "2007-06-24" "2007-06-25"
>
> so that the the days are adjacent...
>
> -Lauri
>
> 2008/8/20, Moshe Olshansky <[EMAIL PROTECTED]>:
> > How about
> >
> > d[sample(length(d),10)]
> >
> >
> > --- On Wed, 20/8/08, Lauri Nikkinen <[EMAIL PROTECTED]> wrote:
> >
> > > From: Lauri Nikkinen <[EMAIL PROTECTED]>
> > > Subject: [R] Random sequence of days?
> > > To: [EMAIL PROTECTED]
> > > Received: Wednesday, 20 August, 2008, 4:04 PM
> > > Dear list,
> > >
> > > I tried to find a solution for this problem from the
> > > archives but
> > > couldn't find any. I would like sample sequence of ten
> > > days from
> > > vector d
> > >
> > > d <- seq(as.Date("2007-02-12"),
> > > as.Date("2008-08-18"), by="days")
> > >
> > > so that the days follow each other (sample(d, 10) is not
> > > the
> > > appropriate solution). Any ideas?
> > >
> > > Best regards,
> > > Lauri
> > >
> > > __
> > > R-help@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@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] FYI: APL in R

[Jan de Leeuw]

http://idisk.mac.com/jdeleeuw-Public/utilities/apl/apl.R

Dedicated to the IBM 2741.

Implemented for general multidimensional arrays:

drop, take, reshape, shape, rank, select, generalized inner product,  
generalized

outer product, representation, base value, join, expand, reduce, scan,
member of, ravel, compress, tranpose, rotate

Basically, the APL-I part is complete, and after some testing
and clean-up this will be 1.0. In the code there is a small section
explaining the relation between aplTP(), the APL transpose
function, and the (less general but more natural) aperm() from R.

I may decide to add some extensions from APL-X, APL-2, J,
Sharpe APL, APL 2000, Dyalog APL, later on. In any case, the
current code adds a lot of array operations (and even matrix and vector
operations) to R.

Although this is all prefix and no infix, we could easily
recreate some of the infamous APL one-liners that nobody
can possibly understand or reproduce.

===
 Jan de Leeuw, 11667 Steinhoff Rd, Frazier Park, CA 93225
 home 661-245-1725 mobile 661-231-5416 work 310-825-9550
 .mac: jdeleeuw +++  aim: deleeuwjan +++ skype: j_deleeuw
===
 I am I because my little dog knows me.
 Gertrude Stein

__
R-help@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.