Re: [R] help: cannot allocate vector of length 828310236

2006-08-15 Thread Prof Brian Ripley
Does it make any statistical sense to do polr or probit regression (not 
the same thing) on `really huge data'?  There are few regression-like 
problems in which model inadequacy does not swamp estimation uncertainty 
for as few as a 1000 cases.

If you want to do that sort of thing, by all means use SAS to do it.
But if you are not prepared to spend a few $$ on adequate RAM, don't 
expect free technical consultancy, especially not from those whose work 
you are using and not crediting.

- The uncredited author of polr().


On Mon, 14 Aug 2006, T Mu wrote:

 Hi all,
 
 I was trying a probit regression using polr() and got this message,

polr is a strange choice of tool for 'probit regression' as the term is 
usually used.  It does 'ordered probit regression'.

 Error in model.matrix.default(Terms, m, contrasts) :
 cannot allocate vector of length 828310236
 
 The data is about 20M (a few days ago I asked a question about large file,
 thank you for responses, then I use MS Access to select those columns I
 would use).
 
 R is 2.3.1, Windows XP, 512M Ram.
 
 I am going to read some help on memory use in R, but hope anybody can give
 me some quick hints.

Quick hint: read and follow the posting guide BEFORE posting.

 Is it because iphysical memory runs out, or some other things could be wrong
 with data or polr()?
 Does R use virtual memory? If so, what options can I set?
 If not, can R deal with really huge data (except adding RAM according to
 data size)? If this is the case, it is too bad that I have to tell my boss
 to go back to SAS. Now it is not a speed issue yet.
 
 Thank you.
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

-- 
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Help with workaround for: Function '`[`' is not in thederivatives table

2006-08-15 Thread Bill.Venables
Earl F. Glynn asks: 
 -Original Message-
 From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Earl F. Glynn
 Sent: Tuesday, 15 August 2006 8:44 AM
 To: r-help@stat.math.ethz.ch
 Subject: [R] Help with workaround for: Function '`[`' is not in
thederivatives table
 
 # This works fine:
  a - 1
  b - 2
  c - 3
  E - expression(a * exp(b*X) + c)
  X - c(0.5, 1.0, 2.0)
  eval(E)
 [1]  5.718282 10.389056 57.598150
  D(E, b)
 a * (exp(b * X) * X)
  eval(D(E, b))
 [1]   1.359141   7.389056 109.196300
 
 # But if (a,b,c) are replaced with (A[1], A[2], A[3]), how can I get a

 derivative using D?

It's well to note what D can differentiate with respect to.  The
second argument is, to quote the help page, a character string giving
the name of a variable...  'A[1]' is a character string, but it is
not the name of a variable.  When parsed it becomes a call to the
function, literally, `[`.

  A - c(1, 2, 3)
  E - expression(A[1] * exp(A[2]*X) + A[3])
  X - c(0.5, 1.0, 2.0)
  eval(E)
 [1]  5.718282 10.389056 57.598150

No problem, because the evaluator does know how to evaluate `[`, along
with everything else in this expr.

 
 
 
 # Why doesn't this work?  Any workarounds?
  D(E, A[2])
 Error in D(E, A[2]) : Function '`[`' is not in the derivatives table

It fails because 'A[2]' is not a (character string giving the) name of
a variable.  I think the error message could be a bit more
informative, but ... all you need to do is read he help page, really.

  If I want to have a long vector of coefficients, A, (perhaps
 dozens) how can I use D to compute partial derivatives?

Essentially you need to turn calls to '[' into names of variables, do
the derivative business and then turn them back again.  This is not
easy to do in complete generality, but if you are only talking about
singly subscripted arrays, you can get somewhere.  Here is an outline
of what I mean.

 Ident - ([A-z][A-z0-9_.]*)
 Subsc - ([A-z][A-z0-9_.]*|[0-9]+)
 patn - paste(Ident, \\[, Subsc, \\], sep = )
 repl - \\1__\\2
 E - expression(A[1] * exp(A[2]*X) + A[3])
 Es - deparse(E[[1]])
 Es
[1] A[1] * exp(A[2] * X) + A[3]
 Ess - gsub(patn, repl, Es)
 Ess
[1] A__1 * exp(A__2 * X) + A__3
 Ex - parse(text = Ess)[[1]]
 Ex
A__1 * exp(A__2 * X) + A__3

OK, the calls to `[` have been replaced by variables with two
underscores in the middle.  We hope this works - there is a strong
assumption here on just how complicated your indices are, for example.
We are assuming they are either identifiers (not using two successive
underscores in their names) or numbers.  If they are not, you need to
get even craftier.


 Ex1 - D(Ex, A__2)
 Ex1
A__1 * (exp(A__2 * X) * X)
 Ex1s - deparse(Ex1)
 Ex1s
[1] A__1 * (exp(A__2 * X) * X)
 pat1 - paste(Ident, __, Subsc, sep = )
 rep1 - \\1\\[\\2\\]
 Ex1ss - gsub(pat1, rep1, Ex1s)
 Ex1ss
[1] A[1] * (exp(A[2] * X) * X)
 Ex2 - parse(text = Ex1ss)[[1]]
 Ex2
A[1] * (exp(A[2] * X) * X)

Which is the required result.  This is messy and gets messier if you
are looking for some kind of generality, but you need to remember, R
is not, and never will be, a replacement for Maple.

 
 
 Thanks for any help with this.

Best of luck in automating this, but the tools are there.

Bill Venables.

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


[R] how to call forth a class definition buried in a package

2006-08-15 Thread Bálint Czúcz
Dear list,

I would like to use this C call:
tmp - .Call(R_MPinv, covariance, tol, svdmem)

but it gives me the following error:
Error in getClass(LinStatExpectCovarMPinv) :
LinStatExpectCovarMPinv is not a defined class

I think the reason for this is that the C code starts with this call:

PROTECT(ans = NEW_OBJECT(MAKE_CLASS(LinStatExpectCovarMPinv)));


and the Class definition for LinStatExpectCovarMPinv hides deep
buried in package party.

Were it not for the C call, i could avoid this error by specifying a
where= argument to GetClass(). But what can I do now?

Thank you!
Bálint

-- 
Czúcz Bálint
PhD hallgató
BCE KTK Talajtan és Vízgazdálkodás Tanszék
1118 Budapest, Villányi út 29-43.

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


Re: [R] Question on .Options$max.print - print/cat extremely long strings on a screen

2006-08-15 Thread Hans-Joerg Bibiko
Dear list members,

sorry for my incompleteness!



My problem is the following:
(R 2.3.1 on Mac OS X 10.4.7 RAM 1GByte using Mac GUI)

I have a function like this:

foo1 - function()
{
out - NULL
for(i in 1:10010) out - paste(out, i, . line\n, sep=)
return(out) 
}
a - foo1()

Now I want to display 'a' on the screen:

cat(a)

This doesn't work. 'a' is displayed only until '830. line'

print(a)

The same, 'a' is displayed only until '\n754. line\n755'.

cat(a, file='test.txt') # OK
works fine. That means, internally 'a' is fine.


Then I tried this way:

foo2 - function()
{
out - NULL
for(i in 1:10010) out - c(out, paste(i, . line, sep=))
return(out) 
}
a - foo2()

With

cat(a,sep=\n)
I see the complete content of 'a'.

paste(a, collapse = \n)
I only see 'a' until '\n754. line\n755'.

cat(paste(a, collapse =\n))
I only see 'a' until '830. line'.

cat(paste(a, collapse =\n),file='test.txt')
This is OK.


My question now is whether there is an option to specifiy the maximum  
size of a string which is displayed on the screen (running R in a  
GUI)? Or is this fixed?

I read the help page about '.Options' and I found a variable  
'max.print' with the comment that is not yet used in basic R.
I don't know whether this variable is responsible for that. I  
increased 'max.print' but nothing changed.

##
If I try this code on a Windows XP machine with 756 MByte RAM R-GUI  
says after executing

foo1 - function()
{
out - NULL
for(i in 1:10010) out - paste(out, i, . line\n, sep=)
return(out) 
}
a - foo1()
cat(a)

...
7548. lineWarning message:
printing of extremely long output is truncated

At least Windows writes a warning message.

###
If I start a R session without the R-GUI via Mac Terminal typing 'R'

a-foo1()
cat(a)

everything works perfectly!!!


I know the issue of outputting long strings and the way to display  
long strings via foo2() would be ok for me, but I spent some time to  
figure out why my function foo1() didn't work. R, running in GUI on  
Mac, don't give you a warning. Now I know that if I print a long  
string at the Mac-GUI-console the missing final quote character is an  
indicator for a truncated output on a Mac. Maybe it would be nice to  
output a warning(?)

Cheers,

Hans

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


[R] Protection stack overflow

2006-08-15 Thread Paul Koufalas
G'day all.

I'm a new user of R -- but an arms-length user, as I'm running it from
Octave via the ROctave interface that Duncan Temple Lang wrote some
years ago and Stephan Van Der Walt recently updated for use with Octave
2.1.71. I'm using R version 2.1.1. ROctave uses libR.so to provide the
interface. My system is Ubuntu linux 5.10 and I'm using the packages
that come with this distro.

I'm getting a protection stack overflow error when recursively
calculating AR(p) time series models using the arima() function. The
recursion is involved because I calculate a new model with each new time
series data point. (I'm trying to reproduce some results in a research
paper and this is what they're doing.)

I've tried setting the expressions parameter to a higher number using
options(expressions=50) but I'm still getting this problem with
stack overflow. I can get about 400-odd iterations and then the overflow
error appears. I yet haven't tried running the recursion natively in R,
and I realise I should do that.

Your advice would be much appreciated.

Cheers,
Paul.

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


[R] Query the kind of calling a funcion

2006-08-15 Thread Hans-Joerg Bibiko

Dear all,


My question is concerned to the kind how a function is called.

Example A:

  foo(1)

Example B:

  a - foo(1)


Is there any way for the function foo() to recognise whether the  
returned value of foo() is stored in a variable or not, i.e. to  
distinguish between Example A and B?


Any comments are welcomed.

Many thanks in advance

Hans

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


Re: [R] Protection stack overflow

2006-08-15 Thread Prof Brian Ripley
R has a command-line option to set the ppstack size,

 --max-ppsize=NSet max size of protect stack to N

Looks like you need to supply this (and it can be done with embedded R)
if the problem persists with current R.

You could also try arima0 or even ar to do the fitting.

On Tue, 15 Aug 2006, Paul Koufalas wrote:

 G'day all.
 
 I'm a new user of R -- but an arms-length user, as I'm running it from
 Octave via the ROctave interface that Duncan Temple Lang wrote some
 years ago and Stephan Van Der Walt recently updated for use with Octave
 2.1.71. I'm using R version 2.1.1. ROctave uses libR.so to provide the
 interface. My system is Ubuntu linux 5.10 and I'm using the packages
 that come with this distro.

Note that your version of R is well outdated, and the posting guide did 
ask you to update it BEFORE posting.

 I'm getting a protection stack overflow error when recursively
 calculating AR(p) time series models using the arima() function. The
 recursion is involved because I calculate a new model with each new time
 series data point. (I'm trying to reproduce some results in a research
 paper and this is what they're doing.)
 
 I've tried setting the expressions parameter to a higher number using
 options(expressions=50) but I'm still getting this problem with
 stack overflow. I can get about 400-odd iterations and then the overflow
 error appears. I yet haven't tried running the recursion natively in R,
 and I realise I should do that.
 
 Your advice would be much appreciated.
 
 Cheers,
 Paul.
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

-- 
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to reply to a thread if .. R-help mails .. in digest

2006-08-15 Thread Martin Maechler
Thanks a lot, Ted, for the good extensive explanation.
Let me try summarize, confirm and add a bit w/o repeating too much:

- If you get regular mailing list e-mails, and reply to
  postings, any decent mail software will take the 
  'Message-Id:' header of the message you reply to, and
  produce 'In-Reply-To:' and 'References:' headers from it, 
  and will add them to your own e-mail.
  [ Most e-mail softwares keep these headers hidden, and quite a
few pieces of e-mail crapware don't even allow you to see these headers.]

 { And as Ted mentioned; unfortunately, we still have too many
   r-help posters who **misuse** their e-mail-software's `reply'
   and then inadvertently hijack threads... }

- Yes, the threads are *not* built from 'Subject:'s but rather
  using the e-mail headers 'References:' and/or 'In-Reply-To:'

This is all based on international e-mail standards (RFC/..) mostly going
back into the age where most R-help readers have not ever used e-mail..

- With mailman, there are 2 ways to receive digests:
   (1) plain text  --- which is default, since some dumb e-mail
programs can only deal with those
   (2) MIME  --- which uses the MIME standard to send the digest.
   With a good e-mail software, this means that you
   get ``each message as attachment'' {that's how it
   typically looks to the user}
   which you then can open - in your mail software(!) -
   and it will behave like a regular e-mail; it has a
   (typically hidden) 'Message-ID:' etc
  == when replying to such a message you automatically
  get both:
   1) correct Subject
   2) correct In-Reply-To + References === correct thread

*So* : What we (and many) recommend to all  digest  subscribers
 is to activate the MIME option - on their mailing list
 membership info page to where you get from
 https://stat.ethz.ch/mailman/listinfo/r-help 
 after entering your e-mail address at the very bottom of the
 page (and then your list password).

-- but, as Michael Dewey just mentioned,
 unfortunately there are (too many) pieces of e-mail crapware
 around which even don't support the above correctly...

Martin Maechler, 
ETH Zurich, provider of (most of) the mailing list infrastructure for R.

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


Re: [R] split a y-axis to show data on different scales

2006-08-15 Thread Johannes Hüsing
 The pro's and con's of using scale breaks were discussed by
 Cleveland (1985) The Elements of Graphing Data (Wadsworth, pp. 85-91,
 149).  I don't know what Cleveland said about this is the second edition

Spencer Graves:
 but I believe there are times when scale breaks are
 appropriate, but the display should make this nonstandard transition
 very clear;

... in which case you are close to having two graphs
sharing an x-axis and therefore saving on ink (yay!).

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


[R] A plot with a bisector

2006-08-15 Thread Amir Safari
 
  Dear R Users,
  How to make a plot with a bisector is my question. Indeed, I want to plot a 
QQplot. But it doesn't have such a bisector.
  Thanks a lot.
  Amir 


-

[[alternative HTML version deleted]]

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


[R] Aliases for arguments in a function

2006-08-15 Thread [EMAIL PROTECTED]
Hi all. 

I have a function that I would like to use either the argument name as 
originally defined or another name. Within the function (and other functions) 
use the argument name as originally written, so I don't want to simply remove 
the old argument name for the new one, but simply allow the function to treat 
both argument names as equivalent. 

Here is an example:

foo - function(arg1, this)
{
if(this  0) stop(this must be positive)
return(arg1/this)
}

foo(arg1=5, this=10)

But, I also want foo() to work equivalently with the following (where 'this' 
and 'that 'are treated as if they were the same):
foo(arg1=5, that=10)

Any thoughts would be appreciated.

Thanks,
Ken


-

[[alternative HTML version deleted]]

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


Re: [R] Aliases for arguments in a function

2006-08-15 Thread Richard M. Heiberger
foo - function(arg1, this, that) {
  if (missing(this)  !missing(that)) this - that
  if(this  0) stop(this must be positive)
  return(arg1/this)
}

foo(arg1=5, this=10)

foo(arg1=5, that=10)

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


[R] nls

2006-08-15 Thread Xiaodong Jin
  Is there anyway to change any y[i] value (i=2,...6) to make following NLS 
workable? 
   
  x - c(0,5,10,15,20,25,30)
  y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000)
  lm(1/y ~~ x)
  nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE)
   
  #0.0920573 :  1.16122 0.01565 1.0 
#Error in numericDeriv(form[[3]], names(ind), env) : 
#Missing value or an infinity produced when evaluating the model


-


[[alternative HTML version deleted]]

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


Re: [R] nls

2006-08-15 Thread Prof Brian Ripley
You problem is x^c for x = 0.  If you intended only c  1, try a starting 
value meeting that condition (but it seems that the optimal c is about 
0.27 is you increase x slightly).

Why have you used ~~ ?  (Maybe because despite being asked not to, you 
sent HTML mail?)

On Tue, 15 Aug 2006, Xiaodong Jin wrote:

   Is there anyway to change any y[i] value (i=2,...6) to make following NLS 
 workable? 

   x - c(0,5,10,15,20,25,30)
   y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000)
   lm(1/y ~~ x)
   nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE)

   #0.0920573 :  1.16122 0.01565 1.0 
 #Error in numericDeriv(form[[3]], names(ind), env) : 
 #Missing value or an infinity produced when evaluating the model
 
   
 -
 
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

