date:20121001


On Sep 30, 2012, at 6:48 PM, Dan Bebber wrote:

 I have spatial data on a sphere (the Earth) for which I would like to run an 
 gls model assuming that the errors are autcorrelated, i.e. including a 
 corSpatial correlation in the model specification.
 
 In this case the distance metric should be calculated on the sphere, 
 therefore metric =  euclidean in (for example) corSpher would be incorrect.
 
 I would be grateful for help on how to write a new distance metric for the 
 corSpatial function.
 I believe there are several ways that distances on a sphere can be calculated 
 in R, for example the distMeeus function in the geosphere library. However, 
 I have no idea how to write this into a corSpatial function.
 
 The aim is to end up with a metric  = sphere option that calculates great 
 circle distances between points using latitude and longitude.

LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html

-- 

David Winsemius, MD
Alameda, CA, USA

__
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] merge.zoo returns unmatched dates

2012-10-01 Thread Achim Zeileis


On Mon, 1 Oct 2012, Vindoggy ! wrote:



Sorry for the lack of reproducible data, but this seems to be a problem 
inherent to my dataset and I can't figure out where the issue is.

I have several data frames set up as a time series with identical POSIXct date 
formats. If I keep the original data in data frame format and merge them using 
base merge- everything is perfect and everyone is happy.

If I transform the data frames to zoo objects, and then do a merge.zoo- the 
data seem to become uncoupled from the original data. Even more unusual is that 
some dates in the new merged data set  are prior to the original data set. I've 
attempted bellow to show what this looks like, and I hope someone has a 
suggestion as to what may be causing the problem.

Here is one set of data in data.frame format

head(Vup)
Date Velocity_m/s
1 2010-01-21 07:42:00 1.217943
2 2010-01-21 07:43:00 1.624395
3 2010-01-21 07:44:00 1.526379
4 2010-01-21 07:45:00 1.456831
5 2010-01-21 07:46:00 1.245390
6 2010-01-21 07:47:00 1.374330

str(Vup)
'data.frame':7168 obs. of  2 variables:
$ Date: POSIXct, format: 2010-01-21 07:42:00 2010-01-21 07:43:00 ...
$ Velocity_m/s: num  1.22 1.62 1.53 1.46 1.25 ...

And here is a second in data.frame format:

head(PAS)
Date   PAS
1 2010-01-21 05:01:00   0.0013938
2 2010-01-21 05:02:00   0.0015331
3 2010-01-21 05:03:00   0.0016725
4 2010-01-21 05:04:00   0.0016725
5 2010-01-21 05:05:00   0.0012265
6 2010-01-21 05:06:00   0.0015889

str(PAS)
'data.frame':5520 obs. of  2 variables:
$ Date   : POSIXct, format: 2010-01-21 05:01:00 2010-01-21 05:02:00 ...
$ PAS: num  0.00139 0.00153 0.00167 0.00167 0.00123 ...



Using zoo:

PASmin-zoo(as.matrix(PAS[,2]),as.POSIXct(PAS[,1],format=%Y-%m-%d 
%H:%M:%S,tz=UTC))

str(PASmin)
?zoo? series from 2010-01-21 05:01:00 to 2010-01-27 13:01:00
 Data: num [1:5520, 1] 0.00139 0.00153 0.00167 0.00167 0.00123 ...
- attr(*, dimnames)=List of 2
 ..$ : NULL
 ..$ : chr PAS
 Index:  POSIXct[1:5520], format: 2010-01-21 05:01:00 2010-01-21 05:02:00 
2010-01-21 05:03:00 ...




ADP_UPmin-zoo(as.matrix(Vup[,2]),as.POSIXct(Vup[,1], format=%Y-%m-%d 
%H:%M,tz=UTC))

str(ADP_UPmin)
?zoo? series from 2010-01-21 07:42:00 to 2010-01-26 20:12:00
 Data: num [1:7168, 1] 1.22 1.62 1.53 1.46 1.25 ...
- attr(*, dimnames)=List of 2
 ..$ : NULL
 ..$ : chr UP_Velocity_m/s
 Index:  POSIXct[1:7168], format: 2010-01-21 07:42:00 2010-01-21 07:43:00 
2010-01-21 07:44:00 ...


And if I merge the two zoo objects I get this:

M-merge(ADP_UPmin,PASmin)
head(M)

   UP_Velocity_m/s   PAS
2010-01-20 21:01:00  NA 0.0013938
2010-01-20 21:02:00  NA 0.0015331
2010-01-20 21:03:00  NA 0.0016725
2010-01-20 21:04:00  NA 0.0016725
2010-01-20 21:05:00  NA 0.0012265
2010-01-20 21:06:00  NA 0.0015889


?zoo? series from 2010-01-20 21:01:00 to 2010-01-27 05:01:00
 Data: num [1:8499, 1:2] NA NA NA NA NA NA NA NA NA NA ...
- attr(*, dimnames)=List of 2
 ..$ : NULL
 ..$ : chr [1:2] UP_Velocity_m/s PAR
 Index:  POSIXct[1:8499], format: 2010-01-20 21:01:00 2010-01-20 21:02:00 
2010-01-20 21:03:00 ...


For some reason I can not figure out, even though both the PAS data frame and 
PAS zoo object starts at 2010-01-21 05:01:00, once merged the PAS data starts a 
day earlier at 2010-01-20 21:01:00.  The actual numeric data looks good, but 
both variables have no come uncoupled from the time series dates (The Velocity 
data is similarity uncoupled). And as stated before, doing an non-zoo merge on 
the data.frame data works fine.

Anyone got any ideas what's going on?


My guess is that you create both zoo series with time zone UTC but that 
the TZ attribute gets lost upon the merge. Then, the time is displayed in 
your systems time zone (which you haven't told us) which apparently is a 
couple of hours before UTC.


On my system (which is in CET) I can create a series with UTC times

R x - zoo(1:2, as.POSIXct(c(2012-01-01 00:00:00,
+2012-01-01 01:00:00), format = %Y-%m-%d %H:%M:%S, tz = UTC))
R x
2012-01-01 00:00:00 2012-01-01 01:00:00
  1   2

The times are in UTC as requested, but applying the c() method, they get 
dropped. See ?c.POSIXct.


R time(x)
[1] 2012-01-01 00:00:00 UTC 2012-01-01 01:00:00 UTC
R c(time(x))
[1] 2012-01-01 01:00:00 CET 2012-01-01 02:00:00 CET

Hence:

R merge(x, x)
x x
2012-01-01 01:00:00 1 1
2012-01-01 02:00:00 2 2

But you can set the system time in your R session to UTC which gives the 
desired result:


R Sys.setenv(TZ = UTC)
R merge(x, x)
x x
2012-01-01 00:00:00 1 1
2012-01-01 01:00:00 2 2

hth,
Z




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

Re: [R] nlme: spatial autocorrelation on a sphere

2012-10-01 Thread Dan Bebber

Thanks, but the problem is quite specific and not addressed on the Spatial Data 
taskview page.
Quite specifically, I would like to know how to edit corSpatial functions to 
calculate great circle distances.
The Bayesian equivalent, georamps in the ramps package, is able to do this, 
therefore I imagine it must be possible for nlme.

Dan

From: David Winsemius [dwinsem...@comcast.net]
Sent: 01 October 2012 08:38
To: Dan Bebber
Cc: r-help@r-project.org
Subject: Re: [R] nlme: spatial autocorrelation on a sphere

On Sep 30, 2012, at 6:48 PM, Dan Bebber wrote:

 I have spatial data on a sphere (the Earth) for which I would like to run an 
 gls model assuming that the errors are autcorrelated, i.e. including a 
 corSpatial correlation in the model specification.

 In this case the distance metric should be calculated on the sphere, 
 therefore metric =  euclidean in (for example) corSpher would be incorrect.

 I would be grateful for help on how to write a new distance metric for the 
 corSpatial function.
 I believe there are several ways that distances on a sphere can be calculated 
 in R, for example the distMeeus function in the geosphere library. However, 
 I have no idea how to write this into a corSpatial function.

 The aim is to end up with a metric  = sphere option that calculates great 
 circle distances between points using latitude and longitude.

LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html

--

David Winsemius, MD
Alameda, CA, USA



__
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] adjust font in ggplot2 to LaTeX document

2012-10-01 Thread Jonas Stein

Hi,

how can i adjust the font in a ggplot2 qplot so that it will look 
similar to the LaTeX font? 
Computer Modern Sans Serif in the same size would be nice.

My output device is 
ggsave(filename=test.pdf, width=5.5, height=3, dpi=300)
and i will include the graphic with 5.5 inch in LaTeX.

I found some pages about that topic but all solutions that i found have
been very complicate and confusing to me. 
But the articles have been from 2009 too, so i guess there will be
some easy solution today...

Is theme the keyword to find the solution?
Or will i have to include a binary font file?

Can someone give me a link to an example?

Kind regards and thanks a lot,

-- 
Jonas Stein n...@jonasstein.de

__
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] merge.zoo returns unmatched dates

2012-10-01 Thread Vindoggy !

