date:20090119

It would be helpful to have a reproducible dataset to track down what
is happening.

On Mon, Jan 19, 2009 at 3:45 AM, Harsh singhal...@gmail.com wrote:
 Hello R List,
 I seem to have a peculiar problem. When using  time series data, I get
 the following error when running the acf and pacf function.
 Using the function acf(dtxts,plot= TRUE,xaxt = n,col=red,na.action
 = na.omit) (where dtxts is a time series object created with package
 xts ) results in the error below.

 Error in na.omit.ts(as.ts(x)) : time series contains internal NAs

 The above error is seen in R 2.8.0 running on Linux.

 The same function does not yield any error in R 2.8.0 on a Windows system.


 I've also tried na.remove(dtxts) from the tseries package to solve
 this problem but to no avail.

 Thank you.

 Harsh Singhal
 Decision Systems
 Mu Sigma Inc.,

 __
 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.




-- 
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.

Re: [R] time series contains internal NAs error

2009-01-19 Thread Harsh

Very Sorry for the oversight.

The dataset that I have used is:

Sunspots,Datefield
9.5,1/1/1900
2.7,1/1/1901
5,1/1/1902
24.4,1/1/1903
42,1/1/1904
63.5,1/1/1905
53.8,1/1/1906
62,1/1/1907
48.5,1/1/1908
43.9,1/1/1909
18.6,1/1/1910
5.7,1/1/1911
3.6,1/1/1912
1.4,1/1/1913
9.6,1/1/1914
47.4,1/1/1915
57.1,1/1/1916
103.9,1/1/1917
80.6,1/1/1918
63.6,1/1/1919
37.6,1/1/1920
26.1,1/1/1921
14.2,1/1/1922
5.8,1/1/1923
16.7,1/1/1924
44.3,1/1/1925
63.9,1/1/1926
69,1/1/1927
77.8,1/1/1928
64.9,1/1/1929
35.7,1/1/1930
21.2,1/1/1931
11.1,1/1/1932
5.7,1/1/1933
8.7,1/1/1934
36.1,1/1/1935
79.7,1/1/1936
114.4,1/1/1937
109.6,1/1/1938
88.8,1/1/1939
67.8,1/1/1940
47.5,1/1/1941
30.6,1/1/1942
16.3,1/1/1943
9.6,1/1/1944
33.2,1/1/1945
92.6,1/1/1946
151.6,1/1/1947
136.3,1/1/1948
134.7,1/1/1949
83.9,1/1/1950
69.4,1/1/1951
31.5,1/1/1952
13.9,1/1/1953
4.4,1/1/1954
38,1/1/1955
141.7,1/1/1956
190.2,1/1/1957
184.8,1/1/1958
159,1/1/1959
112.3,1/1/1960
53.9,1/1/1961
37.5,1/1/1962
27.9,1/1/1963
10.2,1/1/1964
15.1,1/1/1965
47,1/1/1966
93.8,1/1/1967
105.9,1/1/1968
105.5,1/1/1969
104.5,1/1/1970
66.6,1/1/1971
68.9,1/1/1972
38,1/1/1973
34.5,1/1/1974
15.5,1/1/1975
12.6,1/1/1976
27.5,1/1/1977
92.5,1/1/1978
155.4,1/1/1979
32.27,1/1/1980
54.25,1/1/1981
59.65,1/1/1982
63.62,1/1/1983

I convert the Datefield column to class Date, by using Datefield-
as.Date(Datefield,format=%m/%d/%Y)
Then I create a time series object by using
dtxts- xts(Sunspots,order.by = Datefield)

On using the ACF function
acf(dtxts,plot= TRUE,xaxt = n,col=red,na.action = na.omit)
I get the error
Error in na.omit.ts(as.ts(x)) : time series contains internal NAs

The peculiar issue is that I do not get the above error on
- R 2.7 on Ubuntu in India
- R 2.8 on Windows in India

The error appears when I execute this code on an R 2.8 implementation
running on a Fedora machine in the USA.

Is it a problem with the locale?

Thanks for the help. Really appreciate it.

Harsh Singhal
Decisions Systems
Mu Sigma Inc.
Chicago, USA

[[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] about power.law.fit

2009-01-19 Thread Gábor Csárdi

power.law.fit simply ML fits the 'prob(d) = d^\alpha' model to the
input, where d is positive integer. It seems to work for me:

 data - sample(1:1, prob=(1:1)^-3, rep=TRUE)
 power.law.fit(data)

Call:
mle(minuslogl = mlogl, start = list(alpha = start))

Coefficients:
   alpha
3.017056
 data - sample(1:1, prob=(1:1)^-2, rep=TRUE)
Warning message:
In sample(1:1, prob = (1:1)^-2, rep = TRUE) :
  Walker's alias method used: results are different from R  2.2.0
 power.law.fit(data)

Call:
mle(minuslogl = mlogl, start = list(alpha = start))

Coefficients:
   alpha
2.016645

It returns with an mle object, so you can call 'confint',  'logLik',
etc. on it, see mle-class for details.

 tmp - power.law.fit(data)
 summary(tmp)
Maximum likelihood estimation

Call:
mle(minuslogl = mlogl, start = list(alpha = start))

Coefficients:
  Estimate Std. Error
alpha 2.016645 0.01085921

-2 log L: 32150.62
 confint(tmp)
Profiling...
   2.5 %   97.5 %
1.995522 2.038091

Gabor

ps. there is also an igraph-help mailing list, FYI. Just in case I
miss your questions here

On Sun, Jan 18, 2009 at 4:50 PM, Weijia You weiji...@gmail.com wrote:
 Dear all,

 I'm using igraph for some analysis about the network I have. I have a
 question about the function power.law.fit.

 I wonder if there is any test for checking whether the  power.law.fit is
 good for the input, i.e., under which situation, could we use this function
 to get a reliable result. I'm afraid even I input a random graph without any
 property of power-law characteristics, it will returns an outcome which
 seems to be a fit to available data while it has no meaning to us. Is there
 any index like goodness of fit ?

 Thank you for any comments.

 Best!
 Weijia

[[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.




-- 
Gabor Csardi gabor.csa...@unil.ch UNIL DGM

__
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] reference category for factor in regression

2009-01-19 Thread ONKELINX, Thierry

Dear Jos,

In R you don't need to create you own dummy variables. Just create a
factor variable LABOUR (with two levels) and rerun your model. Then you
should be able to calculate all coefficients.

HTH,

Thierry



ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
thierry.onkel...@inbo.be 
www.inbo.be 

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey

-Oorspronkelijk bericht-
Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
Namens Jos Elkink
Verzonden: maandag 19 januari 2009 15:16
Aan: r-help@r-project.org
Onderwerp: [R] reference category for factor in regression

Hi all,

I am struggling with a strange issue in R that I have not encountered
before and I am not sure how to resolve this.

The model looks like this, with all irrelevant variables left out:

LABOUR - a dummy variable
NONLABOUR = 1 - LABOUR
AGE - a categorical variable / factor
VOTE - a dummy variable

glm(VOTE ~ 0 + LABOUR + NONLABOUR + LABOUR : AGE + NONLABOUR : AGE,
family=binomial(link=logit))

In other words, a standard interaction model, but I want to know the
intercepts and coefficients for each of the two cases (LABOUR and
NONLABOUR), instead of getting coefficients for the differences as in
a normal interaction model.

But the strange thing is, for the two occurances of the AGE variable,
it makes a different choice as to which AGE category to leave out of
the regression. The cross-table of AGE with LABOUR does not have empty
cells.

Anyone any idea what might be going wrong? Or what I could do about
this?

Thanks in advance for any help!

Regards,

Jos

-- 
Johan A. Elkink
Lecturer
School of Politics and International Relations  CHS Graduate School
University College Dublin
Ph. +353 1 716 7026  |  Library Building, Rm 512
http://jaeweb.cantr.net

__
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.

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.

__
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] time series contains internal NAs error

2009-01-19 Thread Achim Zeileis


On Mon, 19 Jan 2009, stephen sefick wrote:


It would be helpful to have a reproducible dataset to track down what
is happening.


True.

Although in this case it's relatively easy to guess what went wrong. 
The user probably has some irregular series, for example daily with 
missing days:

   x - xts(rnorm(10), as.Date(2009-01-19) + c(0:4, 7:11))


I seem to have a peculiar problem. When using  time series data, I get
the following error when running the acf and pacf function.
Using the function acf(dtxts,plot= TRUE,xaxt = n,col=red,na.action
= na.omit) (where dtxts is a time series object created with package
xts ) results in the error below.


acf() just works with regular time series of class ts. Everything else 
is coerced via as.ts(). This is where the problem is created for your data 
as the error message conveyed:



Error in na.omit.ts(as.ts(x)) : time series contains internal NAs


as.ts(x) detects that there is an underlying regularity (1-day steps) if 
you include two internal NAs (the weekend).


na.omit(as.ts(x)) cannot omit internal NAs because the ts class cannot 
represent such an object. Even if it could, acf() couldn't compute the ACF 
with internal mising data.


Some people just ignore the weekend effect (i.e., assume that the Mon-Fri 
correlation is the same as Tue-Mon). If one wants to do this, it can be 
obtained via acf(coredata(x)).



The above error is seen in R 2.8.0 running on Linux.

The same function does not yield any error in R 2.8.0 on a Windows system.


I'm fairly certain it does.
Z

__
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] time series contains internal NAs error

The statement to read in the data is missing from your post but
I suspect that you are representing the data as daily data so its filling
in 364 or 365 NA's between points.  Represent it as the annual data
that it is.  Try this:

Lines - Sunspots,Datefield
9.5,1/1/1900
2.7,1/1/1901
5,1/1/1902
24.4,1/1/1903
42,1/1/1904
63.5,1/1/1905
53.8,1/1/1906
62,1/1/1907
48.5,1/1/1908
43.9,1/1/1909
18.6,1/1/1910
5.7,1/1/1911
3.6,1/1/1912
1.4,1/1/1913
9.6,1/1/1914
47.4,1/1/1915
57.1,1/1/1916
103.9,1/1/1917
80.6,1/1/1918
63.6,1/1/1919
37.6,1/1/1920
26.1,1/1/1921
14.2,1/1/1922
5.8,1/1/1923
16.7,1/1/1924
44.3,1/1/1925
63.9,1/1/1926
69,1/1/1927
77.8,1/1/1928
64.9,1/1/1929
35.7,1/1/1930
21.2,1/1/1931
11.1,1/1/1932
5.7,1/1/1933
8.7,1/1/1934
36.1,1/1/1935
79.7,1/1/1936
114.4,1/1/1937
109.6,1/1/1938
88.8,1/1/1939
67.8,1/1/1940
47.5,1/1/1941
30.6,1/1/1942
16.3,1/1/1943
9.6,1/1/1944
33.2,1/1/1945
92.6,1/1/1946
151.6,1/1/1947
136.3,1/1/1948
134.7,1/1/1949
83.9,1/1/1950
69.4,1/1/1951
31.5,1/1/1952
13.9,1/1/1953
4.4,1/1/1954
38,1/1/1955
141.7,1/1/1956
190.2,1/1/1957
184.8,1/1/1958
159,1/1/1959
112.3,1/1/1960
53.9,1/1/1961
37.5,1/1/1962
27.9,1/1/1963
10.2,1/1/1964
15.1,1/1/1965
47,1/1/1966
93.8,1/1/1967
105.9,1/1/1968
105.5,1/1/1969
104.5,1/1/1970
66.6,1/1/1971
68.9,1/1/1972
38,1/1/1973
34.5,1/1/1974
15.5,1/1/1975
12.6,1/1/1976
27.5,1/1/1977
92.5,1/1/1978
155.4,1/1/1979
32.27,1/1/1980
54.25,1/1/1981
59.65,1/1/1982
63.62,1/1/1983

library(chron)
library(zoo)
z - read.zoo(textConnection(Lines), header = TRUE, sep = ,,
FUN = function(x) as.yearmon(chron(x)), index = 2)
acf(z)


On Mon, Jan 19, 2009 at 9:34 AM, Harsh singhal...@gmail.com wrote:
 Very Sorry for the oversight.

 The dataset that I have used is:

 Sunspots,Datefield
 9.5,1/1/1900
 2.7,1/1/1901
 5,1/1/1902
 24.4,1/1/1903
 42,1/1/1904
 63.5,1/1/1905
 53.8,1/1/1906
 62,1/1/1907
 48.5,1/1/1908
 43.9,1/1/1909
 18.6,1/1/1910
 5.7,1/1/1911
 3.6,1/1/1912
 1.4,1/1/1913
 9.6,1/1/1914
 47.4,1/1/1915
 57.1,1/1/1916
 103.9,1/1/1917
 80.6,1/1/1918
 63.6,1/1/1919
 37.6,1/1/1920
 26.1,1/1/1921
 14.2,1/1/1922
 5.8,1/1/1923
 16.7,1/1/1924
 44.3,1/1/1925
 63.9,1/1/1926
 69,1/1/1927
 77.8,1/1/1928
 64.9,1/1/1929
 35.7,1/1/1930
 21.2,1/1/1931
 11.1,1/1/1932
 5.7,1/1/1933
 8.7,1/1/1934
 36.1,1/1/1935
 79.7,1/1/1936
 114.4,1/1/1937
 109.6,1/1/1938
 88.8,1/1/1939
 67.8,1/1/1940
 47.5,1/1/1941
 30.6,1/1/1942
 16.3,1/1/1943
 9.6,1/1/1944
 33.2,1/1/1945
 92.6,1/1/1946
 151.6,1/1/1947
 136.3,1/1/1948
 134.7,1/1/1949
 83.9,1/1/1950
 69.4,1/1/1951
 31.5,1/1/1952
 13.9,1/1/1953
 4.4,1/1/1954
 38,1/1/1955
 141.7,1/1/1956
 190.2,1/1/1957
 184.8,1/1/1958
 159,1/1/1959
 112.3,1/1/1960
 53.9,1/1/1961
 37.5,1/1/1962
 27.9,1/1/1963
 10.2,1/1/1964
 15.1,1/1/1965
 47,1/1/1966
 93.8,1/1/1967
 105.9,1/1/1968
 105.5,1/1/1969
 104.5,1/1/1970
 66.6,1/1/1971
 68.9,1/1/1972
 38,1/1/1973
 34.5,1/1/1974
 15.5,1/1/1975
 12.6,1/1/1976
 27.5,1/1/1977
 92.5,1/1/1978
 155.4,1/1/1979
 32.27,1/1/1980
 54.25,1/1/1981
 59.65,1/1/1982
 63.62,1/1/1983

 I convert the Datefield column to class Date, by using Datefield-
 as.Date(Datefield,format=%m/%d/%Y)
 Then I create a time series object by using
 dtxts- xts(Sunspots,order.by = Datefield)

 On using the ACF function
 acf(dtxts,plot= TRUE,xaxt = n,col=red,na.action = na.omit)
 I get the error
 Error in na.omit.ts(as.ts(x)) : time series contains internal NAs

 The peculiar issue is that I do not get the above error on
 - R 2.7 on Ubuntu in India
 - R 2.8 on Windows in India

 The error appears when I execute this code on an R 2.8 implementation
 running on a Fedora machine in the USA.

 Is it a problem with the locale?

 Thanks for the help. Really appreciate it.

 Harsh Singhal
 Decisions Systems
 Mu Sigma Inc.
 Chicago, USA