-- 
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Aliases for arguments in a function

2006-08-15 Thread Prof Brian Ripley
foo - function(arg1, this=that, that)
...

works.

On Tue, 15 Aug 2006, [EMAIL PROTECTED] wrote:

 Hi all. 
 
 I have a function that I would like to use either the argument name as 
 originally defined or another name. Within the function (and other functions) 
 use the argument name as originally written, so I don't want to simply remove 
 the old argument name for the new one, but simply allow the function to treat 
 both argument names as equivalent. 
 
 Here is an example:
 
 foo - function(arg1, this)
 {
 if(this  0) stop(this must be positive)
 return(arg1/this)
 }
 
 foo(arg1=5, this=10)
 
 But, I also want foo() to work equivalently with the following (where 'this' 
 and 'that 'are treated as if they were the same):
 foo(arg1=5, that=10)
 
 Any thoughts would be appreciated.
 
 Thanks,
 Ken
 
   
 -
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

-- 
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] ARCH Jump diffusion models

2006-08-15 Thread Hazrachoudhury, Avishek \(RSCH\)
Hi 

I was wondering if you could tell me how would I go abt fitting the following 
model to a data series. It is an ARCH jump diffusion model.


X(t)=k+ ó (t)z(t)+J x Y, and

ó(t)= sqrt[a +b [X(t-1)-k]^2]