You nailed it, It was a tz issue. Thanks!

 Date: Mon, 1 Oct 2012 09:41:35 +0200
 From: achim.zeil...@uibk.ac.at
 To: vindo...@hotmail.com
 CC: r-help@r-project.org
 Subject: Re: [R] merge.zoo returns unmatched dates

 On Mon, 1 Oct 2012, Vindoggy ! wrote:

  Sorry for the lack of reproducible data, but this seems to be a problem 
  inherent to my dataset and I can't figure out where the issue is.

  I have several data frames set up as a time series with identical POSIXct 
  date formats. If I keep the original data in data frame format and merge 
  them using base merge- everything is perfect and everyone is happy.

  If I transform the data frames to zoo objects, and then do a merge.zoo- the 
  data seem to become uncoupled from the original data. Even more unusual is 
  that some dates in the new merged data set  are prior to the original data 
  set. I've attempted bellow to show what this looks like, and I hope someone 
  has a suggestion as to what may be causing the problem.

  Here is one set of data in data.frame format

  head(Vup)
  Date Velocity_m/s
  1 2010-01-21 07:42:00 1.217943
  2 2010-01-21 07:43:00 1.624395
  3 2010-01-21 07:44:00 1.526379
  4 2010-01-21 07:45:00 1.456831
  5 2010-01-21 07:46:00 1.245390
  6 2010-01-21 07:47:00 1.374330

  str(Vup)
  'data.frame':7168 obs. of  2 variables:
  $ Date: POSIXct, format: 2010-01-21 07:42:00 2010-01-21 
  07:43:00 ...
  $ Velocity_m/s: num  1.22 1.62 1.53 1.46 1.25 ...

  And here is a second in data.frame format:

  head(PAS)
  Date   PAS
  1 2010-01-21 05:01:00   0.0013938
  2 2010-01-21 05:02:00   0.0015331
  3 2010-01-21 05:03:00   0.0016725
  4 2010-01-21 05:04:00   0.0016725
  5 2010-01-21 05:05:00   0.0012265
  6 2010-01-21 05:06:00   0.0015889

  str(PAS)
  'data.frame':5520 obs. of  2 variables:
  $ Date   : POSIXct, format: 2010-01-21 05:01:00 2010-01-21 05:02:00 
  ...
  $ PAS: num  0.00139 0.00153 0.00167 0.00167 0.00123 ...

  Using zoo:

  PASmin-zoo(as.matrix(PAS[,2]),as.POSIXct(PAS[,1],format=%Y-%m-%d 
  %H:%M:%S,tz=UTC))

  str(PASmin)
  ?zoo? series from 2010-01-21 05:01:00 to 2010-01-27 13:01:00
   Data: num [1:5520, 1] 0.00139 0.00153 0.00167 0.00167 0.00123 ...
  - attr(*, dimnames)=List of 2
   ..$ : NULL
   ..$ : chr PAS
   Index:  POSIXct[1:5520], format: 2010-01-21 05:01:00 2010-01-21 
  05:02:00 2010-01-21 05:03:00 ...

  ADP_UPmin-zoo(as.matrix(Vup[,2]),as.POSIXct(Vup[,1], format=%Y-%m-%d 
  %H:%M,tz=UTC))

  str(ADP_UPmin)
  ?zoo? series from 2010-01-21 07:42:00 to 2010-01-26 20:12:00
   Data: num [1:7168, 1] 1.22 1.62 1.53 1.46 1.25 ...
  - attr(*, dimnames)=List of 2
   ..$ : NULL
   ..$ : chr UP_Velocity_m/s
   Index:  POSIXct[1:7168], format: 2010-01-21 07:42:00 2010-01-21 
  07:43:00 2010-01-21 07:44:00 ...

  And if I merge the two zoo objects I get this:

  M-merge(ADP_UPmin,PASmin)
  head(M)

 UP_Velocity_m/s   PAS
  2010-01-20 21:01:00  NA 0.0013938
  2010-01-20 21:02:00  NA 0.0015331
  2010-01-20 21:03:00  NA 0.0016725
  2010-01-20 21:04:00  NA 0.0016725
  2010-01-20 21:05:00  NA 0.0012265
  2010-01-20 21:06:00  NA 0.0015889

  ?zoo? series from 2010-01-20 21:01:00 to 2010-01-27 05:01:00
   Data: num [1:8499, 1:2] NA NA NA NA NA NA NA NA NA NA ...
  - attr(*, dimnames)=List of 2
   ..$ : NULL
   ..$ : chr [1:2] UP_Velocity_m/s PAR
   Index:  POSIXct[1:8499], format: 2010-01-20 21:01:00 2010-01-20 
  21:02:00 2010-01-20 21:03:00 ...

  For some reason I can not figure out, even though both the PAS data frame 
  and PAS zoo object starts at 2010-01-21 05:01:00, once merged the PAS data 
  starts a day earlier at 2010-01-20 21:01:00.  The actual numeric data looks 
  good, but both variables have no come uncoupled from the time series dates 
  (The Velocity data is similarity uncoupled). And as stated before, doing an 
  non-zoo merge on the data.frame data works fine.

  Anyone got any ideas what's going on?

 My guess is that you create both zoo series with time zone UTC but that 
 the TZ attribute gets lost upon the merge. Then, the time is displayed in 
 your systems time zone (which you haven't told us) which apparently is a 
 couple of hours before UTC.

 On my system (which is in CET) I can create a series with UTC times

 R x - zoo(1:2, as.POSIXct(c(2012-01-01 00:00:00,
 +2012-01-01 01:00:00), format = %Y-%m-%d %H:%M:%S, tz = UTC))
 R x
 2012-01-01 00:00:00 2012-01-01 01:00:00
1   2

 The times are in UTC as requested, but applying the c() method, they get 
 dropped. See ?c.POSIXct.

 R time(x)
 [1] 2012-01-01 00:00:00 UTC 2012-01-01 01:00:00 UTC
 R c(time(x))
 [1] 2012-01-01 01:00:00 CET 2012-01-01 02:00:00 CET

 Hence:

 R merge(x, x)
  x x
 2012-01-01 01:00:00 1 1
 2012-01-01 02:00:00 2 2

 But you can set

Re: [R] REML - quasipoisson

2012-10-01 Thread Greg Dropkin

hi

the explanation seems to be that fix.family.ls is invoked to define a
saturated log likelihood for quasi families, which means that the F2q term
which I had thought to be undefined, is actually defined in this example
as

F2q--500 * log(phiq)/2

The formula then becomes

F1q-F2q+F3q-F4q

which does coincide with m3$gcv

Note that when phiq=1, the F2q term vanishes.

greg




 hi

 I'm puzzled as to the relation between the REML score computed by gam and
 the formula (4) on p.4 here:
 http://opus.bath.ac.uk/22707/1/Wood_JRSSB_2011_73_1_3.pdf

 I'm ok with this for poisson, or for quasipoisson when phi=1.

 However, when phi differs from 1, I'm stuck.

 #simulate some data
 library(mgcv)
 set.seed(1)
 x1-runif(500)
 x2-rnorm(500)
 linp--0.5+x1+exp(-x2^2/2)*sin(4*x2)
 y-rpois(500,exp(linp))

 ##poisson
 #phi=1
 m1-gam(y~s(x1)+s(x2),family=poisson,method=REML)
 phi-m1$scale

 #formula
 #1st term
 S1-m1$smooth[[1]]$S[[1]]*m1$sp[1]
 S2-m1$smooth[[2]]$S[[1]]*m1$sp[2]
 S-matrix(0,19,19)
 for (i in 2:10)
 {
 for (j in 2:10)
 {
 S[i,j]=S1[i-1,j-1]
 S[i+9,j+9]=S2[i-1,j-1]
 }
 }
 beta-m1$coef
 #penalised deviance
 Dp-m1$dev+t(beta)%*%S%*%beta
 F1-Dp/(2*phi)

 #2nd term
 F2-sum(ifelse(y==0,0,y*log(y)-y-log(factorial(y

 #3rd term
 X-predict(m1,type=lpmatrix)
 W-diag(fitted(m1))
 H-t(X)%*%W%*%X
 ldhs-determinant(H+S,log=TRUE)$modulus[1]
 eigS-eigen(S,only.values=TRUE)$val
 lds-sum(log(eigS[1:16]))
 F3-(ldhs-lds)/2

 #4th term
 Mp=3
 F4-Mp/2*log(2*pi*phi)

 F1-F2+F3-F4
 m1$gcv
 #reml score = formula

 ##quasipoisson with scale = 1
 #fitting is identical to the poisson case
 #F1, F3, and F4 unchanged but F2 is now undefined

 m2-gam(y~s(x1)+s(x2),family=quasipoisson,method=REML,scale=1)
 F1+F3-F4
 m2$gcv
 #reml score = formula with F2 omitted

 ##quasipoisson with unknown scale
 m3-gam(y~s(x1)+s(x2),family=quasipoisson,method=REML,scale=-1)
 phiq-m3$scale

 #1st term
 S1q-m3$smooth[[1]]$S[[1]]*m3$sp[1]
 S2q-m3$smooth[[2]]$S[[1]]*m3$sp[2]
 Sq-matrix(0,19,19)
 for (i in 2:10)
 {
 for (j in 2:10)
 {
 Sq[i,j]=S1q[i-1,j-1]
 Sq[i+9,j+9]=S2q[i-1,j-1]
 }
 }
 betaq-m3$coef
 #penalised deviance
 Dpq-m3$dev+t(betaq)%*%Sq%*%betaq
 F1q-Dpq/(2*phiq)

 #2nd term undefined

 #3rd term
 Xq-predict(m3,type=lpmatrix)
 Wq-diag(fitted(m3))
 Hq-t(Xq)%*%Wq%*%Xq
 ldhsq-determinant(Hq+Sq,log=TRUE)$modulus[1]
 eigSq-eigen(Sq,only.values=TRUE)$val
 ldsq-sum(log(eigSq[1:16]))
 F3q-(ldhsq-ldsq)/2

 #4th term
 Mp=3
 F4q-Mp/2*log(2*pi*phiq)

 F1q+F3q-F4q
 m3$gcv

 #quite different

 #but if phiq is replaced by the Pearson estimate of the scale
 P-sum((y-fitted(m3))^2/fitted(m3))
 phip-P/(500-sum(m3$edf))
 F1p-Dpq/(2*phip)
 F3p-F3q
 #third term independent of scale
 F4p-Mp/2*log(2*pi*phip)
 F1p+F3p-F4p
 m3$gcv

 #closer but still wrong

 #is there a value of phi which makes this work?
 f-function(t) (Dpq/(2*t) + F3q + Mp/2*log(2*pi*t) - m3$gcv)^2
 optimise(f,interval=c(0.5,1.5))$min-phix

 Dpq/(2*phix) + F3q + Mp/2*log(2*pi*phix)
 m3$gcv
 #but what is phix - not the Pearson estimate or the scale returned by gam?

 thanks

 Greg Dropkin



__
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] [Fwd: REML - quasipoisson]

2012-10-01 Thread Simon Wood


Hi Greg,

For quasi families I've used extended quasi-likelihood (see Mccullagh 
and Nelder, Generalized Linear Models 2nd ed, section 9.6) in place of 
the likelihood/quasi-likelihood in the expression for the (RE)ML score. 
I hadn't realised that this was possible before the paper was published.


best,
Simon

ps. sorry for slow reply, the original message slipped through my filter 
for mgcv related stuff



hi

I'm puzzled as to the relation between the REML score computed by gam and
the formula (4) on p.4 here:
http://opus.bath.ac.uk/22707/1/Wood_JRSSB_2011_73_1_3.pdf

I'm ok with this for poisson, or for quasipoisson when phi=1.

However, when phi differs from 1, I'm stuck.

#simulate some data
library(mgcv)
set.seed(1)
x1-runif(500)
x2-rnorm(500)
linp--0.5+x1+exp(-x2^2/2)*sin(4*x2)
y-rpois(500,exp(linp))

##poisson
#phi=1
m1-gam(y~s(x1)+s(x2),family=poisson,method=REML)
phi-m1$scale

#formula
#1st term
S1-m1$smooth[[1]]$S[[1]]*m1$sp[1]
S2-m1$smooth[[2]]$S[[1]]*m1$sp[2]
S-matrix(0,19,19)
for (i in 2:10)
{
for (j in 2:10)
{
S[i,j]=S1[i-1,j-1]
S[i+9,j+9]=S2[i-1,j-1]
}
}
beta-m1$coef
#penalised deviance
Dp-m1$dev+t(beta)%*%S%*%beta
F1-Dp/(2*phi)

#2nd term
F2-sum(ifelse(y==0,0,y*log(y)-y-log(factorial(y

#3rd term
X-predict(m1,type=lpmatrix)
W-diag(fitted(m1))
H-t(X)%*%W%*%X
ldhs-determinant(H+S,log=TRUE)$modulus[1]
eigS-eigen(S,only.values=TRUE)$val
lds-sum(log(eigS[1:16]))
F3-(ldhs-lds)/2

#4th term
Mp=3
F4-Mp/2*log(2*pi*phi)

F1-F2+F3-F4
m1$gcv
#reml score = formula

##quasipoisson with scale = 1
#fitting is identical to the poisson case
#F1, F3, and F4 unchanged but F2 is now undefined

m2-gam(y~s(x1)+s(x2),family=quasipoisson,method=REML,scale=1)
F1+F3-F4
m2$gcv
#reml score = formula with F2 omitted

##quasipoisson with unknown scale
m3-gam(y~s(x1)+s(x2),family=quasipoisson,method=REML,scale=-1)
phiq-m3$scale

#1st term
S1q-m3$smooth[[1]]$S[[1]]*m3$sp[1]
S2q-m3$smooth[[2]]$S[[1]]*m3$sp[2]
Sq-matrix(0,19,19)
for (i in 2:10)
{
for (j in 2:10)
{
Sq[i,j]=S1q[i-1,j-1]
Sq[i+9,j+9]=S2q[i-1,j-1]
}
}
betaq-m3$coef
#penalised deviance
Dpq-m3$dev+t(betaq)%*%Sq%*%betaq
F1q-Dpq/(2*phiq)

#2nd term undefined

#3rd term
Xq-predict(m3,type=lpmatrix)
Wq-diag(fitted(m3))
Hq-t(Xq)%*%Wq%*%Xq
ldhsq-determinant(Hq+Sq,log=TRUE)$modulus[1]
eigSq-eigen(Sq,only.values=TRUE)$val
ldsq-sum(log(eigSq[1:16]))
F3q-(ldhsq-ldsq)/2

#4th term
Mp=3
F4q-Mp/2*log(2*pi*phiq)

F1q+F3q-F4q
m3$gcv

#quite different

#but if phiq is replaced by the Pearson estimate of the scale
P-sum((y-fitted(m3))^2/fitted(m3))
phip-P/(500-sum(m3$edf))
F1p-Dpq/(2*phip)
F3p-F3q
#third term independent of scale
F4p-Mp/2*log(2*pi*phip)
F1p+F3p-F4p
m3$gcv

#closer but still wrong

#is there a value of phi which makes this work?
f-function(t) (Dpq/(2*t) + F3q + Mp/2*log(2*pi*t) - m3$gcv)^2
optimise(f,interval=c(0.5,1.5))$min-phix

Dpq/(2*phix) + F3q + Mp/2*log(2*pi*phix)
m3$gcv
#but what is phix - not the Pearson estimate or the scale returned by gam?

thanks

Greg Dropkin






--
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603   http://people.bath.ac.uk/sw283

__
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] nlme: spatial autocorrelation on a sphere


On Oct 1, 2012, at 12:59 AM, Dan Bebber wrote:

 Thanks, but the problem is quite specific and not addressed on the Spatial 
 Data taskview page.
 Quite specifically, I would like to know how to edit corSpatial functions to 
 calculate great circle distances.
 The Bayesian equivalent, georamps in the ramps package, is able to do this, 
 therefore I imagine it must be possible for nlme.
 

Can't you use the corStruct functions in pkg::ramps? They allow specification 
of the 'haversine' metric. The corR* functions inherit from class corStruct.

-- 
David.

 Dan
 
 From: David Winsemius [dwinsem...@comcast.net]
 Sent: 01 October 2012 08:38
 To: Dan Bebber
 Cc: r-help@r-project.org
 Subject: Re: [R] nlme: spatial autocorrelation on a sphere
 
 On Sep 30, 2012, at 6:48 PM, Dan Bebber wrote:
 
 I have spatial data on a sphere (the Earth) for which I would like to run an 
 gls model assuming that the errors are autcorrelated, i.e. including a 
 corSpatial correlation in the model specification.
 
 In this case the distance metric should be calculated on the sphere, 
 therefore metric =  euclidean in (for example) corSpher would be incorrect.
 
 I would be grateful for help on how to write a new distance metric for the 
 corSpatial function.
 I believe there are several ways that distances on a sphere can be 
 calculated in R, for example the distMeeus function in the geosphere 
 library. However, I have no idea how to write this into a corSpatial 
 function.
 
 The aim is to end up with a metric  = sphere option that calculates great 
 circle distances between points using latitude and longitude.
 
 LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html
 
 --
 
 David Winsemius, MD
 Alameda, CA, USA
 
 

David Winsemius, MD
Alameda, CA, USA

__
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] adjust font in ggplot2 to LaTeX document

2012-10-01 Thread Ben Bolker

Jonas Stein news at jonasstein.de writes:

 
 Hi,
 
 how can i adjust the font in a ggplot2 qplot so that it will look 
 similar to the LaTeX font? 
 Computer Modern Sans Serif in the same size would be nice.
 
 My output device is 
 ggsave(filename=test.pdf, width=5.5, height=3, dpi=300)
 and i will include the graphic with 5.5 inch in LaTeX.
 

  The easiest thing to do is to use the tikzDevice device
for graphics output, and in turn the easiest thing to do is
to use it within knitr ... google 'tikzDevice' ...

__
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] installing a package from Archive

2012-10-01 Thread grazia


Dear all,
I'm trying to install in R the package blme that has been removed from R 
cran but it is still in the Archive:

http://cran.r-project.org/src/contrib/Archive/blme/

I have a Mac OSX and I have been trying with the following command from 
the Shell:


R CMD INSTALL --build 
http://cran.r-project.org/src/contrib/Archive/blme/blme_0.01-4.tar.gz


getting the following answer: 
invalid package ‘0.01-4.tar.gz’


Can somebody help me? Are there other ways?
Thanks in advance,
Grazia

--

M. Grazia Pittau, Ph.D.
Associate Professor
Faculty of Statistics
Sapienza, University of Rome
P.le Aldo Moro, 5
00185 ROMA, ITALY

gra...@stat.columbia.edu
grazia.pit...@uniroma1.it
Phone: +39 06 49910782
Fax:   +39 06 4959241


__
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] installing a package from Archive

2012-10-01 Thread Duncan Murdoch


On 01/10/2012 7:48 AM, gra...@stat.columbia.edu wrote:

Dear all,
I'm trying to install in R the package blme that has been removed from R
cran but it is still in the Archive:
http://cran.r-project.org/src/contrib/Archive/blme/

I have a Mac OSX and I have been trying with the following command from
the Shell:

R CMD INSTALL --build
http://cran.r-project.org/src/contrib/Archive/blme/blme_0.01-4.tar.gz

getting the following answer:
invalid package ‘0.01-4.tar.gz’

Can somebody help me? Are there other ways?


You need a file, not a URL. So first download the file, then install it. 
(And you don't need the --build command line option.)


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] installing a package from Archive

2012-10-01 Thread Prof Brian Ripley


On 01/10/2012 12:48, gra...@stat.columbia.edu wrote:

Dear all,
I'm trying to install in R the package blme that has been removed from R
cran but it is still in the Archive:
http://cran.r-project.org/src/contrib/Archive/blme/

I have a Mac OSX and I have been trying with the following command from
the Shell:

R CMD INSTALL --build
http://cran.r-project.org/src/contrib/Archive/blme/blme_0.01-4.tar.gz

getting the following answer: invalid package ‘0.01-4.tar.gz’

Can somebody help me? Are there other ways?


You have to download the file first.  'INSTALL' does not know about 
URLs, as reading ?INSTALL would have reminded you.


 Both ‘lib’ and the elements of ‘pkgs’ may be absolute or relative
 path names of directories.  ‘pkgs’ may also contain names of
 package archive files: these are then extracted to a temporary
 directory.  These are tarballs containing a single directory,
 optionally compressed by ‘gzip’, ‘bzip2’, ‘xz’ or ‘compress’.
 Finally, binary package archive files (as created by ‘R CMD
 INSTALL --build’) can be supplied.




Thanks in advance,
Grazia



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




--
Brian D. Ripley,  rip...@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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

Re: [R] merge.zoo returns unmatched dates

2012-10-01 Thread arun

HI,

You can also try this:
Vup-read.table(text=
    Date, Velocity_m/s
 2010-01-21 07:42:00,    1.217943
 2010-01-21 07:43:00,    1.624395
 2010-01-21 07:44:00,    1.526379
 2010-01-21 07:45:00,    1.456831
 2010-01-21 07:46:00,    1.245390
 2010-01-21 07:47:00,    1.374330
,sep=,,header=TRUE,stringsAsFactors=FALSE)
 
PAS-read.table(text=
    Date,   PAS
 2010-01-21 05:01:00,  0.0013938
 2010-01-21 05:02:00,  0.0015331
 2010-01-21 05:03:00,  0.0016725
 2010-01-21 05:04:00,  0.0016725
 2010-01-21 05:05:00,  0.0012265
 2010-01-21 05:06:00,  0.0015889
,sep=,,header=TRUE,stringsAsFactors=FALSE)

library(xts)
PAS$Date-as.POSIXct(PAS$Date,format=%Y-%m-%d %H:%M:%S,tz=UTC)
Vup$Date-as.POSIXct(Vup$Date,format=%Y-%m-%d %H:%M:%S,tz=UTC)
 Vupxt-xts(Vup[,2],order.by=Vup[,1],tzone=UTC)
 PASxt-xts(PAS[,2],order.by=PAS[,1],tzone=UTC)
 VUPPASxt- merge(Vupxt,PASxt)
 VUPPASzoo-zoo(VUPPASxt)
VUPPASzoo
#   Vupxt PASxt
#2010-01-21 05:01:00   NA 0.0013938
#2010-01-21 05:02:00   NA 0.0015331
#2010-01-21 05:03:00   NA 0.0016725
#2010-01-21 05:04:00   NA 0.0016725
#2010-01-21 05:05:00   NA 0.0012265
#2010-01-21 05:06:00   NA 0.0015889
#2010-01-21 07:42:00 1.217943    NA
#2010-01-21 07:43:00 1.624395    NA
#2010-01-21 07:44:00 1.526379    NA
#2010-01-21 07:45:00 1.456831    NA
#2010-01-21 07:46:00 1.245390    NA
#2010-01-21 07:47:00 1.374330    NA

str(VUPPASzoo)
#‘zoo’ series from 2010-01-21 05:01:00 to 2010-01-21 07:47:00
 # Data: num [1:12, 1:2] NA NA NA NA NA ...
 #- attr(*, dimnames)=List of 2
  #..$ : chr [1:12] 2010-01-21 05:01:00 2010-01-21 05:02:00 2010-01-21 
05:03:00 2010-01-21 05:04:00 ...
  #..$ : chr [1:2] Vupxt PASxt
  #Index:  POSIXct[1:12], format: 2010-01-21 05:01:00 2010-01-21 05:02:00 
...


A.K.

 




- Original Message -
From: Vindoggy ! vindo...@hotmail.com
To: r-help@r-project.org
Cc: 
Sent: Monday, October 1, 2012 2:29 AM
Subject: [R] merge.zoo returns unmatched dates


Sorry for the lack of reproducible data, but this seems to be a problem 
inherent to my dataset and I can't figure out where the issue is. 

I have several data frames set up as a time series with identical POSIXct date 
formats. If I keep the original data in data frame format and merge them using 
base merge- everything is perfect and everyone is happy.

If I transform the data frames to zoo objects, and then do a merge.zoo- the 
data seem to become uncoupled from the original data. Even more unusual is that 
some dates in the new merged data set  are prior to the original data set. I've 
attempted bellow to show what this looks like, and I hope someone has a 
suggestion as to what may be causing the problem.

Here is one set of data in data.frame format

head(Vup)
                 Date Velocity_m/s
1 2010-01-21 07:42:00     1.217943
2 2010-01-21 07:43:00     1.624395
3 2010-01-21 07:44:00     1.526379
4 2010-01-21 07:45:00     1.456831
5 2010-01-21 07:46:00     1.245390
6 2010-01-21 07:47:00     1.374330

str(Vup)
'data.frame':    7168 obs. of  2 variables:
$ Date        : POSIXct, format: 2010-01-21 07:42:00 2010-01-21 07:43:00 ...
$ Velocity_m/s: num  1.22 1.62 1.53 1.46 1.25 ...

And here is a second in data.frame format:

head(PAS)
                 Date               PAS
1 2010-01-21 05:01:00   0.0013938
2 2010-01-21 05:02:00   0.0015331
3 2010-01-21 05:03:00   0.0016725
4 2010-01-21 05:04:00   0.0016725
5 2010-01-21 05:05:00   0.0012265
6 2010-01-21 05:06:00   0.0015889

str(PAS)
'data.frame':    5520 obs. of  2 variables:
$ Date       : POSIXct, format: 2010-01-21 05:01:00 2010-01-21 05:02:00 ...
$ PAS: num  0.00139 0.00153 0.00167 0.00167 0.00123 ...



Using zoo:

PASmin-zoo(as.matrix(PAS[,2]),as.POSIXct(PAS[,1],format=%Y-%m-%d 
%H:%M:%S,tz=UTC))

str(PASmin)
‘zoo’ series from 2010-01-21 05:01:00 to 2010-01-27 13:01:00
  Data: num [1:5520, 1] 0.00139 0.00153 0.00167 0.00167 0.00123 ...
- attr(*, dimnames)=List of 2
  ..$ : NULL
  ..$ : chr PAS
  Index:  POSIXct[1:5520], format: 2010-01-21 05:01:00 2010-01-21 05:02:00 
2010-01-21 05:03:00 ...




ADP_UPmin-zoo(as.matrix(Vup[,2]),as.POSIXct(Vup[,1], format=%Y-%m-%d 
%H:%M,tz=UTC))

str(ADP_UPmin)
‘zoo’ series from 2010-01-21 07:42:00 to 2010-01-26 20:12:00
  Data: num [1:7168, 1] 1.22 1.62 1.53 1.46 1.25 ...
- attr(*, dimnames)=List of 2
  ..$ : NULL
  ..$ : chr UP_Velocity_m/s
  Index:  POSIXct[1:7168], format: 2010-01-21 07:42:00 2010-01-21 07:43:00 
2010-01-21 07:44:00 ...


And if I merge the two zoo objects I get this:

M-merge(ADP_UPmin,PASmin)
head(M)

                    UP_Velocity_m/s       PAS
2010-01-20 21:01:00              NA 0.0013938
2010-01-20 21:02:00              NA 0.0015331
2010-01-20 21:03:00              NA 0.0016725
2010-01-20 21:04:00              NA 0.0016725
2010-01-20 21:05:00              NA 0.0012265
2010-01-20 21:06:00              NA 0.0015889


‘zoo’ series from 2010-01-20 21:01:00 to 2010-01-27 05:01:00
  Data: num [1:8499, 1:2] NA NA NA NA

Re: [R] lme help configuring random effects

2012-10-01 Thread Ben Bolker

Julie Lee-Yaw julleeyaw at yahoo.ca writes:


 

[snip]
 
 I am trying to run a mixed effects model in R using the lme
 package. My experiment is such that I am interested in the effects
 of Temperature (2 levels) and Species (3 levels) on Growth. I
 collected individuals from three populations within each
 species. Because individuals within a population are potentially
 more similar to each other than individuals among populations, I
 want to include population as a random factor in my model. 

 You may have some practical difficulties incorporating a random
effect with only three populations ...

 
 I would have thought that I would structure the model as follows: 
 z-lme(Growth~Temp*Species, random=~1|Species/Population) 
 But the summary for this model includes NAs (e.g. for two of the species). 
 I've considered a model such as 

  z- lme(Growth~Temp*Species,random=~1|Population) 

I think you need

  z - lme(Growth~Temp*Species, random=~1|Species:Population)

 Your specification (Species/Population) expands to Species+
Species:Population , which ends up including Species as both
a random and a fixed factor.

  You probably don't have the data, but in principle you
should consider  random=~Temp|Species:Population (see a paper
by Schielzeth on incorporating treatment-by-block interactions)

[snip]

 I'm also confused as to the naming of populations. 
 Currently I've got the populations named 1,2,3,4,5
 ...9. Should I be naming them A1,A2,A3, B1,B2,B3, C1,C2, C3 
 (where the letters represent different
 species)? When does that matter? 

  This is the difference between implicit and explicit
nesting.  In general naming them A1, ... C1, ... is better,
because it reduces the probability of making mistakes.

  http://glmm.wikidot.com/faq may be useful.
  Further questions should probably go to the r-sig-mixed-models
mailing list.

__
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] svyboxplot - library (survey)

2012-10-01 Thread Muhuri, Pradip (SAMHSA/CBHSQ)

Hello,

I have used the library (survey) package for boxplots using the following code. 

Could anyone please tell me why I am getting only 1  boxplot instead of 2 
boxplots (1-SPD,  2-No SPD).  

What changes in the following code would be required to get 2 boxplots in the 
same plot frame?

Thanks,

Pradip

###
nhis - svydesign (id=~psu, strat=~stratum, weights=~wt8, 
data=tor, nest=TRUE)

svyboxplot (dthage~xspd2, subset (nhis, mortstat==1), col=gray80,
 varwidth=TRUE, ylab=Age at Death, xlab=SPD Status: 1-SPD, 2=No 
SPD) 


Pradip K. Muhuri
Statistician
Substance Abuse  Mental Health Services Administration
The Center for Behavioral Health Statistics and Quality
Division of Population Surveys
1 Choke Cherry Road, Room 2-1071
Rockville, MD 20857
 
Tel: 240-276-1070
Fax: 240-276-1260
e-mail: pradip.muh...@samhsa.hhs.gov
 
The Center for Behavioral Health Statistics and Quality your feedback.  Please 
click on the following link to complete a brief customer survey:   
http://cbhsqsurvey.samhsa.gov

vide 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] installing a package from Archive

2012-10-01 Thread grazia


On Mon, 1 Oct 2012, Duncan Murdoch wrote:
Thanks a lot.
Following your suggestions, I have done the following:


zercona:Downloads graziapittau$ R CMD INSTALL blme 0.01-4.tar.gz
Warning: invalid package ‘0.01-4.tar.gz’
* installing to library 
‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library’

* installing *source* package ‘blme’ ...
** package ‘blme’ successfully unpacked and MD5 sums checked
** libs
*** arch - i386
sh: make: command not found
ERROR: compilation failed for package ‘blme’
* removing 
‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library/blme’
* restoring previous 
‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library/blme’



Unfortunately, I still get an error message and from Rthe following:

library(blme)
Error in library(blme) : ‘blme’ is not a valid installed package

Where am I doing wrong?
Thanks a lot.
Grazia


On 01/10/2012 7:48 AM, gra...@stat.columbia.edu wrote:

 Dear all,
 I'm trying to install in R the package blme that has been removed from R
 cran but it is still in the Archive:
 http://cran.r-project.org/src/contrib/Archive/blme/

 I have a Mac OSX and I have been trying with the following command from
 the Shell:

 R CMD INSTALL --build
 http://cran.r-project.org/src/contrib/Archive/blme/blme_0.01-4.tar.gz

 getting the following answer:
 invalid package ‘0.01-4.tar.gz’

 Can somebody help me? Are there other ways?


You need a file, not a URL. So first download the file, then install it. (And 
you don't need the --build command line option.)


Duncan Murdoch



--

M. Grazia Pittau, Ph.D.
Associate Professor
Faculty of Statistics
Sapienza, University of Rome
P.le Aldo Moro, 5
00185 ROMA, ITALY

gra...@stat.columbia.edu
grazia.pit...@uniroma1.it
Phone: +39 06 49910782
Fax:   +39 06 4959241


__
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] installing a package from Archive

2012-10-01 Thread Duncan Murdoch


On 01/10/2012 10:02 AM, gra...@stat.columbia.edu wrote:

On Mon, 1 Oct 2012, Duncan Murdoch wrote:
Thanks a lot.
Following your suggestions, I have done the following:


zercona:Downloads graziapittau$ R CMD INSTALL blme 0.01-4.tar.gz
Warning: invalid package ‘0.01-4.tar.gz’


The filename should not have a space in it.  I think the original had an 
underscore; if that's been lost, you need to rename the file to 
something with no spaces.  (You can probably quote the filename, but the 
details of how to do that depend on your shell; it's easier just to 
avoid filenames that contain spaces.)


Duncan Murdoch


* installing to library
‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library’
* installing *source* package ‘blme’ ...
** package ‘blme’ successfully unpacked and MD5 sums checked
** libs
*** arch - i386
sh: make: command not found
ERROR: compilation failed for package ‘blme’
* removing
‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library/blme’
* restoring previous
‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library/blme’


Unfortunately, I still get an error message and from Rthe following:

library(blme)
Error in library(blme) : ‘blme’ is not a valid installed package

Where am I doing wrong?
Thanks a lot.
Grazia

 On 01/10/2012 7:48 AM, gra...@stat.columbia.edu wrote:
  Dear all,
  I'm trying to install in R the package blme that has been removed from R
  cran but it is still in the Archive:
  http://cran.r-project.org/src/contrib/Archive/blme/

  I have a Mac OSX and I have been trying with the following command from
  the Shell:

  R CMD INSTALL --build
  http://cran.r-project.org/src/contrib/Archive/blme/blme_0.01-4.tar.gz

  getting the following answer:
  invalid package ‘0.01-4.tar.gz’

  Can somebody help me? Are there other ways?

 You need a file, not a URL. So first download the file, then install it. (And
 you don't need the --build command line option.)

 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] svyboxplot - library (survey)

2012-10-01 Thread Anthony Damico

using a slight modification of the example shown in ?svyboxplot


# load survey library
library(survey)

# load example data
data(api)

# create an example svydesign
dstrat - svydesign(id = ~1, strata = ~stype, weights = ~pw, data =
apistrat,
fpc = ~fpc)

# set the plot window to display 1 plot x 2 plots
par(mfrow=c(1,2))

# generate two example boxplots
svyboxplot(enroll~stype,dstrat,all.outliers=TRUE)
svyboxplot(enroll~1,dstrat)

# done



# alternative: not as nice

# set the plot window to display 2 plots x 1 plot
par(mfrow=c(2,1))

# generate two example boxplots
svyboxplot(enroll~stype,dstrat,all.outliers=TRUE)
svyboxplot(enroll~1,dstrat)

# done







On Mon, Oct 1, 2012 at 9:50 AM, Muhuri, Pradip (SAMHSA/CBHSQ) 
pradip.muh...@samhsa.hhs.gov wrote:

 Hello,

 I have used the library (survey) package for boxplots using the following
 code.

 Could anyone please tell me why I am getting only 1  boxplot instead of 2
 boxplots (1-SPD,  2-No SPD).

 What changes in the following code would be required to get 2 boxplots in
 the same plot frame?

 Thanks,

 Pradip

 ###
 nhis - svydesign (id=~psu, strat=~stratum, weights=~wt8,
 data=tor, nest=TRUE)

 svyboxplot (dthage~xspd2, subset (nhis, mortstat==1), col=gray80,
  varwidth=TRUE, ylab=Age at Death, xlab=SPD Status: 1-SPD,
 2=No SPD)


 Pradip K. Muhuri
 Statistician
 Substance Abuse  Mental Health Services Administration
 The Center for Behavioral Health Statistics and Quality
 Division of Population Surveys
 1 Choke Cherry Road, Room 2-1071
 Rockville, MD 20857

 Tel: 240-276-1070
 Fax: 240-276-1260
 e-mail: pradip.muh...@samhsa.hhs.gov

 The Center for Behavioral Health Statistics and Quality your feedback.
 Please click on the following link to complete a brief customer survey:
 http://cbhsqsurvey.samhsa.gov

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


[[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] adjust font in ggplot2 to LaTeX document

2012-10-01 Thread Brian Diggs


On 9/30/2012 11:38 AM, Jonas Stein wrote:

Hi,

how can i adjust the font in a ggplot2 qplot so that it will look
similar to the LaTeX font?
Computer Modern Sans Serif in the same size would be nice.

My output device is
ggsave(filename=test.pdf, width=5.5, height=3, dpi=300)
and i will include the graphic with 5.5 inch in LaTeX.

I found some pages about that topic but all solutions that i found have
been very complicate and confusing to me.
But the articles have been from 2009 too, so i guess there will be
some easy solution today...


There is a recent (late last month) blog post on this topic:

http://blog.revolutionanalytics.com/2012/09/how-to-use-your-favorite-fonts-in-r-charts.html


Is theme the keyword to find the solution?


In ggplot2, the font family is controlled by the theme. font family 
will probably get more specific results.



Or will i have to include a binary font file?


Possibly, but probably not since you would be using a font already 
present in the (final) PDF.



Can someone give me a link to an example?

Kind regards and thanks a lot,




--
Brian S. Diggs, PhD
Senior Research Associate, Department of Surgery
Oregon Health  Science University

__
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] Better way of Grouping?

2012-10-01 Thread Charles Determan Jr

Thank you Jeff,

The main purpose I am looking to use these subsets for is to do comparisons
between groups and/or timepoints within groups.  For example the difference
in means between 3 and 4, or the percent difference between group L an D
within group 3.  The code I have provided is almost exactly as used
although David accurately noted I mistakenly typed '3' when I intended to
put '4'.  I am just looking for some general guidance in improving how I go
about my code.  Because there are so many subgroups that can be formed with
three layers of groups, it simply becomes a clutter and I wanted to learn
if there was a 'better' way to organize groups for comparisons.  Does that
help clarify?

Thanks again,
Charles

On Fri, Sep 28, 2012 at 4:25 PM, Jeff Newmiller jdnew...@dcn.davis.ca.uswrote:

 You have not specified the objective function you are trying to optimize
 with your term efficient, or what you do with all of these subsets once
 you have them.

 For notational simplification and completeness of coverage (not
 necessarily computational speedup) you might want to look at tapply or
 ddply/dlply from the plyr package. If you build lists of subsets you can
 index into them according to grouping value. You can use expand.grid to
 build all permutations of grouping values to use as indexes into those
 lists of subsets.

 To reiterate, you have not indicated what you want to do with these
 subsets, so there could be special-purpose functions that do what you want.
  As always, reproducible code leads to reproducible answers. :)
 ---
 Jeff NewmillerThe .   .  Go Live...
 DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live
 Go...
   Live:   OO#.. Dead: OO#..  Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
 ---
 Sent from my phone. Please excuse my brevity.

 Charles Determan Jr deter...@umn.edu wrote:

 Hello R users,
 
 This is more of a convenience question that I hope others might find
 useful
 if there is a better answer.  I work with large datasets that requires
 multiple parsing stages for different analysis.  For example, compare
 group
 3 vs. group 4.  A more complicated comparison would be time B in group
 3 of
 group L with B in group 4 of group L.  I normally subset each group
 with
 the following type of code.
 
 data=read(...)
 
 #L v D
 L=data[LvD %in% c(L),]
 D=data[LvD %in% c(D),]
 
 #Groups 3 and 4 within L and D
 group3L=L[group %in% c(3),]
 group4L=L[group %in% c(3),]
 
 group3D=D[group %in% c(3),]
 group4D=D[group %in% c(3),]
 
 #Times B, S45, FR2, FR8
 you get the idea
 
 
 Is there a more efficient way to subset groups?  Thanks for any
 insight.
 
 Regards,
 Charles
 