[[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] lattice question: independent per-row or per-column scaling?

2009-01-19 Thread Chuck Cleland

On 1/19/2009 8:51 AM, hadley wickham wrote:
 On Thu, Jan 8, 2009 at 11:25 AM, René J.V. Bertin rjvber...@gmail.com wrote:
 Hello - and happy newyear to all of you!

 I've got some data that I'm plotting with bwplot, a 3x2x3 design where
 the observable decreases with the principle independent factor, but at
 different rates.

 I'd like to get lattice to impose not a single set of axes ranges
 identical for all panels, but ranges that are identical for each panel
 row or each column. Effects will stand out much better like that.

 I've looked through the documentation of the latest lattice version,
 but I don't see a way to achieve this with a simple argument passed to
 bwplot. Can it be done otherwise and if so, how?

  The argument for xlim or ylim can be a list.  Here is the key part of
the help page for xyplot:

xlim could also be a list, with as many components as the number of
panels (recycled if necessary), with each component as described above.
This is meaningful only when scales$x$relation is free or sliced, in
which case these are treated as if they were the corresponding limit
components returned by prepanel calculations.

  Here is a little example:

library(lattice)

mdf - data.frame(X1 = rep(LETTERS[1:3], each = 100*2*3),
  X2 = rep(c(J,K), 900),
  X3 = rep(LETTERS[24:26], 100*3*2),
  Y = c(runif(600, min=.01,max=.32),
runif(600, min=.33,max=.65),
runif(600, min=.66,max=.99)))

bwplot(Y ~ X3 | X2*X1, data = mdf,
   layout=c(2,3,1),
   ylim=as.data.frame(matrix(c(.01,.32,
   .01,.32,
   .33,.65,
   .33,.65,
   .66,.99,
   .66,.99), nrow=2)),
   scales=list(y=list(relation=free)))

 It's not lattice, but you can do this with ggplot2 - see the examples
 for http://had.co.nz/ggplot2/facet_grid.html
 
 Regards,
 
 Hadley 

-- 
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894

__
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] time series contains internal NAs error

2009-01-19 Thread Achim Zeileis


On Mon, 19 Jan 2009, Gabor Grothendieck wrote:


The statement to read in the data is missing from your post but
I suspect that you are representing the data as daily data so its filling
in 364 or 365 NA's between points.  Represent it as the annual data
that it is.


One further pointer: Also look at
  help(sunspots)


Try this:

Lines - Sunspots,Datefield
9.5,1/1/1900
2.7,1/1/1901
5,1/1/1902
24.4,1/1/1903
42,1/1/1904
63.5,1/1/1905
53.8,1/1/1906
62,1/1/1907
48.5,1/1/1908
43.9,1/1/1909
18.6,1/1/1910
5.7,1/1/1911
3.6,1/1/1912
1.4,1/1/1913
9.6,1/1/1914
47.4,1/1/1915
57.1,1/1/1916
103.9,1/1/1917
80.6,1/1/1918
63.6,1/1/1919
37.6,1/1/1920
26.1,1/1/1921
14.2,1/1/1922
5.8,1/1/1923
16.7,1/1/1924
44.3,1/1/1925
63.9,1/1/1926
69,1/1/1927
77.8,1/1/1928
64.9,1/1/1929
35.7,1/1/1930
21.2,1/1/1931
11.1,1/1/1932
5.7,1/1/1933
8.7,1/1/1934
36.1,1/1/1935
79.7,1/1/1936
114.4,1/1/1937
109.6,1/1/1938
88.8,1/1/1939
67.8,1/1/1940
47.5,1/1/1941
30.6,1/1/1942
16.3,1/1/1943
9.6,1/1/1944
33.2,1/1/1945
92.6,1/1/1946
151.6,1/1/1947
136.3,1/1/1948
134.7,1/1/1949
83.9,1/1/1950
69.4,1/1/1951
31.5,1/1/1952
13.9,1/1/1953
4.4,1/1/1954
38,1/1/1955
141.7,1/1/1956
190.2,1/1/1957
184.8,1/1/1958
159,1/1/1959
112.3,1/1/1960
53.9,1/1/1961
37.5,1/1/1962
27.9,1/1/1963
10.2,1/1/1964
15.1,1/1/1965
47,1/1/1966
93.8,1/1/1967
105.9,1/1/1968
105.5,1/1/1969
104.5,1/1/1970
66.6,1/1/1971
68.9,1/1/1972
38,1/1/1973
34.5,1/1/1974
15.5,1/1/1975
12.6,1/1/1976
27.5,1/1/1977
92.5,1/1/1978
155.4,1/1/1979
32.27,1/1/1980
54.25,1/1/1981
59.65,1/1/1982
63.62,1/1/1983

library(chron)
library(zoo)
z - read.zoo(textConnection(Lines), header = TRUE, sep = ,,
FUN = function(x) as.yearmon(chron(x)), index = 2)
acf(z)


On Mon, Jan 19, 2009 at 9:34 AM, Harsh singhal...@gmail.com wrote:

Very Sorry for the oversight.

The dataset that I have used is:

Sunspots,Datefield
9.5,1/1/1900
2.7,1/1/1901
5,1/1/1902
24.4,1/1/1903
42,1/1/1904
63.5,1/1/1905
53.8,1/1/1906
62,1/1/1907
48.5,1/1/1908
43.9,1/1/1909
18.6,1/1/1910
5.7,1/1/1911
3.6,1/1/1912
1.4,1/1/1913
9.6,1/1/1914
47.4,1/1/1915
57.1,1/1/1916
103.9,1/1/1917
80.6,1/1/1918
63.6,1/1/1919
37.6,1/1/1920
26.1,1/1/1921
14.2,1/1/1922
5.8,1/1/1923
16.7,1/1/1924
44.3,1/1/1925
63.9,1/1/1926
69,1/1/1927
77.8,1/1/1928
64.9,1/1/1929
35.7,1/1/1930
21.2,1/1/1931
11.1,1/1/1932
5.7,1/1/1933
8.7,1/1/1934
36.1,1/1/1935
79.7,1/1/1936
114.4,1/1/1937
109.6,1/1/1938
88.8,1/1/1939
67.8,1/1/1940
47.5,1/1/1941
30.6,1/1/1942
16.3,1/1/1943
9.6,1/1/1944
33.2,1/1/1945
92.6,1/1/1946
151.6,1/1/1947
136.3,1/1/1948
134.7,1/1/1949
83.9,1/1/1950
69.4,1/1/1951
31.5,1/1/1952
13.9,1/1/1953
4.4,1/1/1954
38,1/1/1955
141.7,1/1/1956
190.2,1/1/1957
184.8,1/1/1958
159,1/1/1959
112.3,1/1/1960
53.9,1/1/1961
37.5,1/1/1962
27.9,1/1/1963
10.2,1/1/1964
15.1,1/1/1965
47,1/1/1966
93.8,1/1/1967
105.9,1/1/1968
105.5,1/1/1969
104.5,1/1/1970
66.6,1/1/1971
68.9,1/1/1972
38,1/1/1973
34.5,1/1/1974
15.5,1/1/1975
12.6,1/1/1976
27.5,1/1/1977
92.5,1/1/1978
155.4,1/1/1979
32.27,1/1/1980
54.25,1/1/1981
59.65,1/1/1982
63.62,1/1/1983

I convert the Datefield column to class Date, by using Datefield-
as.Date(Datefield,format=%m/%d/%Y)
Then I create a time series object by using
dtxts- xts(Sunspots,order.by = Datefield)

On using the ACF function
acf(dtxts,plot= TRUE,xaxt = n,col=red,na.action = na.omit)
I get the error
Error in na.omit.ts(as.ts(x)) : time series contains internal NAs

The peculiar issue is that I do not get the above error on
- R 2.7 on Ubuntu in India
- R 2.8 on Windows in India

The error appears when I execute this code on an R 2.8 implementation
running on a Fedora machine in the USA.

Is it a problem with the locale?

Thanks for the help. Really appreciate it.

Harsh Singhal
Decisions Systems
Mu Sigma Inc.
Chicago, USA

   [[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-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] termplot

2009-01-19 Thread John Fox

Dear Bob,

Take a look at the effects package, described in
http://www.jstatsoft.org/v08/i15/paper.

I hope this helps,
 John

--
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On
 Behalf Of Robert Michael Inman
 Sent: January-19-09 8:17 AM
 To: r-help@r-project.org
 Subject: [R] termplot
 
 I have used glm and stepAIC to choose a best model.  I can use termplot to
 assess the contribution of each explanatory variable in the glm.  However
 the final model after running stepAIC includes interaction terms, and when
I
 do termplot I get Error in `[.data.frame`(mf, , i) : undefined columns
 selected.  I also see the termplot detail saying Nothing sensible
happens
 for interaction terms.
 
 Is it possilbe to determine explanatory power of final set of variables
when
 interactions are included?
 
 Thanks
 
 - Bob
 
 __
 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] Sweave encoding problem

2009-01-19 Thread Gerrit Voigt


Hello,
Sweave seems to have trouble processing german letters in R.
For example, my noweb R-input looks like this.
=
Oberflächenfehler = c(4, 11, 6, 2, 7, 9)
@
If I send it through Sweave, I get the following error message.

error:  chunk 1
Error in parse(text = chunk) : unexpected input in OberflÃ¤
extra: Warning message:
In readLines(f[1]) :
  underfull last line in C:\

(my R is in german, so I needed to translate the error message myself.)

I got the impression, that this is an encoding issue of Sweave, since  
the input typed into R directly works just fine. The encoding I use in  
my noweb document is utf8.


Thanks in advance
Gerrit

__
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] plot data with a colour scale - more details!

2009-01-19 Thread Rayhan Ahmed

Hi Jim,

This is clearly the function I need to use. However, I am new to Rplot
and coding.  I know that my output image is not what the data
indicate. My input file is a matrix, where the first row is the
residue number (10 values), and the second row is the hydrophathy
values (10 corresponding values ranging from 0.3-0.9).

At this point, the discontinous shading would be a polishing step, and
I would be just content with the color scale corresponding to the
range of hydropahty values.

Thanks.

-Ray



On Mon, Jan 19, 2009 at 8:55 AM, Ahmed, Rayhan raah...@utmb.edu wrote:


 
 From: Jim Lemon [...@bitwrit.com.au]
 Sent: Saturday, January 17, 2009 5:58 PM
 To: Ahmed, Rayhan
 Cc: r-help@r-project.org
 Subject: Re: [R] plot data with a colour scale - more details!

 Ahmed, Rayhan wrote:
 Hello again,

 I am trying to create a color scale plotting the hydroathy index (y-axis) 
 versus residue (x-axis) -  each residue (1-100) has a value between 0-1. 
 I've been trying to create a scale where:
 0-0.499: increasing intensity of red
 0.5- yellow
 0.51 - 1 increasing intensity of green.

 THis is in lieu of a line graph.

 I'm not sure which type of plot to begin with - I attempted to modify a heat 
 plot, but it did not work.


 Hi Rayhan,
 I think that color2Dmatplot might do what you want. Check the last
 example that calls color.scale twice to see how to get discontinuous
 ranges of color for the values.

 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.

Re: [R] reference category for factor in regression

2009-01-19 Thread Jos Elkink

Hi Thierry,

Thanks for your quick answer. The problem is not so much the LABOUR
variable, however, but the AGE variable, which consists of about 5
categories for which I do indeed not create separate dummy variables.
But R does not behave as expected when deciding on which dummy to use
as reference category ...

Jos

On Mon, Jan 19, 2009 at 2:37 PM, ONKELINX, Thierry
thierry.onkel...@inbo.be wrote:
 Dear Jos,

 In R you don't need to create you own dummy variables. Just create a
 factor variable LABOUR (with two levels) and rerun your model. Then you
 should be able to calculate all coefficients.

 HTH,

 Thierry

 
 
 ir. Thierry Onkelinx
 Instituut voor natuur- en bosonderzoek / Research Institute for Nature
 and Forest
 Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
 methodology and quality assurance
 Gaverstraat 4
 9500 Geraardsbergen
 Belgium
 tel. + 32 54/436 185
 thierry.onkel...@inbo.be
 www.inbo.be

 To call in the statistician after the experiment is done may be no more
 than asking him to perform a post-mortem examination: he may be able to
 say what the experiment died of.
 ~ Sir Ronald Aylmer Fisher

 The plural of anecdote is not data.
 ~ Roger Brinner

 The combination of some data and an aching desire for an answer does not
 ensure that a reasonable answer can be extracted from a given body of
 data.
 ~ John Tukey

 -Oorspronkelijk bericht-
 Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
 Namens Jos Elkink
 Verzonden: maandag 19 januari 2009 15:16
 Aan: r-help@r-project.org
 Onderwerp: [R] reference category for factor in regression

 Hi all,

 I am struggling with a strange issue in R that I have not encountered
 before and I am not sure how to resolve this.

 The model looks like this, with all irrelevant variables left out:

 LABOUR - a dummy variable
 NONLABOUR = 1 - LABOUR
 AGE - a categorical variable / factor
 VOTE - a dummy variable

 glm(VOTE ~ 0 + LABOUR + NONLABOUR + LABOUR : AGE + NONLABOUR : AGE,
 family=binomial(link=logit))

 In other words, a standard interaction model, but I want to know the
 intercepts and coefficients for each of the two cases (LABOUR and
 NONLABOUR), instead of getting coefficients for the differences as in
 a normal interaction model.

 But the strange thing is, for the two occurances of the AGE variable,
 it makes a different choice as to which AGE category to leave out of
 the regression. The cross-table of AGE with LABOUR does not have empty
 cells.

 Anyone any idea what might be going wrong? Or what I could do about
 this?

 Thanks in advance for any help!

 Regards,

 Jos

 --
 Johan A. Elkink
 Lecturer
 School of Politics and International Relations  CHS Graduate School
 University College Dublin
 Ph. +353 1 716 7026  |  Library Building, Rm 512
 http://jaeweb.cantr.net

 __
 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.

 Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer
 en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd 
 is
 door een geldig ondertekend document. The views expressed in  this message
 and any annex are purely those of the writer and may not be regarded as 
 stating
 an official position of INBO, as long as the message is not confirmed by a 
 duly
 signed document.




