Re: [R] getting the results of tapply into a single matrix

[arun]

Hi,
Try:
example <- read.table(text="ID Sex Location CL
1   F    lake1   40
1   F    lake1    
1   F    lake1 43
2   M    lake1    30
3   M    lake2    22
4   F    lake2 25
4   F    lake2 27",sep="",header=TRUE,stringsAsFactors=FALSE,fill=TRUE) 


aggregate(CL~.,example,mean,na.rm=TRUE)

#or
library(plyr)
ddply(example,.(ID,Sex,Location),summarize,CL=mean(CL,na.rm=TRUE))


#or from `result`
 res <- 
setNames(cbind(expand.grid(dimnames(result),stringsAsFactors=FALSE),as.data.frame(as.matrix(result))),c("Location","Sex","ID","CL"))
 res[!is.na(res$CL),c(3,2,1,4)]

###It would be better to store the results in a data.frame than in a matrix for 
this case.

A.K.



I have a table with three categorical columns (ID, Sex, Location) and a 
measurement column (CL). These are measurements from individuals in a 
study. Some individuals were measured more than once and thus have 
multiple rows. For those individuals, I need to take an average of all 
of their measurements so that I can run statistical tests without having 
pseudoreplication. Below is an example table with the code that I am 
using, the results that I am getting, and the results that I want (I am 
calling the table "example"). 
  
ID Sex Location CL 
1   F    lake1     40 
1   F    lake1     
1   F    lake1     43 
2   M    lake1    30 
3   M    lake2    22 
4   F    lake2     25 
4   F    lake2     27 

> result <- with(example, tapply(CL, list(Location, Sex, ID), mean, na.rm=T)) 

this almost does what I want. I takes the mean for each ID, and 
retains its relationship to the categorical variables, the problem is 
the output looks like this: 

, , 1 

             F M 
lake1 41.7   
lake2           

, , 2 

      F        M 
lake1   30.0 
lake2           

, , 3 

      F        M 
lake1           
lake2   22.0 

, , 4 

             F M 
lake1           
lake2 26.0   

How do I get those results back into a single table like I 
originally had. In other words, I want a table that looks like this:   

ID Sex Location CL 
1   F    lake1     41.7 
2   M    lake1    30 
3   M    lake2    22 
4   F    lake2     26.0 
  

Thanks for the help!

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

[R] Error in MuMIn "models are not all fitted to the same data"

[Lilly Dethier]

I'm pretty new to GLMMs and model averaging, but think I'm getting some
understanding of it all through lots of reading. However, I keep receiving
an error message when trying to average models that I don't understand and
can't find any resources about. I'm doing science education research trying
to evaluate population demographic factors that predict biology student
math performance. I have a lot of factors and so I tested a lot of models.
6 of my models had pretty similar AIC values (and evidence ratios of less
than 2.7) so I'm trying to average them. I keep receiving an error message
that says the models are not fitted to the same data, but I have no idea
how this is possible because all the models are from the same set of data
(same file and same variables)...strangely it seems to work when I try to
average MEx7, MEx10, & MEx22 only OR MEx24, MEx29, and MEx47 only. My code
is below. Any ideas? Thanks for any advice you can offer!!

library(MuMIn)
MEx7=lmer(cbind(c.score, w.score) ~ year + transfer + gender + p.math +
(1|section) + (1|quarter), family=binomial, data=survey.full, REML=F)
MEx10=lmer(cbind(c.score, w.score) ~ transfer + gender + p.math + Pmajor +
(1|section) + (1|quarter), family=binomial, data=survey.full, REML=F)
MEx22=lmer(cbind(c.score, w.score) ~ year + transfer + p.math + (1|section)
+ (1|quarter), family=binomial, data=survey.full, REML=F)
MEx24=lmer(cbind(c.score, w.score) ~ transfer + gender + p.math +
(1|section) + (1|quarter), family=binomial, data=survey.full, REML=F)
MEx29=lmer(cbind(c.score, w.score) ~ transfer + p.math + Pmajor +
(1|section) + (1|quarter), family=binomial, data=survey.full, REML=F)
MEx47=lmer(cbind(c.score, w.score) ~ transfer + p.math + (1|section) +
(1|quarter), family=binomial, data=survey.full, REML=F)
MExAvg=model.avg(rank=AIC, MEx24, MEx7, MEx10, MEx47, MEx29, MEx22)

Error in model.avg.default(rank = AIC, MEx24, MEx7, MEx10, MEx47, MEx29,  :
  models are not all fitted to the same data
Lilly Dethier

[[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] Find the cutoff correlation value for Pearson correlation test

[Jim Lemon]

On 11/15/2013 12:53 PM, jpm miao wrote:

Hi,

I find a few Pearson correlation test functions like

fBasics::correlationTest or stats::cor.test

which give the p-value of the test result. Is there a function that
calculate the cutoff correlation value for a specific p-value , e.g., p =
0.05?

I have a plot for the cross correlations between two time series, and I
would like to add a horizontal line that marks the significance of the
correlations.


Hi Miao,
Perhaps you could use the confidence interval.

Jim

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

[R] Find the cutoff correlation value for Pearson correlation test

[jpm miao]

Hi,

   I find a few Pearson correlation test functions like

   fBasics::correlationTest or stats::cor.test

which give the p-value of the test result. Is there a function that
calculate the cutoff correlation value for a specific p-value , e.g., p =
0.05?

I have a plot for the cross correlations between two time series, and I
would like to add a horizontal line that marks the significance of the
correlations.

Thanks,

Miao

[[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] Inconsistent results between caret+kernlab versions

[Andrew Digby]

I'm using caret to assess classifier performance (and it's great!). However, 
I've found that my results differ between R2.* and R3.* - reported accuracies 
are reduced dramatically. I suspect that a code change to kernlab ksvm may be 
responsible (see version 5.16-24 here: 
http://cran.r-project.org/web/packages/caret/news.html). I get very different 
results between caret_5.15-61 + kernlab_0.9-17 and caret_5.17-7 + 
kernlab_0.9-19 (see below).

Can anyone please shed any light on this?

Thanks very much!


### To replicate:

require(repmis)  # For downloading from https
df <- source_data('https://dl.dropboxusercontent.com/u/47973221/data.csv', 
sep=',')
require(caret)
svm.m1 <- 
train(df[,-1],df[,1],method='svmRadial',metric='Kappa',tunelength=5,trControl=trainControl(method='repeatedcv',
 number=10, repeats=10, classProbs=TRUE))
svm.m1
sessionInfo()

### Results - R2.15.2

> svm.m1
1241 samples
   7 predictors
  10 classes: ‘O27479’, ‘O31403’, ‘O32057’, ‘O32059’, ‘O32060’, ‘O32078’, 
‘O32089’, ‘O32663’, ‘O32668’, ‘O32676’ 

No pre-processing
Resampling: Cross-Validation (10 fold, repeated 10 times) 

Summary of sample sizes: 1116, 1116, 1114, 1118, 1118, 1119, ... 

Resampling results across tuning parameters:

  C Accuracy  Kappa  Accuracy SD  Kappa SD
  0.25  0.684 0.63   0.0353   0.0416  
  0.5   0.729 0.685  0.0379   0.0445  
  1 0.756 0.716  0.0357   0.0418  

Tuning parameter ‘sigma’ was held constant at a value of 0.247
Kappa was used to select the optimal model using  the largest value.
The final values used for the model were C = 1 and sigma = 0.247. 
> sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_NZ.UTF-8/en_NZ.UTF-8/en_NZ.UTF-8/C/en_NZ.UTF-8/en_NZ.UTF-8

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

other attached packages:
 [1] e1071_1.6-1 class_7.3-5 kernlab_0.9-17  repmis_0.2.4
caret_5.15-61   reshape2_1.2.2  plyr_1.8lattice_0.20-10 foreach_1.4.0   
cluster_1.14.3 

loaded via a namespace (and not attached):
 [1] codetools_0.2-8 compiler_2.15.2 digest_0.6.0evaluate_0.4.3  
formatR_0.7 grid_2.15.2 httr_0.2iterators_1.0.6 knitr_1.1   
RCurl_1.95-4.1  stringr_0.6.2   tools_2.15.2  

### Results - R3.0.2

> require(caret)
> svm.m1 <- 
> train(df[,-1],df[,1],method=’svmRadial’,metric=’Kappa’,tunelength=5,trControl=trainControl(method=’repeatedcv’,
>  number=10, repeats=10, classProbs=TRUE))
Loading required package: class
Warning messages:
1: closing unused connection 4 
(https://dl.dropboxusercontent.com/u/47973221/df.Rdata) 
2: executing %dopar% sequentially: no parallel backend registered 
> svm.m1
1241 samples
   7 predictors
  10 classes: ‘O27479’, ‘O31403’, ‘O32057’, ‘O32059’, ‘O32060’, ‘O32078’, 
‘O32089’, ‘O32663’, ‘O32668’, ‘O32676’ 

No pre-processing
Resampling: Cross-Validation (10 fold, repeated 10 times) 

Summary of sample sizes: 1118, 1117, 1115, 1117, 1116, 1118, ... 

Resampling results across tuning parameters:

  C Accuracy  Kappa  Accuracy SD  Kappa SD
  0.25  0.372 0.278  0.0330.0371  
  0.5   0.39  0.297  0.0317   0.0358  
  1 0.399 0.307  0.0289   0.0323  

Tuning parameter ‘sigma’ was held constant at a value of 0.2148907
Kappa was used to select the optimal model using  the largest value.
The final values used for the model were C = 1 and sigma = 0.215. 
> sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_NZ.UTF-8/en_NZ.UTF-8/en_NZ.UTF-8/C/en_NZ.UTF-8/en_NZ.UTF-8

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

other attached packages:
 [1] e1071_1.6-1 class_7.3-9 kernlab_0.9-19  repmis_0.2.6.2  
caret_5.17-7reshape2_1.2.2  plyr_1.8lattice_0.20-24 foreach_1.4.1   
cluster_1.14.4 

loaded via a namespace (and not attached):
[1] codetools_0.2-8 compiler_3.0.2  digest_0.6.3grid_3.0.2  httr_0.2
iterators_1.0.6 RCurl_1.95-4.1  stringr_0.6.2   tools_3.0.2  

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] optimization: multiple assignment problem

[Hans W.Borchers]

Jean-Francois Chevalier  bisnode.com> writes:
> 

You have already given the answer yourself. You have binary variables x(j, i),
you need to set up the inequalities, and then apply one of the mixed-integer
linear programming solvers in R, for instance 'lpSolve', 'Rglpk', 'Rsymphony'.

Setting up the inequalities may be slightly involved. You have not provided
reproducible code, so no concrete answer possible.

Hans Werner

> Hello,
> I'm trying to solve a multiple assignment problem.
> I found a package Adagio and its function mknapsack which
> maximize vstar = p(1)*(x(1,1) + ... + x(m,1)) + ... ... + 
>  p(n)*(x(1,n) + ... + x(m,n))
> subject to w(1)*x(i,1) + ... + w(n)*x(i,n) <= k(i) for i=1,...,m
> x(1,j) + ... + x(m,j) <= 1 for j=1,...,n x(i,j) = 0 or 1 
> for i=1,...,m , j=1,...,n ,
> It's close to what I'm trying to do except that
> 1)k(i) = k for any I (not an issue)
> 2)p is dependent of the item AND the knapsack
> 3)each item must be assigned
> maximize vstar = p(1,1)*x(1,1) + ... + p(m,1)*x(m,1) + ... ... + 
> p(1,n)*x(1,n) + ... + p(m,n)*x(m,n)
> with p(j,i) profit of assigning item i to knapsack j
> subject to w(1)*x(i,1) + ... + w(n)*x(i,n) <= k for i=1,...,m
> x(1,j) + ... + x(m,j) = 1 for j=1,...,n x(i,j) = 0 or 1 
>  for i=1,...,m , j=1,...,n ,
> It would be really helpful if you could indicate me any package, 
> function that would solve my problem?
> Thanks in advance,
> Best regards,
> Jean-François

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Double Pareto Log Normal Distribution DPLN

[b. alzahrani]

Thanks Dave for your help on this. Waiting other suggestions from the list on 
this.

Regards
Bander

> On 14 Nov 2013, at 09:29 pm, "David R Forrest"  wrote:
> 
> Hi Bander,
> 
> I'm pushing this discussion back to the list, because I'm not sure of the 
> shape/rate parameters for rpareto and rexp and how they'd be applied across 
> this mix of typo'd papers.
> 
> # Reed Equation 6 http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf 
> exponentiated per end of sec 3:
> rdpln<-function(n,a=1,b=1,t=1,v=0){exp(v+t*rnorm(n,sd=t)+rexp(n,rate=1/a)-rexp(n,rate=1/b))}
> 
> # Reed Equation 10:
> library(VGAM)
> rdpln2<-function(n,a,b,v,t){ rlnorm(n,meanlog=v,sdlog=t)*
>   rpareto(n,location=1,shape=1/a)/
>   rpareto(n,location=1,shape=1/b)}
> 
> boxplot(data.frame(x1= log(rdpln(a=2.5,b=.01,t=0.45,v=6.5,n=10)), x2= 
> log(rdpln2(a=2.5,b=.01,t=0.45,v=6.5,n=10  
> 
> # Reed equation 8 http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf 
> # with S1 errata #1 from 
> http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0048964 
> 
> ddpln <- function(x, a=1, b=1, v=0, t=1){
> # Density of Double Pareto LogNormal distribution
> # from "b. alzahrani"  email of 2013-11-13
> # per formula 8 from http://cs.stanford.edu/people/jure/pubs/dpln-kdd08.pdf
> 
>  c <- (a * b /(a+b))
> 
>  norm1<-pnorm((log(x)-v-(a*t^2))/t)
>  norm2<-pnorm((log(x)-v+(b*t^2))/t)
>  expo1<-  a*v+(a^2*t^2/2)
>  expo2<- -b*v+(b^2*t^2/2)  # edited from the paper's eqn 8  with: s/t/v/
> 
>  z<- (x^(-a-1)) * exp(expo1)*(  norm1)
>  y<- (x^( b-1)) * exp(expo2)*(1-norm2)  # 1-norm is the complementary CDF of 
> N(0,1)
> 
>  c*(z+y)
> }
> 
> 
> Dave
> 
> 
> 
> On Nov 14, 2013, at 9:12 AM, David Forrest 
> wrote:
> 
>> 
>> I think exponentiation of eqn 6 from the Reed paper generates DPLN variates 
>> directly, so maybe:
>> 
>> rdpln=function(n,a=1,b=1,t=1,v=0){exp(v+t*rnorm(n,sd=t)+rexp(n,rate=1/a)-rexp(n,rate=1/b))}
>> 
>> 
>> Dave
>> 
>> 
>> On Nov 13, 2013, at 4:34 PM, "b. alzahrani" 
>> wrote:
>> 
>>> You help is much appreciated. Just one last point if you could help, Now I 
>>> want to pass this curve to a function that can generate random numbers 
>>> distributed according to DPLN ( right curve).
>>> 
>>> I found the package Runuran can do that 
>>> http://cran.r-project.org/web/packages/Runuran/Runuran.pdf  but I do not 
>>> know how to do it I think it would be something similar to page 8 and 9.
>>> 
>>> Regards
>>> **
>>> Bander 
>>> *
>>> 
>>> 
>>> 
>>> From: d...@vims.edu
>>> To: cs_2...@hotmail.com
>>> Subject: Re: [R] Double Pareto Log Normal Distribution DPLN
>>> Date: Wed, 13 Nov 2013 21:13:43 +
>>> 
>>> 
>>> I read the parameters in Fig 4, right as "DPLN[2.5,0.1,0.45,6.5]", so:
>>> 
>>> x<- 10^seq(0,4,by=.1)
>>> plot(x,ddpln(x,a=2.5,b=.01,v=6.5,t=0.45),log='xy',type='l')
>>> 
>>> ... and the attached graph does not look dissimilar the figure--It starts 
>>> at 10^-2, goes through 10^-4 at x=100, and elbows down around 900 and 
>>> passes through 10^-6 at about 2000.
>>> 
>>> The correction of Reed helps -- The uncorrected Reed Eq9 equation suggests 
>>> that the the 't' in Sehshadri Eq9 should be a 'v' , but it doesn't exactly 
>>> make sense with the extra 'a' in there.  If the errata clears that up, then 
>>> your expo2 term looks just like the expo1 term, but with a=-b.
>>> 
>>> 
>>> 
>>> 
>>> 
>>> On Nov 13, 2013, at 3:43 PM, "b. alzahrani" wrote: > Thank you very much 
>>> for the help and the change you suggested in my code, I also found a 
>>> correction on equation 9 that has been published by Reed ( here 
>>> http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf i.e. the original paper on 
>>> Double Pareto Log Normal Distribution ). > > can you please see the 
>>> correction in this 
>>> linkhttp://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0048964
>>>  (in Supporting Information section, appendix S1), does your suggested code 
>>> coincide with the correction on this link? as I can see > > Actually, I am 
>>> interested in the most right curve in figure 4. and if we plot the curve 
>>> with same order of the paper's parameters you suggests: 
>>> plot(x,ddpln(x,a=2.8,b=.01,v=6.5,t=0.45),log='xy',type='l') the curve is 
>>> different from the one in the paper?? > > Thanks > > 
>>> ** > Bander 
>>> Alzahrani, > * > > > !
 > > From: d...@vims.edu > > To: cs_2...@hotmail.com > > CC: 
 > > r-h...@stat.math.ethz.ch > > Subject: Re: [R] Double Pareto Log Normal 
 > > Distribution DPLN > > Date: Wed, 13 Nov 2013 19:09:34 + > > > > 
 > > ...Additionally...your set of parameters match none of the curves in 
 > > figure 4. > > > > I think the ordering of the parameters as listed on the 
 > > graphs is different than in the text of the article. > > > > Th

[R] Warning message during starts up

[Dereje Fentie]

During startup R gives a Warning message: "Setting LC_CTYPE=en_US gnumeric
& failed" What is the error and how can I fix it?

[[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] Windows 7/Rstudio/Rattle - rattle() Causes R Session to abort

[Duncan Murdoch]

On 13-11-14 4:41 PM, Lopez, Dan wrote:

It's happening in the 64bit installation of R I have too (w/o Rstudio). This is 
R version 2.15.1


That's kind of an old R version (from June 2012); there were two more 
releases in the 2.15.x series.  The current release is 3.0.2.  Can you 
upgrade?


Duncan Murdoch



Dan


-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com]
Sent: Thursday, November 14, 2013 1:38 PM
To: Lopez, Dan; R help (r-help@r-project.org)
Subject: Re: [R] Windows 7/Rstudio/Rattle - rattle() Causes R Session to abort

On 13-11-14 3:50 PM, Lopez, Dan wrote:> Windows 7, R 2.15.1 64bit, RStudio 0.97.312, 
Rattle 2.6.26  >  > Hi,  > Please help.
  >
  > I removed rattle then reinstalled then loaded then attempted to open rattle 
with rattle(). But each time I do that My R session in R studio aborts and 
restarts.
  >
  > remove.packages("rattle")
  > install.packages("rattle")
  > library(rattle)
  >
  > How can I fix this? I saw something in rattle-users but that applied to Mac.

Do you have problems running rattle in the plain Rgui.exe?  If not, then this 
sounds like an RStudio problem, and you should use one of their support sites.

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] Windows 7/Rstudio/Rattle - rattle() Causes R Session to abort

[Lopez, Dan]

It's happening in the 64bit installation of R I have too (w/o Rstudio). This is 
R version 2.15.1

Dan


-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com] 
Sent: Thursday, November 14, 2013 1:38 PM
To: Lopez, Dan; R help (r-help@r-project.org)
Subject: Re: [R] Windows 7/Rstudio/Rattle - rattle() Causes R Session to abort

On 13-11-14 3:50 PM, Lopez, Dan wrote:> Windows 7, R 2.15.1 64bit, RStudio 
0.97.312, Rattle 2.6.26  >  > Hi,  > Please help.
 >
 > I removed rattle then reinstalled then loaded then attempted to open rattle 
 > with rattle(). But each time I do that My R session in R studio aborts and 
 > restarts.
 >
 > remove.packages("rattle")
 > install.packages("rattle")
 > library(rattle)
 >
 > How can I fix this? I saw something in rattle-users but that applied to Mac.

Do you have problems running rattle in the plain Rgui.exe?  If not, then this 
sounds like an RStudio problem, and you should use one of their support sites.

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] Windows 7/Rstudio/Rattle - rattle() Causes R Session to abort

[Duncan Murdoch]

On 13-11-14 3:50 PM, Lopez, Dan wrote:> Windows 7, R 2.15.1 64bit, 
RStudio 0.97.312, Rattle 2.6.26

>
> Hi,
> Please help.
>
> I removed rattle then reinstalled then loaded then attempted to open 
rattle with rattle(). But each time I do that My R session in R studio 
aborts and restarts.

>
> remove.packages("rattle")
> install.packages("rattle")
> library(rattle)
>
> How can I fix this? I saw something in rattle-users but that applied 
to Mac.


Do you have problems running rattle in the plain Rgui.exe?  If not, then 
this sounds like an RStudio problem, and you should use one of their 
support sites.


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] Double Pareto Log Normal Distribution DPLN

[David R Forrest]

Hi Bander,

I'm pushing this discussion back to the list, because I'm not sure of the 
shape/rate parameters for rpareto and rexp and how they'd be applied across 
this mix of typo'd papers.

# Reed Equation 6 http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf exponentiated 
per end of sec 3:
rdpln<-function(n,a=1,b=1,t=1,v=0){exp(v+t*rnorm(n,sd=t)+rexp(n,rate=1/a)-rexp(n,rate=1/b))}

# Reed Equation 10:
library(VGAM)
rdpln2<-function(n,a,b,v,t){ rlnorm(n,meanlog=v,sdlog=t)*
   rpareto(n,location=1,shape=1/a)/
   rpareto(n,location=1,shape=1/b)}
 
boxplot(data.frame(x1= log(rdpln(a=2.5,b=.01,t=0.45,v=6.5,n=10)), x2= 
log(rdpln2(a=2.5,b=.01,t=0.45,v=6.5,n=10  

# Reed equation 8 http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf 
# with S1 errata #1 from 
http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0048964 

ddpln <- function(x, a=1, b=1, v=0, t=1){
 # Density of Double Pareto LogNormal distribution
 # from "b. alzahrani"  email of 2013-11-13
 # per formula 8 from http://cs.stanford.edu/people/jure/pubs/dpln-kdd08.pdf

  c <- (a * b /(a+b))

  norm1<-pnorm((log(x)-v-(a*t^2))/t)
  norm2<-pnorm((log(x)-v+(b*t^2))/t)
  expo1<-  a*v+(a^2*t^2/2)
  expo2<- -b*v+(b^2*t^2/2)  # edited from the paper's eqn 8  with: s/t/v/

  z<- (x^(-a-1)) * exp(expo1)*(  norm1)
  y<- (x^( b-1)) * exp(expo2)*(1-norm2)  # 1-norm is the complementary CDF of 
N(0,1)

  c*(z+y)
}


Dave



On Nov 14, 2013, at 9:12 AM, David Forrest 
 wrote:

> 
> I think exponentiation of eqn 6 from the Reed paper generates DPLN variates 
> directly, so maybe:
> 
> rdpln=function(n,a=1,b=1,t=1,v=0){exp(v+t*rnorm(n,sd=t)+rexp(n,rate=1/a)-rexp(n,rate=1/b))}
> 
> 
> Dave
> 
> 
> On Nov 13, 2013, at 4:34 PM, "b. alzahrani" 
> wrote:
> 
>> You help is much appreciated. Just one last point if you could help, Now I 
>> want to pass this curve to a function that can generate random numbers 
>> distributed according to DPLN ( right curve).
>> 
>> I found the package Runuran can do that 
>> http://cran.r-project.org/web/packages/Runuran/Runuran.pdf  but I do not 
>> know how to do it I think it would be something similar to page 8 and 9.
>> 
>> Regards
>> **
>> Bander 
>> *
>> 
>> 
>> 
>> From: d...@vims.edu
>> To: cs_2...@hotmail.com
>> Subject: Re: [R] Double Pareto Log Normal Distribution DPLN
>> Date: Wed, 13 Nov 2013 21:13:43 +
>> 
>> 
>> I read the parameters in Fig 4, right as "DPLN[2.5,0.1,0.45,6.5]", so:
>> 
>> x<- 10^seq(0,4,by=.1)
>> plot(x,ddpln(x,a=2.5,b=.01,v=6.5,t=0.45),log='xy',type='l')
>> 
>> ... and the attached graph does not look dissimilar the figure--It starts at 
>> 10^-2, goes through 10^-4 at x=100, and elbows down around 900 and passes 
>> through 10^-6 at about 2000.
>> 
>> The correction of Reed helps -- The uncorrected Reed Eq9 equation suggests 
>> that the the 't' in Sehshadri Eq9 should be a 'v' , but it doesn't exactly 
>> make sense with the extra 'a' in there.  If the errata clears that up, then 
>> your expo2 term looks just like the expo1 term, but with a=-b.
>> 
>> 
>> 
>> 
>> 
>> On Nov 13, 2013, at 3:43 PM, "b. alzahrani" wrote: > Thank you very much for 
>> the help and the change you suggested in my code, I also found a correction 
>> on equation 9 that has been published by Reed ( here 
>> http://www.math.uvic.ca/faculty/reed/dPlN.3.pdf i.e. the original paper on 
>> Double Pareto Log Normal Distribution ). > > can you please see the 
>> correction in this 
>> linkhttp://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0048964
>>  (in Supporting Information section, appendix S1), does your suggested code 
>> coincide with the correction on this link? as I can see > > Actually, I am 
>> interested in the most right curve in figure 4. and if we plot the curve 
>> with same order of the paper's parameters you suggests: 
>> plot(x,ddpln(x,a=2.8,b=.01,v=6.5,t=0.45),log='xy',type='l') the curve is 
>> different from the one in the paper?? > > Thanks > > 
>> ** > Bander 
>> Alzahrani, > * > > > > > From: 
>> d...@vims.edu > > To: cs_2...@hotmail.com > > CC: r-h...@stat.math.ethz.ch > 
>> > Subject: Re: [R] Double Pareto Log Normal Distribution DPLN > > Date: Wed, 
>> 13 Nov 2013 19:09:34 + > > > > ...Additionally...your set of parameters 
>> match none of the curves in figure 4. > > > > I think the ordering of the 
>> parameters as listed on the graphs is different than in the text of the 
>> article. > > > > The 'v' parameter controls the location of the 'elbow' and 
>> should be near log(x) in each graph, while the 't' parameter is the 
>> sharpness of the elbow. So eyeballing the elbows: > > > > 
>> log(c(80,150,300))= 4.382027 5.010635 5.703782 > > > > These appear to match 
>> positional parameter #4 in the legends, which y

[R] Windows 7/Rstudio/Rattle - rattle() Causes R Session to abort

[Lopez, Dan]

Windows 7, R 2.15.1 64bit, RStudio 0.97.312, Rattle 2.6.26

Hi,
Please help.

I removed rattle then reinstalled then loaded then attempted to open rattle 
with rattle(). But each time I do that My R session in R studio aborts and 
restarts.

remove.packages("rattle")
install.packages("rattle")
library(rattle)

How can I fix this? I saw something in rattle-users but that applied to Mac.
Dan

[[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] Fitting arbitrary curve to 1D data with error bars

[Erkcan Özcan]

Thanks, this was a useful pointer. Since the function I am trying to fit is 
exponential, I decided to use nls. And I was able to reproduce exactly the 
results and the plot in the URL I had posted. For future reference, here is the 
R code I wrote:


require("gplots")
xx <- 1:10
yy <- c(1.56,1.20,1.10,0.74,0.57,0.55,0.31,0.27,0.28,0.11)
dy <- c(0.02,0.02,0.20,0.03,0.03,0.10,0.05,0.02,0.10,0.05)
plotCI(xx,yy,uiw=dy,gap=0)
nlmod <- nls(yy ~ Alpha * exp(Beta * xx), start=list(Alpha=1, Beta=-1))
lines(xx, predict(nlmod), col = "blue")
wnlmod <- nls(yy ~ Alpha * exp(Beta * xx), start=list(Alpha=1, Beta=-1), 
weights=dy^-2)
lines(xx, predict(wnlmod), col = "red")


Thanks to everybody who responded,

e.

On 14 Nov 2013, at 11:34, Suzen, Mehmet wrote:

> Maybe you are after "weights" option given by 'lm' or 'glm'
> 
> See: 
> http://stackoverflow.com/questions/6375650/function-for-weighted-least-squares-estimates
> 
> On 14 November 2013 10:01, Erkcan Özcan  wrote:
>> Thanks, but if you have another closer look to my post, you will see that my 
>> question has nothing to do with drawing error bars on a plot.
>> 
>> What I want is to do a curve fit to a data with error bars.
>> 
>> Best,
>> e.
>> 
>> On 14 Nov 2013, at 04:21, Suzen, Mehmet wrote:
>> 
>>> If you are after adding error bars in a scatter plot; one example is
>>> given below :
>>> 
>>> #some example data
>>> set.seed(42)
>>> df <- data.frame(x = rep(1:10,each=5), y = rnorm(50))
>>> 
>>> #calculate mean, min and max for each x-value
>>> library(plyr)
>>> df2 <- ddply(df,.(x),function(df)
>>> c(mean=mean(df$y),min=min(df$y),max=max(df$y)))
>>> 
>>> #plot error bars
>>> library(Hmisc)
>>> with(df2,errbar(x,mean,max,min))
>>> grid(nx=NA,ny=NULL)
>>> 
>>> (From: 
>>> http://stackoverflow.com/questions/13032777/scatter-plot-with-error-bars)
>>> 
>> 
> 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] polygon circling a graph

[Carl Witthoft]

Please post the packages from which 'barabasi' and 'layout.fruch'
originate (not to mention whatever the plot() method is for whatever class
your 'g' is).  Further, without seeing what your data look like we have no
way of knowing whether you've fed the appropriate elements of "L" to chull.


email mail wrote
> Hi:
> 
> I want to create a polygon encircling a graph. For this i use convex
> hull  to get the coordinate points for polygon.
> 
> g <- barabasi.game(10)
> L<-layout.fruchterman.reingold(g)
> temp1 <- chull(L)
> temp1 <- c(temp1, temp1[1])
> plot(g, layout=layout.fruchterman.reingold)
> 
> 
> But when i plot the polygon with the code below, the polygon dosen't
> encircle the graph.
> 
> polygon(L[temp1, ], col = "#FFAA")
> 
> How can I plot a polygon circling a graph?
> 
> Regards:
> John
> 
> __

> R-help@

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





--
View this message in context: 
http://r.789695.n4.nabble.com/polygon-circling-a-graph-tp4680479p4680484.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] polygon circling a graph

[William Dunlap]

layout.fruchterman.reingold(g) returns a random result, so you
want to call it once and use the one return value.  Also, I think
you need to avoid the rescaling that plot.igraph does.  It looks
like you need to explicitly specify xlim and ylim if you do that,
but I may not have looked long enough at it.

   plot(g, layout=L, rescale=FALSE, xlim=range(L[,1]), ylim=range(L[,2]))
   polygon(L[temp1, ], col = "#FFAA")


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 email
> Sent: Thursday, November 14, 2013 10:04 AM
> To: r-help@r-project.org
> Subject: [R] polygon circling a graph
> 
> Hi:
> 
> I want to create a polygon encircling a graph. For this i use convex
> hull  to get the coordinate points for polygon.
> 
> g <- barabasi.game(10)
> L<-layout.fruchterman.reingold(g)
> temp1 <- chull(L)
> temp1 <- c(temp1, temp1[1])
> plot(g, layout=layout.fruchterman.reingold)
> 
> 
> But when i plot the polygon with the code below, the polygon dosen't
> encircle the graph.
> 
> polygon(L[temp1, ], col = "#FFAA")
> 
> How can I plot a polygon circling a graph?
> 
> Regards:
> John
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Column Name Matching in xts Objects

[Joshua Ulrich]

On Thu, Nov 14, 2013 at 12:06 PM, Pooya Lalehzari
 wrote:
> Hello,
> I noticed an unexpected behavior when using the xts object and I was 
> wondering if anyone knows why that happens. I have a code to create a new 
> column and copy one of the columns to the new column (please see below):
>
> library(xts)
> df = data.frame(stringsAsFactors=FALSE)
> df[1:3,"date"] = c("2001-1-1","2001-1-2","2001-1-3")
> df[1:3,"col1"] = c(1:3)
> df[1:3,"col2"] = c(-1:-3)
> df[1:3,"col3"] = c(101:103)
> rownames(df)= df$date
> df2 = as.xts(df)
> View(df2)
> Col_To_Copy = "col1"
> df2$NewCol = df2[,colnames(df2)==Col_To_Copy]
> View(df2)
>
> This code works, expect, I noticed that when the Col_To_Copy is not found, 
> the last column is renamed to "NewCol". For example:
> library(xts)
> df = data.frame(stringsAsFactors=FALSE)
> df[1:3,"date"] = c("2001-1-1","2001-1-2","2001-1-3")
> df[1:3,"col1"] = c(1:3)
> df[1:3,"col2"] = c(-1:-3)
> df[1:3,"col3"] = c(101:103)
> rownames(df)= df$date
> df2 = as.xts(df)
> View(df2)
> Col_To_Copy = "col5"  # the only line changed
> df2$NewCol = df2[,colnames(df2)==Col_To_Copy]
> View(df2)
>
> In this case, the last column's name is changed to "NewCol". Is that behavior 
> expected for xts objects?
>
xts doesn't have a $<- method, so $<-.zoo is dispatched.  I'm not
familiar with the function, but the behavior you found is explicitly
defined, so that seems to suggest it was intended.

>From "$<-.zoo":

  wi <- match(x, colnames(object))
  if(is.na(wi)) {
object <- cbind(object, value)
if(is.null(dim(object))) dim(object) <- c(length(object), 1)
colnames(object)[NCOL(object)] <- x
  } else {
...

Best,
--
Joshua Ulrich  |  about.me/joshuaulrich
FOSS Trading  |  www.fosstrading.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] optimization: multiple assignment problem

[Simon Zehnder]

It would be more clear if you tell, what you want to do instead of what you do 
not want to do. 

If you start with a usual cost matrix (whatever cost function you have) and you 
have to assign N to N this reduces to the well-known Munkre’s algorithm (see 
for example: http://gallery.rcpp.org/articles/minimal-assignment/). If you have 
to assign M to N, it is an extended version of the Munkre’s algorithm which has 
been covered by Bourgeois and Lassalle (1971). All this is graph theory. 

Best 

Simon
 
On 14 Nov 2013, at 19:24, Jean-Francois Chevalier 
 wrote:

> Hello,
> I'm trying to solve a multiple assignment problem.
> I found a package Adagio and its function mknapsack which
> 
> maximize vstar = p(1)*(x(1,1) + ... + x(m,1)) + ... ... + p(n)*(x(1,n) + ... 
> + x(m,n))
> 
> subject to w(1)*x(i,1) + ... + w(n)*x(i,n) <= k(i) for i=1,...,m
> 
> x(1,j) + ... + x(m,j) <= 1 for j=1,...,n x(i,j) = 0 or 1 for i=1,...,m , 
> j=1,...,n ,
> It's close to what I'm trying to do except that
> 1)k(i) = k for any I (not an issue)
> 2)p is dependent of the item AND the knapsack
> 3)each item must be assigned
> 
> maximize vstar = p(1,1)*x(1,1) + ... + p(m,1)*x(m,1) + ... ... + 
> p(1,n)*x(1,n) + ... + p(m,n)*x(m,n)
> 
> with p(j,i) profit of assigning item i to knapsack j
> 
> subject to w(1)*x(i,1) + ... + w(n)*x(i,n) <= k for i=1,...,m
> 
> x(1,j) + ... + x(m,j) = 1 for j=1,...,n x(i,j) = 0 or 1 for i=1,...,m , 
> j=1,...,n ,
> It would be really helpful if you could indicate me any package, function 
> that would solve my problem?
> Thanks in advance,
> Best regards,
> 
> Jean-François
> 
> 
> 
>  DISCLAIMER \ "This e-mail and any attachments t...{{dropped:13}}
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] polygon circling a graph

[email]

Hi:

I want to create a polygon encircling a graph. For this i use convex
hull  to get the coordinate points for polygon.

g <- barabasi.game(10)
L<-layout.fruchterman.reingold(g)
temp1 <- chull(L)
temp1 <- c(temp1, temp1[1])
plot(g, layout=layout.fruchterman.reingold)


But when i plot the polygon with the code below, the polygon dosen't
encircle the graph.

polygon(L[temp1, ], col = "#FFAA")

How can I plot a polygon circling a graph?

Regards:
John

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] optimization: multiple assignment problem

[Jean-Francois Chevalier]

Hello,
I'm trying to solve a multiple assignment problem.
I found a package Adagio and its function mknapsack which

maximize vstar = p(1)*(x(1,1) + ... + x(m,1)) + ... ... + p(n)*(x(1,n) + ... + 
x(m,n))

subject to w(1)*x(i,1) + ... + w(n)*x(i,n) <= k(i) for i=1,...,m

x(1,j) + ... + x(m,j) <= 1 for j=1,...,n x(i,j) = 0 or 1 for i=1,...,m , 
j=1,...,n ,
It's close to what I'm trying to do except that
1)k(i) = k for any I (not an issue)
2)p is dependent of the item AND the knapsack
3)each item must be assigned