Here k is mean of the series. 
z(t) is iid standard normal error
J is a jump variable (either 1 or 0)
Y is normally distributed with mean m and variance v


Regards
Avishek


If you are not an intended recipient of this e-mail, please notify the sender, 
delete it and do not read, act upon, print, disclose, copy, retain or 
redistribute it. Click here for important additional terms relating to this 
e-mail. http://www.ml.com/email_terms/


[[alternative HTML version deleted]]

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


Re: [R] nls

2006-08-15 Thread Peter Dalgaard
Prof Brian Ripley [EMAIL PROTECTED] writes:

 You problem is x^c for x = 0.  If you intended only c  1, try a starting 
 value meeting that condition (but it seems that the optimal c is about 
 0.27 is you increase x slightly).

Surely you mean c  0.

  nls(1/y ~ a+b*x^exp(c), start=list(a=1.16122,b=0.01565,c=0))
Nonlinear regression model
  model:  1/y ~ a + b * x^exp(c)
   data:  parent.frame()
 a  b  c
 0.9944025  0.1953168 -1.1495206
 residual sum-of-squares:  0.03303464
  nls(1/y ~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=exp(-1.1)))
Nonlinear regression model
  model:  1/y ~ a + b * x^c
   data:  parent.frame()
a b c
0.9944026 0.1953165 0.3167891
 residual sum-of-squares:  0.03303464

(but even setting c=exp(-1) triggers the error; there could be cause
for robustifying the nls algorithm)
 
 Why have you used ~~ ?  (Maybe because despite being asked not to, you 
 sent HTML mail?)
 
 On Tue, 15 Aug 2006, Xiaodong Jin wrote:
 
Is there anyway to change any y[i] value (i=2,...6) to make following NLS 
  workable? 
 
x - c(0,5,10,15,20,25,30)
y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000)
lm(1/y ~~ x)
nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE)
 
#0.0920573 :  1.16122 0.01565 1.0 
  #Error in numericDeriv(form[[3]], names(ind), env) : 
  #Missing value or an infinity produced when evaluating the model
  
  
  -
  
  
  [[alternative HTML version deleted]]
  
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
  
 
 -- 
 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@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

-- 
   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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] nls

2006-08-15 Thread Petr Pikal
Hi

Why do you want to change your variable values? It smells a rat to 
me.

If you just change your a,b,c values nls arrives to some finite 
result (e.g. c=1.5 or c=0.3) . BTW by what magic you obtained such 
precise and wrong estimates for a,b,c?

HTH
Petr





On 15 Aug 2006 at 5:54, Xiaodong Jin wrote:

Date sent:  Tue, 15 Aug 2006 05:54:51 -0700 (PDT)
From:   Xiaodong Jin [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Subject:[R] nls

   Is there anyway to change any y[i] value (i=2,...6) to make
   following NLS workable? 
 
   x - c(0,5,10,15,20,25,30)
   y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000)
   lm(1/y ~~ x) nls(1/y ~~ a+b*x^c,
   start=list(a=1.16122,b=0.01565,c=1), trace=TRUE)
 
   #0.0920573 :  1.16122 0.01565 1.0 
 #Error in numericDeriv(form[[3]], names(ind), env) : 
 #Missing value or an infinity produced when evaluating the
 #model
 
 
 -
 
 
  [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html and provide commented,
 minimal, self-contained, reproducible code.

Petr Pikal
[EMAIL PROTECTED]

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


Re: [R] nls

2006-08-15 Thread Prof Brian Ripley
On Tue, 15 Aug 2006, Peter Dalgaard wrote:

 Prof Brian Ripley [EMAIL PROTECTED] writes:
 
  You problem is x^c for x = 0.  If you intended only c  1, try a starting 
  value meeting that condition (but it seems that the optimal c is about 
  0.27 is you increase x slightly).
 
 Surely you mean c  0.

I did.

   nls(1/y ~ a+b*x^exp(c), start=list(a=1.16122,b=0.01565,c=0))
 Nonlinear regression model
   model:  1/y ~ a + b * x^exp(c)
data:  parent.frame()
  a  b  c
  0.9944025  0.1953168 -1.1495206
  residual sum-of-squares:  0.03303464
   nls(1/y ~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=exp(-1.1)))
 Nonlinear regression model
   model:  1/y ~ a + b * x^c
data:  parent.frame()
 a b c
 0.9944026 0.1953165 0.3167891
  residual sum-of-squares:  0.03303464
 
 (but even setting c=exp(-1) triggers the error; there could be cause
 for robustifying the nls algorithm)

Well, there is an option to use a bounded-region algorithm, e.g.

x - c(0,5,10,15,20,25,30)
y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000)
nls(1/y ~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE,
algorithm=port, lower=c(-Inf, -Inf, 0), upper=rep(Inf, 3))

works.

  
  Why have you used ~~ ?  (Maybe because despite being asked not to, you 
  sent HTML mail?)
  
  On Tue, 15 Aug 2006, Xiaodong Jin wrote:
  
 Is there anyway to change any y[i] value (i=2,...6) to make following 
   NLS workable? 
  
 x - c(0,5,10,15,20,25,30)
 y - c(1.0,0.82000,0.68000,0.64000,0.7,0.68667,0.64000)
 lm(1/y ~~ x)
 nls(1/y ~~ a+b*x^c, start=list(a=1.16122,b=0.01565,c=1), trace=TRUE)
  
 #0.0920573 :  1.16122 0.01565 1.0 
   #Error in numericDeriv(form[[3]], names(ind), env) : 
   #Missing value or an infinity produced when evaluating the model
   
 
   -
   
   
 [[alternative HTML version deleted]]
   
   __
   R-help@stat.math.ethz.ch mailing list
   https://stat.ethz.ch/mailman/listinfo/r-help
   PLEASE do read the posting guide 
   http://www.R-project.org/posting-guide.html
   and provide commented, minimal, self-contained, reproducible code.
   
  
  -- 
  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@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
  
 
 