-- 
Johan A. Elkink
Lecturer
School of Politics and International Relations  CHS Graduate School
University College Dublin
Ph. +353 1 716 7026  |  Library Building, Rm 512
http://jaeweb.cantr.net

__
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] Mac OS X / preview.app / fullrefman.pdf

2009-01-19 Thread Berend Hasselman



mmuurr[AT]gmail.com wrote:
 
 .
 has anyone else experienced this?  i sent in a bug report to Apple,
 but i doubt i'll see any change in preview.app prior to the next OS
 release.  fullrefman.pdf is also the only document that i've ever
 observed this behavior with, and since i can now replicate it on two
 (very) different machines (one a PowerPC G5 PowerMac, one an Intel-
 based MacBook), i'm convinced it is more than a one-time anomaly
 (which is what i originally chalked it up to).
 
 

I have had the same problem. AFAICR, it started with Leopard.
I have also sent bug reports to Apple. Nothing has happened yet.

Nowadays I use Skim to view fullrefman.pdf and that works well
(http://skim-app.sourceforge.net/)

Berend


-- 
View this message in context: 
http://www.nabble.com/Mac-OS-X---preview.app---fullrefman.pdf-tp21535696p21545778.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] Sweave encoding problem

2009-01-19 Thread Rau, Roland

Hi Gerrit,

 -Original Message-
 From: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org] On Behalf Of Gerrit Voigt
 Sent: Monday, January 19, 2009 4:48 PM
 To: r-help@r-project.org
 Subject: [R] Sweave encoding problem

 Hello,
 Sweave seems to have trouble processing german letters in R.
 For example, my noweb R-input looks like this.
 =
 Oberflächenfehler = c(4, 11, 6, 2, 7, 9)
 @
 If I send it through Sweave, I get the following error message.

 error:  chunk 1
 Error in parse(text = chunk) : unexpected input in OberflÃ¤
 extra: Warning message:
 In readLines(f[1]) :
underfull last line in C:\

 (my R is in german, so I needed to translate the error 
 message myself.)

 I got the impression, that this is an encoding issue of 
 Sweave, since  
 the input typed into R directly works just fine. The encoding 
 I use in  
 my noweb document is utf8.

I don't think it has something to do with German letters.
I saved the following text in a file 'sweavy.Snw':
\documentclass{article}

\begin{document}
Hello World!

=
1+1
@ 

=
Oberflächenfehler = c(4, 11, 6, 2, 7, 9)
@
\end{document}

This is what happened in R:
 library(utils)
 Sweave(sweavy.Snw)
Writing to file sweavy.tex
Processing code chunks ...
 1 : echo term verbatim
 2 : echo term verbatim

You can now run LaTeX on 'sweavy.tex'

 sessionInfo()
R version 2.7.0 (2008-04-22) 
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 datasets  methods   base 

And also the dvi looked fine after processing latex sweavy.tex
To make things sure, I did in my editor (GNU Emacs 22.1.50.1)
C-x RET f utf-8
to change 
set-buffer-file-coding-system to utf-8.
Still works fine.

Maybe this helps you further to track down the reason for the problem?!?

Best,
Roland

--
This mail has been sent through the MPI for Demographic Research.  Should you 
receive a mail that is apparently from a MPI user without this text displayed, 
then the address has most likely been faked. If you are uncertain about the 
validity of this message, please check the mail header or ask your system 
administrator for assistance.

__
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] WinBUGS with R

2009-01-19 Thread Lindsay Stirton


Dear UseRs,

I am having some problems using R with WinBUGS using the R2WinBUGS
package. Specifically, when I try to run bugs() I get the following
message.

Error in FUN(X[[1L]], ...) :
   .C(..): 'type' must be real for this format




To give a little more context, my bugs() command (for a multilevel
ordinal logit  similar to Gelman and Hill, Data Analysis Using
Regression and Multilevel/Hierarchical Models p. 383 is:

Wednesbury.data - list (n.judge, n, n.cut, y judge, ct,
ra, lg)

Wednesbury.inits - function(){
 list(C=matrix(0,39,2))
  }
Wednesbury.parameters - c(C, b1, b2, b3)

__
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] reference category for factor in regression

2009-01-19 Thread Berwin A Turlach

G'day Jos,

On Mon, 19 Jan 2009 15:52:00 +
Jos Elkink jos.elk...@ucd.ie wrote:

 Thanks for your quick answer. The problem is not so much the LABOUR
 variable, however, but the AGE variable, which consists of about 5
 categories for which I do indeed not create separate dummy variables.
 But R does not behave as expected when deciding on which dummy to use
 as reference category ...

I guess in this case we need more information.

What is the output of str(AGE)?  That should tell you which category
is the reference category.  And which category do you expect to be the
reference category?

Cheers,

Berwin

=== Full address =
Berwin A TurlachTel.: +65 6515 4416 (secr)
Dept of Statistics and Applied Probability+65 6515 6650 (self)
Faculty of Science  FAX : +65 6872 3919   
National University of Singapore
6 Science Drive 2, Blk S16, Level 7  e-mail: sta...@nus.edu.sg
Singapore 117546http://www.stat.nus.edu.sg/~statba

__
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] reference category for factor in regression

2009-01-19 Thread Stephan Kolassa


Hi Jos,

you can force R to set contrasts for factors the way you like them with 
contrasts(). You seem to be thinking of treatment contrasts, which are 
most easily interpreted, but there are also others.


However: are you sure you want to bin an age variable into categories? 
You will lose power, along with a lot of other unpleasant things:

http://biostat.mc.vanderbilt.edu/twiki/bin/view/Main/CatContinuous

With five categories, you are giving up 4 df. I'd recommend looking into 
splines, where you should be able to get more bang for the buck. Look at 
rcs() in the Design package, and at Frank Harrell's excellent book 
Regression Modeling Strategies.


Of course, if you only have the binned data, all this is irrelevant...

HTH,
Stephan


Jos Elkink schrieb:

Hi Thierry,

Thanks for your quick answer. The problem is not so much the LABOUR
variable, however, but the AGE variable, which consists of about 5
categories for which I do indeed not create separate dummy variables.
But R does not behave as expected when deciding on which dummy to use
as reference category ...

Jos

On Mon, Jan 19, 2009 at 2:37 PM, ONKELINX, Thierry
thierry.onkel...@inbo.be wrote:

Dear Jos,

In R you don't need to create you own dummy variables. Just create a
factor variable LABOUR (with two levels) and rerun your model. Then you
should be able to calculate all coefficients.

HTH,

Thierry



ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
thierry.onkel...@inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey

-Oorspronkelijk bericht-
Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
Namens Jos Elkink
Verzonden: maandag 19 januari 2009 15:16
Aan: r-help@r-project.org
Onderwerp: [R] reference category for factor in regression

Hi all,

I am struggling with a strange issue in R that I have not encountered
before and I am not sure how to resolve this.

The model looks like this, with all irrelevant variables left out:

LABOUR - a dummy variable
NONLABOUR = 1 - LABOUR
AGE - a categorical variable / factor
VOTE - a dummy variable

glm(VOTE ~ 0 + LABOUR + NONLABOUR + LABOUR : AGE + NONLABOUR : AGE,
family=binomial(link=logit))

In other words, a standard interaction model, but I want to know the
intercepts and coefficients for each of the two cases (LABOUR and
NONLABOUR), instead of getting coefficients for the differences as in
a normal interaction model.

But the strange thing is, for the two occurances of the AGE variable,
it makes a different choice as to which AGE category to leave out of
the regression. The cross-table of AGE with LABOUR does not have empty
cells.

Anyone any idea what might be going wrong? Or what I could do about
this?

Thanks in advance for any help!

Regards,

Jos

--
Johan A. Elkink
Lecturer
School of Politics and International Relations  CHS Graduate School
University College Dublin
Ph. +353 1 716 7026  |  Library Building, Rm 512
http://jaeweb.cantr.net

__
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.

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message
and any annex are purely those of the writer and may not be regarded as stating
an official position of INBO, as long as the message is not confirmed by a duly
signed document.







__
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] bootstrapped eigenvector method following prcomp

2009-01-19 Thread Stas Kolenikov

I don't know if there are bugs in the code, but the step 4) does not
compute significance... at least the way statisticians know it. The
fractions above or below 0 are not significance. I don't even know how
to call those... probably cdf of the bootstrap distribution evaluated
at zero.

Let's put the flies and the sausages separately, as a Russian saying
goes. If you bootstrap from your original data, you can get the
confidence intervals, but not the test statistics. What you can do
with your bootstrap from the raw data is accumulate the distribution
of the eigenvectors, along the lines of (assuming you are testing for
the significance of variables in your first component):

function (x, permutations=1000, ...)
{
  pcnull - princomp(x, ... )
  res - pcnull$loadings[,1]
  bsresults = matrix( rep.int(NA, permutations*NROW(res)) ,
nrow=permutations, ncol=NROW(res) )
  N - nrow(x)
  for (i in 1:permutations) {
  pc - princomp(x[sample(N, replace=TRUE), ], ... )
  pred - predict(pc, newdata = x)
  r -  cor(pcnull$scores, pred)
  k - apply(abs(r), 2, which.max)
  reve - sign(diag(r[k,]))
  sol - pc$loadings[ ,k]
  sol - sweep(sol, 2, reve, *)
  bsresults[i,] - t(sol[,1])
  }
  apply( bsresults, 2, quantile, c(0.05, 0.95) )
}

if I am not messing up the dimensions and other stuff too much.
However as a result you will get an approximately degenerate
distribution sitting on the unit sphere since the eigenvectors are
always normalized to have the length of 1. You can still do marginal
confidence intervals with that though, and see if 0 is covered.

The main problem here is I am not entirely sure the bootstrap is
applicable for the problem at hand. In other words, it is not clear
whether the bootstrap estimates are consistent for the true
variability. Eigenproblems are quite difficult and prone to
non-regularities (the above mentioned degeneracy is just one of them,
and probably not the worst one). There are different asymptotics
(Anderson's of fixed p and n \to \infty,
http://www.citeulike.org/user/ctacmo/article/3908837, and Johnstone
joint p and n \to \infty with p/n \to const,
http://www.citeulike.org/user/ctacmo/article/3908846), the
distributions are often non-standard, while the bootstrap works well
when the distributions of the estimates are normal. When you have
something weird, bootstrap may easily break down, and in a lot of
other situations, you need to come up with special schemes. See my
oh-so-favorite paper on the bootstrap,
http://www.citeulike.org/user/ctacmo/article/575126.

One of those special schemes (back to out muttons, or rather flies and
sausages) -- to set up the bootstrap for hypothesis testing and get
the p-values, you need to bootstrap from the distribution that
corresponds to the null. Beran and Srivastava (1985 Annals,
http://www.citeulike.org/user/ctacmo/article/3015345) discuss how to
rotate your data to conform to the null hypothesis of interest for
inference with covariance matrices and their functions (such as
eigenvalues, for instance). Whether you need to go into all this
trouble, I don't really know.

If you have an inferential problem of testing whether a particular
variable contributes to the overall index, and have a pretty good idea
of where each variable goes, may be you need to shift your paradigm
and look at confirmatory factor analysis models instead, estimable in
R with John Fox' sem package.

On 1/19/09, Axel Strauß a.stra...@tu-bs.de wrote:
 G'Day R users!

  Following an ordination using prcomp, I'd like to test which variables
 singnificantly contribute to a principal component. There is a method
 suggested by Peres-Neto and al. 2003. Ecology 84:2347-2363  called
 bootstrapped eigenvector. It was asked for that in this forum in January
 2005 by Jérôme Lemaître:
  1) Resample 1000 times with replacement entire raws from the original data
 sets []
  2) Conduct a PCA on each bootstrapped sample
  3) To prevent axis reflexion and/or axis reordering in the bootstrap, here
 are two more steps for each bootstrapped sample
  3a) calculate correlation matrix between the PCA scores of the original and
 those of the bootstrapped sample
  3b) Examine whether the highest absolute correlation is between the
 corresponding axis for the original and bootstrapped samples. When it is not
 the case, reorder the eigenvectors. This means that if the highest
 correlation is between the first original axis and the second bootstrapped
 axis, the loadings for the second bootstrapped axis and use to estimate the
 confidence interval for the original first PC axis.
  4) Determine the p value for each loading. Obtained as follow: number of
 loadings =0 for loadings that were positive in the original matrix divided
 by the number of boostrap samples (1000) and/or number of loadings =0 for
 loadings that were negative in the original matrix divided by the number of
 boostrap samples (1000).

  (see
 https://stat.ethz.ch/pipermail/r-help/2005-January/065139.html

Re: [R] reference category for factor in regression

2009-01-19 Thread Marc Schwartz

Jos,

See ?relevel for information on how to reorder the levels of a factor,
while being able to specify the reference level.

Basically, the first level of the factor is taken as the reference. If
you want to utilize a different ordering, as an alternative to the
above, simply use:

  AGE - factor(AGE, levels = c(FirstLevel, SecondLevel, ...)

BTW, you might want to review Frank Harrell's page on why categorizing a
continuous variable is not a good idea:

  http://biostat.mc.vanderbilt.edu/twiki/bin/view/Main/CatContinuous

HTH,

Marc Schwartz


on 01/19/2009 09:52 AM Jos Elkink wrote:
 Hi Thierry,
 
 Thanks for your quick answer. The problem is not so much the LABOUR
 variable, however, but the AGE variable, which consists of about 5
 categories for which I do indeed not create separate dummy variables.
 But R does not behave as expected when deciding on which dummy to use
 as reference category ...
 
 Jos
 
 On Mon, Jan 19, 2009 at 2:37 PM, ONKELINX, Thierry
 thierry.onkel...@inbo.be wrote:
 Dear Jos,

 In R you don't need to create you own dummy variables. Just create a
 factor variable LABOUR (with two levels) and rerun your model. Then you
 should be able to calculate all coefficients.

 HTH,

 Thierry

 
 
 ir. Thierry Onkelinx
 Instituut voor natuur- en bosonderzoek / Research Institute for Nature
 and Forest
 Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
 methodology and quality assurance
 Gaverstraat 4
 9500 Geraardsbergen
 Belgium
 tel. + 32 54/436 185
 thierry.onkel...@inbo.be
 www.inbo.be

 To call in the statistician after the experiment is done may be no more
 than asking him to perform a post-mortem examination: he may be able to
 say what the experiment died of.
 ~ Sir Ronald Aylmer Fisher

 The plural of anecdote is not data.
 ~ Roger Brinner

 The combination of some data and an aching desire for an answer does not
 ensure that a reasonable answer can be extracted from a given body of
 data.
 ~ John Tukey

 -Oorspronkelijk bericht-
 Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
 Namens Jos Elkink
 Verzonden: maandag 19 januari 2009 15:16
 Aan: r-help@r-project.org
 Onderwerp: [R] reference category for factor in regression

 Hi all,

 I am struggling with a strange issue in R that I have not encountered
 before and I am not sure how to resolve this.

 The model looks like this, with all irrelevant variables left out:

 LABOUR - a dummy variable
 NONLABOUR = 1 - LABOUR
 AGE - a categorical variable / factor
 VOTE - a dummy variable

 glm(VOTE ~ 0 + LABOUR + NONLABOUR + LABOUR : AGE + NONLABOUR : AGE,
 family=binomial(link=logit))

 In other words, a standard interaction model, but I want to know the
 intercepts and coefficients for each of the two cases (LABOUR and
 NONLABOUR), instead of getting coefficients for the differences as in
 a normal interaction model.

 But the strange thing is, for the two occurances of the AGE variable,
 it makes a different choice as to which AGE category to leave out of
 the regression. The cross-table of AGE with LABOUR does not have empty
 cells.

 Anyone any idea what might be going wrong? Or what I could do about
 this?

 Thanks in advance for any help!

 Regards,

 Jos


__
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] ifelse help?

2009-01-19 Thread rkevinburton

I am having a hard time understanding what is happening with ifelse.

Let me illustrate:

h - numeric(5)
p - 1:5
j - floor(j)
x - 1:1000
 ifelse(h == 0, x[j+2], 1:5)
[1] 2 3 4 5 6

My question is, shouldn't this be retruning 25 numbers? It seems that the 
ifelse should check 5 values of h for zero. For each of the 5 values I am 
thinking it should return an array of 5 (x[j+2] if h[index] == 0). Since the 
dimension of h is 5 that would mean 25 values. But that isn't what is being 
returned.Something about this that I don't understand. Please help my ignorance.

Thank you.

Kevin

__
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] Using apply to generate matrix from rows?