maximize vstar = p(1,1)*x(1,1) + ... + p(m,1)*x(m,1) + ... ... + p(1,n)*x(1,n) 
+ ... + p(m,n)*x(m,n)

with p(j,i) profit of assigning item i to knapsack j

subject to w(1)*x(i,1) + ... + w(n)*x(i,n) <= k for i=1,...,m

x(1,j) + ... + x(m,j) = 1 for j=1,...,n x(i,j) = 0 or 1 for i=1,...,m , 
j=1,...,n ,
It would be really helpful if you could indicate me any package, function that 
would solve my problem?
Thanks in advance,
Best regards,

Jean-François



 DISCLAIMER \ "This e-mail and any attachments t...{{dropped:13}}

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

[R] Column Name Matching in xts Objects

[Pooya Lalehzari]

Hello,
I noticed an unexpected behavior when using the xts object and I was wondering 
if anyone knows why that happens. I have a code to create a new column and copy 
one of the columns to the new column (please see below):

library(xts)
df = data.frame(stringsAsFactors=FALSE)
df[1:3,"date"] = c("2001-1-1","2001-1-2","2001-1-3")
df[1:3,"col1"] = c(1:3)
df[1:3,"col2"] = c(-1:-3)
df[1:3,"col3"] = c(101:103)
rownames(df)= df$date
df2 = as.xts(df)
View(df2)
Col_To_Copy = "col1"
df2$NewCol = df2[,colnames(df2)==Col_To_Copy]
View(df2)