[[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.



[[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] which() with multiple conditions

2012-10-01 Thread pdb

I hope someone can point me in the right direction please.

I have a data frame with a column containing names.  I want to identify the
columns that contain names in a list.

namestofind - c('fred','bill',a long list)

If I only wanted to identify a single name I would use

which(z$name == 'bill')

What syntax would I use to identify all the rows that contain any of the
names in namestofind?

Thanks in advance for the pointer





--
View this message in context: 
http://r.789695.n4.nabble.com/which-with-multiple-conditions-tp4644677.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] can't create png graphics under linux anymore

2012-10-01 Thread Tobias Pilz


Hi,

since my update of Linux from opensuse 11.4 to 12.1 and the update of R 
to 2.15.1 I can't produce graphics with png() or similar (like tiff, 
jpeg etc.) anymore. The following error occurs:


 png()
Error in X11(paste(png::, filename, sep = ), g$width, g$height, 
pointsize,  :

  unable to start Device PNG
Additionally: Warning:
In png() : png isn't supported of this R version

(error message translated from german)

Before the updates everything worked fine and still works under Windows 
7 for R 2.15.1 on the same machine. I couldn't find something useful 
with google. As far as I know R still supports png?! Maybe the problem 
is related to suse?


Thanks for any help,

Tobias

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

2012-10-01 Thread Karan Anand

hi,
I am new to using R ,I have 2 datasets with dates common in them ,how
can i take out the common dates within them.

karan

[[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] calculating probability from the density function

2012-10-01 Thread Eugene Kanshin

Hello,
I have a data x with normal (or very close to normal) distribution, I can
plot a density distribution with density(x,...). My question is is there
any way to calculate an area under this distribution (=probability) for
particular range of x values, let's say for x from 0 to 2? I was not able
to find any kind of simple procedure to do this.
Thanks in advance for your help,
Evgeny.

[[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] which() with multiple conditions

2012-10-01 Thread R. Michael Weylandt

On Monday, October 1, 2012, pdb wrote:

 I hope someone can point me in the right direction please.

 I have a data frame with a column containing names.  I want to identify the
 columns that contain names in a list.

 namestofind - c('fred','bill',a long list)

 If I only wanted to identify a single name I would use

 which(z$name == 'bill')

 What syntax would I use to identify all the rows that contain any of the
 names in namestofind?


%in%

Michael



 Thanks in advance for the pointer





 --
 View this message in context:
 http://r.789695.n4.nabble.com/which-with-multiple-conditions-tp4644677.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org javascript:; mailing list
 https://stat.ethz.ch/mailman/listinfo/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.

Re: [R] enquiry

2012-10-01 Thread R. Michael Weylandt

On Monday, October 1, 2012, Karan Anand wrote:

 hi,
 I am new to using R ,I have 2 datasets with dates common in them ,how
 can i take out the common dates within them.


Depending on what you mean by 'out' either merge() or setdiff()

RMW



 karan

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org javascript:; mailing list
 https://stat.ethz.ch/mailman/listinfo/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] Retrieve hypergeometric results in large scale


I'm going to use 

dhyper(x, m, n, k)

to get a 95% coverage. Let me use an example to explain my problem:

Suppose I have a urn containing 90 red and 10 black balls.
Now I wanna remove 3 from the urn. By the following codes:

m-90;n-10;k-3;
x-0:3
dhyper(x,m,n,k)

I can obtain the probability that 0,1,2,3 red balls will be removed.
 0.000742115 0.025046382 0.247680891 0.726530612

So 95% time, 2 to 3 red balls will be removed and the resultant composition
will be changed to 
87:10 or 88:9, the original percent of red balls will be changed from 90 to
89.69 to 90.72 then.

If now I have 50:50 and again to remove 3 balls, I will obtain the
probability as:
0.1212 0.3788 0.3788 0.1212

To get the resultant range of red balls for 95% time, this time all the
four cases have to consider and so the resultant change of red balls will
become 48.45 to 51.54

So my problem is, is there any convenient built-in function that helps
extract this 95% confidence interval-like data?






--
View this message in context: 
http://r.789695.n4.nabble.com/Retrieve-hypergeometric-results-in-large-scale-tp4644683.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] Retrieve hypergeometric results in large scale

2012-10-01 Thread Jeff Newmiller

Perhaps you should read

?dhyper

and if you have a hard time parsing that, then read

?Distributions

and then go back to

?dhyper
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

jas4710 wata...@post.com wrote:


I'm going to use 

dhyper(x, m, n, k)

to get a 95% coverage. Let me use an example to explain my problem:

Suppose I have a urn containing 90 red and 10 black balls.
Now I wanna remove 3 from the urn. By the following codes:

m-90;n-10;k-3;
x-0:3
dhyper(x,m,n,k)

I can obtain the probability that 0,1,2,3 red balls will be removed.
 0.000742115 0.025046382 0.247680891 0.726530612

So 95% time, 2 to 3 red balls will be removed and the resultant
composition
will be changed to 
87:10 or 88:9, the original percent of red balls will be changed from
90 to
89.69 to 90.72 then.

If now I have 50:50 and again to remove 3 balls, I will obtain the
probability as:
0.1212 0.3788 0.3788 0.1212

To get the resultant range of red balls for 95% time, this time all
the
four cases have to consider and so the resultant change of red balls
will
become 48.45 to 51.54

So my problem is, is there any convenient built-in function that helps
extract this 95% confidence interval-like data?






--
View this message in context:
http://r.789695.n4.nabble.com/Retrieve-hypergeometric-results-in-large-scale-tp4644683.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.

Re: [R] Retrieve hypergeometric results in large scale

Homework? There's a no homework policy on this list.

-- Bert

On Mon, Oct 1, 2012 at 8:10 AM, Jeff Newmiller jdnew...@dcn.davis.ca.us wrote:
 Perhaps you should read

 ?dhyper

 and if you have a hard time parsing that, then read

 ?Distributions

 and then go back to

 ?dhyper
 ---
 Jeff NewmillerThe .   .  Go Live...
 DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
   Live:   OO#.. Dead: OO#..  Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
 ---
 Sent from my phone. Please excuse my brevity.

 jas4710 wata...@post.com wrote:


I'm going to use

dhyper(x, m, n, k)

to get a 95% coverage. Let me use an example to explain my problem:

Suppose I have a urn containing 90 red and 10 black balls.
Now I wanna remove 3 from the urn. By the following codes:

m-90;n-10;k-3;
x-0:3
dhyper(x,m,n,k)

I can obtain the probability that 0,1,2,3 red balls will be removed.
 0.000742115 0.025046382 0.247680891 0.726530612

So 95% time, 2 to 3 red balls will be removed and the resultant
composition
will be changed to
87:10 or 88:9, the original percent of red balls will be changed from
90 to
89.69 to 90.72 then.

If now I have 50:50 and again to remove 3 balls, I will obtain the
probability as:
0.1212 0.3788 0.3788 0.1212

To get the resultant range of red balls for 95% time, this time all
the
four cases have to consider and so the resultant change of red balls
will
become 48.45 to 51.54

So my problem is, is there any convenient built-in function that helps
extract this 95% confidence interval-like data?






--
View this message in context:
http://r.789695.n4.nabble.com/Retrieve-hypergeometric-results-in-large-scale-tp4644683.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.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

__
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] Problem with nls regression fit

2012-10-01 Thread Gyanendra Pokharel

Hi all,
I got following problem in fitting the data.
Any kind of suggestions are welcome


 beta - 3.5
 d - seq(0.1,62.5,0.1)
 y - exp(-beta*d)
 y1 - y
 x - read.table(epidist.txt, header = TRUE)
 data.nls - as.data.frame(cbind(y1,x))
 #attach(data.nls)
 nls.fit - nls(y1~dist,data.nls)
Error in cll[[1L]] : object of type 'symbol' is not subsettable


Best

[[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] svyboxplot - library (survey)

2012-10-01 Thread Muhuri, Pradip (SAMHSA/CBHSQ)

Dear Anthony,

Yes, I can follow the example code you have given.  But, do you know from the 
code shown below (following Thomas Lumley's Complex Surveys) why I am getting 
the boxplot of dthage for just xspd=1, not xspd2=2?

My intent is the make this code work so that I can generate similar plots on 
other continuous variable.

Any help will be appreciated.

Thanks,

Pradip




nhis - svydesign (id=~psu, strat=~stratum, weights=~wt8,
data=tor, nest=TRUE)

svyboxplot (dthage~xspd2, subset (nhis, mortstat==1), col=gray80,
 varwidth=TRUE, ylab=Age at Death, xlab=SPD Status: 1-SPD, 2=No 
SPD)



Pradip K. Muhuri, PhD
Statistician
Substance Abuse  Mental Health Services Administration
The Center for Behavioral Health Statistics and Quality
Division of Population Surveys
1 Choke Cherry Road, Room 2-1071
Rockville, MD 20857

Tel: 240-276-1070
Fax: 240-276-1260
e-mail: pradip.muh...@samhsa.hhs.govmailto:pradip.muh...@samhsa.hhs.gov

The Center for Behavioral Health Statistics and Quality your feedback.  Please 
click on the following link to complete a brief customer survey:   
http://cbhsqsurvey.samhsa.govhttp://cbhsqsurvey.samhsa.gov/

From: Anthony Damico [mailto:ajdam...@gmail.com]
Sent: Monday, October 01, 2012 10:07 AM
To: Muhuri, Pradip (SAMHSA/CBHSQ)
Cc: R help
Subject: Re: [R] svyboxplot - library (survey)

using a slight modification of the example shown in ?svyboxplot


# load survey library
library(survey)

# load example data
data(api)

# create an example svydesign
dstrat - svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat,
fpc = ~fpc)

# set the plot window to display 1 plot x 2 plots
par(mfrow=c(1,2))

# generate two example boxplots
svyboxplot(enroll~stype,dstrat,all.outliers=TRUE)
svyboxplot(enroll~1,dstrat)

# done



# alternative: not as nice

# set the plot window to display 2 plots x 1 plot
par(mfrow=c(2,1))

# generate two example boxplots
svyboxplot(enroll~stype,dstrat,all.outliers=TRUE)
svyboxplot(enroll~1,dstrat)

# done






On Mon, Oct 1, 2012 at 9:50 AM, Muhuri, Pradip (SAMHSA/CBHSQ) 
pradip.muh...@samhsa.hhs.govmailto:pradip.muh...@samhsa.hhs.gov wrote:
Hello,

I have used the library (survey) package for boxplots using the following code.

Could anyone please tell me why I am getting only 1  boxplot instead of 2 
boxplots (1-SPD,  2-No SPD).

What changes in the following code would be required to get 2 boxplots in the 
same plot frame?

Thanks,

Pradip

###
nhis - svydesign (id=~psu, strat=~stratum, weights=~wt8,
data=tor, nest=TRUE)

svyboxplot (dthage~xspd2, subset (nhis, mortstat==1), col=gray80,
 varwidth=TRUE, ylab=Age at Death, xlab=SPD Status: 1-SPD, 2=No 
SPD)


Pradip K. Muhuri
Statistician
Substance Abuse  Mental Health Services Administration
The Center for Behavioral Health Statistics and Quality
Division of Population Surveys
1 Choke Cherry Road, Room 2-1071
Rockville, MD 20857

Tel: 240-276-1070
Fax: 240-276-1260
e-mail: pradip.muh...@samhsa.hhs.govmailto:pradip.muh...@samhsa.hhs.gov

The Center for Behavioral Health Statistics and Quality your feedback.  Please 
click on the following link to complete a brief customer survey:   
http://cbhsqsurvey.samhsa.gov

vide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.orgmailto: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.

Re: [R] calculating probability from the density function

Forget density(). Smooth the ecdf instead.

?lowess
?smooth.spline
?monotone.smooth (in package fda, for monotone smoothing, which may be
preferable)
?locfit (in package locfit)

... and many many others

-- Bert

On Mon, Oct 1, 2012 at 6:46 AM, Eugene Kanshin kanshin...@gmail.com wrote:
 Hello,
 I have a data x with normal (or very close to normal) distribution, I can
 plot a density distribution with density(x,...). My question is is there
 any way to calculate an area under this distribution (=probability) for
 particular range of x values, let's say for x from 0 to 2? I was not able
 to find any kind of simple procedure to do this.
 Thanks in advance for your help,
 Evgeny.

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



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

__
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] can't create png graphics under linux anymore

2012-10-01 Thread Prof Brian Ripley


See the 'R Installation and Administration Manual'.

And yes, it is related to the actions of whoever installed R, so please 
discuss it with them.


On 01/10/2012 14:54, Tobias Pilz wrote:

Hi,

since my update of Linux from opensuse 11.4 to 12.1 and the update of R
to 2.15.1 I can't produce graphics with png() or similar (like tiff,
jpeg etc.) anymore. The following error occurs:

  png()
Error in X11(paste(png::, filename, sep = ), g$width, g$height,
pointsize,  :
   unable to start Device PNG
Additionally: Warning:
In png() : png isn't supported of this R version

(error message translated from german)

Before the updates everything worked fine and still works under Windows
7 for R 2.15.1 on the same machine. I couldn't find something useful
with google. As far as I know R still supports png?! Maybe the problem
is related to suse?

Thanks for any help,

Tobias

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



--
Brian D. Ripley,  rip...@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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

Re: [R] calculating probability from the density function


On Oct 1, 2012, at 6:46 AM, Eugene Kanshin wrote:

 Hello,
 I have a data x with normal (or very close to normal) distribution, I can
 plot a density distribution with density(x,...). My question is is there
 any way to calculate an area under this distribution (=probability) for
 particular range of x values, let's say for x from 0 to 2? I was not able
 to find any kind of simple procedure to do this.

? pnorm

-- 

David Winsemius, MD
Alameda, CA, USA

__
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] Hmisc describe error

2012-10-01 Thread Bond, Stephen

Describe fails for me with a message similar to what was an issue in 2008 and 
got fixed according to posts.

R version 2.15.0 (2012-03-30)
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-pc-mingw32/i386 (32-bit)

# output truncated

 options(chmhelp = FALSE, help_type = text)
 .help.ESS - help
 options(STERM='iESS', editor='gnuclient.exe')
 load('prostate.sav')
 library(rms)
Loading required package: Hmisc
Loading required package: survival
Loading required package: splines
Hmisc library by Frank E Harrell Jr

Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview')
to see overall documentation.

NOTE:Hmisc no longer redefines [.factor to drop unused levels when
subsetting.  To get the old behavior of Hmisc type dropUnusedLevels().


Attaching package: 'Hmisc'

The following object(s) are masked from 'package:survival':

untangle.specials

The following object(s) are masked from 'package:base':

format.pval, round.POSIXt, trunc.POSIXt, units


Attaching package: 'rms'

The following object(s) are masked from 'package:survival':

Surv

 describe(prostate)
Error in format(dates(x)) : could not find function dates


Thanks everybody.
Stephen 

__
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] nlme: spatial autocorrelation on a sphere


On Oct 1, 2012, at 8:34 AM, Spencer Graves wrote:

 On 10/1/2012 12:38 AM, David Winsemius wrote:
 snip
 
 LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html
 
 
 
  What's LMCTVIFY?

Please accept my apologies for attempting to be cute and I also apologize to Mr 
Bebber for reading his request far too quickly. I was trying to use (C)RAN 
(T)ask (V)iews as a verb.

  LMGTFY for LMCTVIY produced only LMGTFY, at least for me. There's space 
 for it on Wikipedia (RTFM) after LMGTFY (which I found using LMGTFY 
 Wikipedia).

I don't think it warrants being enshrined. I might also  have written:  
require(sos); findFn(metric spherical latitude longitude)  might produce. In 
this instance that approach found the ramps::corRSpher function, which has a 
haversine metric, much  more quickly than did my subsequent efforts with 
RSeek, which I conducted after I realized that Bebber had already made a good 
faith effort at identifying resources.

-- 


David Winsemius, MD
Alameda, CA, USA

__
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] Retrieve hypergeometric results in large scale

Thanks Jeff 
The documentation pages, if I haven't missed any crucial points, illustrate
how to get probability and cumulative probability values.

I can first retrieve the data structures and use Perl (I don't know how to
use R...) to sort the derived ratios and sum the probability values until
the cumulative probability exceeds 95%. Well, I just don't know whether such
seemingly routine procedures have already been implemented...

Thanks again!



--
View this message in context: 
http://r.789695.n4.nabble.com/Retrieve-95-coverage-of-results-from-a-hypergeometric-distribution-tp4644683p4644701.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] Error messages when attempting to calculate polychoric correlation matrices

2012-10-01 Thread Kevin Cheung

Dear R users,

I am a psychology postgraduate student who is relatively new to using R. I am 
currently developing a psychometric scale and have run into a few problems when 
using R to calculate a polychoric correlation matrix for my dataset. I am 
trying to produce a polychoric correlation matrix for calculating ordinal 
reliability estimates (eg. Alpha, omega).The set consists of 439 observations 
of 18 variables. Each of the 18 variables was measured using a 6 point Likert 
scale but some variables only feature responses across 5 ordinal steps because 
there were no observed responses for one extreme of the item. The dataset has 
no missing data following multiple imputation and preliminary analyses were 
conducted using SPSS. I have included a command that will reproduce the dataset 
below, apologies for the size of the command:

data.matrix18 -
structure(list(SABAS02 = structure(c(6, 2, 5, 6, 4, 6, 5, 3,
5, 1, 4, 4, 6, 5, 5, 5, 4, 4, 4, 3, 4, 6, 5, 2, 5, 4, 5, 4, 4,
5, 5, 5, 2, 5, 4, 5, 5, 5, 6, 3, 6, 4, 5, 5, 5, 5, 4, 4, 5, 5,
5, 4, 5, 5, 4, 6, 6, 5, 6, 4, 4, 5, 4, 5, 5, 5, 6, 4, 4, 5, 4,
6, 5, 6, 5, 6, 6, 5, 4, 6, 6, 5, 3, 4, 2, 4, 5, 5, 6, 2, 6, 5,
4, 6, 2, 5, 3, 5, 4, 4, 5, 4, 4, 3, 4, 4, 6, 5, 6, 6, 4, 4, 5,
5, 5, 6, 5, 5, 5, 5, 4, 5, 5, 5, 4, 4, 5, 6, 4, 5, 5, 4, 5, 5,
4, 5, 5, 4, 5, 5, 5, 4, 2, 5, 4, 5, 5, 5, 5, 5, 5, 4, 5, 4, 4,
5, 4, 6, 2, 4, 6, 4, 4, 5, 4, 4, 6, 6, 5, 3, 3, 2, 4, 5, 5, 5,
5, 6, 2, 5, 5, 4, 6, 4, 5, 6, 5, 4, 5, 5, 4, 4, 5, 4, 5, 5, 3,
5, 4, 3, 5, 4, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 5, 4, 4, 5, 5, 6,
4, 3, 5, 6, 5, 4, 5, 5, 2, 5, 5, 4, 2, 4, 5, 5, 5, 5, 5, 3, 4,
5, 2, 4, 3, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 5, 4, 5, 6, 4, 5, 3,
5, 6, 5, 4, 5, 5, 5, 5, 5, 4, 4, 5, 6, 5, 4, 5, 5, 6, 6, 6, 4,
3, 4, 4, 5, 6, 4, 5, 5, 4, 6, 5, 5, 4, 5, 4, 4, 4, 4, 5, 3, 5,
5, 4, 3, 4, 5, 4, 4, 5, 4, 5, 5, 3, 5, 6, 5, 5, 5, 5, 5, 4, 3,
5, 4, 4, 5, 5, 6, 4, 5, 4, 2, 6, 5, 5, 5, 5, 5, 5, 5, 6, 6, 5,
5, 5, 4, 5, 5, 6, 5, 4, 4, 5, 4, 4, 5, 5, 4, 6, 4, 4, 5, 6, 5,
5, 4, 5, 4, 5, 4, 6, 5, 5, 5, 4, 5, 6, 5, 5, 5, 5, 5, 4, 3, 4,
5, 4, 4, 6, 4, 4, 4, 2, 5, 4, 4, 6, 2, 6, 4, 5, 5, 6, 5, 6, 4,
4, 3, 4, 5, 5, 6, 5, 5, 5, 5, 4, 5, 3, 5, 5, 5, 5, 6, 3, 5, 6,
5, 6, 4, 5, 4, 4, 5, 6, 2, 5, 6), value.labels = structure(c(6,
5, 4, 3, 2, 1), .Names = c(Strongly Agree, Agree, Slightly Agree,
Slightly Disagree, Disagree, Strongly Disagree))), SABAS06 = 
structure(c(6,
6, 6, 6, 6, 6, 5, 5, 6, 3, 6, 6, 6, 5, 5, 6, 5, 6, 6, 5, 5, 6,
6, 6, 6, 5, 6, 3, 6, 5, 5, 5, 5, 6, 6, 6, 3, 6, 6, 5, 5, 4, 5,
6, 6, 5, 6, 5, 6, 6, 6, 6, 6, 6, 5, 5, 6, 6, 6, 6, 6, 6, 5, 6,
6, 6, 6, 6, 4, 6, 5, 5, 6, 6, 6, 6, 6, 6, 3, 5, 6, 6, 6, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 6, 6, 6, 5, 5, 5, 6, 5, 6,
5, 6, 6, 6, 6, 5, 5, 6, 5, 5, 6, 2, 5, 5, 6, 5, 5, 5, 6, 5, 6,
5, 5, 6, 6, 6, 5, 5, 5, 6, 6, 5, 6, 6, 6, 4, 6, 6, 5, 4, 6, 5,
5, 6, 5, 5, 5, 6, 6, 6, 6, 5, 5, 6, 6, 5, 6, 4, 6, 6, 6, 6, 5,
5, 6, 4, 5, 6, 5, 6, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5,
6, 6, 5, 5, 6, 6, 5, 6, 5, 6, 5, 4, 6, 5, 6, 6, 5, 6, 5, 6, 5,
5, 5, 5, 6, 6, 5, 6, 5, 4, 6, 5, 6, 6, 5, 5, 6, 5, 5, 6, 6, 6,
6, 5, 5, 6, 6, 5, 5, 6, 6, 4, 5, 5, 6, 6, 6, 6, 4, 6, 6, 6, 6,
6, 5, 6, 6, 6, 6, 5, 6, 6, 6, 5, 5, 6, 5, 6, 6, 6, 5, 2, 6, 6,
5, 5, 6, 5, 6, 6, 5, 6, 5, 6, 5, 6, 6, 6, 5, 5, 4, 5, 6, 6, 5,
5, 4, 4, 6, 6, 5, 4, 5, 5, 6, 6, 6, 5, 4, 5, 5, 3, 5, 5, 6, 6,
5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 4, 5, 5, 6,
5, 6, 6, 6, 5, 6, 5, 6, 6, 6, 5, 5, 6, 5, 6, 6, 6, 6, 5, 6, 6,
5, 6, 5, 6, 6, 5, 5, 5, 6, 5, 5, 5, 5, 6, 5, 6, 6, 5, 5, 6, 5,
5, 6, 6, 4, 3, 3, 4, 3, 5, 4, 5, 5, 6, 6, 6, 6, 6, 6, 2, 5, 6,
5, 6, 6, 6, 6, 6, 6, 5, 6, 5, 6, 6, 6, 6, 5, 6, 5, 5, 5, 6, 6,
6, 6, 5, 6, 5, 6, 5, 6, 6, 6, 5, 4, 6, 6, 6, 6, 6, 5), value.labels = 
structure(c(6,
5, 4, 3, 2, 1), .Names = c(Strongly Agree, Agree, Slightly Agree,
Slightly Disagree, Disagree, Strongly Disagree))), SABAS09 = 
structure(c(5,
3, 5, 5, 4, 5, 4, 2, 5, 1, 4, 5, 5, 5, 3, 5, 4, 4, 6, 4, 4, 6,
5, 6, 4, 4, 4, 4, 4, 4, 5, 5, 3, 4, 5, 5, 2, 2, 5, 4, 5, 4, 6,
5, 6, 4, 3, 6, 5, 5, 5, 5, 6, 5, 4, 5, 6, 6, 5, 5, 5, 5, 5, 5,
4, 5, 6, 5, 5, 4, 4, 2, 5, 5, 5, 5, 6, 1, 4, 5, 5, 5, 4, 4, 5,
5, 5, 5, 5, 5, 6, 5, 4, 5, 6, 4, 5, 5, 5, 5, 5, 4, 4, 5, 2, 4,
5, 6, 4, 6, 5, 5, 5, 5, 4, 4, 4, 3, 5, 3, 1, 4, 4, 4, 4, 5, 5,
5, 4, 5, 4, 4, 5, 3, 5, 6, 5, 3, 4, 5, 5, 3, 3, 5, 5, 5, 5, 5,
5, 5, 5, 2, 5, 4, 4, 4, 5, 5, 4, 5, 5, 2, 5, 3, 2, 5, 4, 5, 4,
3, 4, 2, 4, 5, 5, 4, 4, 5, 5, 5, 4, 5, 5, 5, 4, 5, 5, 5, 5, 4,
4, 4, 6, 4, 5, 5, 2, 4, 3, 4, 5, 4, 5, 5, 5, 5, 4, 3, 5, 3, 5,
5, 5, 4, 5, 4, 4, 6, 4, 5, 3, 5, 5, 5, 5, 2, 5, 3, 4, 4, 3, 5,
6, 5, 4, 5, 4, 4, 4, 4, 5, 2, 3, 5, 5, 5, 5, 4, 3, 5, 5, 5, 6,
5, 4, 4, 5, 4, 5, 4, 5, 6, 5, 3, 3, 5, 5, 4, 4, 4, 3, 4, 6, 5,
5, 4, 5, 5, 5, 6, 5, 4, 6, 4, 4, 4, 4, 5, 5, 2, 4, 6, 6, 6, 4,
3, 3, 5, 4, 6, 3, 5, 5, 3, 3, 6, 5, 3, 4, 5, 4, 4, 4, 4, 4, 5,
5, 4, 5, 4, 6, 4, 4, 5, 5, 5, 4, 5, 2, 5, 2, 3, 2, 5, 5, 4, 5,
5, 5, 5, 5, 5, 5, 4, 5, 5, 6, 3, 4, 5, 4, 5, 4, 4, 3, 4, 5, 4,
4, 5, 3, 5, 6, 1, 5, 4, 5, 5, 3, 4, 4, 5, 2, 4,

Re: [R] Retrieve hypergeometric results in large scale

Hi Bert. This is not a homework. If I can do some basic programming in R like
Perl, then I'll have a better chance to accomplish this task but the matrix
concept is not quickly comprehensible...



--
View this message in context: 
http://r.789695.n4.nabble.com/Retrieve-95-coverage-of-results-from-a-hypergeometric-distribution-tp4644683p4644703.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] Transform pairwise observations into a table

2012-10-01 Thread AHJ

Hi,

I have a table of pairs of individuals and a coefficient that belongs to the
pair:

ind1ind2coef
1   1   1
1   2   0.25
1   3   0.125
1   4   0.5
2   2   1
2   1   0.25
2   3   0.125
2   4   0.5
3   3   1
3   1   0.125
3   2   0.125
3   4   0.5
4   4   1
4   1   0.5
4   2   0.5
4   3   0.5

And I want to convert it to a matrix where each individual is both a row and
a column and at the intersection of each pair is the coefficient that
belongs to that pair:

1  2   3 4
1   1  0.250.1250.5
2   0.25  10.1250.5
3   0.1250.1251 0.5
4   0.50.5 0.5  1

If table() would allow me to specify something other than frequencies to
fill the table with, it would be what I need. I tried a few different
combinations of t() and unique() but none of it made enough sense to post as
my starting code... I am just lost. Any help would be greatly appreciated.

Thank you,
AHJ



--
View this message in context: 
http://r.789695.n4.nabble.com/Transform-pairwise-observations-into-a-table-tp4644706.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] Hmisc describe error


On Oct 1, 2012, at 9:33 AM, Bond, Stephen wrote:

 Describe fails for me with a message similar to what was an issue in 2008 and 
 got fixed according to posts.
 
 R version 2.15.0 (2012-03-30)
 Copyright (C) 2012 The R Foundation for Statistical Computing
 ISBN 3-900051-07-0
 Platform: i386-pc-mingw32/i386 (32-bit)
 
 # output truncated
 
 options(chmhelp = FALSE, help_type = text)
 .help.ESS - help
 options(STERM='iESS', editor='gnuclient.exe')
 load('prostate.sav')
 library(rms)
 Loading required package: Hmisc
 Loading required package: survival
 Loading required package: splines
 Hmisc library by Frank E Harrell Jr
 
 Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview')
 to see overall documentation.
 
 NOTE:Hmisc no longer redefines [.factor to drop unused levels when
 subsetting.  To get the old behavior of Hmisc type dropUnusedLevels().
 
 
 Attaching package: 'Hmisc'
 
 The following object(s) are masked from 'package:survival':
 
untangle.specials
 
 The following object(s) are masked from 'package:base':
 
format.pval, round.POSIXt, trunc.POSIXt, units
 
 
 Attaching package: 'rms'
 
 The following object(s) are masked from 'package:survival':
 
Surv
 
 describe(prostate)
 Error in format(dates(x)) : could not find function dates
 
 
Iget anerror when I use that code.
Error in describe(prostate) : object 'prostate' not found

And when I attempt loading it get:
data(prostate)

So you must have a different set of objects or packages loaded on your 
installation than I do.


David Winsemius, MD
Alameda, CA, USA

__
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] Transform pairwise observations into a table

2012-10-01 Thread Hadjixenofontos, Athena

Thank you. I had looked at xtabs but misunderstood the syntax. This is great. :)

AHJ



On Oct 1, 2012, at 12:53 PM, arun smartpink...@yahoo.com wrote:

 Hi,
 Try this:
 
 dat1-read.table(text= 
 ind1 ind2 coef 
 1 1 1 
 1 2 0.25 
 1 3 0.125 
 1 4 0.5 
 2 2 1 
 2 1 0.25 
 2 3 0.125 
 2 4 0.5 
 3 3 1 
 3 1 0.125 
 3 2 0.125 
 3 4 0.5 
 4 4 1 
 4 1 0.5 
 4 2 0.5 
 4 3 0.5 
 ,sep=,header=TRUE) 
 mat1-as.matrix(xtabs(coef~ind1+ind2,data=dat1)) 
 
 #ind2 
 #ind1 1 2 3 4 
   # 1 1.000 0.250 0.125 0.500 
#2 0.250 1.000 0.125 0.500 
#3 0.125 0.125 1.000 0.500 
#4 0.500 0.500 0.500 1.000 
 
 A.K.
 
 
 
 - Original Message -
 From: AHJ ahadjixenofon...@med.miami.edu
 To: r-help@r-project.org
 Cc: 
 Sent: Monday, October 1, 2012 12:17 PM
 Subject: [R] Transform pairwise observations into a table
 
 Hi,
 
 I have a table of pairs of individuals and a coefficient that belongs to the
 pair:
 
 ind1ind2coef
 111
 120.25
 130.125
 140.5
 221
 210.25
 230.125
 240.5
 331
 310.125
 320.125
 340.5
 441
 410.5
 420.5
 430.5
 
 And I want to convert it to a matrix where each individual is both a row and
 a column and at the intersection of each pair is the coefficient that
 belongs to that pair:
 
 1   2   3 4
 11   0.25   0.1250.5
 20.25  1   0.1250.5
 30.1250.12510.5
 40.5   0.5   0.51
 
 If table() would allow me to specify something other than frequencies to
 fill the table with, it would be what I need. I tried a few different
 combinations of t() and unique() but none of it made enough sense to post as
 my starting code... I am just lost. Any help would be greatly appreciated.
 
 Thank you,
 AHJ
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Transform-pairwise-observations-into-a-table-tp4644706.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.

Re: [R] Retrieve hypergeometric results in large scale

If you have not already done so, stop what you are doing and work
through the Introduction to R tutorial that ships with R (or other R
tutorial on the web that you may prefer).

The tutorials are written to help you climb the R learning curve much
more efficiently than the fooling around that you appear to be doing
now.

-- Bert

On Mon, Oct 1, 2012 at 8:31 AM, jas4710 wata...@post.com wrote:
 Hi Bert. This is not a homework. If I can do some basic programming in R like
 Perl, then I'll have a better chance to accomplish this task but the matrix
 concept is not quickly comprehensible...



 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Retrieve-95-coverage-of-results-from-a-hypergeometric-distribution-tp4644683p4644703.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.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

__
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] Hmisc describe error

2012-10-01 Thread stephenb

The error is just a misleading error message. loading the data produced
column sdate as 
 str(prostate)
$ sdate :Classes 'labelled', 'dates'  atomic [1:502] 2778 2820 2933 2999
3002 ...
  .. ..- attr(*, format)= symbol ddmmmyy
  .. ..- attr(*, label)= chr Date on study

prostate$sdate - as.Date(prostate$sdate,origin=1970-01-01)  # fixes the
problem and 'describe' runs

just as a comment: maybe the R dataset should have sdate as a Date class.



--
View this message in context: 
http://r.789695.n4.nabble.com/Hmisc-describe-error-tp4644708p4644711.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] Transform pairwise observations into a table

2012-10-01 Thread arun

Hi,
Try this:

dat1-read.table(text= 
ind1 ind2 coef 
1 1 1 
1 2 0.25 
1 3 0.125 
1 4 0.5 
2 2 1 
2 1 0.25 
2 3 0.125 
2 4 0.5 
3 3 1 
3 1 0.125 
3 2 0.125 
3 4 0.5 
4 4 1 
4 1 0.5 
4 2 0.5 
4 3 0.5 
,sep=,header=TRUE) 
mat1-as.matrix(xtabs(coef~ind1+ind2,data=dat1)) 

#    ind2 
#ind1     1     2     3     4 
  # 1 1.000 0.250 0.125 0.500 
   #2 0.250 1.000 0.125 0.500 
   #3 0.125 0.125 1.000 0.500 
   #4 0.500 0.500 0.500 1.000 

A.K.



- Original Message -
From: AHJ ahadjixenofon...@med.miami.edu
To: r-help@r-project.org
Cc: 
Sent: Monday, October 1, 2012 12:17 PM
Subject: [R] Transform pairwise observations into a table

Hi,

I have a table of pairs of individuals and a coefficient that belongs to the
pair:

ind1    ind2    coef
1    1    1
1    2    0.25
1    3    0.125
1    4    0.5
2    2    1
2    1    0.25
2    3    0.125
2    4    0.5
3    3    1
3    1    0.125
3    2    0.125
3    4    0.5
4    4    1
4    1    0.5
4    2    0.5
4    3    0.5

And I want to convert it to a matrix where each individual is both a row and
a column and at the intersection of each pair is the coefficient that
belongs to that pair:

    1           2           3             4
1    1           0.25           0.125    0.5
2    0.25      1           0.125    0.5
3    0.125    0.125    1            0.5
4    0.5           0.5           0.5            1

If table() would allow me to specify something other than frequencies to
fill the table with, it would be what I need. I tried a few different
combinations of t() and unique() but none of it made enough sense to post as
my starting code... I am just lost. Any help would be greatly appreciated.

Thank you,
AHJ



--
View this message in context: 
http://r.789695.n4.nabble.com/Transform-pairwise-observations-into-a-table-tp4644706.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.

Re: [R] Retrieve hypergeometric results in large scale

Thanks Jeff~~~

In fact I do not know how to combine and extract vectors in R.

ans-sort(dhyper(x, m, n, k),decreasing=TRUE)
rbind(ans,cumsum(ans)

will show the first point that exceeds 95% threshold. The problem is:
*information is lost*
I can no longer identify where are the first few elements from. e.g. for 10
numbers, maybe they are from 4,5,6,7 or for 100 numbers, from 45 to 68

So to append ID's to the data for later retrieval? rbind appears to do the
job but not so exactly...



--
View this message in context: 
http://r.789695.n4.nabble.com/Retrieve-95-coverage-of-results-from-a-hypergeometric-distribution-tp4644683p4644715.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] svyboxplot - library (survey)

2012-10-01 Thread Anthony Damico

try

svyboxplot(
dthage ~ factor( xspd2 ) ,
subset(nhis, mortstat==1) ,
col=gray80,
varwidth=TRUE ,
ylab=Age at Death ,
xlab=SPD Status: 1-SPD, 2=No SPD
)




On Mon, Oct 1, 2012 at 11:28 AM, Muhuri, Pradip (SAMHSA/CBHSQ) 
pradip.muh...@samhsa.hhs.gov wrote:

 Dear Anthony,

 ** **

 Yes, I can follow the example code you have given.  But, do you know from
 the code shown below (following Thomas Lumleys Complex Surveys) why I am
 getting the boxplot of dthage for just xspd=1, not xspd2=2?

 ** **

 My intent is the make this code work so that I can generate similar plots
 on other continuous variable.

 ** **

 Any help will be appreciated.

 ** **

 Thanks,

 ** **

 Pradip

 ** **

 ** **

 ** **

 

 nhis - svydesign (id=~psu, strat=~stratum, weights=~wt8,
 data=tor, nest=TRUE)

 svyboxplot (dthage~xspd2, subset (nhis, mortstat==1), col=gray80,
  varwidth=TRUE, ylab=Age at Death, xlab=SPD Status: 1-SPD,
 2=No SPD)


 

 ** **

 Pradip K. Muhuri, PhD

 Statistician

 Substance Abuse  Mental Health Services Administration

 The Center for Behavioral Health Statistics and Quality

 Division of Population Surveys

 1 Choke Cherry Road, Room 2-1071

 Rockville, MD 20857

  

 Tel: 240-276-1070

 Fax: 240-276-1260

 e-mail: pradip.muh...@samhsa.hhs.gov

  

 The Center for Behavioral Health Statistics and Quality your feedback.
 Please click on the following link to complete a brief customer survey:
 http://cbhsqsurvey.samhsa.gov

 ** **

 *From:* Anthony Damico [mailto:ajdam...@gmail.com]
 *Sent:* Monday, October 01, 2012 10:07 AM
 *To:* Muhuri, Pradip (SAMHSA/CBHSQ)
 *Cc:* R help
 *Subject:* Re: [R] svyboxplot - library (survey)

 ** **

 using a slight modification of the example shown in ?svyboxplot


 # load survey library
 library(survey)

 # load example data
 data(api)

 # create an example svydesign
 dstrat - svydesign(id = ~1, strata = ~stype, weights = ~pw, data =
 apistrat,
 fpc = ~fpc)

 # set the plot window to display 1 plot x 2 plots
 par(mfrow=c(1,2))

 # generate two example boxplots
 svyboxplot(enroll~stype,dstrat,all.outliers=TRUE)
 svyboxplot(enroll~1,dstrat)

 # done



 # alternative: not as nice

 # set the plot window to display 2 plots x 1 plot
 par(mfrow=c(2,1))

 # generate two example boxplots
 svyboxplot(enroll~stype,dstrat,all.outliers=TRUE)
 svyboxplot(enroll~1,dstrat)

 # done






 

 On Mon, Oct 1, 2012 at 9:50 AM, Muhuri, Pradip (SAMHSA/CBHSQ) 
 pradip.muh...@samhsa.hhs.gov wrote:

 Hello,

 I have used the library (survey) package for boxplots using the following
 code.

 Could anyone please tell me why I am getting only 1  boxplot instead of 2
 boxplots (1-SPD,  2-No SPD).

 What changes in the following code would be required to get 2 boxplots in
 the same plot frame?

 Thanks,

 Pradip

 ###
 nhis - svydesign (id=~psu, strat=~stratum, weights=~wt8,
 data=tor, nest=TRUE)

 svyboxplot (dthage~xspd2, subset (nhis, mortstat==1), col=gray80,
  varwidth=TRUE, ylab=Age at Death, xlab=SPD Status: 1-SPD,
 2=No SPD)


 Pradip K. Muhuri
 Statistician
 Substance Abuse  Mental Health Services Administration
 The Center for Behavioral Health Statistics and Quality
 Division of Population Surveys
 1 Choke Cherry Road, Room 2-1071
 Rockville, MD 20857

 Tel: 240-276-1070
 Fax: 240-276-1260
 e-mail: pradip.muh...@samhsa.hhs.gov

 The Center for Behavioral Health Statistics and Quality your feedback.
 Please click on the following link to complete a brief customer survey:
 http://cbhsqsurvey.samhsa.gov

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

 ** **


[[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] Transform pairwise observations into a table

You might be amused by this alternative:

with(dat1,matrix(coef[order(ind2,ind1)],nr=4,nc=4))

It works because matrices are vectors stored in column major order.

Cheers,
Bert


On Mon, Oct 1, 2012 at 9:53 AM, arun smartpink...@yahoo.com wrote:
 Hi,
 Try this:

 dat1-read.table(text=
 ind1 ind2 coef
 1 1 1
 1 2 0.25
 1 3 0.125
 1 4 0.5
 2 2 1
 2 1 0.25
 2 3 0.125
 2 4 0.5
 3 3 1
 3 1 0.125
 3 2 0.125
 3 4 0.5
 4 4 1
 4 1 0.5
 4 2 0.5
 4 3 0.5
 ,sep=,header=TRUE)
 mat1-as.matrix(xtabs(coef~ind1+ind2,data=dat1))

 #ind2
 #ind1 1 2 3 4
   # 1 1.000 0.250 0.125 0.500
#2 0.250 1.000 0.125 0.500
#3 0.125 0.125 1.000 0.500
#4 0.500 0.500 0.500 1.000

 A.K.



 - Original Message -
 From: AHJ ahadjixenofon...@med.miami.edu
 To: r-help@r-project.org
 Cc:
 Sent: Monday, October 1, 2012 12:17 PM
 Subject: [R] Transform pairwise observations into a table

 Hi,

 I have a table of pairs of individuals and a coefficient that belongs to the
 pair:

 ind1ind2coef
 111
 120.25
 130.125
 140.5
 221
 210.25
 230.125
 240.5
 331
 310.125
 320.125
 340.5
 441
 410.5
 420.5
 430.5

 And I want to convert it to a matrix where each individual is both a row and
 a column and at the intersection of each pair is the coefficient that
 belongs to that pair:

 1   2   3 4
 11   0.25   0.1250.5
 20.25  1   0.1250.5
 30.1250.12510.5
 40.5   0.5   0.51

 If table() would allow me to specify something other than frequencies to
 fill the table with, it would be what I need. I tried a few different
 combinations of t() and unique() but none of it made enough sense to post as
 my starting code... I am just lost. Any help would be greatly appreciated.

 Thank you,
 AHJ



 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Transform-pairwise-observations-into-a-table-tp4644706.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.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

__
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] Hmisc describe error


On Oct 1, 2012, at 9:48 AM, stephenb wrote:

 The error is just a misleading error message. loading the data produced
 column sdate as 
 str(prostate)
 $ sdate :Classes 'labelled', 'dates'  atomic [1:502] 2778 2820 2933 2999
 3002 ...
  .. ..- attr(*, format)= symbol ddmmmyy
  .. ..- attr(*, label)= chr Date on study
 
 prostate$sdate - as.Date(prostate$sdate,origin=1970-01-01)  # fixes the
 problem and 'describe' runs
 
 just as a comment: maybe the R dataset should have sdate as a Date class.
 

What R dataset? From what package?

And it's already a dates-classed variable. 

-- 
David Winsemius, MD
Alameda, CA, USA

__
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] Hmisc describe error

2012-10-01 Thread Bond, Stephen

Just as a clarification: I downloaded 'prostate.sav' from F. Harrell's website. 
For some reason

 data(prostate)
Warning message:
In data(prostate) : data set 'prostate' not found

I don't have the prostate data set as is.

-Original Message-
From: David Winsemius [mailto:dwinsem...@comcast.net] 
Sent: Monday, October 01, 2012 12:53 PM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] Hmisc describe error


On Oct 1, 2012, at 9:33 AM, Bond, Stephen wrote:

 Describe fails for me with a message similar to what was an issue in 2008 and 
 got fixed according to posts.
 
 R version 2.15.0 (2012-03-30)
 Copyright (C) 2012 The R Foundation for Statistical Computing
 ISBN 3-900051-07-0
 Platform: i386-pc-mingw32/i386 (32-bit)
 
 # output truncated
 
 options(chmhelp = FALSE, help_type = text)
 .help.ESS - help
 options(STERM='iESS', editor='gnuclient.exe')
 load('prostate.sav')
 library(rms)
 Loading required package: Hmisc
 Loading required package: survival
 Loading required package: splines
 Hmisc library by Frank E Harrell Jr
 
 Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview')
 to see overall documentation.
 
 NOTE:Hmisc no longer redefines [.factor to drop unused levels when
 subsetting.  To get the old behavior of Hmisc type dropUnusedLevels().
 
 
 Attaching package: 'Hmisc'
 
 The following object(s) are masked from 'package:survival':
 
untangle.specials
 
 The following object(s) are masked from 'package:base':
 
format.pval, round.POSIXt, trunc.POSIXt, units
 
 
 Attaching package: 'rms'
 
 The following object(s) are masked from 'package:survival':
 
Surv
 
 describe(prostate)
 Error in format(dates(x)) : could not find function dates
 
 
Iget anerror when I use that code.
Error in describe(prostate) : object 'prostate' not found

And when I attempt loading it get:
data(prostate)

So you must have a different set of objects or packages loaded on your 
installation than I do.


David Winsemius, MD
Alameda, CA, USA

__
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] Hmisc describe error


On Oct 1, 2012, at 10:27 AM, David Winsemius wrote:

 
 On Oct 1, 2012, at 9:48 AM, stephenb wrote:
 
 The error is just a misleading error message. loading the data produced
 column sdate as 
 str(prostate)
 $ sdate :Classes 'labelled', 'dates'  atomic [1:502] 2778 2820 2933 2999
 3002 ...
 .. ..- attr(*, format)= symbol ddmmmyy
 .. ..- attr(*, label)= chr Date on study
 
 prostate$sdate - as.Date(prostate$sdate,origin=1970-01-01)  # fixes the
 problem and 'describe' runs
 
 just as a comment: maybe the R dataset should have sdate as a Date class.
 
 
 What R dataset? From what package?
 
 And it's already a dates-classed variable. 

Not only that, but there is a message at the top of Harrell's Datasets page 
that says:

For R users of the prostate dataset, put library(chron) into effect to handle 
date variables. A simpler approach is to just convert the one date variable to 
the built-in R format by running the command prostate$sdate - 
as.Date(prostate$sdate).

-- 

David Winsemius, MD
Alameda, CA, USA

__
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] Transform pairwise observations into a table

2012-10-01 Thread Rui Barradas


Hello,

Try the following.


dat - read.table(text=
ind1 ind2 coef
1 1 1
1 2 0.25
1 3 0.125
1 4 0.5
2 2 1
2 1 0.25
2 3 0.125
2 4 0.5
3 3 1
3 1 0.125
3 2 0.125
3 4 0.5
4 4 1
4 1 0.5
4 2 0.5
4 3 0.5
, header=TRUE)
dat

reshape(dat, v.names = coef, idvar = ind1, timevar = ind2,
direction = wide)

Hope this helps,

Rui Barradas
Em 01-10-2012 17:17, AHJ escreveu:

Hi,

I have a table of pairs of individuals and a coefficient that belongs to the
pair:

ind1ind2coef
1   1   1
1   2   0.25
1   3   0.125
1   4   0.5
2   2   1
2   1   0.25
2   3   0.125
2   4   0.5
3   3   1
3   1   0.125
3   2   0.125
3   4   0.5
4   4   1
4   1   0.5
4   2   0.5
4   3   0.5

And I want to convert it to a matrix where each individual is both a row and
a column and at the intersection of each pair is the coefficient that
belongs to that pair:

1  2   3 4
1   1  0.250.1250.5
2   0.25  10.1250.5
3   0.1250.1251 0.5
4   0.50.5 0.5  1

If table() would allow me to specify something other than frequencies to
fill the table with, it would be what I need. I tried a few different
combinations of t() and unique() but none of it made enough sense to post as
my starting code... I am just lost. Any help would be greatly appreciated.

Thank you,
AHJ



--
View this message in context: 
http://r.789695.n4.nabble.com/Transform-pairwise-observations-into-a-table-tp4644706.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.

Re: [R] Retrieve hypergeometric results in large scale

2012-10-01 Thread Rui Barradas

Hello,

See the differences.

k - 3
p - 0.95

m - 90; n - 10
dhyper(0:k, m, n, k) # Prob(X = x), with x = 0:k
phyper(0:k, m, n, k) # Prob(X = x)
# quantiles, what you want
qhyper(p, m, n, k) # inverse of phyper

m - 50; n - 50
dhyper(0:k, m, n, k)
phyper(0:k, m, n, k)
qhyper(p, m, n, k)

In your original post (op) you gave the output of dhyper. I would find
phyper much better, it shows that 0.95 is between the last two values of
phyper. Apparently you want its inverse, qhyper, but are assuming there
are such things as decimals. The hypergeometric is a discrete
distribution, so, from the help page,

The quantile is defined as the smallest value /x/ such that /F(x) ? p/,
where /F/ is the distribution function.

And, R is surely more competent at distribution functions than Perl.
Stick to it and read An Introduction to R, file R-intro.pdf in the doc
directory of your R installation, Chapter 8. All computer languages have
some learning time and as far as statistics is concerned, learning R pays.

Hope this helps,

Rui Barradas
Em 01-10-2012 16:28, jas4710 escreveu:
Thanks Jeff
The documentation pages, if I haven't missed any crucial points, illustrate
how to get probability and cumulative probability values.

I can first retrieve the data structures and use Perl (I don't know how to
use R...) to sort the derived ratios and sum the probability values until
the cumulative probability exceeds 95%. Well, I just don't know whether such
seemingly routine procedures have already been implemented...

Thanks again!

--
View this message in context:
http://r.789695.n4.nabble.com/Retrieve-95-coverage-of-results-from-a-hypergeometric-distribution-tp4644683p4644701.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.

[[alternative HTML version deleted]]

[R] ffbase, help with %in%

2012-10-01 Thread Lucas Chaparro

Hello to everyone.
I'm trying to use the %in% to match to vectors in ff format.
a-as.ff(data[,1]) %in% fire$fecha

 aff (open) logical length=3653 (3653)
   [1][2][3][4][5][6][7][8][3646]
 FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  :  FALSE
[3647] [3648] [3649] [3650] [3651] [3652] [3653]
 FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE


Here you see a part of the data:

data[1:20,]  (just a sample, data has 3653 obs)

fecha juliano altura  UTM.E   UTM.N
1  1990-07-01 182 15 248500 6239500
2  1990-07-02 183 15 248500 6239500
3  1990-07-03 184 15 248500 6239500
4  1990-07-04 185 15 248500 6239500
5  1990-07-05 186 15 248500 6239500
6  1990-07-06 187 15 248500 6239500
7  1990-07-07 188 15 248500 6239500
8  1990-07-08 189 15 248500 6239500
9  1990-07-09 190 15 248500 6239500
10 1990-07-10 191 15 248500 6239500
11 1990-07-11 192 15 248500 6239500
12 1990-07-12 193 15 248500 6239500
13 1990-07-13 194 15 248500 6239500
14 1990-07-14 195 15 248500 6239500
15 1990-07-15 196 15 248500 6239500
16 1990-07-16 197 15 248500 6239500
17 1990-07-17 198 15 248500 6239500
18 1990-07-18 199 15 248500 6239500
19 1990-07-19 200 15 248500 6239500
20 1990-07-20 201 15 248500 6239500


 fire$fecha[1:20,] [1] 1984-11-08 1984-11-08 1984-11-09 1984-11-09 
 1984-11-09
 [6] 1984-11-10 1984-11-10 1984-11-11 1984-11-11 1984-11-11
[11] 1984-11-11 1984-11-11 1984-11-11 1984-11-12 1984-11-12
[16] 1984-11-13 1984-11-13 1984-11-13 1984-11-14 1984-11-14


to see if a got any match:


table.ff(a)


FALSE  TRUE
 1687  1966 Mensajes de aviso perdidosIn if (useNA == no) c(NA, NaN) :
  la condición tiene longitud  1 y sólo el primer elemento será usado


in a regular data.frame I use data[a,] to extract the rows that a ==
TRUE, but when i do this in a ffdf i get this error:


 data[a,]Error: vmode(index) == integer is not TRUE


I'm just learning how to use the ff package so, obviously I'm missing something


If any of you knows how to solve this, please teach me.

Thank you so much.


Lucas.

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

2012-10-01 Thread Rui Barradas


Hello,

Or maybe %in%.

x - Sys.Date() + 1:10
y - Sys.Date() + 7:15

x %in% y # logical index into 'x'
x[x %in% y] # common
x[!x %in% y] # not common

setdiff(x, y) # doesn't keep class Date


Hope this helps,

Rui Barradas

Em 01-10-2012 15:57, R. Michael Weylandt escreveu:

On Monday, October 1, 2012, Karan Anand wrote:


hi,
 I am new to using R ,I have 2 datasets with dates common in them ,how
can i take out the common dates within them.



Depending on what you mean by 'out' either merge() or setdiff()

RMW




karan

 [[alternative HTML version deleted]]

__
R-help@r-project.org javascript:; mailing list
https://stat.ethz.ch/mailman/listinfo/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-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] modified Ward method

2012-10-01 Thread David L Carlson

You asked this question before. I think you need to tell us what modified
ward method you are talking about and what you are trying to accomplish. I
did a quick Google Scholar search and turned up only a few papers using that
phrase, but they were not defining it in the same way. 

--
David L Carlson
Associate Professor of Anthropology
Texas AM University
College Station, TX 77843-4352

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of eliza botto
 Sent: Saturday, September 29, 2012 9:29 PM
 To: r-help@r-project.org
 Subject: [R] modified Ward method
 
 
 
 
 
 Dear useRs,Is there a command or package in R to do heirarchical
 clustering by modified ward method??regardseliza
 
   [[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] calculating probability from the density function

2012-10-01 Thread Greg Snow

You may find it easier to use the logspline density fits (logspline
package) rather than the kernel density estimators for this.

On Mon, Oct 1, 2012 at 7:46 AM, Eugene Kanshin kanshin...@gmail.com wrote:
 Hello,
 I have a data x with normal (or very close to normal) distribution, I can
 plot a density distribution with density(x,...). My question is is there
 any way to calculate an area under this distribution (=probability) for
 particular range of x values, let's say for x from 0 to 2? I was not able
 to find any kind of simple procedure to do this.
 Thanks in advance for your help,
 Evgeny.

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



-- 
Gregory (Greg) L. Snow Ph.D.
538...@gmail.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] (senza oggetto)

2012-10-01 Thread MYRIAM TABASSO

hi, can I help me to cancel my name in this list?
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.

Re: [R] (senza oggetto)

2012-10-01 Thread peter dalgaard


On Oct 1, 2012, at 20:28 , MYRIAM TABASSO wrote:

 hi, can I help me to cancel my name in this list?

Sure, just go to the mailman site listed in the footer (look near bottom of 
page).

 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.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.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] ffbase, help with %in%

2012-10-01 Thread Christiaan Pauw

Hi Lucas

I don't know the ff package very well but here is what I found. Maybe
there is a clue in here

 z   - as.Date(1970-01-01)+1:10
 zff - as.ff(z)
 z %in% zff

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

#but
 z %in% zff[1:length(zff),]
 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# and also
 z %in% zff[1:10]
 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

That seems to be because an ff object is a list

 str(zff)
 list()
 - attr(*, physical)=Class 'ff_pointer' externalptr
  ..- attr(*, vmode)= chr double
  ..- attr(*, maxlength)= int 10
  ..- attr(*, pattern)= chr clone
  ..- attr(*, filename)= chr
/private/var/folders/ub/ubvWLUkKHf8WAywv5rmtcE+++TI/-Tmp-/RtmpmMKHRx/clone406828b135cc.ff
  ..- attr(*, pagesize)= int 65536
  ..- attr(*, finalizer)= chr close
  ..- attr(*, finonexit)= logi TRUE
  ..- attr(*, readonly)= logi FALSE
  ..- attr(*, caching)= chr mmnoflush
 - attr(*, virtual)= list()
  ..- attr(*, Length)= int 10
  ..- attr(*, Symmetric)= logi FALSE
  ..- attr(*, ramclass)= chr Date
 - attr(*, class) =  chr [1:2] ff_vector ff

If you make a data frame e.g
df - data.frame(date=z,letter=letters[1:10])
# this doesn't work
 as.ff(df[1:10,1])%in% z
logical(0)
# but this does
 as.ff(df[1:10,1]%in% z)
ff (open) logical length=10 (10)
 [1]  [2]  [3]  [4]  [5]  [6]  [7]  [8]  [9] [10]
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

The question is: In your example, do you want a to be a ff object or
the first column of data to be an ff object (I could get this to work
- propably for the reasons above :
 df - data.frame(date=as.ff(z),letter=letters[1:10])
Error in as.data.frame.default(x[[i]], optional = TRUE,
stringsAsFactors = stringsAsFactors) :
  cannot coerce class 'c(ff_vector, ff)' into a data.frame

On 1 October 2012 20:07, Lucas Chaparro lpchaparro...@gmail.com wrote:
 Hello to everyone.
 I'm trying to use the %in% to match to vectors in ff format.
 a-as.ff(data[,1]) %in% fire$fecha

 aff (open) logical length=3653 (3653)
[1][2][3][4][5][6][7][8][3646]
  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  :  FALSE
 [3647] [3648] [3649] [3650] [3651] [3652] [3653]
  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE


 Here you see a part of the data:

 data[1:20,]  (just a sample, data has 3653 obs)

 fecha juliano altura  UTM.E   UTM.N
 1  1990-07-01 182 15 248500 6239500
 2  1990-07-02 183 15 248500 6239500
 3  1990-07-03 184 15 248500 6239500
 4  1990-07-04 185 15 248500 6239500
 5  1990-07-05 186 15 248500 6239500
 6  1990-07-06 187 15 248500 6239500
 7  1990-07-07 188 15 248500 6239500
 8  1990-07-08 189 15 248500 6239500
 9  1990-07-09 190 15 248500 6239500
 10 1990-07-10 191 15 248500 6239500
 11 1990-07-11 192 15 248500 6239500
 12 1990-07-12 193 15 248500 6239500
 13 1990-07-13 194 15 248500 6239500
 14 1990-07-14 195 15 248500 6239500
 15 1990-07-15 196 15 248500 6239500
 16 1990-07-16 197 15 248500 6239500
 17 1990-07-17 198 15 248500 6239500
 18 1990-07-18 199 15 248500 6239500
 19 1990-07-19 200 15 248500 6239500
 20 1990-07-20 201 15 248500 6239500


 fire$fecha[1:20,] [1] 1984-11-08 1984-11-08 1984-11-09 1984-11-09 
 1984-11-09
  [6] 1984-11-10 1984-11-10 1984-11-11 1984-11-11 1984-11-11
 [11] 1984-11-11 1984-11-11 1984-11-11 1984-11-12 1984-11-12
 [16] 1984-11-13 1984-11-13 1984-11-13 1984-11-14 1984-11-14


 to see if a got any match:


 table.ff(a)


 FALSE  TRUE
  1687  1966 Mensajes de aviso perdidosIn if (useNA == no) c(NA, NaN) :
   la condición tiene longitud  1 y sólo el primer elemento será usado


 in a regular data.frame I use data[a,] to extract the rows that a ==
 TRUE, but when i do this in a ffdf i get this error:


 data[a,]Error: vmode(index) == integer is not TRUE


 I'm just learning how to use the ff package so, obviously I'm missing 
 something


 If any of you knows how to solve this, please teach me.

 Thank you so much.


 Lucas.

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




-- 
Christiaan Pauw
Nova Institute
www.nova.org.za

__
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] Problem with nls regression fit

2012-10-01 Thread Jean V Adams

You have not specified a nonlinear formula.  There are no parameters to 
estimate in the formula you provide, y1~dist.  What is the nonlinear 
relation you are trying to fit?  Look at the help file for nls to see some 
examples worked.

?nls

Jean


Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote on 10/01/2012 
10:27:23 AM:
 
 Hi all,
 I got following problem in fitting the data.
 Any kind of suggestions are welcome
 
 
  beta - 3.5
  d - seq(0.1,62.5,0.1)
  y - exp(-beta*d)
  y1 - y
  x - read.table(epidist.txt, header = TRUE)
  data.nls - as.data.frame(cbind(y1,x))
  #attach(data.nls)
  nls.fit - nls(y1~dist,data.nls)
 Error in cll[[1L]] : object of type 'symbol' is not subsettable
 
 
 Best

[[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] clustering spline-based models

2012-10-01 Thread Wyatt McMahon

Hello playeRs!

 

I'm working on a project for a client.  She's modeling hormone levels
periodically, and trying to develop a model and fit her data to that
model, and subsequently she's trying to cluster individuals based on how
well each fits the model.

 

I've been looking at grofit for this, but it doesn't appear that it can do
the sort of post-hoc analysis I'm looking for, although it can use a
model-free spline fit nicely (reportedly).  

 

Does anyone have any packages that can accomplish what I'm looking to
accomplish?

 

Thanks in advance for any help you can offer,

 

Wyatt


[[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] Retrieve hypergeometric results in large scale

2012-10-01 Thread William Dunlap

order() is usually a lot more useful than sort(), since, as you noticed,
sort() drops information about where each element in its output came
from.

Your example was incomplete so I made up one which I
think is similar.
   n - 10 ; p - 0.7 ; k - 0:n ; d - dbinom(k, n, p)
   plot(k, d) # density of binomial over its domain
If you want the indices of the largest density values whose
cumulative sum is less than 0.95 you
   ord - order(d, decreasing=TRUE) # indices such that d[ord] is in 
decreasing order
   big - ord[cumsum(d[ord])  0.95]
   data.frame(big, d=d[big], cumsum=cumsum(d[big]))
big dcumsum
  1   8 0.2668279 0.2668279
  2   9 0.2334744 0.5003024
  3   7 0.2001209 0.7004233
  4  10 0.1210608 0.8214841
  5   6 0.1029193 0.9244035
  points(cex=2, k[big], d[big])

If you want to include the index of the density value that puts
you just over 0.95 first find the complement of the desired indices
and use setdiff to compute its complement.  E.g.,
   ord - order(d)
   little - ord[cumsum(d[ord])  0.05]
   big - setdiff(seq_along(d), little) # difference of two sets of numbers
   big
  [1]  5  6  7  8  9 10

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
 Behalf
 Of jas4710
 Sent: Monday, October 01, 2012 9:59 AM
 To: r-help@r-project.org
 Subject: Re: [R] Retrieve hypergeometric results in large scale
 
 Thanks Jeff~~~
 
 In fact I do not know how to combine and extract vectors in R.
 
 ans-sort(dhyper(x, m, n, k),decreasing=TRUE)
 rbind(ans,cumsum(ans)
 
 will show the first point that exceeds 95% threshold. The problem is:
 *information is lost*
 I can no longer identify where are the first few elements from. e.g. for 10
 numbers, maybe they are from 4,5,6,7 or for 100 numbers, from 45 to 68
 
 So to append ID's to the data for later retrieval? rbind appears to do the
 job but not so exactly...
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Retrieve-95-coverage-of-
 results-from-a-hypergeometric-distribution-tp4644683p4644715.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.

Re: [R] svyboxplot - library (survey)

2012-10-01 Thread Thomas Lumley

The documentation says The grouping variable in svyboxplot, if
present, must be a factor

  -thomas

On Tue, Oct 2, 2012 at 4:28 AM, Muhuri, Pradip (SAMHSA/CBHSQ)
pradip.muh...@samhsa.hhs.gov wrote:
 Dear Anthony,

 Yes, I can follow the example code you have given.  But, do you know from the 
 code shown below (following Thomas Lumley's Complex Surveys) why I am 
 getting the boxplot of dthage for just xspd=1, not xspd2=2?

 My intent is the make this code work so that I can generate similar plots on 
 other continuous variable.

 Any help will be appreciated.

 Thanks,

 Pradip



 
 nhis - svydesign (id=~psu, strat=~stratum, weights=~wt8,
 data=tor, nest=TRUE)

 svyboxplot (dthage~xspd2, subset (nhis, mortstat==1), col=gray80,
  varwidth=TRUE, ylab=Age at Death, xlab=SPD Status: 1-SPD, 
 2=No SPD)



 Pradip K. Muhuri, PhD
 Statistician
 Substance Abuse  Mental Health Services Administration
 The Center for Behavioral Health Statistics and Quality
 Division of Population Surveys
 1 Choke Cherry Road, Room 2-1071
 Rockville, MD 20857

 Tel: 240-276-1070
 Fax: 240-276-1260
 e-mail: pradip.muh...@samhsa.hhs.govmailto:pradip.muh...@samhsa.hhs.gov

 The Center for Behavioral Health Statistics and Quality your feedback.  
 Please click on the following link to complete a brief customer survey:   
 http://cbhsqsurvey.samhsa.govhttp://cbhsqsurvey.samhsa.gov/

 From: Anthony Damico [mailto:ajdam...@gmail.com]
 Sent: Monday, October 01, 2012 10:07 AM
 To: Muhuri, Pradip (SAMHSA/CBHSQ)
 Cc: R help
 Subject: Re: [R] svyboxplot - library (survey)

 using a slight modification of the example shown in ?svyboxplot


 # load survey library
 library(survey)

 # load example data
 data(api)

 # create an example svydesign
 dstrat - svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat,
 fpc = ~fpc)

 # set the plot window to display 1 plot x 2 plots
 par(mfrow=c(1,2))

 # generate two example boxplots
 svyboxplot(enroll~stype,dstrat,all.outliers=TRUE)
 svyboxplot(enroll~1,dstrat)

 # done



 # alternative: not as nice

 # set the plot window to display 2 plots x 1 plot
 par(mfrow=c(2,1))

 # generate two example boxplots
 svyboxplot(enroll~stype,dstrat,all.outliers=TRUE)
 svyboxplot(enroll~1,dstrat)

 # done






 On Mon, Oct 1, 2012 at 9:50 AM, Muhuri, Pradip (SAMHSA/CBHSQ) 
 pradip.muh...@samhsa.hhs.govmailto:pradip.muh...@samhsa.hhs.gov wrote:
 Hello,

 I have used the library (survey) package for boxplots using the following 
 code.

 Could anyone please tell me why I am getting only 1  boxplot instead of 2 
 boxplots (1-SPD,  2-No SPD).

 What changes in the following code would be required to get 2 boxplots in the 
 same plot frame?

 Thanks,

 Pradip

 ###
 nhis - svydesign (id=~psu, strat=~stratum, weights=~wt8,
 data=tor, nest=TRUE)

 svyboxplot (dthage~xspd2, subset (nhis, mortstat==1), col=gray80,
  varwidth=TRUE, ylab=Age at Death, xlab=SPD Status: 1-SPD, 
 2=No SPD)


 Pradip K. Muhuri
 Statistician
 Substance Abuse  Mental Health Services Administration
 The Center for Behavioral Health Statistics and Quality
 Division of Population Surveys
 1 Choke Cherry Road, Room 2-1071
 Rockville, MD 20857

 Tel: 240-276-1070
 Fax: 240-276-1260
 e-mail: pradip.muh...@samhsa.hhs.govmailto:pradip.muh...@samhsa.hhs.gov

 The Center for Behavioral Health Statistics and Quality your feedback.  
 Please click on the following link to complete a brief customer survey:   
 http://cbhsqsurvey.samhsa.gov

 vide commented, minimal, self-contained, reproducible code.
 __
 R-help@r-project.orgmailto: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.



-- 
Thomas Lumley
Professor of Biostatistics
University of Auckland

__
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] Error messages when attempting to calculate polychoric correlation matrices

2012-10-01 Thread John Fox

Dear Kevin,

From ?polychor (in the polycor package):

polychor(x, y, ML = FALSE, control = list(), std.err = FALSE, maxcor=.)

Arguments

x   
a contingency table of counts or an ordered categorical variable; the latter
can be numeric, logical, a factor, or an ordered factor, but if a factor,
its levels should be in proper order.

y   
if x is a variable, a second ordered categorical variable.

So polychor doesn't take a data frame as its argument, and computes only one
polychoric correlation at a time. (It generally helps in R, and even more
generally, to read the documentation.)

Thus, e.g.,

 with(data.matrix18, polychor(SABAS02, SABAS06))
[1] 0.1811532

 with(data.matrix18, polychor(SABAS02, SABAS06, ML=TRUE))
[1] 0.1817071

If you want to compute a matrix of polychoric correlations using the polycor
package, then you can use the hetcor function -- see ?hetcor -- but you'll
then have to change your variables from numeric to factors or ordered
factors, or hetcor will compute Pearson pairwise correlations among them.
You could use, for example

for (j in 1:ncol(data.matrix18)) data.matrix18[, j] -
as.factor(data.matrix18[, j])

I hope this helps,
 John

---
John Fox
Senator McMaster Professor of Social Statistics
Department of Sociology
McMaster University
Hamilton, Ontario, Canada




 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
 On Behalf Of Kevin Cheung
 Sent: Monday, October 01, 2012 12:19 PM
 To: 'r-help@r-project.org'
 Subject: [R] Error messages when attempting to calculate polychoric
 correlation matrices
 
 Dear R users,
 
 I am a psychology postgraduate student who is relatively new to using R.
 I am currently developing a psychometric scale and have run into a few
 problems when using R to calculate a polychoric correlation matrix for
 my dataset. I am trying to produce a polychoric correlation matrix for
 calculating ordinal reliability estimates (eg. Alpha, omega).The set
 consists of 439 observations of 18 variables. Each of the 18 variables
 was measured using a 6 point Likert scale but some variables only
 feature responses across 5 ordinal steps because there were no observed
 responses for one extreme of the item. The dataset has no missing data
 following multiple imputation and preliminary analyses were conducted
 using SPSS. I have included a command that will reproduce the dataset
 below, apologies for the size of the command:
 
 data.matrix18 -
 structure(list(SABAS02 = structure(c(6, 2, 5, 6, 4, 6, 5, 3,
 5, 1, 4, 4, 6, 5, 5, 5, 4, 4, 4, 3, 4, 6, 5, 2, 5, 4, 5, 4, 4,
 5, 5, 5, 2, 5, 4, 5, 5, 5, 6, 3, 6, 4, 5, 5, 5, 5, 4, 4, 5, 5,
 5, 4, 5, 5, 4, 6, 6, 5, 6, 4, 4, 5, 4, 5, 5, 5, 6, 4, 4, 5, 4,
 6, 5, 6, 5, 6, 6, 5, 4, 6, 6, 5, 3, 4, 2, 4, 5, 5, 6, 2, 6, 5,
 4, 6, 2, 5, 3, 5, 4, 4, 5, 4, 4, 3, 4, 4, 6, 5, 6, 6, 4, 4, 5,
 5, 5, 6, 5, 5, 5, 5, 4, 5, 5, 5, 4, 4, 5, 6, 4, 5, 5, 4, 5, 5,
 4, 5, 5, 4, 5, 5, 5, 4, 2, 5, 4, 5, 5, 5, 5, 5, 5, 4, 5, 4, 4,
 5, 4, 6, 2, 4, 6, 4, 4, 5, 4, 4, 6, 6, 5, 3, 3, 2, 4, 5, 5, 5,
 5, 6, 2, 5, 5, 4, 6, 4, 5, 6, 5, 4, 5, 5, 4, 4, 5, 4, 5, 5, 3,
 5, 4, 3, 5, 4, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 5, 4, 4, 5, 5, 6,
 4, 3, 5, 6, 5, 4, 5, 5, 2, 5, 5, 4, 2, 4, 5, 5, 5, 5, 5, 3, 4,
 5, 2, 4, 3, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 5, 4, 5, 6, 4, 5, 3,
 5, 6, 5, 4, 5, 5, 5, 5, 5, 4, 4, 5, 6, 5, 4, 5, 5, 6, 6, 6, 4,
 3, 4, 4, 5, 6, 4, 5, 5, 4, 6, 5, 5, 4, 5, 4, 4, 4, 4, 5, 3, 5,
 5, 4, 3, 4, 5, 4, 4, 5, 4, 5, 5, 3, 5, 6, 5, 5, 5, 5, 5, 4, 3,
 5, 4, 4, 5, 5, 6, 4, 5, 4, 2, 6, 5, 5, 5, 5, 5, 5, 5, 6, 6, 5,
 5, 5, 4, 5, 5, 6, 5, 4, 4, 5, 4, 4, 5, 5, 4, 6, 4, 4, 5, 6, 5,
 5, 4, 5, 4, 5, 4, 6, 5, 5, 5, 4, 5, 6, 5, 5, 5, 5, 5, 4, 3, 4,
 5, 4, 4, 6, 4, 4, 4, 2, 5, 4, 4, 6, 2, 6, 4, 5, 5, 6, 5, 6, 4,
 4, 3, 4, 5, 5, 6, 5, 5, 5, 5, 4, 5, 3, 5, 5, 5, 5, 6, 3, 5, 6,
 5, 6, 4, 5, 4, 4, 5, 6, 2, 5, 6), value.labels = structure(c(6,
 5, 4, 3, 2, 1), .Names = c(Strongly Agree, Agree, Slightly Agree,
 Slightly Disagree, Disagree, Strongly Disagree))), SABAS06 =
 structure(c(6,
 6, 6, 6, 6, 6, 5, 5, 6, 3, 6, 6, 6, 5, 5, 6, 5, 6, 6, 5, 5, 6,
 6, 6, 6, 5, 6, 3, 6, 5, 5, 5, 5, 6, 6, 6, 3, 6, 6, 5, 5, 4, 5,
 6, 6, 5, 6, 5, 6, 6, 6, 6, 6, 6, 5, 5, 6, 6, 6, 6, 6, 6, 5, 6,
 6, 6, 6, 6, 4, 6, 5, 5, 6, 6, 6, 6, 6, 6, 3, 5, 6, 6, 6, 6, 6,
 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 6, 6, 6, 5, 5, 5, 6, 5, 6,
 5, 6, 6, 6, 6, 5, 5, 6, 5, 5, 6, 2, 5, 5, 6, 5, 5, 5, 6, 5, 6,
 5, 5, 6, 6, 6, 5, 5, 5, 6, 6, 5, 6, 6, 6, 4, 6, 6, 5, 4, 6, 5,
 5, 6, 5, 5, 5, 6, 6, 6, 6, 5, 5, 6, 6, 5, 6, 4, 6, 6, 6, 6, 5,
 5, 6, 4, 5, 6, 5, 6, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5,
 6, 6, 5, 5, 6, 6, 5, 6, 5, 6, 5, 4, 6, 5, 6, 6, 5, 6, 5, 6, 5,
 5, 5, 5, 6, 6, 5, 6, 5, 4, 6, 5, 6, 6, 5, 5, 6, 5, 5, 6, 6, 6,
 6, 5, 5, 6, 6, 5, 5, 6, 6, 4, 5, 5, 6, 6, 6, 6, 4, 6, 6, 6, 6,
 6, 5, 6, 6, 6, 6, 5, 6, 6, 6, 5, 5, 6, 5, 6, 6, 6, 5, 2, 6, 6,
 5, 5, 6, 5, 6, 6, 5, 6, 5, 6, 5, 6, 6, 6, 5, 5, 4, 5, 6, 6, 5,
 5, 4, 4, 6, 6,

Re: [R] Transform pairwise observations into a table

2012-10-01 Thread arun

HI AHJ,
No problem.

One more way in addition to reshape() (Rui's suggestion) to get the same result.
library(reshape)

as.matrix(cast(melt(dat1,id=c(ind1,ind2)),ind1~ind2,value=value))
#  1 2 3   4
#1 1.000 0.250 0.125 0.5
#2 0.250 1.000 0.125 0.5
#3 0.125 0.125 1.000 0.5
#4 0.500 0.500 0.500 1.0
A.K.



- Original Message -
From: Hadjixenofontos, Athena ahadjixenofon...@med.miami.edu
To: arun smartpink...@yahoo.com
Cc: R help r-help@r-project.org
Sent: Monday, October 1, 2012 12:59 PM
Subject: Re: [R] Transform pairwise observations into a table

Thank you. I had looked at xtabs but misunderstood the syntax. This is great. :)

AHJ



On Oct 1, 2012, at 12:53 PM, arun smartpink...@yahoo.com wrote:

 Hi,
 Try this:
 
 dat1-read.table(text= 
 ind1 ind2 coef 
 1 1 1 
 1 2 0.25 
 1 3 0.125 
 1 4 0.5 
 2 2 1 
 2 1 0.25 
 2 3 0.125 
 2 4 0.5 
 3 3 1 
 3 1 0.125 
 3 2 0.125 
 3 4 0.5 
 4 4 1 
 4 1 0.5 
 4 2 0.5 
 4 3 0.5 
 ,sep=,header=TRUE) 
 mat1-as.matrix(xtabs(coef~ind1+ind2,data=dat1)) 
 
 #    ind2 
 #ind1     1     2     3     4 
   # 1 1.000 0.250 0.125 0.500 
    #2 0.250 1.000 0.125 0.500 
    #3 0.125 0.125 1.000 0.500 
    #4 0.500 0.500 0.500 1.000 
 
 A.K.
 
 
 
 - Original Message -
 From: AHJ ahadjixenofon...@med.miami.edu
 To: r-help@r-project.org
 Cc: 
 Sent: Monday, October 1, 2012 12:17 PM
 Subject: [R] Transform pairwise observations into a table
 
 Hi,
 
 I have a table of pairs of individuals and a coefficient that belongs to the
 pair:
 
 ind1    ind2    coef
 1    1    1
 1    2    0.25
 1    3    0.125
 1    4    0.5
 2    2    1
 2    1    0.25
 2    3    0.125
 2    4    0.5
 3    3    1
 3    1    0.125
 3    2    0.125
 3    4    0.5
 4    4    1
 4    1    0.5
 4    2    0.5
 4    3    0.5
 
 And I want to convert it to a matrix where each individual is both a row and
 a column and at the intersection of each pair is the coefficient that
 belongs to that pair:
 
     1           2           3             4
 1    1           0.25           0.125    0.5
 2    0.25      1           0.125    0.5
 3    0.125    0.125    1            0.5
 4    0.5           0.5           0.5            1
 
 If table() would allow me to specify something other than frequencies to
 fill the table with, it would be what I need. I tried a few different
 combinations of t() and unique() but none of it made enough sense to post as
 my starting code... I am just lost. Any help would be greatly appreciated.
 
 Thank you,
 AHJ
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Transform-pairwise-observations-into-a-table-tp4644706.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 


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

[R] Reading labels for very large heatmaps

2012-10-01 Thread JIMonroe

Hello,
I have a large (919 X 919), hierarchically clustered heatmap that I would
like to read the labels off of.  I have tried saving the figure in pdf
format and enlarging it until I can see the labels, but if I make the labels
small enough to read (so that they don't overlap) using cexRow and cexCol,
they do not appear in the pdf.  The limit seems to be anything below cexRow
or Col = 0.06.  Is there a way to have smaller labels visible in the pdf? 
Is this an issue with resolution?  I could get by with just a quarter of the
heatmap (i.e. half of a row and half of a column) so that the labels don't
have to be so small, but the heatmap must be clustered before this is done. 
I tried to cut the dendrogram at just the halfway point of the heatmap, but
was not successful.  Any suggestions?
Thanks,
Jacob



--
View this message in context: 
http://r.789695.n4.nabble.com/Reading-labels-for-very-large-heatmaps-tp4644739.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] (no subject)

2012-10-01 Thread bibek sharma

Hello,
I am a new R -user and request your help for the following problem.
I need to merge two dataset of longitudinal study which has two column
(id and respose) common. when I used merge option to join the datas
side be side,  because of the repeated subject id, I got larger data
set which is not accurate.

 I would like to connect twi data sets by id and response in such a
way that data are connected by same id and response type  and if the
same subject has less data point in one data set, then produce NA.
Sample data sets is attached.

Bibek
__
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] mlogit and model-based recursive partitioning

2012-10-01 Thread tudor

Hello:

Has anyone tried to model-based recursive partition (using mob from package
party; thanks Achim and colleagues) a data set based on a multinomial logit
model (using mlogit from package mlogit; thanks Yves)? 

I attempted to do so, but there are at least two reasons why I could not. 
First, in mob I am not quite sure that a model of class StatModel exists for
mlogit models.  Second, as mlogit uses the pipe character | to specify the
model, I wonder how this would interact with mob which uses pipe to
differentiate between explanatory and segmentation variables.

An example (not working) of what I would like to accomplish follows below.

Thanks a lot.
 
Tudor

library(party)
library(mlogit)
data(Fishing, package = mlogit)
Fish - mlogit.data(Fishing, varying = c(2:9), shape = wide, choice =
mode)
# FIT AN mlogit MODEL
m1 - mlogit(mode ~ price + catch, data=Fish)
# THE DESIRED END RESULT:  SEGMENT m1 BASED ON INCOME AND/OR OTHER POSSIBLE
COVARIATES





--
View this message in context: 
http://r.789695.n4.nabble.com/mlogit-and-model-based-recursive-partitioning-tp4644743.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] glmmPQL and spatial autocorrelation

2012-10-01 Thread Luis Huckstadt

Hi all,

I am analyzing data on habitat utilization of seals in the Southern Ocean.
My data show spatial autocorrelation, which I'm interested in incorporating
into my model. I am trying to model the presence of dives (versus simulated
pseudo-absences) using a binomial generalized binomial model (glmmPQL),
since I can incorporate the autocorrelation structure to the model using
that package. However, I'm trying to come up with ways to evaluate which
semivariogram is the most adequate for my model, but not having AICs
complicated things a bit for me. Any suggestions on how to best achieve
this?

Thanks in advance

Luis A. Huckstadt, Ph.D.
Department of Ecology and Evolutionary Biology
University of California Santa Cruz
Long Marine Lab
100 Shaffer Road
Santa Cruz, CA  95060

[[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] (no subject)

2012-10-01 Thread Bhupendrasinh Thakre

It will be great if you can reproduce some code and  data you are working on.


Bhupendrasinh Thakre





On Oct 1, 2012, at 4:25 PM, bibek sharma mbhpat...@gmail.com wrote:

 Hello,
 I am a new R -user and request your help for the following problem.
 I need to merge two dataset of longitudinal study which has two column
 (id and respose) common. when I used merge option to join the datas
 side be side,  because of the repeated subject id, I got larger data
 set which is not accurate.
 
 I would like to connect twi data sets by id and response in such a
 way that data are connected by same id and response type  and if the
 same subject has less data point in one data set, then produce NA.
 Sample data sets is attached.
 
 Bibek
 __
 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] Ifelse Execution

2012-10-01 Thread Bhupendrasinh Thakre

Hi Everyone,

I am trying to run a time based query and need some of your help.
Not much of data or packages.

Just a simple one.
Query I am trying to execute.

ifelse ((as.numeric(as.POSIXct(2012-10-01 20:38:00))), 
(rnorm(1,2,1)),(Sys.sleep()))

Note. Why I am using as.numeric is as I have a list of time at which I wanted 
to run the command.

Something like 

1349142243
1349138667

Question.

1. The query is running before the time reaches in both R-Studio and R-Terminal 
in Mac based System.
2. Sys.sleep () is ineffective in putting R on sleep till the time comes. 


Bhupendrasinh Thakre






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

2012-10-01 Thread Jeff Newmiller

Ifelse is a vector function and is absolutely inappropriate for that use.
Use if instead.
Also, read the help for Sys.sleep... you need to tell it how long you want to 
sleep. You should compute how long that is from now and sleep that long.
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

Bhupendrasinh Thakre vickytha...@gmail.com wrote:

Hi Everyone,

I am trying to run a time based query and need some of your help.
Not much of data or packages.

Just a simple one.
Query I am trying to execute.

ifelse ((as.numeric(as.POSIXct(2012-10-01 20:38:00))),
(rnorm(1,2,1)),(Sys.sleep()))

Note. Why I am using as.numeric is as I have a list of time at which I
wanted to run the command.

Something like 

1349142243
1349138667

Question.

1. The query is running before the time reaches in both R-Studio and
R-Terminal in Mac based System.
2. Sys.sleep () is ineffective in putting R on sleep till the time
comes. 


Bhupendrasinh Thakre






   [[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] Transform pairwise observations into a table