2009-01-19 Thread Stephan Lindner

Dear all,


I have a simple question which I unfortunately do not seem to be able
to solve myself. I have a (NxK) matrix and want to generate a new
matrix by multiplying each row with itself such that the new matrix
has dimension ((N*K)xK) (or better, generate an array with dimension
(K,K,N)). I tried apply, but that did not work. Any suggestions? 

Thanks! 


Stephan



## Here is a simple example: 

u - matrix(c(1,2,3,4,5,6,7,8,9,10),nrow=2)

## What I want to obtain 

u[1,]%*%t(u[1,])
u[2,]%*%t(u[2,])

## stacked together -- 10x5 matrix


## This does not work 

sq - function(x)x%*%t(x)
apply(u,1,function(y)sq(y))






-- 
---
Stephan Lindner
University of Michigan

__
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] lazy evaluation question

2009-01-19 Thread Peter Dalgaard


Gabor Grothendieck wrote:

Note that

rm(i)
for(j in 1:4) F(j)

raises an error due to scoping issues.


Yes. This has nothing to do with lazy evaluation, and everything to do 
with scoping: f is not defined in the scope of F, so does not know about 
its variables (nor those in the implicit loop of lapply()).


Notice also that in

lapply(1:4,function(i) F(i))

it would be pretty weird if lapply would behave differently depending on 
the name of formal arguments of the function, i.e. if


lapply(1:4,function(meep) F(meep))

gave a different result. And f() depends on looking for a variable i 
outside of the function.



On Sun, Jan 18, 2009 at 10:02 PM,  markle...@verizon.net wrote:

I've been going back to old difficult R-list evaluation emails that I save
 in order to understand evaluation better and below still confuses me. Could
someone explain why A) works and B) doesn't. A variant of below is in the
 Pat's Inferno book also but I'm still not clear on what is happening.
Thanks.

f - function() {
 # FORCING i here doesn't help
 i*i
}

F - function(i) {
 force(i)
 print(f())
 }

A) THIS WORKS
for ( i in 1:4 ) {
 F(i)
}

B) THIS DOESN'T
lapply(1:4,function(i) F(i))



--
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - (p.dalga...@biostat.ku.dk)  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] lazy evaluation question

On Mon, Jan 19, 2009 at 12:40 PM, Peter Dalgaard
p.dalga...@biostat.ku.dk wrote:
 Gabor Grothendieck wrote:

 Note that

 rm(i)
 for(j in 1:4) F(j)

 raises an error due to scoping issues.

 Yes. This has nothing to do with lazy evaluation, and everything to do with
 scoping: f is not defined in the scope of F, so does not know about its
 variables (nor those in the implicit loop of lapply()).

 Notice also that in

 lapply(1:4,function(i) F(i))

 it would be pretty weird if lapply would behave differently depending on the
 name of formal arguments of the function, i.e. if

 lapply(1:4,function(meep) F(meep))

 gave a different result. And f() depends on looking for a variable i outside
 of the function.


I believe that is what I was already referring to when I referred to
scoping issues.
Since not works was never defined in the original post its hard to
know what the
problem being asked is but it might have been confusion over the
different uses of i.

__
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] Using apply to generate matrix from rows?

2009-01-19 Thread Henrique Dallazuanna

Try this:

matrix(apply(u, 1, tcrossprod), nr = nrow(u)*ncol(u), byrow = T)

On Mon, Jan 19, 2009 at 3:39 PM, Stephan Lindner lindn...@umich.edu wrote:

 Dear all,


 I have a simple question which I unfortunately do not seem to be able
 to solve myself. I have a (NxK) matrix and want to generate a new
 matrix by multiplying each row with itself such that the new matrix
 has dimension ((N*K)xK) (or better, generate an array with dimension
 (K,K,N)). I tried apply, but that did not work. Any suggestions?

 Thanks!


Stephan



 ## Here is a simple example:

 u - matrix(c(1,2,3,4,5,6,7,8,9,10),nrow=2)

 ## What I want to obtain

 u[1,]%*%t(u[1,])
 u[2,]%*%t(u[2,])

 ## stacked together -- 10x5 matrix


 ## This does not work

 sq - function(x)x%*%t(x)
 apply(u,1,function(y)sq(y))






 --
 ---
 Stephan Lindner
 University of Michigan

 __
 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