This code works, expect, I noticed that when the Col_To_Copy is not found, the 
last column is renamed to "NewCol". For example:
library(xts)
df = data.frame(stringsAsFactors=FALSE)
df[1:3,"date"] = c("2001-1-1","2001-1-2","2001-1-3")
df[1:3,"col1"] = c(1:3)
df[1:3,"col2"] = c(-1:-3)
df[1:3,"col3"] = c(101:103)
rownames(df)= df$date
df2 = as.xts(df)
View(df2)
Col_To_Copy = "col5"  # the only line changed
df2$NewCol = df2[,colnames(df2)==Col_To_Copy]
View(df2)

In this case, the last column's name is changed to "NewCol". Is that behavior 
expected for xts objects?


Thank you.



THIS E-MAIL IS FOR THE SOLE USE OF THE INTENDED RECIPIENT(S) AND MAY CONTAIN
CONFIDENTIAL AND PRIVILEGED INFORMATION.ANY UNAUTHORIZED REVIEW, USE, DISCLOSURE
OR DISTRIBUTION IS PROHIBITED. IF YOU ARE NOT THE INTENDED RECIPIENT, PLEASE
CONTACT THE SENDER BY REPLY E-MAIL AND DESTROY ALL COPIES OF THE ORIGINAL 
E-MAIL.
[[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 aggregated data to individual data

[arun]

HI,

A more general form would be:
 D[rep(1:nrow(D),D[,3]),-3]

A.K.


On Thursday, November 14, 2013 12:27 PM, arun  wrote:


Hi,
Try:
 D1 <- D[rep(row.names(D),D[,3]),-3] ##assuming rownames(D) are from 1:nrow(D)
 row.names(D1) <- 1:nrow(D1)
A.K.


On Thursday, November 14, 2013 5:32 AM, peron  wrote:
Hello 



I have data in following form : 100 ind describe by two variables x and y.



D<-data .frame( x=rnorm(3), y=rnorm(3), size=c(50,10,40))



I want data for individual, i.e, 100 observations for my 100 ind.



Thank for your help



Olivier Peron


    [[alternative HTML version deleted]]

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



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

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

Re: [R] From list to dataframe

[Hermann Norpois]

do.call (rbind, testlist) and Reduce (rbind, testlist) are perfect.
Thanks
hermann


2013/11/14 arun 

>
>
> Hi Hermann,
>
> You may try:
> do.call(rbind,testlist)
> #or
> Reduce(rbind,testlist)
> #or
> library(plyr)
>  ldply(testlist)
> A.K.
>
>
>
>
> On Thursday, November 14, 2013 8:27 AM, Hermann Norpois <
> hnorp...@gmail.com> wrote:
> Hello,
>
> having a list like testlist I would like to transform it in dataframe. How
> does it work?
>
> Thanks
> Hermann
>
> > testlist
> [[1]]
>   BP_A SNP_A BP_B SNP_B   R2
> 2 27001689 rs4822747 27002392 rs4820690 0.695642
> 3 27001689 rs4822747 27004298 rs5761627 0.695642
> 4 27001689 rs4822747 27004902 rs5752355 0.687861
> 5 27001689 rs4822747 27004964 rs4822748 0.682715
> 6 27001689 rs4822747 27005122 rs4820691 0.695642
> 7 27001689 rs4822747 27005158 rs4822749 0.695642
>
> [[2]]
>BP_A SNP_A BP_B SNP_B   R2
> 9  27002392 rs4820690 27001689 rs4822747 0.695642
> 11 27002392 rs4820690 27004298 rs5761627 1.00
> 12 27002392 rs4820690 27004902 rs5752355 0.988681
> 13 27002392 rs4820690 27004964 rs4822748 0.988655
> 14 27002392 rs4820690 27005122 rs4820691 1.00
> 15 27002392 rs4820690 27005158 rs4822749 1.00
> 16 27002392 rs4820690 27005194 rs4822750 0.988655
> 17 27002392 rs4820690 27005470 rs5761629 0.695642
>
> [[3]]
>BP_A SNP_A BP_B SNP_B   R2
> 19 27004298 rs5761627 27001689 rs4822747 0.695642
> 20 27004298 rs5761627 27002392 rs4820690 1.00
> 22 27004298 rs5761627 27004902 rs5752355 0.988681
> 23 27004298 rs5761627 27004964 rs4822748 0.988655
> 24 27004298 rs5761627 27005122 rs4820691 1.00
> 25 27004298 rs5761627 27005158 rs4822749 1.00
>
> #The aim is the dataframe
>BP_A SNP_A BP_B SNP_B   R2
> 2  27001689 rs4822747 27002392 rs4820690 0.695642
> 3  27001689 rs4822747 27004298 rs5761627 0.695642
> 4  27001689 rs4822747 27004902 rs5752355 0.687861
> 5  27001689 rs4822747 27004964 rs4822748 0.682715
> 6  27001689 rs4822747 27005122 rs4820691 0.695642
> 7  27001689 rs4822747 27005158 rs4822749 0.695642
> 9  27002392 rs4820690 27001689 rs4822747 0.695642
> 11 27002392 rs4820690 27004298 rs5761627 1.00
> 12 27002392 rs4820690 27004902 rs5752355 0.988681
> 13 27002392 rs4820690 27004964 rs4822748 0.988655
> 14 27002392 rs4820690 27005122 rs4820691 1.00
> 15 27002392 rs4820690 27005158 rs4822749 1.00
> 16 27002392 rs4820690 27005194 rs4822750 0.988655
> 17 27002392 rs4820690 27005470 rs5761629 0.695642
> 19 27004298 rs5761627 27001689 rs4822747 0.695642
> 20 27004298 rs5761627 27002392 rs4820690 1.00
> 22 27004298 rs5761627 27004902 rs5752355 0.988681
> 23 27004298 rs5761627 27004964 rs4822748 0.988655
> 24 27004298 rs5761627 27005122 rs4820691 1.00
> 25 27004298 rs5761627 27005158 rs4822749 1.00
> >
> > dput (testlist)
> list(structure(list(BP_A = c(27001689L, 27001689L, 27001689L,
> 27001689L, 27001689L, 27001689L), SNP_A = c("rs4822747", "rs4822747",
> "rs4822747", "rs4822747", "rs4822747", "rs4822747"), BP_B = c(27002392L,
> 27004298L, 27004902L, 27004964L, 27005122L, 27005158L), SNP_B =
> c("rs4820690",
> "rs5761627", "rs5752355", "rs4822748", "rs4820691", "rs4822749"
> ), R2 = c(0.695642, 0.695642, 0.687861, 0.682715, 0.695642, 0.695642
> )), .Names = c("BP_A", "SNP_A", "BP_B", "SNP_B", "R2"), row.names = 2:7,
> class = "data.frame"),
> structure(list(BP_A = c(27002392L, 27002392L, 27002392L,
> 27002392L, 27002392L, 27002392L, 27002392L, 27002392L), SNP_A =
> c("rs4820690",
> "rs4820690", "rs4820690", "rs4820690", "rs4820690", "rs4820690",
> "rs4820690", "rs4820690"), BP_B = c(27001689L, 27004298L,
> 27004902L, 27004964L, 27005122L, 27005158L, 27005194L, 27005470L
> ), SNP_B = c("rs4822747", "rs5761627", "rs5752355", "rs4822748",
> "rs4820691", "rs4822749", "rs4822750", "rs5761629"), R2 = c(0.695642,
> 1, 0.988681, 0.988655, 1, 1, 0.988655, 0.695642)), .Names = c("BP_A",
> "SNP_A", "BP_B", "SNP_B", "R2"), row.names = c(9L, 11L, 12L,
> 13L, 14L, 15L, 16L, 17L), class = "data.frame"), structure(list(
> BP_A = c(27004298L, 27004298L, 27004298L, 27004298L,
> 27004298L, 27004298L), SNP_A = c("rs5761627", "rs5761627",
> "rs5761627", "rs5761627", "rs5761627", "rs5761627"),
> BP_B = c(27001689L, 27002392L, 27004902L, 27004964L,
> 27005122L, 27005158L), SNP_B = c("rs4822747", "rs4820690",
> "rs5752355", "rs4822748", "rs4820691", "rs4822749"),
> R2 = c(0.695642, 1, 0.988681, 0.988655, 1, 1)), .Names = c("BP_A",
> "SNP_A", "BP_B", "SNP_B", "R2"), row.names = c(19L, 20L,
> 22L, 23L, 24L, 25L), class = "data.frame"))
>
> [[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] 2SLS for panel data, re

[Chanita Holmes]

Hi,

I am trying to estimate a 2sls using panel data (random effect model). I
tried the same estimation in STATA using the ivtreg2 command. However STATA
and R are giving me two different results. I figure there is something with
my R code:

iv=plm(formula=wecon~fdistockgdp +trade + polrightsreversed +lnrgdpch +
execleft + muslim2+c100rat +c111rat +yeardum| polrightsreversed+lnrgdpch+
execleft+muslim2+c100rat+c111rat+yeardum
+lnpop+lnarea+devcountrycomlanguage+bitcum,
data = women, index = c("country", "year"), random.method = c("swar"),
inst.method = c("bvk"), model="random")
summary(iv)

Coefficients :
Estimate Std. Error t-value  Pr(>|t|)
(Intercept)   -0.2258528  0.2951301 -0.7653 0.4441954
fdistockgdp   -0.0067207  0.0077315 -0.8693 0.3847993
trade  0.0068462  0.0023687  2.8903 0.0038863 **
polrightsreversed  0.0092366  0.0106174  0.8699 0.3844229
lnrgdpch   0.1246679  0.0389043  3.2045 0.0013724 **
execleft   0.1118046  0.0340817  3.2805 0.0010524 **
muslim2   -0.0044742  0.0012433 -3.5986 0.0003270 ***
c100rat0.0226208  0.0595134  0.3801 0.7039114
c111rat0.0165951  0.0618339  0.2684 0.7884310
yeardum19820.1479947  0.0588824  2.5134 0.0120282 *
yeardum19830.1783255  0.0606153  2.9419 0.0032958 **
yeardum19840.0344572  0.0597167  0.5770 0.5639907
yeardum19850.2206961  0.0610344  3.6159 0.0003060 ***
yeardum19860.2428015  0.0649779  3.7367 0.0001912 ***
yeardum19870.0489043  0.0615708  0.7943 0.4271186
yeardum19880.2243599  0.0605343  3.7063 0.0002155 ***
yeardum19890.2215060  0.0624042  3.5495 0.0003940 ***
yeardum19900.0688333  0.0607056  1.1339 0.2569648
yeardum19910.1370871  0.0638830  2.1459 0.0319892 *
yeardum19920.1851857  0.0630868  2.9354 0.0033655 **
yeardum19930.0904620  0.0698526  1.2950 0.1954420
yeardum19940.1003735  0.0737431  1.3611 0.1736137
yeardum19950.1164818  0.0721240  1.6150 0.1064494
yeardum19960.0482520  0.0787232  0.6129 0.5399837
yeardum19970.1049161  0.0895001  1.1722 0.2412247
yeardum19980.2191887  0.1109757  1.9751 0.0483807 *
yeardum19990.1573342  0.1397150  1.1261 0.2602422
yeardum20000.1532796  0.1627206  0.9420 0.3463059
---
 However STATA gives me

---
wecon |  Coef.   Std. Err.  zP>|z| [95% Conf.
Interval]
--+-
---
trade |   .0093915   .0027483 3.42   0.001  .004005
.014778
  fdistockgdp |  -.0169171   .0092405-1.83   0.067-.0350281
.0011938
polrightsreversed |   .0165855   .0119176 1.39   0.164-.0067726
.0399436
 lnrgdpch |   .1045675   .0431179 2.43   0.015 .0200579
.189077
 execleft |   .1373652   .0384442 3.57   0.000 .0620159
.2127145
  muslim2 |  -.0043645   .0013551-3.22   0.001-.0070205
-.0017085
  c100rat |   .0480539   .0657304 0.73   0.465-.0807752
.1768831
  c111rat |   .0170048   .0676272 0.25   0.801-.1155421
.1495516

Really would appreciate any help explaining why the results are so different

[[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] Substring and extract a certain number of characters

[arun]

Also,
library(stringr)
str_sub(s2,-5,-1)

str_extract(s1,"[[:alpha:]]+")

A.K.





On Thursday, November 14, 2013 10:50 AM, arun  wrote:
Try:
s1 <- "hello.world"
s2 <- c("GLEm0045", "GLEn0042", "GLEz0048")

 substr(s1,1,5)
#or
gsub("\\..*","",s1)

n <- 5
substr(s2,nchar(s2)-n +1,nchar(s2))


A.K.



If I have a column that has "hello.world"  how   could I extract hello? 
 In addition in a column how can I extract the last 5 characters?    For 
example from character  G1Em0017  I wish to extract m0017 .

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Substring and extract a certain number of characters

[arun]

Try:
s1 <- "hello.world"
s2 <- c("GLEm0045", "GLEn0042", "GLEz0048")

 substr(s1,1,5)
#or
gsub("\\..*","",s1)

n <- 5
substr(s2,nchar(s2)-n +1,nchar(s2))


A.K.



If I have a column that has "hello.world"  how   could I extract hello? 
 In addition in a column how can I extract the last 5 characters?    For 
example from character  G1Em0017  I wish to extract m0017 .

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Replace NA's with value in the next row

[arun]

Hi,

I think you used a column that doesn't exist in the dataset.

Targetstation <- read.table(text="V1 V2 V3 V4 V5 V6 V7  
0 0 0 1.2 0 0 0.259
0 0 12.8 0 23.7 0 8.495
 6 0 81.7 0.2 0 20 19.937
 0 1.5 60.9 0 0 15.5 13.900
 1 13 56.8 17.5 32.8 6.4 27.654
  4 3 66.4 2 0.3 NA 17.145",sep="",header=TRUE)


within(Targetstation, V6 <- replace(V6,is.na(V6),V7[is.na(V6)]))
  V1   V2   V3   V4   V5 V6 V7
1  0  0.0  0.0  1.2  0.0  0.000  0.259
2  0  0.0 12.8  0.0 23.7  0.000  8.495
3  6  0.0 81.7  0.2  0.0 20.000 19.937
4  0  1.5 60.9  0.0  0.0 15.500 13.900
5  1 13.0 56.8 17.5 32.8  6.400 27.654
6  4  3.0 66.4  2.0  0.3 17.145 17.145

#if you use:
 !is.na(Targetstation$v6) #'v6' and 'V6' are different
logical(0)
Warning message:
In is.na(Targetstation$v6) :
  is.na() applied to non-(list or vector) of type 'NULL'



A.K.









On Thursday, November 14, 2013 2:26 AM, dila radi  wrote:
Hi all,

I have a data set which treat missing value as NA and now I need to replace
all these NA's by using number in the same row but different column.

Here is the part of my data:
V1 V2 V3 V4 V5 V6 V7  0 0 0 1.2 0 0 0.259  0 0 12.8 0 23.7 0 8.495  6 0
81.7 0.2 0 20 19.937  0 1.5 60.9 0 0 15.5 13.900  1 13 56.8 17.5 32.8 6.4
27.654  4 3 66.4 2 0.3 NA 17.145


I want to replace (V6, 6) with (V7, 6). I have about 1000 NA's in V6 which
I want to replace  with the same row in V7. The other values in V6, I want
to keep remain the same.

How to achieve this? Assuming my data is called "Targetstation",  I have
tried this:

Targetstation <- within(Targetstation, v6 <- replace(v6, is.na(v6), v7[is.na
(v6)]))

But R gives me this:

Warning messages:

1: In is.na(v6) : is.na() applied to non-(list or vector) of type 'NULL'

2: In is.na(v6) : is.na() applied to non-(list or vector) of type 'NULL'


How to solve this?

Thank you in advance.

Regards,

Dila.

    [[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 aggregated data to individual data

[arun]

Hi,
Try:
 D1 <- D[rep(row.names(D),D[,3]),-3] ##assuming rownames(D) are from 1:nrow(D)
 row.names(D1) <- 1:nrow(D1)
A.K.


On Thursday, November 14, 2013 5:32 AM, peron  wrote:
Hello 



I have data in following form : 100 ind describe by two variables x and y.



D<-data .frame( x=rnorm(3), y=rnorm(3), size=c(50,10,40))



I want data for individual, i.e, 100 observations for my 100 ind.



Thank for your help



Olivier Peron


    [[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] From list to dataframe

[arun]

Hi Hermann,

You may try:
do.call(rbind,testlist)
#or
Reduce(rbind,testlist)
#or
library(plyr)
 ldply(testlist)
A.K.




On Thursday, November 14, 2013 8:27 AM, Hermann Norpois  
wrote:
Hello,

having a list like testlist I would like to transform it in dataframe. How
does it work?

Thanks
Hermann

> testlist
[[1]]
      BP_A     SNP_A     BP_B     SNP_B       R2
2 27001689 rs4822747 27002392 rs4820690 0.695642
3 27001689 rs4822747 27004298 rs5761627 0.695642
4 27001689 rs4822747 27004902 rs5752355 0.687861
5 27001689 rs4822747 27004964 rs4822748 0.682715
6 27001689 rs4822747 27005122 rs4820691 0.695642
7 27001689 rs4822747 27005158 rs4822749 0.695642

[[2]]
       BP_A     SNP_A     BP_B     SNP_B       R2
9  27002392 rs4820690 27001689 rs4822747 0.695642
11 27002392 rs4820690 27004298 rs5761627 1.00
12 27002392 rs4820690 27004902 rs5752355 0.988681
13 27002392 rs4820690 27004964 rs4822748 0.988655
14 27002392 rs4820690 27005122 rs4820691 1.00
15 27002392 rs4820690 27005158 rs4822749 1.00
16 27002392 rs4820690 27005194 rs4822750 0.988655
17 27002392 rs4820690 27005470 rs5761629 0.695642

[[3]]
       BP_A     SNP_A     BP_B     SNP_B       R2
19 27004298 rs5761627 27001689 rs4822747 0.695642
20 27004298 rs5761627 27002392 rs4820690 1.00
22 27004298 rs5761627 27004902 rs5752355 0.988681
23 27004298 rs5761627 27004964 rs4822748 0.988655
24 27004298 rs5761627 27005122 rs4820691 1.00
25 27004298 rs5761627 27005158 rs4822749 1.00

#The aim is the dataframe
       BP_A     SNP_A     BP_B     SNP_B       R2
2  27001689 rs4822747 27002392 rs4820690 0.695642
3  27001689 rs4822747 27004298 rs5761627 0.695642
4  27001689 rs4822747 27004902 rs5752355 0.687861
5  27001689 rs4822747 27004964 rs4822748 0.682715
6  27001689 rs4822747 27005122 rs4820691 0.695642
7  27001689 rs4822747 27005158 rs4822749 0.695642
9  27002392 rs4820690 27001689 rs4822747 0.695642
11 27002392 rs4820690 27004298 rs5761627 1.00
12 27002392 rs4820690 27004902 rs5752355 0.988681
13 27002392 rs4820690 27004964 rs4822748 0.988655
14 27002392 rs4820690 27005122 rs4820691 1.00
15 27002392 rs4820690 27005158 rs4822749 1.00
16 27002392 rs4820690 27005194 rs4822750 0.988655
17 27002392 rs4820690 27005470 rs5761629 0.695642
19 27004298 rs5761627 27001689 rs4822747 0.695642
20 27004298 rs5761627 27002392 rs4820690 1.00
22 27004298 rs5761627 27004902 rs5752355 0.988681
23 27004298 rs5761627 27004964 rs4822748 0.988655
24 27004298 rs5761627 27005122 rs4820691 1.00
25 27004298 rs5761627 27005158 rs4822749 1.00
>
> dput (testlist)
list(structure(list(BP_A = c(27001689L, 27001689L, 27001689L,
27001689L, 27001689L, 27001689L), SNP_A = c("rs4822747", "rs4822747",
"rs4822747", "rs4822747", "rs4822747", "rs4822747"), BP_B = c(27002392L,
27004298L, 27004902L, 27004964L, 27005122L, 27005158L), SNP_B =
c("rs4820690",
"rs5761627", "rs5752355", "rs4822748", "rs4820691", "rs4822749"
), R2 = c(0.695642, 0.695642, 0.687861, 0.682715, 0.695642, 0.695642
)), .Names = c("BP_A", "SNP_A", "BP_B", "SNP_B", "R2"), row.names = 2:7,
class = "data.frame"),
    structure(list(BP_A = c(27002392L, 27002392L, 27002392L,
    27002392L, 27002392L, 27002392L, 27002392L, 27002392L), SNP_A =
c("rs4820690",
    "rs4820690", "rs4820690", "rs4820690", "rs4820690", "rs4820690",
    "rs4820690", "rs4820690"), BP_B = c(27001689L, 27004298L,
    27004902L, 27004964L, 27005122L, 27005158L, 27005194L, 27005470L
    ), SNP_B = c("rs4822747", "rs5761627", "rs5752355", "rs4822748",
    "rs4820691", "rs4822749", "rs4822750", "rs5761629"), R2 = c(0.695642,
    1, 0.988681, 0.988655, 1, 1, 0.988655, 0.695642)), .Names = c("BP_A",
    "SNP_A", "BP_B", "SNP_B", "R2"), row.names = c(9L, 11L, 12L,
    13L, 14L, 15L, 16L, 17L), class = "data.frame"), structure(list(
        BP_A = c(27004298L, 27004298L, 27004298L, 27004298L,
        27004298L, 27004298L), SNP_A = c("rs5761627", "rs5761627",
        "rs5761627", "rs5761627", "rs5761627", "rs5761627"),
        BP_B = c(27001689L, 27002392L, 27004902L, 27004964L,
        27005122L, 27005158L), SNP_B = c("rs4822747", "rs4820690",
        "rs5752355", "rs4822748", "rs4820691", "rs4822749"),
        R2 = c(0.695642, 1, 0.988681, 0.988655, 1, 1)), .Names = c("BP_A",
    "SNP_A", "BP_B", "SNP_B", "R2"), row.names = c(19L, 20L,
    22L, 23L, 24L, 25L), class = "data.frame"))

    [[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] lapply?

[Toth, Denes]

Hi,

the output of lapply() is a list; see ?lapply and ?sapply.

# if you know the length of your list in advance,
# this definition is better:
uu <- vector("list", 2)

# list elements
uu[[1]] <- c(1,2,3)
uu[[2]] <- c(3,4,5)


# some options to achieve what you want:
matrix(unlist(uu), 2, 3, T)
do.call(rbind, uu)
t(sapply(uu, I))


HTH,
  Denes


> Hi,
>
> I was trying to use lapply to create a matrix from a list:
>
> uu <- list()
> uu[[1]] <- c(1,2,3)
> uu[[2]] <- c(3,4,5)
>
> The output I desire is a matrix with 2 rows and 3 columns, so I try:
>
> xx <- lapply(uu,rbind)
>
> Obviously, I'm not doing something right, but what!?
>
>   [[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] volume of ellipsoid

[yuanzhi]

Carl Witthoft wrote
> Well, given that an upvoted question on math.stackexchange got no answers,
> I'd say you're asking a very difficult question.   Perhaps this paper  
> http://www.geometrictools.com/Documentation/IntersectionOfEllipsoids.pdf
> 
> will be of some help.   It's possible that you could do:
> 1) find the ellipse of intersection of the two ellipsoids
> 2) find some formula for the equivalent of the "spherical cap" volume
> 3) apply that formula to each ellipsoid, using the intersection ellipse as
> the base of the caps.
> 
> All this, btw, is math,  not really related to R or to software
> programming in general.
> yuanzhi wrote
>> Hi, everyone!
>> 
>> I have a matrix X(n*p) which is n samples from a p-dimensional normal  
>> distribution. Now I want to make the ellipsoid containing (1-α)% of  
>> the probability in the distribution based on Mahalanobis distance:
>> 
>> μ:(x-μ)'Σ^(-1)(x-μ)≤χ2p(α)
>> 
>> where x and Σ is the mean and variance-covariance matrix of X. Are  
>> there some functions in R which can plot the ellipsoid and calculate  
>> the volume (area for p=2, hypervolume for p>3) of the ellipsoid? I  
>> only find the "ellipsoidhull" function in "cluster" package, but it  
>> deal with ellipsoid hull containing X, which obvious not the one I want.
>> 
>> After that, a more difficult is the following. If we have a series of  
>> matrix X1, X2, X3,… Xm which follow different p-dimensional normal  
>> distributions. And we make m ellipsoids (E1, E2, … Em) for each matrix  
>> like the before. How can we calculate the volume of union of the m  
>> ellipsoids? One possible solution for this problem is the  
>> inclusion-exclusion principle:
>> 
>> V(⋃Ei)(1≤i≤m)=
>> V(Ei)(1≤i≤m)-∑V(Ei⋂Ej)(1≤i> 
>> But the problem is how to calculate the volume of intersection between  
>> 2, 3 or more ellipsoids. Are there some functions which can calculate  
>> the volume of intersection between two region or functions which  
>> directly calculate the volume of a union of two region(the region here  
>> is ellipsoid). OR yo you have any good ideas solving this problem in  
>> R? Thank you all in advance!
>> 
>> Yuanzhi
>> 
>> __
>> R-help@

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

Hi, Carl Witthoft 

yes, it looks like a mathematical question. I will try based on your
suggestion to calculate the volume of the intersection. But I still want to
know whether there are some functions in R which can calculate the volume of
an ellipsoid(area for p=2, hypervolume for p>3) containing X, just like the
"convhulln" function in "geometry" package which can calculate the volume of
convex hull containing X.
thanks!



--
View this message in context: 
http://r.789695.n4.nabble.com/volume-of-ellipsoid-tp4680409p4680446.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] volume of ellipsoid

[yuanzhi]

Hi, Carl Witthoft

yes, it looks like a mathematical question. I will try based on your
suggestion to calculate the volume of the intersection. But I still want to
know whether there are some functions in R which can calculate the volume of
an ellipsoid(area for p=2, hypervolume for p>3) containing X, just like the
"convhulln" function in "geometry" package which can calculate the volume of
convex hull containing X. 



--
View this message in context: 
http://r.789695.n4.nabble.com/volume-of-ellipsoid-tp4680409p4680445.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] Survival analysis with truncated data

[Nicolas Palix]

Hi,

Thanks for your response.

On Thu, Nov 14, 2013 at 4:44 PM, Terry Therneau  wrote:
> I think that your data is censored, not truncated.
>   For a fault introduced 1/2005 and erased 2/2006, duration = 13 months
>   For a fault introduced 4/2010 and still in existence at the last
> observation 12/2010, duration> 8 months.
>   For a fault introduced before 2004, erased  3/2005, in a machine installed
> 2/1998, the duration is somewhere between 15 and 87 months.
>   For a fault introduced before 2004, smachine installed 5/2000, still
> present 11/2010 at last check, the duration is > 126 months.
>
> For type=interval2 the data would be (13,13), (8,NA), (15,87), (126, NA).

I have done this that way. My problem is that I have no information
when a fault is introduced before 2004. Indeed, this is about the lifespan
of software faults in the code.

In your example, this means I could not set the upper bound to 87 months.
As I know for sure that the first software release was in 1994. For a fault
which is observed from 2004 up to 2005 I set the range to (12, 120+12). That is
12 observed + 10 years from 1994 to 2004. The estimation is almost
similar if I use (12, NA) and gives me an upper bound.
I tried (12, 12) to have the lower bound.


I tried with 5 years instead of 10. This seems to give an
over-estimation too.

Could I use some properties of the data from ]2004;2010] to give
an average extension to these faults ? The average or median
for instance.

Thanks in advance.


>
> Terry T.
>
>
> On 11/14/2013 05:00 AM, r-help-requ...@r-project.org wrote:
>>
>> Hi,
>>
>> I would like to know how to handle truncated data.
>> My intend is to have the survival curve of a software fault in order
>> to have some information
>> about fault lifespan.
>>
>> I have some observations of a software system between 2004 and 2010.
>> The system was first released in 1994.
>> The event considered is the disappearance of a software fault. The
>> faults can have been
>> introduced at any time, between 1994 and 2010. But for fault
>> introduced before 2004,
>> there is not mean to know their age.
>>
>> I used the Surv and survfit functions with type interval2.
>> For the faults that are first observed in 2004, I set the lower bound
>> to the lifespan
>> observed between 2004 and 2010.
>>
>> How could I set the upper bound ? Using 1994 as a starting point to not
>> seems
>> to be meaningful. Neither is using only the lower bound.
>>
>> Should I consider another survival estimator ?
>>
>> Thanks in advance.



-- 
Nicolas Palix
Tel: +33 4 76 51 46 27
http://membres-liglab.imag.fr/palix/

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] making a barplot with table of experimental conditions underneath (preferably ggplot2)

[John Kane]

Hi Nina, 
I think the following code does what you want (thanks to Jim Lemon for showing 
me what you wanted in terms of x-axis tick labels)

However I  think that barcharts are generally evil  so I changed your geom_bar 
to geom_point.  Feel free to change it back if your discipline requires it but 
I think it shows your data better

Also you were doing some unnecessary extraction of dat from the data.frame so I 
just included a "data =" statement in qplot and changed the variable names in 
it to the original df names. This makes for cleaner and more readable code.

df <- data.frame (experiment=c("E1","E2","E3","E4"), mean = c(3,4,5,6),
stdev=c(0.1,0.1,0.05,0.2), method = c("STD","STD", "FP", "FP"), enzyme =c
("T","T/L","T","T/L"), denaturation=c("U","U","0.05%RG", "0.1%RG"))

df$labs  <-  paste(df[,4],"\n ",df[,5], "\n ",df[,6]) # create labels

df.plot <- qplot(experiment,mean,data = df, xlab="", ylab="# peptides 
identified")+
  geom_point(fill="grey")+
  geom_errorbar(aes(x=experiment, ymin=mean-stdev, ymax=mean+stdev), width=0.25)

p + scale_x_discrete( labels = df$labs)

I hop
John Kane
Kingston ON Canada


> -Original Message-
> From: n.hub...@ncmls.ru.nl
> Sent: Wed, 13 Nov 2013 14:24:28 +
> To: r-help@r-project.org
> Subject: [R] making a barplot with table of experimental conditions
> underneath (preferably ggplot2)
> 
> Dear all,
>
> my data looks the following:
> 
> df <- data.frame (experiment=c("E1","E2","E3","E4"), mean = c(3,4,5,6),
> stdev=c(0.1,0.1,0.05,0.2), method = c("STD","STD", "FP", "FP"), enzyme =c
> ("T","T/L","T","T/L"), denaturation=c("U","U","0.05%RG", "0.1%RG"))

> I would like to make a bar plot with standard deviation which I solved
> the following way:
> 
> x <- df$experiment
> y <- df$mean
> sd <- df$stdev
> 
> df.plot <- qplot(x,y,xlab="", ylab="# peptides identified")+
>   geom_bar(colour="black", fill="darkgrey")+
>   geom_errorbar(aes(x=x, ymin=y-sd, ymax=y+sd), width=0.25)
> 
> df.plot
> 

> However, as the labels for the x-axis (the bars) I do not want the
> experiment number, as now, but instead a table containing the other
> columns of my data.frame (method, enzyme, denaturation) with the
> description in the front and the certain 'value' below the bars.
> 
> I am looking forward to your suggestions!
>
> With best wishes,
> 
> Nina
> __
> 
> Dr. Nina C. Hubner
> scientist quantitative proteomics
> 
> Dep. of Molecular Biology, Radboud University Nijmegen, The Netherlands
> e-mail: n.hub...@ncmls.ru.nl
> tel: +31-24-3613655
> 
> Visiting address:
> NCMLS, Dept of Molecular Biology (route 274)
> Geert Grooteplein 26/28
> 6525 GA Nijmegen
> The Netherlands
> 
> 
> 
> 
> The Radboud University Medical Centre is listed in the Commercial
> Register of the Chamber of Commerce under file number 41055629.
> 
> 
>   [[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.


FREE 3D EARTH SCREENSAVER - Watch the Earth right on your desktop!

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

[Simon Zehnder]

Thanks Martin, for making this clear to me, I thought of Pearson’s 
Method-of-Moments.
On 14 Nov 2013, at 16:08, Martin Maechler  wrote:

>> "SZ" == Simon Zehnder 
>>on Thu, 14 Nov 2013 08:52:16 +0100 writes:
> 
>SZ> Check the gmm package with a weighting matrix equal to
>SZ> the identity.  Best
> 
>SZ> Simon
> 
> No."MM" is ambigous: The  gmm package is about
>   "Generalized Method of Moments"
> 
> but Izhak is really asking about MM estimation in the field of
> robust statistics, where an MM estimator is a special kind of 
> "M-Estimator" (and these, introduced by Peter Huber (1964), Annals,
>  fortunately are well known unambigously,
>  as a generalization of ML estimators (= MLE)),
> namely an M-estimator with redescending psi (==> non-convex
> problem ==> solution typically depends on starting value),
> started with a high-breakdown point initial estimate.
> 
> Izhak mentioned  nlrob() from package robustbase, and he is
> right that this does not provide an MM estimator, but only an M
> estimate started by Least Squares.
> 
> Fortunately, we have had contributions to the robustbase
> package from Eduardo Conceição since last summer.
> He provided alternative versions of  nlrob(),
> notably a  nlrob.MM()  
> currently *hidden* {and undocumented}.
> 
> I.e. you need
> 
> install.packages("robustbase", repos="http://R-Forge.R-project.org";)
> 
> library(robustbase)
> 
> and then use  robustbase:::nlrob.MM()
> 
> The reason for all this is that I plan to have one nlrob()
> function  but with  a new argument  'method'
> and so eventually
> 
>nlrob.MM(, method="MM")
> 
> will correspond to current to the R-forge version
> 
>robustbase:::nlrob.MM()
> 
> ---
> Please for all more, use the dedicated mailing list
> R-SIG-robust only (i.e. do *not* just reply-all to this
> e-mail!).
> 
> Martin Maechler,
> ETH Zurich
> 
> 
>SZ> On 14 Nov 2013, at 04:37, IZHAK shabsogh
>SZ>  wrote:
> 
>>> hi
>>> 
>>> I have a nonlinear regression model with two parameter, i
>>> have estimated using nls in R and i want to also find the
>>> estimate using MM, someone refer me to this function
>>> nlrob but this is main for only M estimate. can you
>>> please help me findout a function in R for MM nonlinear
>>> regression
>>> 
>>> 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.
> 
>SZ> __
>SZ> R-help@r-project.org mailing list
>SZ> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do
>SZ> read the posting guide
>SZ> http://www.R-project.org/posting-guide.html and provide
>SZ> 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] lapply?

[Brian Smith]

Thanks all! So many ways


On Thu, Nov 14, 2013 at 10:35 AM, Rui Barradas  wrote:

> Hello,
>
> You are applying rbind to each element of the list, not rbinding it with
> the others. Try instead
>
> do.call(rbind, uu)
>
> Hope this helps,
>
> Rui Barradas
>
> Em 14-11-2013 15:20, Brian Smith escreveu:
>
>> Hi,
>>
>> I was trying to use lapply to create a matrix from a list:
>>
>> uu <- list()
>> uu[[1]] <- c(1,2,3)
>> uu[[2]] <- c(3,4,5)
>>
>> The output I desire is a matrix with 2 rows and 3 columns, so I try:
>>
>> xx <- lapply(uu,rbind)
>>
>> Obviously, I'm not doing something right, but what!?
>>
>> [[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.

Re: [R] Generating bootstrap samples from a panel data frame

[Dereje Fentie]

With a sample data frame:
d = data.frame(id = rep(c(1,2,3,4),5), x = round(rexp(20), digits=2), y =
round(runif(20), digits=2))

I would like to generate 100 bootstrap data with replacement and save each
bootstrap data as b1, b2, ..., b100. I attempted the sample function but
could not make it work



On Mon, Nov 11, 2013 at 6:11 PM, Dereje Fentie  wrote:

>
>   With a data frame (call it *d*) composed of 2000 individuals and 
> *n*observations for each individual (thus
> *2000n* observations in total), I would like to generate *k* bootstrap
> samples with replacement from *d*. Amongst other variables, *d* has a
> numeric variable *id* taking on identical value for observations
> belonging to the same individual.
>
> Taking into consideration the panel nature of the data, I want to generate
> many bootstrap samples with replacement and store each bootstrap sample
> data frame for further use. Sampling (or selection into the bootstrap
> sample) shall be based on individuals (on unique values of *id*) such
> that if an individual is in a particular bootstrap sample, so will all
> observations belonging to that individual.
>
> How can I do this in r?
>

[[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] Survival analysis with truncated data

[Terry Therneau]

I think that your data is censored, not truncated.
  For a fault introduced 1/2005 and erased 2/2006, duration = 13 months
  For a fault introduced 4/2010 and still in existence at the last observation 12/2010, 
duration> 8 months.
  For a fault introduced before 2004, erased  3/2005, in a machine installed 2/1998, the 
duration is somewhere between 15 and 87 months.
  For a fault introduced before 2004, smachine installed 5/2000, still present 11/2010 at 
last check, the duration is > 126 months.


For type=interval2 the data would be (13,13), (8,NA), (15,87), (126, NA).

Terry T.


On 11/14/2013 05:00 AM, r-help-requ...@r-project.org wrote:

Hi,

I would like to know how to handle truncated data.
My intend is to have the survival curve of a software fault in order
to have some information
about fault lifespan.

I have some observations of a software system between 2004 and 2010.
The system was first released in 1994.
The event considered is the disappearance of a software fault. The
faults can have been
introduced at any time, between 1994 and 2010. But for fault
introduced before 2004,
there is not mean to know their age.

I used the Surv and survfit functions with type interval2.
For the faults that are first observed in 2004, I set the lower bound
to the lifespan
observed between 2004 and 2010.

How could I set the upper bound ? Using 1994 as a starting point to not seems
to be meaningful. Neither is using only the lower bound.

Should I consider another survival estimator ?

Thanks in advance.


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

Re: [R] lapply?

[Rui Barradas]

Hello,

You are applying rbind to each element of the list, not rbinding it with 
the others. Try instead


do.call(rbind, uu)

Hope this helps,

Rui Barradas

Em 14-11-2013 15:20, Brian Smith escreveu:

Hi,

I was trying to use lapply to create a matrix from a list:

uu <- list()
uu[[1]] <- c(1,2,3)
uu[[2]] <- c(3,4,5)

The output I desire is a matrix with 2 rows and 3 columns, so I try:

xx <- lapply(uu,rbind)

Obviously, I'm not doing something right, but what!?

[[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] lapply?

[Berend Hasselman]

On 14-11-2013, at 16:20, Brian Smith  wrote:

> Hi,
> 
> I was trying to use lapply to create a matrix from a list:
> 
> uu <- list()
> uu[[1]] <- c(1,2,3)
> uu[[2]] <- c(3,4,5)
> 
> The output I desire is a matrix with 2 rows and 3 columns, so I try:
> 
> xx <- lapply(uu,rbind)
> 
> Obviously, I'm not doing something right, but what!?


do.call(rbind,uu) 

Berend

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] issues with calling predict.coxph.penal (survival) inside a function

[Terry Therneau]

Thanks for the reproducable example.  I can confirm that it fails on my machine using 
survival 2-37.5, the next soon-to-be-released version,


The issue is with NextMethod, and my assumption that the called routine inherited 
everything from the parent, including the environment chain.  A simple test this AM showed 
me that the assumption is false.  It might have been true for Splus.  Working this out may 
take some time -- every other one of my wrestling matches with predict inside a function 
has -- and there is a reasonable chance that it won't make this already overdue release.


In the meantime, here is a workaround that I have sometimes used in other situations. 
Inside your function do the following: fit a new coxph model with fixed coefficients, and 
do prediction on that.


myfun <- function(oldfit, subset) {
   newX <- model.matrix(oldfit)[subset,]
   newY <- oldfit$y[subset]
   newfit <- coxph(newY ~ newX, iter=0, init=coef(oldfit))
   newfit$var <- oldfit$var
   predict(newfit)
   }

If the subset is all of a particular strata, as you indicated, then all of the predictions 
will be correct.  If not, then those that make use of the the baseline hazard (type= 
expect) will be incorrect but all others are ok.


Terry Therneau


On 11/14/2013 05:00 AM, r-help-requ...@r-project.org wrote:

Hello everyone,



I got an issue with calling predict.coxph.penal inside a function.



Regarding the context: My original problem is that I wrote a function that
uses predict.coxph and survfit(model) to predict

a lot of survival-curves using only the basis-curves for the strata (as
delivered by survfit(model) ) and then adapts them with

the predicted risk-scores. Because there are cases where my new data has
strata which didn't exist in the original model I exclude

them, using a Boolean vector inside the function.

I end up with a call like this: predict (coxph_model,
newdata[subscript_vector,] )



This works fine for coxph.model, but when I fit a model with a spline
(class coxph.penal), I get an error:

"Error in `[.data.frame`(newdata, [subscript_vector, ) : object
'[subscript_vector ' not found"



I suppose this is because of NextMethod, but I am not sure how to work
around it. I also read a little bit about all those
matching-and-frame-issues,

But must confess I am not really into it.



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

[Brian Smith]

Hi,

I was trying to use lapply to create a matrix from a list:

uu <- list()
uu[[1]] <- c(1,2,3)
uu[[2]] <- c(3,4,5)

The output I desire is a matrix with 2 rows and 3 columns, so I try:

xx <- lapply(uu,rbind)

Obviously, I'm not doing something right, but what!?

[[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] MM estimator

[Martin Maechler]

> "SZ" == Simon Zehnder 
> on Thu, 14 Nov 2013 08:52:16 +0100 writes:

SZ> Check the gmm package with a weighting matrix equal to
SZ> the identity.  Best

SZ> Simon

No."MM" is ambigous: The  gmm package is about
   "Generalized Method of Moments"

but Izhak is really asking about MM estimation in the field of
robust statistics, where an MM estimator is a special kind of 
"M-Estimator" (and these, introduced by Peter Huber (1964), Annals,
   fortunately are well known unambigously,
   as a generalization of ML estimators (= MLE)),
namely an M-estimator with redescending psi (==> non-convex
problem ==> solution typically depends on starting value),
started with a high-breakdown point initial estimate.

Izhak mentioned  nlrob() from package robustbase, and he is
right that this does not provide an MM estimator, but only an M
estimate started by Least Squares.

Fortunately, we have had contributions to the robustbase
package from Eduardo Conceição since last summer.
He provided alternative versions of  nlrob(),
notably a  nlrob.MM()  
currently *hidden* {and undocumented}.

I.e. you need

 install.packages("robustbase", repos="http://R-Forge.R-project.org";)

 library(robustbase)

and then use  robustbase:::nlrob.MM()

The reason for all this is that I plan to have one nlrob()
function  but with  a new argument  'method'
and so eventually

nlrob.MM(, method="MM")

will correspond to current to the R-forge version

robustbase:::nlrob.MM()

---
Please for all more, use the dedicated mailing list
R-SIG-robust only (i.e. do *not* just reply-all to this
e-mail!).

Martin Maechler,
ETH Zurich

 
SZ> On 14 Nov 2013, at 04:37, IZHAK shabsogh
SZ>  wrote:

>> hi
>> 
>> I have a nonlinear regression model with two parameter, i
>> have estimated using nls in R and i want to also find the
>> estimate using MM, someone refer me to this function
>> nlrob but this is main for only M estimate. can you
>> please help me findout a function in R for MM nonlinear
>> regression
>> 
>> 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.

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

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

Re: [R] R Beginner - Need Perhaps 5 - 10 Minutes of R User Time to Learn Few Basics

[John Kane]

This is not an R question per se.  It really seems like an RStudio question.  
They have their own help forum and it is probably best to ask there.  

My first thought was that you had just closed the output window in RStudio but 
in my version, 0.97.449 under Ubuntu 13.10 the window automatically opens when 
you send a command to R.

You have not said what operating system you are using. If in Windows, IIRC, you 
can just open an RGui or RTerminal and work there.  Just open the script in a 
text editor and paste it (preferably line-by-line) and see what happens.  
Actually just type in something like 1 + 1 and see what happens

In Linux open a terminal window, start R and do the same. No idea about Macs

This way you will get an idea if the R installation or the RStudio installation 
is likely the culprit.

It sounds like a faulty RStudio installation to me but who knows. An RStudio 
session should look more or less like this http://www.rstudio.com/ide/ . If you 
don't have those four panels then you have an RStudio problem.

John Kane
Kingston ON Canada


> -Original Message-
> From: zfeinst...@isgmn.com
> Sent: Wed, 13 Nov 2013 13:57:17 +
> To: r-help@r-project.org
> Subject: [R] R Beginner - Need Perhaps 5 - 10 Minutes of R User Time to
> Learn Few Basics
> 
> I have finally decided that I will learn R and learn it very well. For
> now I am using a program that a friend of mine developed to do some
> advanced statistical analyses. I downloaded RStudio to my machine.
> [Perhaps RStudio is not the best platform to work from -  I have heard
> that Rattle is sort of the new standard.] I have so far been able to
> highlight the rows of the code that I wish to run, but then I somehow
> turned off seeing the output. I also cannot find where I would locate the
> output window. Yes, frustrated.
> 
> Would any kind soul be interested in helping kickstart my R learning? I
> have JoinMe installed on my machine so I figure we can do it
> interactively. It should not take more than a few minutes. I am already
> very experienced with both C and VBA languages as well as SPSS syntax so
> there is not much need to worry about me being too much of a novice.
> 
> Thank you very much in advance.
> 
> Zach Feinstein
> zfeinst...@isgmn.com
> (952) 277-0162
> (612) 590-4813 (mobile)
> 
> 
>   [[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.


FREE 3D EARTH SCREENSAVER - Watch the Earth right on your desktop!

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Replace NA's with value in the next row

[Bert Gunter]

On Thu, Nov 14, 2013 at 1:13 AM, Frans Marcelissen
 wrote:
> R purists forbid the use of the for loop,

That is utter nonsense. Please do not make such statements when you
have no idea what you're talking about. It just promulgates confusion.

-- Bert

> simple solution:
>
> for (i in 1:(length(V6)-1)) if(is.na(V6[i])) V6[i]<-V6[i+1]
>
>
>
>
> 2013/11/14 Jim Lemon 
>
>> On 11/14/2013 04:02 PM, dila radi wrote:
>>
>>> Hi all,
>>>
>>> I have a data set which treat missing value as NA and now I need to
>>> replace
>>> all these NA's by using number in the same row but different column.
>>>
>>> Here is the part of my data:
>>>   V1 V2 V3 V4 V5 V6 V7  0 0 0 1.2 0 0 0.259  0 0 12.8 0 23.7 0 8.495  6 0
>>> 81.7 0.2 0 20 19.937  0 1.5 60.9 0 0 15.5 13.900  1 13 56.8 17.5 32.8 6.4
>>> 27.654  4 3 66.4 2 0.3 NA 17.145
>>>
>>>
>>> I want to replace (V6, 6) with (V7, 6). I have about 1000 NA's in V6 which
>>> I want to replace  with the same row in V7. The other values in V6, I want
>>> to keep remain the same.
>>>
>>> How to achieve this? Assuming my data is called "Targetstation",  I have
>>> tried this:
>>>
>>> Targetstation<- within(Targetstation, v6<- replace(v6, is.na(v6), v7[
>>> is.na
>>> (v6)]))
>>>
>>> But R gives me this:
>>>
>>> Warning messages:
>>>
>>> 1: In is.na(v6) : is.na() applied to non-(list or vector) of type 'NULL'
>>>
>>> 2: In is.na(v6) : is.na() applied to non-(list or vector) of type 'NULL'
>>>
>>>
>>>  Hi Dila,
>> You could do this like this:
>>
>> V6NA<-is.na(Targetstation$V6)
>> Targetstation$V6[V6NA]<-Targetstation$V7[V6NA]
>>
>> but if you want to use the above, I think you will have to replace the
>> is.na(V6) with is.na(Targetstation$V6) or use the method above.
>>
>> Jim
>>
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/
>> posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
>
> --
>
>
> ---
> dr F.H.G. (Frans) Marcelissen
> DigiPsy (www.DigiPsy.nl )
> Pomperschans 26
> 5595 AV Leende
> tel: 040 2065030/06 2325 06 53
> skype adres: frans.marcelissen
> email: frans.marcelis...@digipsy.nl
>
> [[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

(650) 467-7374

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

[Joshua Ulrich]

On Thu, Nov 14, 2013 at 2:26 AM, Tstudent  wrote:
> I have the following two xts object:
>
> https://dl.dropboxusercontent.com/u/102669/series1.rdata
>
> https://dl.dropboxusercontent.com/u/102669/series2.rdata
>
> With str i see that they are both xts objects
>
> I can't understand why it's impossible to compare each element.
>
> For example: series1 > series2
> Why i don't get an xts object with a sequence of true and false?
>
Because series1 has a POSIXct index and series2 has a yearmon index.
The two objects have no index values in common, so there's nothing to
compare.

Convert series1's index to yearmon, and the comparison works.

> index(series1) <- as.yearmon(index(series1))
> tail(series1 > series2)
   GSPC
2013-06-01 TRUE
2013-07-01 TRUE
2013-08-01 TRUE
2013-09-01 TRUE
2013-10-01 TRUE
2013-11-01 TRUE

Best,
--
Joshua Ulrich  |  about.me/joshuaulrich
FOSS Trading  |  www.fosstrading.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] How to sum a function over a specific range in R?

[Adams, Jean]

You should cc r-help on all correspondence, so others can follow the thread.

Here is a very simple example of a home-made smooth function.  Perhaps you
can modify it to suit your needs.

Jean


# simple smooth function, using a weighted mean
smoothfunction <- function(allxs, allys, centerindex, halfwidth) {
 distfromcenter <- abs(allxs - allxs[centerindex])
weight <- ifelse(distfromcenter <= halfwidth, 1/(distfromcenter+1), 0)
 weighted.mean(allys, weight)
}

# fake data
x <- runif(1000)
y <- 3*x^2 + x + rnorm(1000)

# set desired halfwidth of window for smooth
halfwid <- 0.1

# determine the smoothed value of y for each value of x
L <- length(x)
smoothedy <- numeric(L)
for(i in 1:L) {
smoothedy[i] <- smoothfunction(x, y, i, halfwid)
}

# plot the results
plot(x, y)
lines(x[order(x)], smoothedy[order(x)], lwd=3)



On Wed, Nov 13, 2013 at 5:59 PM, umair durrani wrote:

> *Thanks Jeanat-least someone replied. I have gone through the link you
> provided but  the real problem is that it is more complex to understand
> than the R documentation. Actually, I don't have any background in noise
> reduction / smoothing of data. Can you guide me how I could just apply the
> equation included in my question?*
>
> *best regards,*
> *Umair Durrani*
> *email: umairdurr...@outlook.com *
>
>
> --
> From: jvad...@usgs.gov
> Date: Wed, 13 Nov 2013 15:10:07 -0600
> Subject: Re: [R] How to sum a function over a specific range in R?
> To: umairdurr...@outlook.com
> CC: r-help@r-project.org
>
>
>
> On Tue, Nov 12, 2013 at 11:45 AM, umair durrani 
> wrote:
>
> I am new to R and have already posted this question on stack overflow. The
> problem is that I did not understand the answers as the R documentation
> about the discussed functions (e.g. 'convolve') is quite complicated for a
> newbie like me. Here's the question:
> I have a big text file with more than 3 million rows. The following is the
> example of the three columns I want to use:
> indxvehID   LocalY
> 1   2   35.381
> 2   2   39.381
> 3   2   43.381
> 4   2   47.38
> 5   2   51.381
> 6   2   55.381
> 7   2   59.381
> 8   2   63.379
> 9   2   67.383
> 10  2   71.398
> where,indx = IndexvehID = Vehicle ID (Here only '2' is shown but infact
> there are 2169 vehicle IDs and each one repeats several times because the
> data was collected at every 0.1 seconds)LocalY = The y coordinate of the
> vehicle at a particular time (The time column is not shown here)
> What I want to do is to create a new column of 'SmoothedY' using the
> following formula:
> SmoothedY = 1/Z * Summation from (i-15) to (i+15) (LocalY *
> exp(-abs(i-k))/5))
> where,i = indxZ = Summation from (k =i-15) to (k = i+15) (
> exp(-abs(i-k))/5))
> How can I apply this formula to create the new column 'SmoothedY'? This is
> actually a data smoothing problem but default smoothing algorithms in R are
> not suitable for my data and I have to use this custom formula.
> Thanks in advance.
>
> Umair Durrani
>
>
> I have never tried this myself, but it appears as if you can define your
> own smoothing function using Simon Wood's mgcv package.  Check out
> http://www.maths.bath.ac.uk/~sw283/talks/snw-R-talk.pdf for more
> information.
>
> Jean
>

[[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] From list to dataframe

[andrija djurovic]

Hi. Here is one way:
l <- list(structure(list(BP_A = c(27001689L, 27001689L, 27001689L,
27001689L, 27001689L, 27001689L), SNP_A = c("rs4822747", "rs4822747",
"rs4822747", "rs4822747", "rs4822747", "rs4822747"), BP_B = c(27002392L,
27004298L, 27004902L, 27004964L, 27005122L, 27005158L), SNP_B =
c("rs4820690",
"rs5761627", "rs5752355", "rs4822748", "rs4820691", "rs4822749"
), R2 = c(0.695642, 0.695642, 0.687861, 0.682715, 0.695642, 0.695642
)), .Names = c("BP_A", "SNP_A", "BP_B", "SNP_B", "R2"), row.names = 2:7,
class = "data.frame"),
structure(list(BP_A = c(27002392L, 27002392L, 27002392L,
27002392L, 27002392L, 27002392L, 27002392L, 27002392L), SNP_A =
c("rs4820690",
"rs4820690", "rs4820690", "rs4820690", "rs4820690", "rs4820690",
"rs4820690", "rs4820690"), BP_B = c(27001689L, 27004298L,
27004902L, 27004964L, 27005122L, 27005158L, 27005194L, 27005470L
), SNP_B = c("rs4822747", "rs5761627", "rs5752355", "rs4822748",
"rs4820691", "rs4822749", "rs4822750", "rs5761629"), R2 = c(0.695642,
1, 0.988681, 0.988655, 1, 1, 0.988655, 0.695642)), .Names = c("BP_A",
"SNP_A", "BP_B", "SNP_B", "R2"), row.names = c(9L, 11L, 12L,
13L, 14L, 15L, 16L, 17L), class = "data.frame"), structure(list(
BP_A = c(27004298L, 27004298L, 27004298L, 27004298L,
27004298L, 27004298L), SNP_A = c("rs5761627", "rs5761627",
"rs5761627", "rs5761627", "rs5761627", "rs5761627"),
BP_B = c(27001689L, 27002392L, 27004902L, 27004964L,
27005122L, 27005158L), SNP_B = c("rs4822747", "rs4820690",
"rs5752355", "rs4822748", "rs4820691", "rs4822749"),
R2 = c(0.695642, 1, 0.988681, 0.988655, 1, 1)), .Names = c("BP_A",
"SNP_A", "BP_B", "SNP_B", "R2"), row.names = c(19L, 20L,
22L, 23L, 24L, 25L), class = "data.frame"))

do.call("rbind", l)


Andrija


On Thu, Nov 14, 2013 at 9:41 AM, Hermann Norpois  wrote:

> Hello,
>
> having a list like testlist I would like to transform it in dataframe. How
> does it work?
>
> Thanks
> Hermann
>
> > testlist
> [[1]]
>   BP_A SNP_A BP_B SNP_B   R2
> 2 27001689 rs4822747 27002392 rs4820690 0.695642
> 3 27001689 rs4822747 27004298 rs5761627 0.695642
> 4 27001689 rs4822747 27004902 rs5752355 0.687861
> 5 27001689 rs4822747 27004964 rs4822748 0.682715
> 6 27001689 rs4822747 27005122 rs4820691 0.695642
> 7 27001689 rs4822747 27005158 rs4822749 0.695642
>
> [[2]]
>BP_A SNP_A BP_B SNP_B   R2
> 9  27002392 rs4820690 27001689 rs4822747 0.695642
> 11 27002392 rs4820690 27004298 rs5761627 1.00
> 12 27002392 rs4820690 27004902 rs5752355 0.988681
> 13 27002392 rs4820690 27004964 rs4822748 0.988655
> 14 27002392 rs4820690 27005122 rs4820691 1.00
> 15 27002392 rs4820690 27005158 rs4822749 1.00
> 16 27002392 rs4820690 27005194 rs4822750 0.988655
> 17 27002392 rs4820690 27005470 rs5761629 0.695642
>
> [[3]]
>BP_A SNP_A BP_B SNP_B   R2
> 19 27004298 rs5761627 27001689 rs4822747 0.695642
> 20 27004298 rs5761627 27002392 rs4820690 1.00
> 22 27004298 rs5761627 27004902 rs5752355 0.988681
> 23 27004298 rs5761627 27004964 rs4822748 0.988655
> 24 27004298 rs5761627 27005122 rs4820691 1.00
> 25 27004298 rs5761627 27005158 rs4822749 1.00
>
> #The aim is the dataframe
>BP_A SNP_A BP_B SNP_B   R2
> 2  27001689 rs4822747 27002392 rs4820690 0.695642
> 3  27001689 rs4822747 27004298 rs5761627 0.695642
> 4  27001689 rs4822747 27004902 rs5752355 0.687861
> 5  27001689 rs4822747 27004964 rs4822748 0.682715
> 6  27001689 rs4822747 27005122 rs4820691 0.695642
> 7  27001689 rs4822747 27005158 rs4822749 0.695642
> 9  27002392 rs4820690 27001689 rs4822747 0.695642
> 11 27002392 rs4820690 27004298 rs5761627 1.00
> 12 27002392 rs4820690 27004902 rs5752355 0.988681
> 13 27002392 rs4820690 27004964 rs4822748 0.988655
> 14 27002392 rs4820690 27005122 rs4820691 1.00
> 15 27002392 rs4820690 27005158 rs4822749 1.00
> 16 27002392 rs4820690 27005194 rs4822750 0.988655
> 17 27002392 rs4820690 27005470 rs5761629 0.695642
> 19 27004298 rs5761627 27001689 rs4822747 0.695642
> 20 27004298 rs5761627 27002392 rs4820690 1.00
> 22 27004298 rs5761627 27004902 rs5752355 0.988681
> 23 27004298 rs5761627 27004964 rs4822748 0.988655
> 24 27004298 rs5761627 27005122 rs4820691 1.00
> 25 27004298 rs5761627 27005158 rs4822749 1.00
> >
> > dput (testlist)
> list(structure(list(BP_A = c(27001689L, 27001689L, 27001689L,
> 27001689L, 27001689L, 27001689L), SNP_A = c("rs4822747", "rs4822747",
> "rs4822747", "rs4822747", "rs4822747", "rs4822747"), BP_B = c(27002392L,
> 27004298L, 27004902L, 27004964L, 27005122L, 27005158L), SNP_B =
> c("rs4820690",
> "rs5761627", "rs5752355", "rs4822748", "rs4820691", "rs4822749"
> ), R2 = c(0.695642, 0.695642, 0.687861, 0.682715, 0.695642, 0.695642
> )), .Names = c("BP_A", "SNP_A", "BP_B", "SNP_B", "R2"), row.names = 2:7,
> class = "data.frame"),
> structure(list(BP_A = c(27002392L, 27002392L, 27002392L,
> 27002392L

[R] beta package for 3D PDF output

[Michail Vidiassov]

Dear All,

recent desktop versions of Adobe Acrobat and Adobe Reader have
built-in support for 3D models in PDF.

You can take a look at http://www.pdf3d.com/gallery.php for a gallery
of professional results achieved with a commercial tool.

For an overview of what 3D PDF can be used for I (naturally) suggest the article
http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0069446 .

To bring together 3D PDF and R, I have enhanced "rgl" package with
necessary output routine, so you draw something in 3D using rgl and
after that say "save it as 3D PDF".
The resulting 3D PDF model can be used standalone or in your LaTeX document.

The source code of the package and demo output are at
http://www2.iaas.msu.ru/tmp/u3d/rgl/

I have successfully compiled and run my package on Mac OS 10.9 64bit,
Windows XP 32 bit, Ubuntu 13.04 i686  and 13.10 ppc, Fedora 18 64bit,
but received a crash report in undisclosed Linux environment on some
CRAN test platform.

Testers welcome, feel free to e-mail me if necessary (by the time
search finds you this post I may be off the list).

Sincerely, Michail

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

[Hermann Norpois]

Hello,

having a list like testlist I would like to transform it in dataframe. How
does it work?

Thanks
Hermann

> testlist
[[1]]
  BP_A SNP_A BP_B SNP_B   R2
2 27001689 rs4822747 27002392 rs4820690 0.695642
3 27001689 rs4822747 27004298 rs5761627 0.695642
4 27001689 rs4822747 27004902 rs5752355 0.687861
5 27001689 rs4822747 27004964 rs4822748 0.682715
6 27001689 rs4822747 27005122 rs4820691 0.695642
7 27001689 rs4822747 27005158 rs4822749 0.695642

[[2]]
   BP_A SNP_A BP_B SNP_B   R2
9  27002392 rs4820690 27001689 rs4822747 0.695642
11 27002392 rs4820690 27004298 rs5761627 1.00
12 27002392 rs4820690 27004902 rs5752355 0.988681
13 27002392 rs4820690 27004964 rs4822748 0.988655
14 27002392 rs4820690 27005122 rs4820691 1.00
15 27002392 rs4820690 27005158 rs4822749 1.00
16 27002392 rs4820690 27005194 rs4822750 0.988655
17 27002392 rs4820690 27005470 rs5761629 0.695642

[[3]]
   BP_A SNP_A BP_B SNP_B   R2
19 27004298 rs5761627 27001689 rs4822747 0.695642
20 27004298 rs5761627 27002392 rs4820690 1.00
22 27004298 rs5761627 27004902 rs5752355 0.988681
23 27004298 rs5761627 27004964 rs4822748 0.988655
24 27004298 rs5761627 27005122 rs4820691 1.00
25 27004298 rs5761627 27005158 rs4822749 1.00

#The aim is the dataframe
   BP_A SNP_A BP_B SNP_B   R2
2  27001689 rs4822747 27002392 rs4820690 0.695642
3  27001689 rs4822747 27004298 rs5761627 0.695642
4  27001689 rs4822747 27004902 rs5752355 0.687861
5  27001689 rs4822747 27004964 rs4822748 0.682715
6  27001689 rs4822747 27005122 rs4820691 0.695642
7  27001689 rs4822747 27005158 rs4822749 0.695642
9  27002392 rs4820690 27001689 rs4822747 0.695642
11 27002392 rs4820690 27004298 rs5761627 1.00
12 27002392 rs4820690 27004902 rs5752355 0.988681
13 27002392 rs4820690 27004964 rs4822748 0.988655
14 27002392 rs4820690 27005122 rs4820691 1.00
15 27002392 rs4820690 27005158 rs4822749 1.00
16 27002392 rs4820690 27005194 rs4822750 0.988655
17 27002392 rs4820690 27005470 rs5761629 0.695642
19 27004298 rs5761627 27001689 rs4822747 0.695642
20 27004298 rs5761627 27002392 rs4820690 1.00
22 27004298 rs5761627 27004902 rs5752355 0.988681
23 27004298 rs5761627 27004964 rs4822748 0.988655
24 27004298 rs5761627 27005122 rs4820691 1.00
25 27004298 rs5761627 27005158 rs4822749 1.00
>
> dput (testlist)
list(structure(list(BP_A = c(27001689L, 27001689L, 27001689L,
27001689L, 27001689L, 27001689L), SNP_A = c("rs4822747", "rs4822747",
"rs4822747", "rs4822747", "rs4822747", "rs4822747"), BP_B = c(27002392L,
27004298L, 27004902L, 27004964L, 27005122L, 27005158L), SNP_B =
c("rs4820690",
"rs5761627", "rs5752355", "rs4822748", "rs4820691", "rs4822749"
), R2 = c(0.695642, 0.695642, 0.687861, 0.682715, 0.695642, 0.695642
)), .Names = c("BP_A", "SNP_A", "BP_B", "SNP_B", "R2"), row.names = 2:7,
class = "data.frame"),
structure(list(BP_A = c(27002392L, 27002392L, 27002392L,
27002392L, 27002392L, 27002392L, 27002392L, 27002392L), SNP_A =
c("rs4820690",
"rs4820690", "rs4820690", "rs4820690", "rs4820690", "rs4820690",
"rs4820690", "rs4820690"), BP_B = c(27001689L, 27004298L,
27004902L, 27004964L, 27005122L, 27005158L, 27005194L, 27005470L
), SNP_B = c("rs4822747", "rs5761627", "rs5752355", "rs4822748",
"rs4820691", "rs4822749", "rs4822750", "rs5761629"), R2 = c(0.695642,
1, 0.988681, 0.988655, 1, 1, 0.988655, 0.695642)), .Names = c("BP_A",
"SNP_A", "BP_B", "SNP_B", "R2"), row.names = c(9L, 11L, 12L,
13L, 14L, 15L, 16L, 17L), class = "data.frame"), structure(list(
BP_A = c(27004298L, 27004298L, 27004298L, 27004298L,
27004298L, 27004298L), SNP_A = c("rs5761627", "rs5761627",
"rs5761627", "rs5761627", "rs5761627", "rs5761627"),
BP_B = c(27001689L, 27002392L, 27004902L, 27004964L,
27005122L, 27005158L), SNP_B = c("rs4822747", "rs4820690",
"rs5752355", "rs4822748", "rs4820691", "rs4822749"),
R2 = c(0.695642, 1, 0.988681, 0.988655, 1, 1)), .Names = c("BP_A",
"SNP_A", "BP_B", "SNP_B", "R2"), row.names = c(19L, 20L,
22L, 23L, 24L, 25L), class = "data.frame"))

[[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] xts objects comparison

[Tstudent]

I have the following two xts object:

https://dl.dropboxusercontent.com/u/102669/series1.rdata

https://dl.dropboxusercontent.com/u/102669/series2.rdata

With str i see that they are both xts objects

I can't understand why it's impossible to compare each element.

For example: series1 > series2
Why i don't get an xts object with a sequence of true and false?

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Fw: Lower 'Pi' 3.1415926... and Exact 'Pi' and Squaring of Circle

[sarvajannnadha reddy]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Replace NA's with value in the next row

[Frans Marcelissen]

R purists forbid the use of the for loop, but I am afraid this is the most
simple solution:

for (i in 1:(length(V6)-1)) if(is.na(V6[i])) V6[i]<-V6[i+1]




2013/11/14 Jim Lemon 

> On 11/14/2013 04:02 PM, dila radi wrote:
>
>> Hi all,
>>
>> I have a data set which treat missing value as NA and now I need to
>> replace
>> all these NA's by using number in the same row but different column.
>>
>> Here is the part of my data:
>>   V1 V2 V3 V4 V5 V6 V7  0 0 0 1.2 0 0 0.259  0 0 12.8 0 23.7 0 8.495  6 0
>> 81.7 0.2 0 20 19.937  0 1.5 60.9 0 0 15.5 13.900  1 13 56.8 17.5 32.8 6.4
>> 27.654  4 3 66.4 2 0.3 NA 17.145
>>
>>
>> I want to replace (V6, 6) with (V7, 6). I have about 1000 NA's in V6 which
>> I want to replace  with the same row in V7. The other values in V6, I want
>> to keep remain the same.
>>
>> How to achieve this? Assuming my data is called "Targetstation",  I have
>> tried this:
>>
>> Targetstation<- within(Targetstation, v6<- replace(v6, is.na(v6), v7[
>> is.na
>> (v6)]))
>>
>> But R gives me this:
>>
>> Warning messages:
>>
>> 1: In is.na(v6) : is.na() applied to non-(list or vector) of type 'NULL'
>>
>> 2: In is.na(v6) : is.na() applied to non-(list or vector) of type 'NULL'
>>
>>
>>  Hi Dila,
> You could do this like this:
>
> V6NA<-is.na(Targetstation$V6)
> Targetstation$V6[V6NA]<-Targetstation$V7[V6NA]
>
> but if you want to use the above, I think you will have to replace the
> is.na(V6) with is.na(Targetstation$V6) or use the method above.
>
> Jim
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/
> posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 


---
dr F.H.G. (Frans) Marcelissen
DigiPsy (www.DigiPsy.nl )
Pomperschans 26
5595 AV Leende
tel: 040 2065030/06 2325 06 53
skype adres: frans.marcelissen
email: frans.marcelis...@digipsy.nl

[[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] volume of ellipsoid

[Carl Witthoft]

Well, given that an upvoted question on math.stackexchange got no answers,
I'd say you're asking a very difficult question.   Perhaps this paper  
http://www.geometrictools.com/Documentation/IntersectionOfEllipsoids.pdf

will be of some help.   It's possible that you could do:
1) find the ellipse of intersection of the two ellipsoids
2) find some formula for the equivalent of the "spherical cap" volume
3) apply that formula to each ellipsoid, using the intersection ellipse as
the base of the caps.

All this, btw, is math,  not really related to R or to software programming
in general.


yuanzhi wrote
> Hi, everyone!
> 
> I have a matrix X(n*p) which is n samples from a p-dimensional normal  
> distribution. Now I want to make the ellipsoid containing (1-α)% of  
> the probability in the distribution based on Mahalanobis distance:
> 
> μ:(x-μ)'Σ^(-1)(x-μ)≤χ2p(α)
> 
> where x and Σ is the mean and variance-covariance matrix of X. Are  
> there some functions in R which can plot the ellipsoid and calculate  
> the volume (area for p=2, hypervolume for p>3) of the ellipsoid? I  
> only find the "ellipsoidhull" function in "cluster" package, but it  
> deal with ellipsoid hull containing X, which obvious not the one I want.
> 
> After that, a more difficult is the following. If we have a series of  
> matrix X1, X2, X3,… Xm which follow different p-dimensional normal  
> distributions. And we make m ellipsoids (E1, E2, … Em) for each matrix  
> like the before. How can we calculate the volume of union of the m  
> ellipsoids? One possible solution for this problem is the  
> inclusion-exclusion principle:
> 
> V(⋃Ei)(1≤i≤m)=
> V(Ei)(1≤i≤m)-∑V(Ei⋂Ej)(1≤i 
> But the problem is how to calculate the volume of intersection between  
> 2, 3 or more ellipsoids. Are there some functions which can calculate  
> the volume of intersection between two region or functions which  
> directly calculate the volume of a union of two region(the region here  
> is ellipsoid). OR yo you have any good ideas solving this problem in  
> R? Thank you all in advance!
> 
> Yuanzhi
> 
> __
> R-help@

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





--
View this message in context: 
http://r.789695.n4.nabble.com/volume-of-ellipsoid-tp4680409p4680435.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] Printing from windows device; was: (sin asunto)

[Uwe Ligges]

On 13.11.2013 22:45, Elisa Frutos Bernal wrote:

Hi!

I need to print a graph that I have in a window. Previously I used:

try(win.print(), silent = TRUE)
 if (geterrmessage() != "Error in win.print() : unable to start
device devWindows\n") {
 plotFunctiond(screen = FALSE)
 dev.off()

but now it is not possible.
Can you help me?


No, since I do not know what
plotFunctiond()
does.

Please use a sensible subject line.

Best,
Uwe Ligges

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

[Jim Lemon]

On 11/14/2013 09:28 PM, peron wrote:

Hello



I have data in following form : 100 ind describe by two variables x and y.



D<-data .frame( x=rnorm(3), y=rnorm(3), size=c(50,10,40))



I want data for individual, i.e, 100 observations for my 100 ind.


Hi Olivier,
From the above, it seems you have three values for the three groups. If 
you try:


df100<-data.frame(x=rep(D$x,D$size),y=rep(D$y,D$size))

You will get what you request, if the values are the same for all 
members of each group. I suspect that those three values may be means or 
some other summary measure.


Jim

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

[R] Transform aggregated data to individual data

[peron]

Hello 

 

I have data in following form : 100 ind describe by two variables x and y.

 

D<-data .frame( x=rnorm(3), y=rnorm(3), size=c(50,10,40))

 

I want data for individual, i.e, 100 observations for my 100 ind.

 

Thank for your help

 

Olivier Peron


[[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] Fitting arbitrary curve to 1D data with error bars

[Suzen, Mehmet]

Maybe you are after "weights" option given by 'lm' or 'glm'

See: 
http://stackoverflow.com/questions/6375650/function-for-weighted-least-squares-estimates

On 14 November 2013 10:01, Erkcan Özcan  wrote:
> Thanks, but if you have another closer look to my post, you will see that my 
> question has nothing to do with drawing error bars on a plot.
>
> What I want is to do a curve fit to a data with error bars.
>
> Best,
> e.
>
> On 14 Nov 2013, at 04:21, Suzen, Mehmet wrote:
>
>> If you are after adding error bars in a scatter plot; one example is
>> given below :
>>
>> #some example data
>> set.seed(42)
>> df <- data.frame(x = rep(1:10,each=5), y = rnorm(50))
>>
>> #calculate mean, min and max for each x-value
>> library(plyr)
>> df2 <- ddply(df,.(x),function(df)
>> c(mean=mean(df$y),min=min(df$y),max=max(df$y)))
>>
>> #plot error bars
>> library(Hmisc)
>> with(df2,errbar(x,mean,max,min))
>> grid(nx=NA,ny=NULL)
>>
>> (From: 
>> http://stackoverflow.com/questions/13032777/scatter-plot-with-error-bars)
>>
>

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

Re: [R] Fitting arbitrary curve to 1D data with error bars

[Erkcan Özcan]

Thanks, but if you have another closer look to my post, you will see that my 
question has nothing to do with drawing error bars on a plot.

What I want is to do a curve fit to a data with error bars.

Best,
e.

On 14 Nov 2013, at 04:21, Suzen, Mehmet wrote:

> If you are after adding error bars in a scatter plot; one example is
> given below :
> 
> #some example data
> set.seed(42)
> df <- data.frame(x = rep(1:10,each=5), y = rnorm(50))
> 
> #calculate mean, min and max for each x-value
> library(plyr)
> df2 <- ddply(df,.(x),function(df)
> c(mean=mean(df$y),min=min(df$y),max=max(df$y)))
> 
> #plot error bars
> library(Hmisc)
> with(df2,errbar(x,mean,max,min))
> grid(nx=NA,ny=NULL)
> 
> (From: 
> http://stackoverflow.com/questions/13032777/scatter-plot-with-error-bars)
> 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Replace NA's with value in the next row

[Jim Lemon]

On 11/14/2013 04:02 PM, dila radi wrote:

Hi all,

I have a data set which treat missing value as NA and now I need to replace
all these NA's by using number in the same row but different column.

Here is the part of my data:
  V1 V2 V3 V4 V5 V6 V7  0 0 0 1.2 0 0 0.259  0 0 12.8 0 23.7 0 8.495  6 0
81.7 0.2 0 20 19.937  0 1.5 60.9 0 0 15.5 13.900  1 13 56.8 17.5 32.8 6.4
27.654  4 3 66.4 2 0.3 NA 17.145


I want to replace (V6, 6) with (V7, 6). I have about 1000 NA's in V6 which
I want to replace  with the same row in V7. The other values in V6, I want
to keep remain the same.

How to achieve this? Assuming my data is called "Targetstation",  I have
tried this:

Targetstation<- within(Targetstation, v6<- replace(v6, is.na(v6), v7[is.na
(v6)]))

But R gives me this:

Warning messages:

1: In is.na(v6) : is.na() applied to non-(list or vector) of type 'NULL'

2: In is.na(v6) : is.na() applied to non-(list or vector) of type 'NULL'



Hi Dila,
You could do this like this:

V6NA<-is.na(Targetstation$V6)
Targetstation$V6[V6NA]<-Targetstation$V7[V6NA]

but if you want to use the above, I think you will have to replace the 
is.na(V6) with is.na(Targetstation$V6) or use the method above.


Jim

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