-- 
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] question re: summarry.lm and NA values

2006-08-15 Thread r user
Is there a way to get the following code to include
NA values where the coefficients are “NA”?

((summary(reg))$coefficients)

explanation:

Using a loop, I am running regressions on several
“subsets” of “data1”.

“reg - ( lm(lm(data1[,1] ~., data1[,2:l])) )”

My regression has 10 independent variables, and I
therefore expect 11 coefficients.
After each regression, I wish to save the coefficients
and standard errors of the coefficients in a table
with 22 columns.

I successfully extract the coefficients using the
following code:
“reg$coefficients”

I attempt to extract the standard errors using :

aperm((summary(reg))$coefficients)[2,]

((summary(reg))$coefficients)

My problem:
For some of my subsets, I am missing data for one or
more of the independent variables.  This of course
causes the coefficients and standard erros for this
variable to be “NA”.

Is there a way to include the NA standard errors, so
that I have the same number of standard erros and
coefficients for each regression, and can then store
the coefficients and standard erros in my table of 22
columns?

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


[R] Presentation of multiple models in one table using xtable

2006-08-15 Thread Ajay Narottam Shah
Consider this situation:
 x1 - runif(100); x2 - runif(100); y - 2 + 3*x1 - 4*x2 + rnorm(100)
 m1 - summary(lm(y ~ x1))
 m2 - summary(lm(y ~ x2))
 m3 - summary(lm(y ~ x1 + x2))

Now you have estimated 3 different competing models, and suppose you
want to present the set of models in one table. xtable(m1) is cool,
but doing that thrice would give us 3 different tables.

What I want is this one table:

---
M1 M2  M3
---
Intercept 0.0816 3.6292 2.2272
 (0.5533)   (0.2316)***(0.2385)***

x12.81512.7606
 (0.5533)***   (0.3193)***

x2  -4.2899-4.2580
(0.401)*** (0.3031)***

$\sigma_e$1.538  1.175  0.8873
$R^2$ 0.2089 0.5385 0.7393
---

How would one set about doing this? I am hoping that it's possible to
write a function xtable.multi.lm where one would say
xtable.multi.lm(m1,m2,m3) and get the above table.

My sense is there are three challenges:

1. How to write a general R function which eats a unpredictable number
   of summary(lm()) objects, and fill out a matrix with results such
   as the above.

2. How to get a good xtable(), with decimal alignment and with the ***
   stuff (actually, $^{***}$). Will there be any catch in dropping
   into mathmode for $R^2$? After each pair of lines, I'd like to have
   a \\[2mm] so as to get a nice spacing in the table.

3. This style of presentation seems relevant to a whole host of models
   - whether ARCH models or survival models - not just OLS
   regressions. It would be very nice if one supported all manner of
   model objects and not just what comes out of lm().

I'm happy to take a crack at writing this function, which should
ideally go back into the xtable library. But I don't have an idea on
how to go about these two questions. If you will guide me, I am happy
to work on it. :-)

-- 
Ajay Shah  http://www.mayin.org/ajayshah  
[EMAIL PROTECTED] http://ajayshahblog.blogspot.com
*(:-? - wizard who doesn't know the answer.

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


[R] fMultivar OLS - how to do dynamic regression?

2006-08-15 Thread Kerpel, John
Hi folks!

 

Does anybody know how to use the OLS function in fMultivar to do dynamic
regression?  I've tried specifying lags in OLS using a data series
created in fSeries and it doesn't seem to work.  I've done dynamic
regression using dyn$lm and I was wondering how to accomplish the same
thing using the OLS function from fMultivar.  Thanks!

 

John


[[alternative HTML version deleted]]

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


Re: [R] question re: summarry.lm and NA values

2006-08-15 Thread Petr Pikal
Hi

On 15 Aug 2006 at 7:01, r user wrote:

Date sent:  Tue, 15 Aug 2006 07:01:13 -0700 (PDT)
From:   r user [EMAIL PROTECTED]
To: rhelp r-help@stat.math.ethz.ch
Subject:[R] question re: summarry.lm and NA values

 Is there a way to get the following code to include
 NA values where the coefficients are  NA ?
 
 ((summary(reg))$coefficients)

better
coef(reg)

 
 explanation:
 
 Using a loop, I am running regressions on several
  subsets  of  data1 .
 
  reg - ( lm(lm(data1[,1] ~., data1[,2:l])) ) 
 
 My regression has 10 independent variables, and I
 therefore expect 11 coefficients.
 After each regression, I wish to save the coefficients
 and standard errors of the coefficients in a table
 with 22 columns.
 
 I successfully extract the coefficients using the
 following code:
  reg$coefficients 
 
 I attempt to extract the standard errors using :
 
 aperm((summary(reg))$coefficients)[2,]
 
 ((summary(reg))$coefficients)
 
 My problem:
 For some of my subsets, I am missing data for one or
 more of the independent variables.  This of course
 causes the coefficients and standard erros for this
 variable to be  NA .

??%^*^??

What version? My lm behaves in accordance with na.action and it 
throws an error in case na.fail, computes a value in case of na.omit 
or na.exclude and again throws an error if the variable consist 
exclusively from NA values. 

The only way how to get NA in coeficient is when a variable is either 
constant or linear combination of other variable(s). Then
coef(reg) 
will give you correctly NA in the variable which appears constant and 
in this case you could use it for setting standard error also as NA 
let say by using ifelse statement and matching of names.

HTH
Petr

 
 Is there a way to include the NA standard errors, so
 that I have the same number of standard erros and
 coefficients for each regression, and can then store
 the coefficients and standard erros in my table of 22
 columns?
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html and provide commented,
 minimal, self-contained, reproducible code.

Petr Pikal
[EMAIL PROTECTED]

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


Re: [R] question re: summarry.lm and NA values

2006-08-15 Thread Petr Pikal
Hi

just as a quick workaround you probably can use aliased value from 
summary

fff-rep(summary(reg)$aliased,4)
dim(fff)-c(no.of.your.variables,4)
fff[which(fff)]-NA
fff[which(!fff)]-coef(summary(reg))

to get coef matrix with NA values

HTH
Petr


On 15 Aug 2006 at 17:15, r user [EMAIL PROTECTED], wrote:

From:   Petr Pikal [EMAIL PROTECTED]
To: r user [EMAIL PROTECTED], rhelp 
r-help@stat.math.ethz.ch
Subject:Re: [R] question re: summarry.lm and NA values
Date sent:  Tue, 15 Aug 2006 17:15:01 +0200

 Hi
 
 On 15 Aug 2006 at 7:01, r user wrote:
 
 Date sent:Tue, 15 Aug 2006 07:01:13 -0700 (PDT)
 From: r user [EMAIL PROTECTED]
 To:   rhelp r-help@stat.math.ethz.ch
 Subject:  [R] question re: summarry.lm and NA values
 
  Is there a way to get the following code to include
  NA values where the coefficients are  NA ?
  
  ((summary(reg))$coefficients)
 
 better
 coef(reg)
 
  
  explanation:
  
  Using a loop, I am running regressions on several
   subsets  of  data1 .
  
   reg - ( lm(lm(data1[,1] ~., data1[,2:l])) ) 
  
  My regression has 10 independent variables, and I
  therefore expect 11 coefficients.
  After each regression, I wish to save the coefficients
  and standard errors of the coefficients in a table
  with 22 columns.
  
  I successfully extract the coefficients using the
  following code:
   reg$coefficients 
  
  I attempt to extract the standard errors using :
  
  aperm((summary(reg))$coefficients)[2,]
  
  ((summary(reg))$coefficients)
  
  My problem:
  For some of my subsets, I am missing data for one or
  more of the independent variables.  This of course
  causes the coefficients and standard erros for this
  variable to be  NA .
 
 ??%^*^??
 
 What version? My lm behaves in accordance with na.action and it 
 throws an error in case na.fail, computes a value in case of na.omit
 or na.exclude and again throws an error if the variable consist
 exclusively from NA values. 
 
 The only way how to get NA in coeficient is when a variable is either
 constant or linear combination of other variable(s). Then coef(reg)
 will give you correctly NA in the variable which appears constant and
 in this case you could use it for setting standard error also as NA
 let say by using ifelse statement and matching of names.
 
 HTH
 Petr
 
  
  Is there a way to include the NA standard errors, so
  that I have the same number of standard erros and
  coefficients for each regression, and can then store
  the coefficients and standard erros in my table of 22
  columns?
  
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html and provide commented,
  minimal, self-contained, reproducible code.
 
 

Petr Pikal
[EMAIL PROTECTED]

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


[R] Looking for info on the Regression Modeling Strategies in R course in DC area

2006-08-15 Thread Tim McDonald
Hello list,
   
  A colleague of mine mentioned a great course on  Regression Modeling 
Strategies in R. Anyone knows if this course is offered as public course in DC 
area?
   
  Thanks a bunch - TM


-


[[alternative HTML version deleted]]

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


[R] REML with random slopes and random intercepts giving strange results

2006-08-15 Thread Simon Pickett
Hi everyone,
I have been using REML to derive intercepts and coeficients for each
individual in a growth study. So the code is
m2 - lmer(change.wt ~ newwt+(newwt|id), data = grow)

Calling coef(model.lmer) gives a matrix with this information which is
what I want. However, as a test I looked at each individual on its own and
used a simple linear regression to obtain the same information, then I
compared the results. It looks like the REML method doesnt seem to
approximate the two parameters as well as using the simple linear
regression on each individual separately, as judged by looking at graphs.
Indeed, why do the results differ at all?
Excuse my naivety if this is a silly question.
Thanks to everyone for replying to my previous questions, very much
appreciated.
Simon Pickett
PhD student
Centre For Ecology and Conservation
Tremough Campus
University of Exeter in Cornwall
TR109EZ
Tel 01326371852

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


Re: [R] question re: summarry.lm and NA values

2006-08-15 Thread Prof Brian Ripley
On Tue, 15 Aug 2006, Petr Pikal wrote:

 Hi
 
 On 15 Aug 2006 at 7:01, r user wrote:
 
 Date sent:Tue, 15 Aug 2006 07:01:13 -0700 (PDT)
 From: r user [EMAIL PROTECTED]
 To:   rhelp r-help@stat.math.ethz.ch
 Subject:  [R] question re: summarry.lm and NA values
 
  Is there a way to get the following code to include
  NA values where the coefficients are  NA ?
  
  ((summary(reg))$coefficients)
 
 better
 coef(reg)

coef(summary(reg)), perhaps.

  explanation:
  
  Using a loop, I am running regressions on several
   subsets  of  data1 .
  
   reg - ( lm(lm(data1[,1] ~., data1[,2:l])) ) 
  
  My regression has 10 independent variables, and I
  therefore expect 11 coefficients.
  After each regression, I wish to save the coefficients
  and standard errors of the coefficients in a table
  with 22 columns.
  
  I successfully extract the coefficients using the
  following code:
   reg$coefficients 
  
  I attempt to extract the standard errors using :
  
  aperm((summary(reg))$coefficients)[2,]
  
  ((summary(reg))$coefficients)
  
  My problem:
  For some of my subsets, I am missing data for one or
  more of the independent variables.  This of course
  causes the coefficients and standard erros for this
  variable to be  NA .
 
 ??%^*^??
 
 What version? My lm behaves in accordance with na.action and it 
 throws an error in case na.fail, computes a value in case of na.omit 
 or na.exclude and again throws an error if the variable consist 
 exclusively from NA values. 
 
 The only way how to get NA in coeficient is when a variable is either 
 constant or linear combination of other variable(s). Then
 coef(reg) 
 will give you correctly NA in the variable which appears constant and 
 in this case you could use it for setting standard error also as NA 
 let say by using ifelse statement and matching of names.

That happens in the print method, stats:::print.summary.lm contains

coefs - x$coefficients
if (!is.null(aliased - x$aliased)  any(aliased)) {
cn - names(aliased)
coefs - matrix(NA, length(aliased), 4, dimnames = list(cn,
colnames(coefs)))
coefs[!aliased, ] - x$coefficients
}

so the code is already available

  
  Is there a way to include the NA standard errors, so
  that I have the same number of standard erros and
  coefficients for each regression, and can then store
  the coefficients and standard erros in my table of 22
  columns?
  
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html and provide commented,
  minimal, self-contained, reproducible code.
 
 Petr Pikal
 [EMAIL PROTECTED]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

-- 
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] question re: summarry.lm and NA values