[[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] Using apply to generate matrix from rows?

2009-01-19 Thread Jorge Ivan Velez

Dear Stephan,
Try this:

do.call(rbind,lapply(1:2,function(x) matrix(u[x,]%*%t(u[x,]),ncol=ncol(u

HTH,

Jorge


On Mon, Jan 19, 2009 at 12:39 PM, Stephan Lindner lindn...@umich.eduwrote:

 Dear all,


 I have a simple question which I unfortunately do not seem to be able
 to solve myself. I have a (NxK) matrix and want to generate a new
 matrix by multiplying each row with itself such that the new matrix
 has dimension ((N*K)xK) (or better, generate an array with dimension
 (K,K,N)). I tried apply, but that did not work. Any suggestions?

 Thanks!


Stephan



 ## Here is a simple example:

 u - matrix(c(1,2,3,4,5,6,7,8,9,10),nrow=2)

 ## What I want to obtain

 u[1,]%*%t(u[1,])
 u[2,]%*%t(u[2,])

 ## stacked together -- 10x5 matrix


 ## This does not work

 sq - function(x)x%*%t(x)
 apply(u,1,function(y)sq(y))






 --
 ---
 Stephan Lindner
 University of Michigan

 __
 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.


[[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] download/retain text file structure with RCurl/getURL()

2009-01-19 Thread zack holden


Dear list,
 
I'm trying to download a text file directly from the internet using the RCurl 
package and the command getURL. Duncan Lang graciously helped me solve the 
first step in this problem using the following command:

# 
txtfile - 
getURL('ftp://ftp.wcc.nrcs.usda.gov/data/snow/snow_course/table/history/idaho/13e19.txt',
 
ftp.use.epsv = FALSE)
#
 
This brings the text file into R in a single long character string. I've spent 
many hours now trying to bring this text file into R into a sensible form. I've 
tried every variant of different commands in getURL help file, as well as 
different 
strsplit() commands to try to break this character string into a sensible rows 
and columns, to no avail. 
 
Can anyone suggest a solution for doing this? I suspect there is a getURL 
command I'm missing. Alternatively, do I really have to break this long 
character string into rows and columns that I can then assemble into a table? 
 
I'd be grateful for any advice.
 
Thanks in advance, 
 
Zack

 
__
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] ifelse help?

2009-01-19 Thread Charles C. Berry


On Mon, 19 Jan 2009, rkevinbur...@charter.net wrote:


I am having a hard time understanding what is happening with ifelse.

Let me illustrate:

h - numeric(5)
p - 1:5
j - floor(j)


And j is 0:4 + epsilon , where 0 = epsilon  1, evidently.


x - 1:1000
ifelse(h == 0, x[j+2], 1:5)
[1] 2 3 4 5 6

My question is, shouldn't this be retruning 25 numbers?



No. It should be

A vector of the same length and attributes (including class) as
test and data values from the values of yes or no

according to ?ifelse, and 'test' is what you have as 'h'

Consider


z - as.data.frame( diag(2 ) )
 ifelse(z == 0, letters , 1:5)

 V1  V2
[1,] 1 c
[2,] b 4




all args have different lengths and classes. But the result has those of 
the 'test' arg.



It seems that the ifelse should check 5 values of h for zero. For each 
of the 5 values I am thinking it should return an array of 5 (x[j+2] if 
h[index] == 0). Since the dimension of h is 5 that would mean 25 values.


h has no attributes, therefore no 'dimension'

HTH,

Chuck

But that isn't what is being returned.Something about this that I don't 
understand. Please help my ignorance.


Thank you.

Kevin

__
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
Dept of Family/Preventive Medicine
E mailto:cbe...@tajo.ucsd.edu   UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

__
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] download/retain text file structure with RCurl/getURL()

If you are having problems with the default download.file method you
can try method = wget:

f -
ftp://ftp.wcc.nrcs.usda.gov/data/snow/snow_course/table/history/idaho/13e19.txt;
download.file(f, basename(f), method = wget)

On Mon, Jan 19, 2009 at 1:26 PM, zack holden zack_hol...@hotmail.com wrote:

Dear list,

I'm trying to download a text file directly from the internet using the RCurl
package and the command getURL. Duncan Lang graciously helped me solve the
first step in this problem using the following command:

#
txtfile -
getURL('ftp://ftp.wcc.nrcs.usda.gov/data/snow/snow_course/table/history/idaho/13e19.txt',
ftp.use.epsv = FALSE)
#

This brings the text file into R in a single long character string. I've
spent many hours now trying to bring this text file into R into a sensible
form. I've tried every variant of different commands in getURL help file, as
well as different
strsplit() commands to try to break this character string into a sensible
rows and columns, to no avail.

Can anyone suggest a solution for doing this? I suspect there is a getURL
command I'm missing. Alternatively, do I really have to break this long
character string into rows and columns that I can then assemble into a table?

I'd be grateful for any advice.

Thanks in advance,

Zack

__
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] gz netCDF files

2009-01-19 Thread Magdalena Lucini

Hello,

I am trying to access several netCDF files that are also zipped (via
gzip I guess) (and stored in a directory that I only have reading
permit)

I tried to unzip them using gzfile, gzcon, etc, and then open them
with open.ncdf (from ncdf package). Everything was unsuccesful.

I had no problem using the ncdf package with normal netCDF files

I would really appreciate any suggestion

Thanks

Magdalena Lucini
Dept.  of Physics
U of Toronto

__
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] lazy evaluation question

2009-01-19 Thread Wacek Kusnierczyk

Peter Dalgaard wrote:

snip
 Notice also that in

 lapply(1:4,function(i) F(i))

 it would be pretty weird if lapply would behave differently depending
 on the name of formal arguments of the function, i.e. if

 lapply(1:4,function(meep) F(meep))

 gave a different result. And f() depends on looking for a variable i
 outside of the function.

here's one example:

d = data.frame(a=1, b=2)

lapply(3, function(a) subset(d, select=a))
lapply(3, function(b) subset(d, select=b))
lapply(3, function(c) subset(d, select=c))

vQ

__
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] candisc

2009-01-19 Thread Pete Shepard

Hello,

I have a question regarding the candisc package. My data are:

speciesthreefive
12.956.63
12.537.79
13.575.65
13.165.47
22.584.46
22.166.22
23.273.52

I put these in a table and then a linear model
 newdata - lm(cbind(three, five) ~ species, data=rawdata)

and then do a candisc on them
 candata-candisc(newdata)

Here are my scores;
candata$scores

  species   Can1
1   1 -2.3769280
2   1 -2.7049437
3   1 -3.4748309
4   1 -0.9599825
5   2  4.2293774
6   2  2.6052193
7   2  2.6820884

and here are my coefficients
 candata$coeffs.raw
   Can1
three -5.185380
five  -2.160237
 candata$coeffs.std
   Can1
three -2.530843
five  -2.586620


My question is, what is the precise equation that gives the candata$scores?

Thanks

[[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] Error while adding legends to xyplot

2009-01-19 Thread V Prasanth

Dear All:

 

Greetings!

 

I am able to produce an xyplot in R; But I am not able to put multiple
legends on it! So for that matter, I have saved the xyplot and reproduced
the same using the simple plot option. Then using the legend option I
successfully placed the require text and its corresponding value over it.
Now, when I try to copy it as a Meta file, it shows error and R terminates
abruptly. So I will be grateful if somebody could guide me to overcome this
since I am interested in saving the graphs as a Meta file than any other
file!

 

Below furnished is the R Code that I have used for the same.

 

CorrCoeff - 0

Corr - cor.test(FCV, FCV_Fitted)

CorrCoeff - as.numeric(round(Corr$estimate,4))

myPlot - xyplot(FCV + FCV_Fitted ~ Year, data=Tobacco,
col=c(orange,blue4), freq=TRUE, 

main=FCV Tobacco Production  Forecasted Figures Over The Years
: India, font=2, type=b,

auto.key = list (type = , points=F, lines = F, border = F,
text = c(Observed Production Figures, Fitted Values), font=2,

col=c(orange,blue4), x = 0, y = 0.95), xlab=list(Year,
font=2), ylab=list('000 Metric Tones, font=2))

old.prompt - grid::grid.prompt(TRUE)

 

detach(package:lattice)

 

plot.new()

par(new=TRUE, font=4)

plot(myPlot, xlab = , ylab = , ann=FALSE, xaxt=n, yaxt=n)

legend(x=bottomright, bty=n, lty=c(0), col=c(black),

title= (Correlation Coefficient), legend=c(CorrCoeff))

 

 

Many Thanks in advance,

Prasanth VP,

Global Manager - Biometrics, 

Delta Technology  Management Services Pvt Ltd, 

Plot No: 13/2, Sector - I, 

Third Floor, HUDA Techno Enclave, 

Madhapur, Hyderabad - 500 033. 

Office: +91-40-3024 9508. 
Mobile : +91-9848 290025 

 http://www.deltaintech.com/ www.deltaintech.com

 


**

'The information contained in this email is confidential and may contain
proprietary information. It is meant solely for the intended recipient.
Access to this email by anyone else is unauthorized. If you are not the
intended recipient, any disclosure, copying, distribution or any action
taken or omitted in reliance on this, is prohibited and may be unlawful. No
liability or responsibility is accepted if information or data is, for
whatever reason corrupted or does not reach its intended recipient. No
warranty is given that this email is free of viruses. The views expressed in
this email are, unless otherwise stated, those of the author and not those
of DELTA Technology and Management Services pvt ltd or its management. DELTA
Technology and Management Services pvt ltd reserves the right to monitor
intercept and block emails addressed to its users or take any other action
in accordance with its email use policy' 

Thank you in advance for your cooperation.


**

P Please don't print this e-mail unless you really need to.

 


[[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] Perl-R bridge

2009-01-19 Thread Adam Witney



you could take a look at this:

http://www.omegahat.org/RSPerl/

but I'm not sure how well maintained it is currently.

adam

On 19 Jan 2009, at 02:00, ANJAN PURKAYASTHA wrote:


Hi,
I'm planning to access R from my perl scripts.
The only noteworthy bridge seems to be
Statistics-R-0.03http://search.cpan.org/%7Ectbrown/Statistics-R/lib/Statistics/R.pm 
.

Would anyone like to share their experience with this Perl-R bridge?
I'd like to install it in a Mac OS X.
Suggestions on alternate solutions will be appreciated.
Thanks in advance,
Anjan

--
=
anjan purkayastha, phd
bioinformatics analyst
whitehead institute for biomedical research
nine cambridge center
cambridge, ma 02142

purkayas [at] wi [dot] mit [dot] edu
703.740.6939

[[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] predict.tree

2009-01-19 Thread H.H.H.B.M.van_Haagen

Dear Brian,

 

I looked on the internet for this warning message Warning message:

'newdata' had X rows but variable(s) found have Y rows that happens
when using some kind of fitting and prediting. I still can't figure out
what is going wrong by reading the documentation or the posts in de
forum. So I made a small example

 

pos = rep(1,times=10)

neg = rep(0,times=10)

 

y = c(pos,neg)

x = c(pos*30,neg*2)

 

result1 - glm(y~x, family=binomial)

 

p = predict.glm(result1,as.data.frame(x),type='respons')  

 

 

This works fine. But when I change the argument of the last row into
as.data.frame(x[1:10])  (or some new data) I get that warning again.
Why? Can you help me out? 

 

Best regards

Herman van Haagen

(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] using R how to read a one column alone from a database table from MySQL

2009-01-19 Thread sankar82

hello dieter,

sorry for the late reply...

yes my homework needs to be done in R and thats y i have these questions
posted in this forum...hope u reply me the solutions and also thank you for
spending ur precious time in my problem...

BR
sankar.

Dieter Menne wrote:

sankar82 sankar.arughadhoss at tkk.fi writes:

i have a created a database table in MYSQL consisting of 11 columns.
throught RMYSQL i managed to read the entire table in R. but i have few
qureries which i need solutions...here they are:

1. Using R how to read a one column alone from a database table from
MYSQL.
2. Using R how to print on screen those column value.
3. Using R how to print one particular row (in this case row is X)value
alone.
4. Using R how to read all the column (11) and print (X) value alone on
screen.
5. Using R with logic print those (X) value which is/between say 20 to 30
degrees.

Have you tried the examples in Rmysql-package?
For 5., you could use something like:

dbSendQuery(con, select * from WL where width\_nm between 0.5 and 1)

Ok, here we let SQL do the job, and you homework want it to be done in R.

Dieter

__
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.

--
View this message in context:
http://www.nabble.com/using-R-how-to-read-a-one-column-alone-from-a-database-table-from-MySQL-tp21478288p21543180.html
Sent from the R help mailing list archive at Nabble.com.

Re: [R] lattice question: independent per-row or per-column scaling?

2009-01-19 Thread René J.V. Bertin

Thanks for all the answers. I'll have a look at ggplot2. I'd seen the
possibility to set panel-specific limits via ylim, but I was in fact
looking for a switch to achieve non-global automatic scaling.

Given the fact that there is no built in provision for that, I take it
it's a functionality that isn't (considered) very interesting for most
lattice users, but out of curiosity, how difficult would it be to
implement it?

René

INRETS - LEPSIS

__
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] download/retain text file structure with RCurl/getURL()

2009-01-19 Thread David Winsemius

It's a fixed width format, with irregular entries, perhaps something  
along the lines of:


read.fwf(textConnection(txtfile), skip = 8, # skips the header
 widths = column widths vector,
 colnames= colnames ,
 nrows=48 )#drops the trailing summary text

perhaps :

widths = c(2, -1, 1, -1 ,4, -1, 3    the rest  # the -col  
entries drop the white-space
names = c(year,card, Jan.date, Jan.dep  .  
the rest


Just the first few columns seem to come in acceptably, although the  
lines with all NA's will need to be deleted:

 read.fwf(textConnection(txtfile), skip = 8, # skips the header
+ widths = c(2, -1, 1, -1 ,4, -1, 3),  # the -col entries drop  
the white-space
+ col.names = c(year,card, Jan.date, Jan.dep), 
nrows=48 )

   year card Jan.date Jan.dep
1611 E/ST  NA
2621 E/ST  NA
3631 K/31  15
4641 K/30  12
5NA   NA NA  NA
6651 E/ST  NA
7661 1/07  17
8671 E/ST  NA
9681 K/28  12
10   691 K/31  22
11   NA   NA NA  NA
12   701 K/30  16
13   711 K/29  28
14   721 K/28  32
15   731 1/02  16
snip
--
David Winsemius

On Jan 19, 2009, at 1:26 PM, zack holden wrote:



Dear list,

I'm trying to download a text file directly from the internet using  
the RCurl package and the command getURL. Duncan Lang graciously  
helped me solve the first step in this problem using the following  
command:


#
txtfile - getURL('ftp://ftp.wcc.nrcs.usda.gov/data/snow/snow_course/table/history/idaho/13e19.txt' 
,

ftp.use.epsv = FALSE)
#

This brings the text file into R in a single long character string.  
I've spent many hours now trying to bring this text file into R into  
a sensible form. I've tried every variant of different commands in  
getURL help file, as well as different
strsplit() commands to try to break this character string into a  
sensible rows and columns, to no avail.


Can anyone suggest a solution for doing this? I suspect there is a  
getURL command I'm missing. Alternatively, do I really have to break  
this long character string into rows and columns that I can then  
assemble into a table?


I'd be grateful for any advice.

Thanks in advance,

Zack


__
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] lattice question: independent per-row or per-column scaling?

2009-01-19 Thread René J.V. Bertin

Thanks for all the answers. I'll have a look at ggplot2. I'd seen the
possibility to set panel-specific limits via ylim, but I was in fact
looking for a switch to achieve non-global automatic scaling.

Given the fact that there is no built in provision for that, I take it
it's a functionality that isn't (considered) very interesting for most
lattice users, but out of curiosity, how difficult would it be to
implement it?

René

INRETS - LEPSIS

__
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] Trend.spatial function in geoR

2009-01-19 Thread Julia L. Angstmann

I am having difficulty getting the trend.spatial function in geoR to
work properly.  After creating a trend.spatial object with a covariate,
I try to add the command into my likfit() function as follows:

 

trend1.trend.spatial - trend.spatial(1st, trend1.geodata)

trend1.spatial.EC0.1.reml - likfit(spatial.geodata,
trend1.trend.spatial, ini.cov.pars = spatial.EC0.1.eyefit, fix.kappa =
FALSE, fix.psiA = FALSE, fix.psiR = FALSE,  

+ lik.method = REML, components = FALSE)

 

I receive an error stating: trend matrix has dimension incompatible with
the data

 

My data files have an equal number of rows or spatial locations.  Any
suggestions would be useful.

 

Also, is there a way to add multiple covariates into the trend.spatial
function?

 

Thanks,

 

*

Julia L. Angstmann

Department of Botany

University of Wyoming

1000 E. University Ave.

Laramie, WY 82071

 


[[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] ifelse help?

2009-01-19 Thread rkevinburton

Sorry I didn't give the proper initialization of j. But you are right j should 
also be an array of 5. So x[j + 5] would return 5 values. 

So if the array returned from 'ifelse' is the same dimention as test (h), then 
are all the values of h being tested? So since h as you say has no dimensions 
is the test only testing h[1]? Again it seems that if all of the elements of h 
are tested (there are 5 elements) and each element produces an array of 5 the 
resulting array should be 25. 

Kevin

 Charles C. Berry cbe...@tajo.ucsd.edu wrote: 
 On Mon, 19 Jan 2009, rkevinbur...@charter.net wrote:
 
  I am having a hard time understanding what is happening with ifelse.
 
  Let me illustrate:
 
  h - numeric(5)
  p - 1:5
  j - floor(j)
 
 And j is 0:4 + epsilon , where 0 = epsilon  1, evidently.
 
  x - 1:1000
  ifelse(h == 0, x[j+2], 1:5)
  [1] 2 3 4 5 6
 
  My question is, shouldn't this be retruning 25 numbers?
 
 
 No. It should be
 
   A vector of the same length and attributes (including class) as
   test and data values from the values of yes or no
 
 according to ?ifelse, and 'test' is what you have as 'h'
 
 Consider
 
  z - as.data.frame( diag(2 ) )
   ifelse(z == 0, letters , 1:5)
   V1  V2
 [1,] 1 c
 [2,] b 4
 
 
 all args have different lengths and classes. But the result has those of 
 the 'test' arg.
 
 
  It seems that the ifelse should check 5 values of h for zero. For each 
  of the 5 values I am thinking it should return an array of 5 (x[j+2] if 
  h[index] == 0). Since the dimension of h is 5 that would mean 25 values.
 
 h has no attributes, therefore no 'dimension'
 
 HTH,
 
 Chuck
 
  But that isn't what is being returned.Something about this that I don't 
  understand. Please help my ignorance.
 
  Thank you.
 
  Kevin
 
  __
  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
  Dept of Family/Preventive 
 Medicine
 E mailto:cbe...@tajo.ucsd.edu UC San Diego
 http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901
 


__
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] plotting arrows with different colors and varying head size

2009-01-19 Thread Héctor Villalobos

Dear list,

I would like to plot arrows with different colors according to arrow length, 
and also (if
possible) with head size proportional to arrow length. The idea is to make a 
quiver-like plot of
matlab with wind speed data.

So far, I´ve been able to use different colors, but I need to find a more 
efficient way to recode
arrow length intervals into colors. On the contrary, I can't define different 
head sizes, because
the length argument in the arrows() function seems to control the head size 
of all the
arrows at once.

Any help will be greatly appreciated.



## Generate random data
 set.seed(1)
 xcoor - matrix(runif(20), ncol=5)
 ycoor - matrix(runif(20), ncol=5)
   u - matrix(rnorm(20)/10, ncol=5)
   v - matrix(rnorm(20)/10, ncol=5)

## calculate arrows length and look histogram
 arrowlen - sqrt(u^2 + v^2)
 hist(arrowlen)

## recode arrow lengths (by interval) into colors
## sorry, I don't know how to do it without the car  package
 library(car)
 arrowL - recode(as.vector(arrowlen),
 0='grey90';
 0:0.05='grey50';
   0.05:0.1='grey';
   0.1:0.15='cyan';
   0.15:0.2='blue';
   0.2:0.25='red')

length=0.1

par.uin - function()
  # determine scale of inches/userunits in x and y
  # from http://tolstoy.newcastle.edu.au/R/help/01c/2714.html
  # Brian Ripley Tue 20 Nov 2001 - 20:13:52 EST
 {
u - par(usr)
p - par(pin)
c(p[1]/(u[2] - u[1]), p[2]/(u[4] - u[3]))
  }


## plot arrows
 plot(as.vector(xcoor), as.vector(ycoor), type=p, pch=., xlim=c(-0.2, 1.3),
  ylim=c(-0.2, 1.3))
   arrows( xcoor, ycoor, xcoor + u, ycoor + v,
  length = length*as.vector(arrowlen)*min(par.uin()), col=arrowL)



  sessionInfo()
R version 2.8.1 (2008-12-22)
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 datasets  methods
[7] base

other attached packages:
[1] car_1.2-9  hdf5_1.6.7



--
Héctor Villalobos hvill...@ipn.mx
 CICIMAR - IPN
  La Paz, Baja California Sur, MÉXICO


[[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] reference category for factor in regression

2009-01-19 Thread Jos Elkink

Hi all,

Thanks for the advice.

 See ?relevel for information on how to reorder the levels of a factor,
 while being able to specify the reference level.
 Basically, the first level of the factor is taken as the reference.

Yes, that is how I always used it. But the problem is, in this
particular regression R does *not* take the first level as reference.
In fact, AGE appears twice in the same regression (two different
interactions) and in one case it selects the 1st category and in
another case a different one.

 BTW, you might want to review Frank Harrell's page on why categorizing a
 continuous variable is not a good idea:

I most certainly agree, but the categorisation has been imposed in the
survey itself, so it is all the data I have. I did not design the
questions :-) ... Thanks for this reference, though, as it is
certainly interesting to inform my teaching.

 str(AGE)
 Factor w/ 5 levels 65+,18-24,..: 5 5 1 4 5 5 2 4 1 3 ...

So I expect 65+ to be the reference category, but it is not.

Here is a little bit more R code to show the problem:

 str(AGE)
 Factor w/ 5 levels 65+,18-24,..: 5 5 1 4 5 5 2 4 1 3 ...
 table(LABOUR)
LABOUR
   01
 692 1409
 NONLABOUR - 1 - LABOUR
 m - glm(NOVOTE ~ 0 + LABOUR + NONLABOUR + AGE : LABOUR + AGE : NONLABOUR, 
 family=binomial)
 m

Call:  glm(formula = NOVOTE ~ 0 + LABOUR + NONLABOUR + AGE:LABOUR +
  AGE:NONLABOUR, family = binomial)

Coefficients:
LABOUR   NONLABOUR   LABOUR:AGE65+ LABOUR:AGE18-24
  -0.35110-0.30486-0.11890-0.66444
   LABOUR:AGE25-34 LABOUR:AGE35-49 LABOUR:AGE50-64  NONLABOUR:AGE18-24
  -0.23893-0.15860  NA-0.65655
NONLABOUR:AGE25-34  NONLABOUR:AGE35-49  NONLABOUR:AGE50-64
  -0.72815 0.04951 0.17481

As you can see, 65+ is taken as reference category in the interaction
with NONLABOUR, but not in the interaction with LABOUR.

I know glm(NOVOTE ~ LABOUR * AGE, family=binomial) would be a more
conventional specification, but the above should be equivalent and
should give me the coefficients and standard errors for the two groups
(LABOUR and NONLABOUR) separately, rather than for the difference /
interaction term).

Perhaps the NA in the above output (which I only notice now) is a hint
at the problem, but I am not sure why that occurs.

 table(m$model$AGE, m$model$LABOUR, m$model$NOVOTE)
, ,  = 0


  0   1
  65+   137  24
  18-24  68 127
  25-34  59 267
  35-49  71 298
  50-64  82 179

, ,  = 1


  0   1
  65+   101  15
  18-24  26  46
  25-34  21 148
  35-49  55 179
  50-64  72 126

Anyone any idea? So there must be a reason R decides *not* to use 65+
as reference in that particular scenario, and I am missing why.

Jos

__
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] Compare matrices

2009-01-19 Thread Andrej Kastrin


Dear all,

Suppose that I have a matrix A

   A - matrix(c(3,3,3,3,3,3,3,3,3),3,3)

and a logical matrix B

   B - matrix(c(T,T,T,F,T,T,F,T,F),3,3)

The result  matrix should be

   C - matrix(c(3,3,3,NA,3,3,NA,3,NA),3,3)

Is there any simple tip or trick to perform this without looping?

Thanks in advance for any suggestion.

Best regards, Andrej

__
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] Compare matrices

2009-01-19 Thread Dimitris Rizopoulos


try this:

A - matrix(c(3,3,3,3,3,3,3,3,3),3,3)
B - matrix(c(T,T,T,F,T,T,F,T,F),3,3)

C - A
C[!B] - NA
C


I hope it helps.

Best,
Dimitris


Andrej Kastrin wrote:

Dear all,

Suppose that I have a matrix A

   A - matrix(c(3,3,3,3,3,3,3,3,3),3,3)

and a logical matrix B

   B - matrix(c(T,T,T,F,T,T,F,T,F),3,3)

The result  matrix should be

   C - matrix(c(3,3,3,NA,3,3,NA,3,NA),3,3)

Is there any simple tip or trick to perform this without looping?

Thanks in advance for any suggestion.

Best regards, Andrej

__
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.



--
Dimitris Rizopoulos
Assistant Professor
Department of Biostatistics
Erasmus Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014

__
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] reference category for factor in regression

2009-01-19 Thread Stephan Kolassa


Hi Jos,

does explicitly recoding AGE help?

AGE - factor(c(65+,18-24,18-24,25-34))
str(AGE)
AGE - 
factor(c(65+,18-24,18-24,25-34),levels=c(65+,18-24,25-34))

str(AGE)

Best,
Stephan


Jos Elkink schrieb:

Hi all,

Thanks for the advice.


See ?relevel for information on how to reorder the levels of a factor,
while being able to specify the reference level.
Basically, the first level of the factor is taken as the reference.


Yes, that is how I always used it. But the problem is, in this
particular regression R does *not* take the first level as reference.
In fact, AGE appears twice in the same regression (two different
interactions) and in one case it selects the 1st category and in
another case a different one.


BTW, you might want to review Frank Harrell's page on why categorizing a
continuous variable is not a good idea:


I most certainly agree, but the categorisation has been imposed in the
survey itself, so it is all the data I have. I did not design the
questions :-) ... Thanks for this reference, though, as it is
certainly interesting to inform my teaching.


str(AGE)

 Factor w/ 5 levels 65+,18-24,..: 5 5 1 4 5 5 2 4 1 3 ...

So I expect 65+ to be the reference category, but it is not.

Here is a little bit more R code to show the problem:


str(AGE)

 Factor w/ 5 levels 65+,18-24,..: 5 5 1 4 5 5 2 4 1 3 ...

table(LABOUR)

LABOUR
   01
 692 1409

NONLABOUR - 1 - LABOUR
m - glm(NOVOTE ~ 0 + LABOUR + NONLABOUR + AGE : LABOUR + AGE : NONLABOUR, 
family=binomial)
m


Call:  glm(formula = NOVOTE ~ 0 + LABOUR + NONLABOUR + AGE:LABOUR +
  AGE:NONLABOUR, family = binomial)

Coefficients:
LABOUR   NONLABOUR   LABOUR:AGE65+ LABOUR:AGE18-24
  -0.35110-0.30486-0.11890-0.66444
   LABOUR:AGE25-34 LABOUR:AGE35-49 LABOUR:AGE50-64  NONLABOUR:AGE18-24
  -0.23893-0.15860  NA-0.65655
NONLABOUR:AGE25-34  NONLABOUR:AGE35-49  NONLABOUR:AGE50-64
  -0.72815 0.04951 0.17481

As you can see, 65+ is taken as reference category in the interaction
with NONLABOUR, but not in the interaction with LABOUR.

I know glm(NOVOTE ~ LABOUR * AGE, family=binomial) would be a more
conventional specification, but the above should be equivalent and
should give me the coefficients and standard errors for the two groups
(LABOUR and NONLABOUR) separately, rather than for the difference /
interaction term).

Perhaps the NA in the above output (which I only notice now) is a hint
at the problem, but I am not sure why that occurs.


table(m$model$AGE, m$model$LABOUR, m$model$NOVOTE)

, ,  = 0


  0   1
  65+   137  24
  18-24  68 127
  25-34  59 267
  35-49  71 298
  50-64  82 179

, ,  = 1


  0   1
  65+   101  15
  18-24  26  46
  25-34  21 148
  35-49  55 179
  50-64  72 126

Anyone any idea? So there must be a reason R decides *not* to use 65+
as reference in that particular scenario, and I am missing why.

Jos

__
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] Compare matrices



On 20/01/2009, at 9:48 AM, Andrej Kastrin wrote:


Dear all,

Suppose that I have a matrix A

A - matrix(c(3,3,3,3,3,3,3,3,3),3,3)

and a logical matrix B

B - matrix(c(T,T,T,F,T,T,F,T,F),3,3)

The result  matrix should be

C - matrix(c(3,3,3,NA,3,3,NA,3,NA),3,3)

Is there any simple tip or trick to perform this without looping?


(a) It is helpful to state a clear question rather than expecting the
reader to infer from an example just what your question is.

(b) It is advisable to use TRUE and FALSE as the values of a logical
variate, rather than T and F.  The former are reserved words, the latter
are not, which can cause all hell to break loose.

(c) You are wasting a good many key strokes.  You could simply say

A - matrix(3,3,3)

(d) I infer that what you want is for C to equal A where B is TRUE
and NA otherwise.  One way to accomplish this is:

C - matrix(,3,3)
C[B] - A[B]

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] Compare matrices

2009-01-19 Thread Steven McKinney


Use the is.na() function to assign
NA values:


 is.na(A) - !B
 A
 [,1] [,2] [,3]
[1,]3   NA   NA
[2,]333
[3,]33   NA
 C - matrix(c(3,3,3,NA,3,3,NA,3,NA),3,3)
 all.equal(A, C)
[1] TRUE





Steven McKinney

Statistician
Molecular Oncology and Breast Cancer Program
British Columbia Cancer Research Centre

email: smckinney +at+ bccrc +dot+ ca

tel: 604-675-8000 x7561

BCCRC
Molecular Oncology
675 West 10th Ave, Floor 4
Vancouver B.C. 
V5Z 1L3
Canada




-Original Message-
From: r-help-boun...@r-project.org on behalf of Dimitris Rizopoulos
Sent: Mon 1/19/2009 12:54 PM
To: Andrej Kastrin
Cc: r-help@r-project.org
Subject: Re: [R] Compare matrices
 
try this:

A - matrix(c(3,3,3,3,3,3,3,3,3),3,3)
B - matrix(c(T,T,T,F,T,T,F,T,F),3,3)

C - A
C[!B] - NA
C


I hope it helps.

Best,
Dimitris


Andrej Kastrin wrote:
 Dear all,
 
 Suppose that I have a matrix A
 
A - matrix(c(3,3,3,3,3,3,3,3,3),3,3)
 
 and a logical matrix B
 
B - matrix(c(T,T,T,F,T,T,F,T,F),3,3)
 
 The result  matrix should be
 
C - matrix(c(3,3,3,NA,3,3,NA,3,NA),3,3)
 
 Is there any simple tip or trick to perform this without looping?
 
 Thanks in advance for any suggestion.
 
 Best regards, Andrej
 
 __
 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.
 

-- 
Dimitris Rizopoulos
Assistant Professor
Department of Biostatistics
Erasmus Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014

__
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] Compare matrices

2009-01-19 Thread David Winsemius



On Jan 19, 2009, at 3:54 PM, Dimitris Rizopoulos wrote:


try this:

A - matrix(c(3,3,3,3,3,3,3,3,3),3,3)
B - matrix(c(T,T,T,F,T,T,F,T,F),3,3)

C - A
C[!B] - NA
C



Very elegant. Another, perhaps less elegant, effort:

B[which(B == FALSE)] - NA
 B
 [,1] [,2] [,3]
[1,] TRUE   NA   NA
[2,] TRUE TRUE TRUE
[3,] TRUE TRUE   NA
 C - matrix(A * B, 3,3)   # A * B is *not* matrix multiplication


 C
 [,1] [,2] [,3]
[1,]3   NA   NA
[2,]333
[3,]33   NA

--
David Winsemius


I hope it helps.

Best,
Dimitris


Andrej Kastrin wrote:

Dear all,
Suppose that I have a matrix A
  A - matrix(c(3,3,3,3,3,3,3,3,3),3,3)
and a logical matrix B
  B - matrix(c(T,T,T,F,T,T,F,T,F),3,3)
The result  matrix should be
  C - matrix(c(3,3,3,NA,3,3,NA,3,NA),3,3)
Is there any simple tip or trick to perform this without looping?
Thanks in advance for any suggestion.
Best regards, Andrej
__
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.


--
Dimitris Rizopoulos
Assistant Professor
Department of Biostatistics
Erasmus Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014

__
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] Fitting of lognormal distribution to lower tail experimental data

2009-01-19 Thread Mattias Brännström

Thank you very much, Göran!
I had to install R 2.8.1 since it did not work with 2.4.1.
This is exactly what I wanted, now I can move on with my analysis! (And
learn more about cencoring...)

Best regards,
Mattias

__
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] maptools, sunriset, POSIX timezones

2009-01-19 Thread Phil Taylor


Hi ...

I wonder if anyone can provide some insight into why the first three 
examples using the sunriset function (appended below, with results) give 
the correct answer, but the fourth generates and error.


The first two use ISOdatetime with and without a time zone attribute, 
and the sunriset function returns the correct sunset time.


The third and fourth adds 10 seconds to the ISOdatetime (with and 
without the time zone attribute) but the function only works when the 
time zone is specified (example 3).


When I look at the objects (print or str) they appear the same, and when 
I check to see if they are equivalent (e.g.)


 time.1 - ISOdatetime(1970, 1, 1, 10, 0, 0) + 10
 time.2 - ISOdatetime(1970, 1, 1, 10, 0, 0, tz=GMT) + 10
 time.1 == time.2
[1] TRUE

they appear to be the same.

I wonder if I am either missing something important, doing something 
improperly, or if there is a small bug somewhere.


I'm using windows Vista and R 2.8.1

Thanks,

Phil Taylor
ptay...@resalliance.org

 require(maptools)
 Sys.setenv(TZ = GMT)
 location - matrix(c(-80.1,42.5), nrow=1)

 sunriset(location, ISOdatetime(1970, 1, 1, 10, 0, 0, tz=GMT), 
direction=sunset, POSIXct.out=TRUE)

  day_fractime