2006-08-15 Thread Berton Gunter
Is there a way to... always has the answer yes in R (or C or any
language for that matter). The question is: Is there a GOOD way...? where
good depends on the specifics of the situation. So after that polemic,
below is an effort to answer, (adding to what Petr Pikal already said):

-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
The business of the statistician is to catalyze the scientific learning
process.  - George E. P. Box
 
 

 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of r user
 Sent: Tuesday, August 15, 2006 7:01 AM
 To: rhelp
 Subject: [R] question re: summarry.lm and NA values
 
 Is there a way to get the following code to include
 NA values where the coefficients are NA?
 
 ((summary(reg))$coefficients)
BAAAD! Don't so this. Use the extractor on the object: coef(reg) 
This suggests that you haven't read the documentation carefully, which tends
to arouse the ire of would-be helpers.

 
 explanation:
 
 Using a loop, I am running regressions on several
 subsets of data1.
 
 reg - ( lm(lm(data1[,1] ~., data1[,2:l])) )
??? There's an error here I think. Do you mean update()? Do you have your
subscripting correct?

 
 My regression has 10 independent variables, and I
 therefore expect 11 coefficients.
 After each regression, I wish to save the coefficients
 and standard errors of the coefficients in a table
 with 22 columns.
 
 I successfully extract the coefficients using the
 following code:
 reg$coefficients
Use the extractor, coef()

 
 I attempt to extract the standard errors using :
 
 aperm((summary(reg))$coefficients)[2,]

BAAAD! Use the extractor vcov(): sqrt(diag(vcov(reg)))
 
 ((summary(reg))$coefficients)
 
 My problem:
 For some of my subsets, I am missing data for one or
 more of the independent variables.  This of course
 causes the coefficients and standard erros for this
 variable to be NA.
Not it doesn't, as Petr said.