1 0.915226 1970-01-01 21:57:55

 sunriset(location, ISOdatetime(1970, 1, 1, 10, 0, 0), 
direction=sunset, POSIXct.out=TRUE)

  day_fractime
1 0.915226 1970-01-01 21:57:55

 sunriset(location, ISOdatetime(1970, 1, 1, 10, 0, 0, tz=GMT) + 10, 
direction=sunset, POSIXct.out=TRUE)

  day_fractime
1 0.915226 1970-01-01 21:57:55

 sunriset(location, ISOdatetime(1970, 1, 1, 10, 0, 0) + 10, 
direction=sunset, POSIXct.out=TRUE)
Error in structure(.Internal(as.POSIXct(x, tz)), class = c(POSIXt, 
POSIXct),  :

invalid 'tz' value

__
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] Compare matrices



On 20/01/2009, at 10:05 AM, David Winsemius (who should know better)  
wrote:


snip


B[which(B == FALSE)] - NA


snip

This sort of syntax drives me, and all right-thinking people,  
subclinically neurotic

(or as a psychiatrist would say, stark staring bonkers).

Admittedly it's not (quite) as bad as using ``Y==TRUE'' where ``Y''
is what is wanted, but it's pretty bad.

One should say ``B[!B] - NA''.

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] plotting arrows with different colors and varying head size

2009-01-19 Thread Greg Snow

Look at the my.symbols function in the TeachingDemos package (along with the 
ms.arrows function in the same package), that may do what you want.

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of Héctor Villalobos
 Sent: Monday, January 19, 2009 12:51 PM
 To: r-help@r-project.org
 Subject: [R] plotting arrows with different colors and varying head
 size
 
 Dear list,
 
 I would like to plot arrows with different colors according to arrow
 length, and also (if
 possible) with head size proportional to arrow length. The idea is to
 make a quiver-like plot of
 matlab with wind speed data.
 
 So far, I´ve been able to use different colors, but I need to find a
 more efficient way to recode
 arrow length intervals into colors. On the contrary, I can't define
 different head sizes, because
 the length argument in the arrows() function seems to control the
 head size of all the
 arrows at once.
 
 Any help will be greatly appreciated.
 
 
 
 ## Generate random data
  set.seed(1)
  xcoor - matrix(runif(20), ncol=5)
  ycoor - matrix(runif(20), ncol=5)
u - matrix(rnorm(20)/10, ncol=5)
v - matrix(rnorm(20)/10, ncol=5)
 
 ## calculate arrows length and look histogram
  arrowlen - sqrt(u^2 + v^2)
  hist(arrowlen)
 
 ## recode arrow lengths (by interval) into colors
 ## sorry, I don't know how to do it without the car  package
  library(car)
  arrowL - recode(as.vector(arrowlen),
  0='grey90';
  0:0.05='grey50';
0.05:0.1='grey';
0.1:0.15='cyan';
0.15:0.2='blue';
0.2:0.25='red')
 
 length=0.1
 
 par.uin - function()
   # determine scale of inches/userunits in x and y
   # from http://tolstoy.newcastle.edu.au/R/help/01c/2714.html
   # Brian Ripley Tue 20 Nov 2001 - 20:13:52 EST
  {
 u - par(usr)
 p - par(pin)
 c(p[1]/(u[2] - u[1]), p[2]/(u[4] - u[3]))
   }
 
 
 ## plot arrows
  plot(as.vector(xcoor), as.vector(ycoor), type=p, pch=., xlim=c(-
 0.2, 1.3),
   ylim=c(-0.2, 1.3))
arrows( xcoor, ycoor, xcoor + u, ycoor + v,
   length = length*as.vector(arrowlen)*min(par.uin()),
 col=arrowL)
 
 
 
   sessionInfo()
 R version 2.8.1 (2008-12-22)
 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 datasets  methods
 [7] base
 
 other attached packages:
 [1] car_1.2-9  hdf5_1.6.7
 
 
 
 --
 Héctor Villalobos hvill...@ipn.mx
  CICIMAR - IPN
   La Paz, Baja California Sur, MÉXICO
 
 
   [[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] XYplot in Lattice Package

2009-01-19 Thread jimdare


Thanks for your help David.   I managed to track down this solution in regard
to the second question:

http://www.nabble.com/lattice-xyplot-with-bty%3D%22l%22-tt12486052.html#a12489170

Regards,
James






David Winsemius wrote:
 
 
 On Jan 15, 2009, at 9:27 PM, jimdare wrote:
 

 Dear R-Users

 I have 2 questions to do with XYplot.

 1)

 I am trying to use the XYplot function to generate multiple line  
 graphs with
 the legend outside the plot.
 I am using the following loop for each graph:

 library(lattice)

 for (i in x.sp){
 xyplot(Catch~Year, df, groups = Stock, type=a,auto.key =
   list(space = top, points = FALSE, lines = TRUE,columns = 4))
 }
   From the help on package lattice:
 Note
 High level Lattice functions (like xyplot) are different from  
 conventional R graphics functions because they don't actually draw  
 anything. Instead, they return an object of class trellis which has  
 to be then printed or plotted to create the actual plot. This is  
 normally done automatically, but not when the high level functions are  
 called inside another function (most often source) or other contexts  
 where automatic printing is suppressed (e.g. for or while loops). In  
 such situations, an explicit call to print or plot is required.
 
 
 
 Or the FAQ:
 
 http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-do-lattice_002ftrellis-graphics-not-work_003f
 
 
 


 When I run the script I don't get any output graphs, however if I  
 change
 'XYplot' to 'plot', or generate the plots manually using XYplot, It  
 seems to
 work.  Is there a bug of some sort or is it me?
 
 It's you.
 


 2)

 How do I remove the top and right axis from a plot?  If I add  
 'axis(side
 = c(bottom, left)' to the xyplot call it comes up with the  
 message:

 Error in axis(side = c(bottom, left)) :
  plot.new has not been called yet
 
 see at the top of the axis help page:
 
 axis {graphics} so axis is not part of lattice
 
 Had you continued reading the next sentence in the lattice intro help  
 page, you would have seen:
 
 Lattice plots are highly customizable via user-modifiable settings.  
 However, these are completely unrelated to base graphics settings; in  
 particular, changing par() settings usually have no effect on lattice  
 plots.
 
 
 Try:
 http://stat.ethz.ch/R-manual/R-patched/library/lattice/html/axis.default.html
 
 And chapter 8 of:
 
 http://lmdvr.r-forge.r-project.org/figures/figures.html
 
 Best of luck and buy Sarkar's book.
 David Winsemius
 
 
 


 Any help is much appreciated :)


 -- 
 View this message in context:
 http://www.nabble.com/XYplot-in-Lattice-Package-tp21491296p21491296.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.
 
 

-- 
View this message in context: 
http://www.nabble.com/XYplot-in-Lattice-Package-tp21491296p21552650.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] Concave Hull

2009-01-19 Thread Greg Snow

I don't know if it is the same algorithm or not, but there is the function 
chull that finds the convex hull.

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of Michael Kubovy
 Sent: Saturday, January 17, 2009 9:49 AM
 To: r-help
 Subject: [R] Concave Hull
 
 Dear Friends,
 
 Here is an algorithm for finding concave hulls:
 http://get.dsi.uminho.pt/local/
 
 Has anyone implemented such an algorithm in R?
 
 RSiteSearch('concave hull') didn't reveal one (I think).
 
 _
 Professor Michael Kubovy
 University of Virginia
 Department of Psychology
 Postal Address:
   P.O.Box 400400, Charlottesville, VA 22904-4400
 Express Parcels Address:
   Gilmer Hall, Room 102, McCormick Road, Charlottesville, VA 22903
 Office:B011;  Phone: +1-434-982-4729
 Lab:B019; Phone: +1-434-982-4751
 WWW:http://www.people.virginia.edu/~mk9y/
 Skype name: polyurinsane
 
 
 
 
 
   [[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] Bar Plot ggplot2 Filling bars with cross hatching

2009-01-19 Thread Greg Snow

I think the fact that the grid package does not support cross-hatching is a 
feature not a bug (or deficiency), and I hope that this is not fixed.  
Tufte's book (The Visual Display of Quantitative Information) has a section on 
why cross-hatching should be avoided (unless of course your goal is to induce 
nausea in the observer rather than convey information).

I would edit Hadley's statement below to say fortunately there's no way to do 
this in ggplot2.

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of hadley wickham
 Sent: Thursday, January 15, 2009 10:55 AM
 To: stephen sefick
 Cc: R-help
 Subject: Re: [R] Bar Plot ggplot2 Filling bars with cross hatching
 
 Hi Stephen,
 
  #I am putting a test together for an introductory biology class and I
  would like to put different cross hatching inside of each bar for the
  bar plot below
 
 ggplot2 uses the grid package to do all the drawing, and currently
 grid doesn't support cross-hatching, so unfortunately there's no way
 to do this in ggplot2.
 
 Regards,
 
 Hadley
 
 --
 http://had.co.nz/
 
 __
 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] MH algorithm syntax help

2009-01-19 Thread ekwaters


Well spotted, the b is a muck up (what happens when you are basing stuff on
someone else's code). The bit you though might be a dimension mismatch seems
to work ok, but the bit that I was worried about from the start doesn't; the
reason I have sqrt of var y/ var x is because my posterior is for a reduced
major axis regression, that is the slope, but I somehow thought R might spit
that back at me. I might have to go back to the drawing board, and find a
different way of doing the posterior.

I get this, having fixed the algorithm.

 vx=sd(x)
 vy=sd(y)
 s2y=matrix(vy,m)
 s2x=matrix(vx,m)

 #begin MH sampling
 for(i in 2:m){
+ for(j in vy){s2y[i,j]=s2y[i-1,j]+rnorm(1,mean=0,sd=s2yscale[j]);acc=1}
+
if((post(y,x,s2ey[i-1],s2x,s2y[i])-post(y,x,s2ey[i-1],s2x,s2y[i-1]))log(runif(1,min=0,max=1))){s2y[i,j]=s2y[i-1,j];
acc=0}
+ accrate[i,j]=(accrate[i-1,j]*(i-1)+acc)/i}
Error in x * s2y/s2x : non-conformable arrays

Ned


David Winsemius wrote:
 
 The comma *before*  acc=1 ?
 
 I also wondered whether (further up)  this should work:
 s2y[i,]=s2y[i-1]  # would think this to result in a dimension mismatch
 
 This looks sketchy as well:
 s2y[i,j] = s2y[b[i-1,j] + rnorm(1,mean=0, sd=s2yscale[j])
^ ^ ^  # unmatched sqr-brackets
 
 It would be easier to run through a paren matching editor if you gave  
 the original loop code as well as the error output.
 -- 
 David Winsemius
 
 
 

-- 
View this message in context: 
http://www.nabble.com/MH-algorithm-syntax-help-tp21534889p21553281.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] bootstrapped eigenvector method following prcomp

2009-01-19 Thread Axel Strauß


@ Stas

Thanks for the extensive answer! I squeezed my data in your function but 
still need to mull over it and your comments for some time.


спасибо,
Axel




Stas Kolenikov schrieb:

I don't know if there are bugs in the code, but the step 4) does not
compute significance... at least the way statisticians know it. The
fractions above or below 0 are not significance. I don't even know how
to call those... probably cdf of the bootstrap distribution evaluated
at zero.

Let's put the flies and the sausages separately, as a Russian saying
goes. If you bootstrap from your original data, you can get the
confidence intervals, but not the test statistics. What you can do
with your bootstrap from the raw data is accumulate the distribution
of the eigenvectors, along the lines of (assuming you are testing for
the significance of variables in your first component):

function (x, permutations=1000, ...)
{
  pcnull - princomp(x, ... )
  res - pcnull$loadings[,1]
  bsresults = matrix( rep.int(NA, permutations*NROW(res)) ,
nrow=permutations, ncol=NROW(res) )
  N - nrow(x)
  for (i in 1:permutations) {
  pc - princomp(x[sample(N, replace=TRUE), ], ... )
  pred - predict(pc, newdata = x)
  r -  cor(pcnull$scores, pred)
  k - apply(abs(r), 2, which.max)
  reve - sign(diag(r[k,]))
  sol - pc$loadings[ ,k]
  sol - sweep(sol, 2, reve, *)
  bsresults[i,] - t(sol[,1])
  }
  apply( bsresults, 2, quantile, c(0.05, 0.95) )
}

if I am not messing up the dimensions and other stuff too much.
However as a result you will get an approximately degenerate
distribution sitting on the unit sphere since the eigenvectors are
always normalized to have the length of 1. You can still do marginal
confidence intervals with that though, and see if 0 is covered.

The main problem here is I am not entirely sure the bootstrap is
applicable for the problem at hand. In other words, it is not clear
whether the bootstrap estimates are consistent for the true
variability. Eigenproblems are quite difficult and prone to
non-regularities (the above mentioned degeneracy is just one of them,
and probably not the worst one). There are different asymptotics
(Anderson's of fixed p and n \to \infty,
http://www.citeulike.org/user/ctacmo/article/3908837, and Johnstone
joint p and n \to \infty with p/n \to const,
http://www.citeulike.org/user/ctacmo/article/3908846), the
distributions are often non-standard, while the bootstrap works well
when the distributions of the estimates are normal. When you have
something weird, bootstrap may easily break down, and in a lot of
other situations, you need to come up with special schemes. See my
oh-so-favorite paper on the bootstrap,
http://www.citeulike.org/user/ctacmo/article/575126.

One of those special schemes (back to out muttons, or rather flies and
sausages) -- to set up the bootstrap for hypothesis testing and get
the p-values, you need to bootstrap from the distribution that
corresponds to the null. Beran and Srivastava (1985 Annals,
http://www.citeulike.org/user/ctacmo/article/3015345) discuss how to
rotate your data to conform to the null hypothesis of interest for
inference with covariance matrices and their functions (such as
eigenvalues, for instance). Whether you need to go into all this
trouble, I don't really know.

If you have an inferential problem of testing whether a particular
variable contributes to the overall index, and have a pretty good idea
of where each variable goes, may be you need to shift your paradigm
and look at confirmatory factor analysis models instead, estimable in
R with John Fox' sem package.

On 1/19/09, Axel Strauß a.stra...@tu-bs.de wrote:
  

G'Day R users!

 Following an ordination using prcomp, I'd like to test which variables
singnificantly contribute to a principal component. There is a method
suggested by Peres-Neto and al. 2003. Ecology 84:2347-2363  called
bootstrapped eigenvector. It was asked for that in this forum in January
2005 by Jérôme Lemaître:
 1) Resample 1000 times with replacement entire raws from the original data
sets []
 2) Conduct a PCA on each bootstrapped sample
 3) To prevent axis reflexion and/or axis reordering in the bootstrap, here
are two more steps for each bootstrapped sample
 3a) calculate correlation matrix between the PCA scores of the original and
those of the bootstrapped sample
 3b) Examine whether the highest absolute correlation is between the
corresponding axis for the original and bootstrapped samples. When it is not
the case, reorder the eigenvectors. This means that if the highest
correlation is between the first original axis and the second bootstrapped
axis, the loadings for the second bootstrapped axis and use to estimate the
confidence interval for the original first PC axis.
 4) Determine the p value for each loading. Obtained as follow: number of
loadings =0 for loadings that were positive in the original matrix divided
by the number of boostrap samples (1000) and/or number of loadings =0 for

[R] Month tick marks on a plot()

2009-01-19 Thread glenn

Hi All,

I have a small dataframe [dates, values) I am plotting with
plot(df,type=²l²)

And the date date covers a year. The graph only have marks at 2008¹ and
2009¹.

How do I get the months labeled at the bottom please

Thanks as always

Glenn

[[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] Month tick marks on a plot()

2009-01-19 Thread jim holtman

?axis

On Mon, Jan 19, 2009 at 6:39 PM, glenn g1enn.robe...@btinternet.com wrote:
 Hi All,

 I have a small dataframe [dates, values) I am plotting with
 plot(df,type=²l²)

 And the date date covers a year. The graph only have marks at Œ2008¹ and
 Œ2009¹.

 How do I get the months labeled at the bottom please

 Thanks as always

 Glenn

[[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.

[R] easiest way to integrate own functions on startup

2009-01-19 Thread Jörg Groß


Hi,

I am currently writing some own functions that I frequently need.

So, it would be perfect if I could load these functions at the  
beginning of each R-session with a small command.



I tried to generate a R-package and install it that way.

But it seems that it is not so easy to add new functions to an  
existing R-package.

So I am not so flexible by that.


Is there a way to just load an .R-file which contains the function- 
definitions with a small command?

So that the functions are useable in the current session?


I tried load() but I get an error 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] easiest way to integrate own functions on startup

2009-01-19 Thread Duncan Murdoch


On 19/01/2009 7:13 PM, Jörg Groß wrote:

Hi,

I am currently writing some own functions that I frequently need.

So, it would be perfect if I could load these functions at the  
beginning of each R-session with a small command.



I tried to generate a R-package and install it that way.

But it seems that it is not so easy to add new functions to an  
existing R-package.

So I am not so flexible by that.


Is there a way to just load an .R-file which contains the function- 
definitions with a small command?

So that the functions are useable in the current session?


I tried load() but I get an error message...


load() loads a binary image that was saved by save().  You probably want 
source(), which reads and executes a file of R source.


Duncan Murdoch

__
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] Bar Plot ggplot2 Filling bars with cross hatching

what is your suggestion for distinguishing between many bars without
color?  I have grown up in the time of standarized tests - good or bad
I never felt nauseous.

Stephen

On Mon, Jan 19, 2009 at 5:20 PM, Greg Snow greg.s...@imail.org wrote:
 I think the fact that the grid package does not support cross-hatching is a 
 feature not a bug (or deficiency), and I hope that this is not fixed.  
 Tufte's book (The Visual Display of Quantitative Information) has a section 
 on why cross-hatching should be avoided (unless of course your goal is to 
 induce nausea in the observer rather than convey information).

 I would edit Hadley's statement below to say fortunately there's no way to 
 do this in ggplot2.

 --
 Gregory (Greg) L. Snow Ph.D.
 Statistical Data Center
 Intermountain Healthcare
 greg.s...@imail.org
 801.408.8111


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of hadley wickham
 Sent: Thursday, January 15, 2009 10:55 AM
 To: stephen sefick
 Cc: R-help
 Subject: Re: [R] Bar Plot ggplot2 Filling bars with cross hatching

 Hi Stephen,

  #I am putting a test together for an introductory biology class and I
  would like to put different cross hatching inside of each bar for the
  bar plot below

 ggplot2 uses the grid package to do all the drawing, and currently
 grid doesn't support cross-hatching, so unfortunately there's no way
 to do this in ggplot2.

 Regards,

 Hadley

 --
 http://had.co.nz/

 __
 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.




-- 
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.

Re: [R] Bar Plot ggplot2 Filling bars with cross hatching

2009-01-19 Thread Duncan Murdoch


On 19/01/2009 7:36 PM, stephen sefick wrote:

what is your suggestion for distinguishing between many bars without
color?  I have grown up in the time of standarized tests - good or bad
I never felt nauseous.


Use gray levels or labels.  If many is bigger than 5, it's not going 
to be easy, whatever method you are using.


Duncan Murdoch



Stephen

On Mon, Jan 19, 2009 at 5:20 PM, Greg Snow greg.s...@imail.org wrote:

I think the fact that the grid package does not support cross-hatching is a feature not a 
bug (or deficiency), and I hope that this is not fixed.  Tufte's book (The 
Visual Display of Quantitative Information) has a section on why cross-hatching should be 
avoided (unless of course your goal is to induce nausea in the observer rather than 
convey information).

I would edit Hadley's statement below to say fortunately there's no way to do this 
in ggplot2.

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111



-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
project.org] On Behalf Of hadley wickham
Sent: Thursday, January 15, 2009 10:55 AM
To: stephen sefick
Cc: R-help
Subject: Re: [R] Bar Plot ggplot2 Filling bars with cross hatching

Hi Stephen,


#I am putting a test together for an introductory biology class and I
would like to put different cross hatching inside of each bar for the
bar plot below

ggplot2 uses the grid package to do all the drawing, and currently
grid doesn't support cross-hatching, so unfortunately there's no way
to do this in ggplot2.

Regards,

Hadley

--
http://had.co.nz/

__
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] Bar Plot ggplot2 Filling bars with cross hatching



On 20/01/2009, at 1:36 PM, stephen sefick wrote:


what is your suggestion for distinguishing between many bars without
color?


Exactly.  Sometimes colour printing can be expensive.


I have grown up in the time of standarized tests - good or bad
I never felt nauseous.


Ni moi non plus.

cheers,

Rolf Turner


Stephen

On Mon, Jan 19, 2009 at 5:20 PM, Greg Snow greg.s...@imail.org  
wrote:
I think the fact that the grid package does not support cross- 
hatching is a feature not a bug (or deficiency), and I hope that  
this is not fixed.  Tufte's book (The Visual Display of  
Quantitative Information) has a section on why cross-hatching  
should be avoided (unless of course your goal is to induce nausea  
in the observer rather than convey information).


I would edit Hadley's statement below to say fortunately there's  
no way to do this in ggplot2.


--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111



-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
project.org] On Behalf Of hadley wickham
Sent: Thursday, January 15, 2009 10:55 AM
To: stephen sefick
Cc: R-help
Subject: Re: [R] Bar Plot ggplot2 Filling bars with cross hatching

Hi Stephen,

#I am putting a test together for an introductory biology class  
and I
would like to put different cross hatching inside of each bar  
for the

bar plot below


ggplot2 uses the grid package to do all the drawing, and currently
grid doesn't support cross-hatching, so unfortunately there's no way
to do this in ggplot2.

Regards,

Hadley

--
http://had.co.nz/

__
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.






--
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.



##
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] Bar Plot ggplot2 Filling bars with cross hatching

If classic graphics is ok try this which uses hatches and different
shades of grey:

barplot(lizards, names.arg = color, col = grey(c(.2, .5, 1)), density = 20,
angle = c(45, -45, 0), legend = color)




On Wed, Jan 14, 2009 at 10:18 PM, stephen sefick ssef...@gmail.com wrote:
 #I am putting a test together for an introductory biology class and I
 would like to put different cross hatching inside of each bar for the
 bar plot below

 color - c(Brightly Colored, Dull, Neither)
 lizards - c(277, 70, 3)
 liz.col - data.frame(color, lizards)
 qplot(color, lizards, data=liz.col, geom=bar, ylab=Observed
 Matings, main=Counts Out of 350 Aquariums, ylim=c(0,400),
 fill=color)+scale_y_continuous(breaks=c(0, 70, 277, 350))


 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-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] Bar Plot ggplot2 Filling bars with cross hatching



On 20/01/2009, at 1:46 PM, Duncan Murdoch wrote:


On 19/01/2009 7:36 PM, stephen sefick wrote:

what is your suggestion for distinguishing between many bars without
color?  I have grown up in the time of standarized tests - good or  
bad

I never felt nauseous.


Use gray levels or labels.  If many is bigger than 5, it's not going
to be easy, whatever method you are using.


I disagree.  Grey levels suck; labels are a kludge.  It is an issue
for ``many'' == 2, for which crosshatching works perfectly.

cheers,

Rolf


Duncan Murdoch



Stephen

On Mon, Jan 19, 2009 at 5:20 PM, Greg Snow greg.s...@imail.org  
wrote:
I think the fact that the grid package does not support cross- 
hatching is a feature not a bug (or deficiency), and I hope that  
this is not fixed.  Tufte's book (The Visual Display of  
Quantitative Information) has a section on why cross-hatching  
should be avoided (unless of course your goal is to induce nausea  
in the observer rather than convey information).


I would edit Hadley's statement below to say fortunately there's  
no way to do this in ggplot2.


--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111



-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
project.org] On Behalf Of hadley wickham
Sent: Thursday, January 15, 2009 10:55 AM
To: stephen sefick
Cc: R-help
Subject: Re: [R] Bar Plot ggplot2 Filling bars with cross hatching

Hi Stephen,

#I am putting a test together for an introductory biology class  
and I
would like to put different cross hatching inside of each bar  
for the

bar plot below

ggplot2 uses the grid package to do all the drawing, and currently
grid doesn't support cross-hatching, so unfortunately there's no  
way

to do this in ggplot2.

Regards,

Hadley

--
http://had.co.nz/

__
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.



##
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] easiest way to integrate own functions on startup

2009-01-19 Thread Jörg Groß


Is there a way to execute this command on every R-start?


I tried to add something in the Startup.h file - but that didn't work.

(working on a mac)




Thanks!

Am 20.01.2009 um 01:25 schrieb Duncan Murdoch:


On 19/01/2009 7:13 PM, Jörg Groß wrote:

Hi,
I am currently writing some own functions that I frequently need.
So, it would be perfect if I could load these functions at the   
beginning of each R-session with a small command.

I tried to generate a R-package and install it that way.
But it seems that it is not so easy to add new functions to an   
existing R-package.

So I am not so flexible by that.
Is there a way to just load an .R-file which contains the function-  
definitions with a small command?

So that the functions are useable in the current session?
I tried load() but I get an error message...


load() loads a binary image that was saved by save().  You probably  
want source(), which reads and executes a file of R source.


Duncan Murdoch


__
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] easiest way to integrate own functions on startup

?.Rprofile

On Mon, Jan 19, 2009 at 8:03 PM, Jörg Groß jo...@licht-malerei.de wrote:
 Is there a way to execute this command on every R-start?


 I tried to add something in the Startup.h file - but that didn't work.

 (working on a mac)




 Thanks!

 Am 20.01.2009 um 01:25 schrieb Duncan Murdoch:

 On 19/01/2009 7:13 PM, Jörg Groß wrote:

 Hi,
 I am currently writing some own functions that I frequently need.
 So, it would be perfect if I could load these functions at the  beginning
 of each R-session with a small command.
 I tried to generate a R-package and install it that way.
 But it seems that it is not so easy to add new functions to an  existing
 R-package.
 So I am not so flexible by that.
 Is there a way to just load an .R-file which contains the function-
 definitions with a small command?
 So that the functions are useable in the current session?
 I tried load() but I get an error message...

 load() loads a binary image that was saved by save().  You probably want
 source(), which reads and executes a file of R source.

 Duncan Murdoch

 __
 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.




-- 
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.

Re: [R] Month tick marks on a plot()

The zoo faq has an example:

library(zoo)
vignette(zoo-faq) # see question 8

On Mon, Jan 19, 2009 at 6:39 PM, glenn g1enn.robe...@btinternet.com wrote:
 Hi All,

 I have a small dataframe [dates, values) I am plotting with
 plot(df,type=²l²)

 And the date date covers a year. The graph only have marks at Œ2008¹ and
 Œ2009¹.

 How do I get the months labeled at the bottom please

 Thanks as always

 Glenn

[[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] easiest way to integrate own functions on startup

2009-01-19 Thread Jörg Groß


Hi,

I tried to generate a .Rprofile file.
But R does not load it automatically.

Is there a tutorial on the web on generating such a file?
(haven't found anything that helped me)

And where do I have to put this .Rprofile-file?
In the working directory?

Does R generate a .Rprofile file when R is installed?
And if yes; where can I find this file on a mac?


thanks for any help!

__
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] easiest way to integrate own functions on startup

2009-01-19 Thread Bernardo Rangel Tura

On Tue, 2009-01-20 at 01:13 +0100, Jörg Groß wrote:
 Hi,
 
 I am currently writing some own functions that I frequently need.
 
 So, it would be perfect if I could load these functions at the  
 beginning of each R-session with a small command.
 
 
 I tried to generate a R-package and install it that way.
 
 But it seems that it is not so easy to add new functions to an  
 existing R-package.
 So I am not so flexible by that.
 
 
 Is there a way to just load an .R-file which contains the function- 
 definitions with a small command?
 So that the functions are useable in the current session?
 
 
 I tried load() but I get an error message...


Hi Jörg,

I do you not use .Fisrt.

Something like this:

.First - function(){
source(MyFunction1.txt)
source(MyFunction2.txt)
source(MyFunction3.txt)
...
source(MyFunctionn.txt)
}

-- 
Bernardo Rangel Tura, M.D,MPH,Ph.D
National Institute of Cardiology
Brazil

__
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] easiest way to integrate own functions on startup