One possible approach: Assuming that a variable is actually missing (all
NA's), note that coef(reg) is a named vector, so that the character string
names of the regressors actually used are available. You can thus check for
what's missing and add them as NA's at each return. Though I confess that I
see no reason to put things ina matrix rather than just using a list. But
that's a matter of personal taste I suppose.

 
 Is there a way to include the NA standard errors, so
 that I have the same number of standard erros and
 coefficients for each regression, and can then store
 the coefficients and standard erros in my table of 22
 columns?
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


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


Re: [R] REML with random slopes and random intercepts giving strange results

2006-08-15 Thread Doran, Harold
I don't this is because you are using REML. The BLUPs from a mixed model
experience some shrinkage whereas the OLS estimates would not.  

 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of Simon Pickett
 Sent: Tuesday, August 15, 2006 11:34 AM
 To: r-help@stat.math.ethz.ch
 Subject: [R] REML with random slopes and random intercepts 
 giving strange results
 
 Hi everyone,
 I have been using REML to derive intercepts and coeficients 
 for each individual in a growth study. So the code is
 m2 - lmer(change.wt ~ newwt+(newwt|id), data = grow)
 
 Calling coef(model.lmer) gives a matrix with this information 
 which is what I want. However, as a test I looked at each 
 individual on its own and used a simple linear regression to 
 obtain the same information, then I compared the results. It 
 looks like the REML method doesnt seem to approximate the two 
 parameters as well as using the simple linear regression on 
 each individual separately, as judged by looking at graphs.
 Indeed, why do the results differ at all?
 Excuse my naivety if this is a silly question.
 Thanks to everyone for replying to my previous questions, 
 very much appreciated.
 Simon Pickett
 PhD student
 Centre For Ecology and Conservation
 Tremough Campus
 University of Exeter in Cornwall
 TR109EZ
 Tel 01326371852
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


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


[R] nls convergence problem

2006-08-15 Thread Earl F. Glynn
I'm having problems getting nls to agree that convergence has occurred in a 
toy problem.

nls.out never gets defined when there is an error in nls.  Reaching the 
maximum number of iterations is alway an error, so nls.out never gets 
defined when the maximum number of iterations is reched.

From ?nls.control:
  tol: A positive numeric value specifying the tolerance level for
  the relative offset convergence criterion.

From some S-Plus documentation I found online via Google:
http://www.uni-muenster.de/ZIV/Mitarbeiter/BennoSueselbeck/s-html/helpfiles/nls.control.html

tolerance:
tolerance for the relative offset convergence criterion in the algorithm. 
Default 0.001. Note that the convergence test used in nls() is strictly 
relative. Therefore if the solution to a problem turned out to be a perfect 
fit (unlikely except in artificial examples), convergence is not guaranteed 
to be recognized by the algorithm.


Here's my toy problem:
 ?nls.control
 ?nls
 # Method 2
 X - 0:15
 Y - 9.452 * exp(-0.109*X) + 5.111   # Toy problem

 nls.out - nls(Y ~ a*exp(b*X)+c,
+start=list(a=6,b=-0.5,c=1),
+control=nls.control(maxiter=15, tol=0.01),  # nothing makes 
sense here
+trace=TRUE)
1016.507 :   6.0 -0.5  1.0
143.5807 :   6.1680290 -0.1506021  4.4013020
7.306365 :   9.10093164 -0.09114858  5.44620298
0.0342703 :   9.2801070 -0.1109063  5.2795803
3.001985e-05 :   9.4506654 -0.1089749  5.1122982
1.918531e-14 :   9.452 -0.109  5.111
6.894644e-28 :   9.452 -0.109  5.111
2.208811e-29 :   9.452 -0.109  5.111
7.888609e-30 :   9.452 -0.109  5.111
7.888609e-30 :   9.452 -0.109  5.111
7.099748e-30 :   9.452 -0.109  5.111
3.155444e-30 :   9.452 -0.109  5.111
3.155444e-30 :   9.452 -0.109  5.111
3.155444e-30 :   9.452 -0.109  5.111
3.155444e-30 :   9.452 -0.109  5.111
3.155444e-30 :   9.452 -0.109  5.111
Error in nls(Y ~ a * exp(b * X) + c, start = list(a = 6, b = -0.5, c = 1), 
:
number of iterations exceeded maximum of 15

There is near-perfect convergence after 12 iterations, but I cannot figure 
out a way for R to recognize it.

What does relative offset convergence criterion mean?  How can nls.control 
be used to say this problem has converged, and have nls exit without an 
error?

Should the R documentation be modified to explain what relative offset 
convergence criterion means?  Should the R documentation be expanded to 
include the comment from the S-Plus: Therefore if the solution to a problem 
turned out to be a perfect fit (unlikely except in artificial examples), 
convergence is not guaranteed to be recognized by the algorithm.  If this 
is true, this seems like a suboptimal design.

Thanks for any insight about this.

efg

Earl F. Glynn
Scientific Programmer
Stowers Institute for Medical Research

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


[R] How to show classes of all columns of a data frame?

2006-08-15 Thread T Mu
Hi all,

Suppose I have a data frame myDF, col A is factor, col B is numeric, col C
is character. I can get their classes by

 class(myDF$A)

but is there a quick way to show what classes of all columns are? Thank you.

Tian

[[alternative HTML version deleted]]

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


Re: [R] How to show classes of all columns of a data frame?

2006-08-15 Thread Marc Schwartz (via MN)
On Tue, 2006-08-15 at 13:10 -0400, T Mu wrote:
 Hi all,
 
 Suppose I have a data frame myDF, col A is factor, col B is numeric, col C
 is character. I can get their classes by
 
  class(myDF$A)
 
 but is there a quick way to show what classes of all columns are? Thank you.
 
 Tian

Depending upon the output format you desire:

 lapply(iris, class)
$Sepal.Length
[1] numeric

$Sepal.Width
[1] numeric

$Petal.Length
[1] numeric

$Petal.Width
[1] numeric

$Species
[1] factor


or

 sapply(iris, class)
Sepal.Length  Sepal.Width Petal.Length  Petal.Width  Species
   numericnumericnumericnumeric factor


See ?lapply and ?sapply


HTH,

Marc Schwartz

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


[R] zlim not working in persp3d

2006-08-15 Thread Ken Nussear
Hi list,

I'm trying to limit the z axis plotted using persp3d {rgl} to a given  
range. I'm trying the following statement
persp3d(g0,pa,D, xlim=range(g0),ylim=range(pa), zlim=range 
(zbounds),col=lightblue)

where g0  and pa are both a range of 20 numbers from 0.05 to 1, and D  
is a 20 x 20 matrix with values ranging from 0 to 3. I'm trying  
to set a cap on the Z displayed in the surface to 1000 so I can see  
the lower levels in greater detail, but the persp3d statement doesn't  
appear to respond to any zlim parameter that I have tried.

Can anyone suggest a method?

I'm using R Gui for OSX ver 1.16, R version 2.3.1, and rgl is version  
0.67-2

Thanks

  Ken



[[alternative HTML version deleted]]

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


[R] rexp question

2006-08-15 Thread Spencer Jones
I am using rexp to generate several exponential distributions. I am passing
rexp a vector of rates , r. I am wanting to simulate a sample of size 200
for each rate so the code looks like: rexp(n=200*length(r),rate=r) this
gives me a vector of the random exponential variables, but they are all
disjointed b/c rexp goes through and simulates an exponential variable for
each rate and it does that 200  times, is there any way to get it to do it
the oppisite way,i.e., 200 times in a row for each rate. I have already
tried using a for loop but it is slow.


thanks,

Spencer

[[alternative HTML version deleted]]

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


[R] Grasper model error

2006-08-15 Thread Ken Nussear


I tried this over a the grasp users yahoo group and got no  
responseSo I wonder if anyone here knows about grasper

I keep getting this error when trying to run a model.


Error in smooth.construct.tp.smooth.spec(object, data, knots) :
Too many knots for t.p.r.s term: see `gam.control' to increase limit,  
or use a
different
basis, or see large data set help for `gam'.


I'm using R for OSX Gui Version 1.16 R-version 2.3.1, running grasper  
version
0.4-4


My dataset has ~7500 rows of input, and 21 cols of input.

Does anyone have an idea how to troubleshoot this?

Thanks


Ken



[[alternative HTML version deleted]]

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


Re: [R] rexp question

2006-08-15 Thread Greg Snow
Try:

 rexp(n=200*length(r), rate=rep(r, each=200) )

Hope this helps, 


-- 
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 Spencer Jones
Sent: Tuesday, August 15, 2006 12:15 PM
To: R-Help
Subject: [R] rexp question

I am using rexp to generate several exponential distributions. I am
passing rexp a vector of rates , r. I am wanting to simulate a sample of
size 200 for each rate so the code looks like:
rexp(n=200*length(r),rate=r) this gives me a vector of the random
exponential variables, but they are all disjointed b/c rexp goes through
and simulates an exponential variable for each rate and it does that 200
times, is there any way to get it to do it the oppisite way,i.e., 200
times in a row for each rate. I have already tried using a for loop but
it is slow.


thanks,

Spencer

[[alternative HTML version deleted]]

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

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


Re: [R] nls convergence problem

2006-08-15 Thread Dieter Menne
Earl F. Glynn efg at stowers-institute.org writes:

 
 Here's my toy problem:
  ?nls.control
  ?nls
  # Method 2
  X - 0:15
  Y - 9.452 * exp(-0.109*X) + 5.111   # Toy problem
 
  nls.out - nls(Y ~ a*exp(b*X)+c,
 +start=list(a=6,b=-0.5,c=1),
 +control=nls.control(maxiter=15, tol=0.01),  # nothing makes 
 sense here
 +trace=TRUE)

This toy problem is exactly what the warning is for:

Warning
Do not use nls on artificial zero-residual data. 

Add some noise and try again.

Dieter Menne

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


[R] Hierarchical clustering

2006-08-15 Thread Davendra Sohal
Hi,

I'm new to R; this is my second email to this forum. I have a dataset that
looks like this:


GeneA B
C D



A 511  6116
15151515

B  164  115
1515494

C 21646

D 1616

E  1461

F

G

H



There are about 4 genes with 3-10 patients per dataset. Gene expression
values fill the matrix.

I want to do a simple hierarchical clustering along BOTH rows and columns
and get a tree dendrogram output for both in one panel.

Is there a simple code for this, leaving out fancy options?



Thank you very much.

-DS.

[[alternative HTML version deleted]]

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


[R] model = F causing error in polr()

2006-08-15 Thread T Mu
Hi all,

I got an error message if I set model =F in polr(), like

 polr(y ~ x1 + x2, data1, model = F, method = probit)

Error in model.frame(formula, rownames, variables, varnames, extras,
extranames,  :
variable lengths differ (found for '(model)')
but

 polr(y ~ x1 + x2, data1, method = probit)

would work.

Why? Thank you,

Tian

[[alternative HTML version deleted]]

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


Re: [R] model = F causing error in polr()

2006-08-15 Thread Rolf Turner

Try ``model = FALSE'' rather than ``model = F'' and see if it makes a
difference.  You make have an unwanted variable named ``F'' lurking
somewhere.

(In general it is a *bad* idea to use ``F'' when ``FALSE'' is
intended.)

cheers,

Rolf Turner
[EMAIL PROTECTED]

===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===
Original message:

 Hi all,
 
 I got an error message if I set model =F in polr(), like
 
  polr(y ~ x1 + x2, data1, model = F, method = probit)
 
 Error in model.frame(formula, rownames, variables, varnames, extras,
 extranames,  :
 variable lengths differ (found for '(model)')
 but
 
  polr(y ~ x1 + x2, data1, method = probit)
 
 would work.
 
 Why? Thank you,
 
 Tian

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


Re: [R] nls convergence problem

2006-08-15 Thread Earl F. Glynn
Dieter Menne [EMAIL PROTECTED] wrote in message 
news:[EMAIL PROTECTED]
 Earl F. Glynn efg at stowers-institute.org writes:
 This toy problem is exactly what the warning is for:

 Warning
 Do not use nls on artificial zero-residual data.

 Add some noise and try again.

Thank you!

I had adapted some code and must confess I had read ?nls.control thoroughly, 
but not ?nls. I had even used debug on nls, traced it through line by line 
to the .Call statement, trying to figure out why nls.out never got defined. 
The source code has no comments at all.

IMHO, the warning should be in the Description at the top of the ?nls 
page, not at the bottom of the page. The warning should also appear on the 
?nls.control page. But, a better way would be to have a software design that 
eliminated the warning.

It's not clear to me why this problem cannot be fixed somehow. You 
shouldn't need to add noise to a problem to solve it. (It's a bit like 
saying addition works, but not for integers without adding some noise.) If 
there can be arbitrary defaults of maxiter=50, and (relative) tol=1e-5 in 
nls.control, there could be another arbitrary (absolute) convergence 
criterion.  Or, maybe there's something I don't understand about the 
algorithm being used.

Just my $0.02 and minority opinion,
efg

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


[R] Hierarchical clustering

2006-08-15 Thread Davendra Sohal
Hi,

I'm new to R; this is my second email to this forum. I have a dataset that
looks like this:

Gene AB
12536.25  2532.2
22527.35  4583.3



There are about 4 genes with 3-10 patients per dataset. Gene expression
values fill the matrix.

I want to do a simple hierarchical clustering along BOTH rows and columns
and get a tree dendrogram output for both in one panel.

Is there a simple code for this, leaving out fancy options?

If a heat map can be generated, that'll be even better.



Thank you very much.

-DS.

[[alternative HTML version deleted]]

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


Re: [R] nls convergence problem

2006-08-15 Thread Berton Gunter
   Or, maybe there's something I don't understand about the 
 algorithm being used.

Indeed! So before making such comments, why don't you try to learn about it?
Doug Bates is a pretty smart guy,  and I think you do him a disservice when
you assume that he somehow overlooked something that he explicitly warned
you about. I am fairly confident that if he could have made the problem go
away, he would have. So I think your vent was a bit inconsiderate and
perhaps even intemperate. The R Core folks have produced a minor miracle
IMO, and we should all be careful before assuming that they have overlooked
easily fixable problems. They're certainly not infallible -- but they're a
lot less fallible than most of the rest of us when it comes to R.

 
 Just my $0.02 and minority opinion,
 efg
 

... and mine.

-- Bert Gunter

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


[R] A model for possibly periodic data with varying amplitude [repost, much edited]

2006-08-15 Thread Andrew Robinson
Hi dear R community,

I have up to 12 measures of a protein for each of 6 patients, taken
every two or three days. The pattern of the protein looks periodic,
but the height of the peaks is highly variable.  I'm testing for
periodicity using a Monte Carlo simulation envelope approach applied
to a cumulative periodogram.  Now I want to predict the location of
the peaks in time.  Of course, the peaks might be occurring on
unmeasured days.

Sadly, an NDA prohibits me from sharing the real data. The data look
something like this:

##

patient - data.frame(
day = c(1, 3, 5, 8, 10, 12, 15, 17, 19, 22, 24, 26),
protein = c(5, 3, 10, 7, 2, 8, 25, 22, 7, 10, 12, 5)
)

plot(patient$day, patient$protein, type=b)

# This is my model:

wave.form -
  deriv3( ~ sin(2*pi*((day-offset)/period + 0.25)) * amplitude + mean,
 c(period, offset, amplitude, mean),
 function(day, period, offset, amplitude, mean){})

curve(wave.form(x, period=7, offset=2, mean=5, amplitude=4), 
from=1, to=30)

require(nlme)

wave.1 - gnls(protein ~ wave.form(day, period, offset, amplitude, mean),
   start=list(period=7, offset=0, amplitude=10, mean=5),
   weights=varPower(), data=patient)

wave.1

#

I emphasize that I don't care about the mean or amplitude, just the
offset and the period. The problem for my model is that spikes in the
data, such as at day 15, seem to make fitting the model quite
unstable, (although it is fine in this artificial case).  The best I
can think of right now is to use the weights=varPower() model in
gnls(), which I expect will down-weight the extreme spikes.  Often
that model fails to converge.  Then, all I can do is fall back on
modelling the log response.  If I do both, then when the weights model
converges the estimates sometimes differ, and sometimes considerably.

I know that I need more data, but more data are not available.  So,
instead I need to do the best I can :)

I have a half-baked idea of conditionally shrinking the peaks, so that
they all look like the same size - perhaps by discretizing the data
into, say, three possible values defined by quantiles of the response
variable or the estimated slope - but I'm concerned that doing so
risks creating the periodicity that I want to detect.  Also the series
is not necessarily stationary, so I wouild have to smooth it somehow
first.

On the other hand, maybe I can test for periodicity in the original
data, and then fit on some transformed data.  If I reject the null
hypothesis of no periodic oscillation, then is it reasonable to
condition subsequent analysis on the belief that periodicity exists,
and try to determine its characteristics?

I wonder if anyone can suggest an analysis technique or model that
more directly accommodates the varying heights of the peaks?  Or is
this pretty much the best I can do?

Any suggestions will be gratefully received.

Cheers,

Andrew
-- 
Andrew Robinson  
Department of Mathematics and StatisticsTel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
Email: [EMAIL PROTECTED] http://www.ms.unimelb.edu.au

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


Re: [R] nls convergence problem

2006-08-15 Thread Earl F. Glynn
Berton Gunter [EMAIL PROTECTED] wrote in message 
news:[EMAIL PROTECTED]
   Or, maybe there's something I don't understand about the
 algorithm being used.

 Indeed! So before making such comments, why don't you try to learn about 
 it?
 Doug Bates is a pretty smart guy,  and I think you do him a disservice 
 when
 you assume that he somehow overlooked something that he explicitly warned
 you about. I am fairly confident that if he could have made the problem go
 away, he would have. So I think your vent was a bit inconsiderate and
 perhaps even intemperate. The R Core folks have produced a minor miracle
 IMO, and we should all be careful before assuming that they have 
 overlooked
 easily fixable problems. They're certainly not infallible -- but they're a
 lot less fallible than most of the rest of us when it comes to R.

I meant no disrespect to Doug Bates or any of the R Core folks. I thought 
what I wrote had a neutral tone and was respectful.  I am sorry if anyone 
was offended by my comments and suggestions.  I am certainly thankful for 
all the hard work that has gone into developing R.

efg

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


[R] coefficients' order in polr()?

2006-08-15 Thread T Mu
Hi all,

I am using polr(). The resulting coefficients of first levels are always 0.

What to do if I wnat to get the coefficients of the last level 0.

For example, suppose x has 3 levels, 1, 2, 3

probit - plor(y ~ x, data1, method='probit')

will get coefficients of level 2, 3 of x, but I want coefficients of level
1, 2

Thank you,
Tian

[[alternative HTML version deleted]]

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


Re: [R] Looking for info on the Regression Modeling Strategies in R course in DC area

2006-08-15 Thread Frank E Harrell Jr
Tim McDonald wrote:
 Hello list,

   A colleague of mine mentioned a great course on  Regression Modeling 
 Strategies in R. Anyone knows if this course is offered as public course in 
 DC area?

   Thanks a bunch - TM

I'll be teaching this course in DC Sept 28-29 for XL Solutions. 
Information about the course may be obtained from 
biostat.mc.vanderbilt.edu/rms

Frank Harrell

 
   
 -
 
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] coefficients' order in polr()?

2006-08-15 Thread ronggui

you can use _relevel_ to Reorder Levels of Factor.
use ?relevel to get more information.

2006/8/16, T Mu [EMAIL PROTECTED]:

Hi all,

I am using polr(). The resulting coefficients of first levels are always 0.

What to do if I wnat to get the coefficients of the last level 0.

For example, suppose x has 3 levels, 1, 2, 3

probit - plor(y ~ x, data1, method='probit')

will get coefficients of level 2, 3 of x, but I want coefficients of level
1, 2

Thank you,
Tian

[[alternative HTML version deleted]]

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




--
黄荣贵
Department of Sociology
Fudan University

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


[R] Specifying Path Model in SEM for CFA

2006-08-15 Thread Rick Bilonick
I'm using specify.model for the sem package. I can't figure out how to
represent the residual errors for the observed variables for a CFA
model. (Once I get this working I need to add some further constraints.)

Here is what I've tried:

model.sa - specify.model()
  F1 - X1,l11, NA
  F1 - X2,l21, NA
  F1 - X3,l31, NA
  F1 - X4,l41, NA
  F1 - X5, NA, 0.20
  F2 - X1,l12, NA
  F2 - X2,l22, NA
  F2 - X3,l32, NA
  F2 - X4,l42, NA
  F2 - X6, NA, 0.25
  F1- F2,g12, 1
  F1- F1,g11, 1
  F2- F2,g22, 1
  X1- X1, NA, 1
  X2- X2, NA, 1
  X3- X3, NA, 1
  X4- X4, NA, 1
  X5- X5, NA, 1
  X6- X6, NA, 1

This at least converges:

 summary(fit.sem)

 Model Chisquare =  2147   Df =  10 Pr(Chisq) = 0
 Chisquare (null model) =  2934   Df =  15
 Goodness-of-fit index =  0.4822
 Adjusted goodness-of-fit index =  -0.087387
 RMSEA index =  0.66107   90 % CI: (NA, NA)
 Bentler-Bonnett NFI =  0.26823
 Tucker-Lewis NNFI =  -0.098156
 Bentler CFI =  0.26790
 BIC =  2085.1

 Normalized Residuals
   Min. 1st Qu.  MedianMean 3rd Qu.Max.
 -5.990  -0.618   0.192   0.165   1.700   3.950

 Parameter Estimates
Estimate  Std Error z value  Pr(|z|)
l11 -0.245981 0.21863   -1.12510 0.26054748 X1 --- F1
l21 -0.308249 0.22573   -1.36555 0.17207875 X2 --- F1
l31  0.202590 0.079102.56118 0.01043175 X3 --- F1
l41 -0.235156 0.21980   -1.06985 0.28468885 X4 --- F1
l12  0.839985 0.219623.82476 0.00013090 X1 --- F2
l22  0.828460 0.225483.67418 0.00023862 X2 --- F2
l32  0.066722 0.083690.79725 0.42530606 X3 --- F2
l42  0.832037 0.218403.80963 0.00013917 X4 --- F2
g12  0.936719 0.643311.45609 0.14536647 F2 -- F1
g11  2.567669 1.256082.04418 0.04093528 F1 -- F1
g22  1.208497 0.550402.19567 0.02811527 F2 -- F2

 Iterations =  59

And it produces the following path diagram:

 path.diagram(fit.sem)
digraph fit.sem {
  rankdir=LR;
  size=8,8;
  node [fontname=Helvetica fontsize=14 shape=box];
  edge [fontname=Helvetica fontsize=10];
  center=1;
  F2 [shape=ellipse]
  F1 [shape=ellipse]
  F1 - X1 [label=l11];
  F1 - X2 [label=l21];
  F1 - X3 [label=l31];
  F1 - X4 [label=l41];
  F1 - X5 [label=];
  F2 - X1 [label=l12];
  F2 - X2 [label=l22];
  F2 - X3 [label=l32];
  F2 - X4 [label=l42];
  F2 - X6 [label=];
}

But I don't see the residual error terms that go into each of the
observed variables X1 - X6. I've tried:

model.sa - specify.model()
  E1 - X1, e1,  1
  E2 - X2, e2,  1
  E3 - X3, e3,  1
  E4 - X4, e4,  1
  E5 - X5, e5,  1
  E6 - X6, e6,  1
  E1- E1, s1, NA
  E2- E2, s2, NA
  E3- E3, s3, NA
  E4- E4, s4, NA
  E5- E5, s5, NA
  E6- E6, s6, NA
  F1 - X1,l11, NA
  F1 - X2,l21, NA
  F1 - X3,l31, NA
  F1 - X4,l41, NA
  F1 - X5, NA,  1
  F2 - X1,l12, NA
  F2 - X2,l22, NA
  F2 - X3,l32, NA
  F2 - X4,l42, NA
  F2 - X6, NA,  1
  F1- F2, NA, 1
  F1- F1, NA, 1
  F2- F2,g22, NA
  X1- X1, NA, 1
  X2- X2, NA, 1
  X3- X3, NA, 1
  X4- X4, NA, 1
  X5- X5, NA, 1
  X6- X6, NA, 1

I'm trying to use E1 - E6 as the residual error terms. But I get warning
messages about no variances for X1-X6 and it doesn't converge. Also, the
associated path diagram:

digraph fit.sem {
  rankdir=LR;
  size=8,8;
  node [fontname=Helvetica fontsize=14 shape=box];
  edge [fontname=Helvetica fontsize=10];
  center=1;
  E1 [shape=ellipse]
  E2 [shape=ellipse]
  E3 [shape=ellipse]
  E4 [shape=ellipse]
  E5 [shape=ellipse]
  E6 [shape=ellipse]
  F2 [shape=ellipse]
  F1 [shape=ellipse]
  E1 - X1 [label=];
  E2 - X2 [label=];
  E3 - X3 [label=];
  E4 - X4 [label=];
  E5 - X5 [label=];
  E6 - X6 [label=];
  F1 - X1 [label=l11];
  F1 - X2 [label=l21];
  F1 - X3 [label=l31];
  F1 - X4 [label=l41];
  F1 - X5 [label=];
  F2 - X1 [label=l12];
  F2 - X2 [label=l22];
  F2 - X3 [label=l32];
  F2 - X4 [label=l42];
  F2 - X6 [label=];
}

Has ellipses around the E1-E6 which I believe indicates they are latent
factors and not residual errors.

If anyone could point in the right direction I would appreciate it.

Rick B.

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