[R] Find peaks in dataset(x,y) and area for each peak

[Ramya]

http://r.789695.n4.nabble.com/file/n3312061/x_and_y_values.txt
x_and_y_values.txt 


I have the absorbance values form HPLC Chromatogram. I need to find the
peaks in datapoints and area under each peak.

I am not sure how how to find the peaks, I tried couple of libraries and
peak function but they return me a list of values which when computed from
area turns out to be huge nnumbers. In the dataset I know there are 4 peaks
but i couldnt find through R.

getPeaks method was kind of promising but i am not sure how to use it.

Any help would be appreciated.

Thanks
Ramya
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Find-peaks-in-dataset-x-y-and-area-for-each-peak-tp3312061p3312061.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] CDF of Sample Quantile

[Bentley Coffey]

Duncan,

I'm not sure how I missed your message. Sorry. What you describe is what I
do when (n-1)p is an integer so that R computes the sample quantile using a
single order statistic. My later posting has that exact binomial expression
in there as a special case. When (n-1)p is not an integer then R uses a
weighted average of 2 order statistics, in which case I'm left with my
standing problem...

On Mon, Feb 14, 2011 at 2:26 PM, Duncan Murdoch wrote:

> On 14/02/2011 9:58 AM, Bentley Coffey wrote:
>
>> I need to calculate the probability that a sample quantile will exceed a
>> threshold given the size of the iid sample and the parameters describing
>> the
>> distribution of each observation (normal, in my case). I can compute the
>> probability with brute force simulation: simulate a size N sample, apply
>> R's
>> quantile() function on it, compare it to the threshold, replicate this
>> MANY
>> times, and count the number of times the sample quantile exceeded the
>> threshold (dividing by the total number of replications yields the
>> probability of interest). The problem is that the number of replications
>> required to get sufficient precision (3 digits say) is so HUGE that this
>> takes FOREVER. I have to perform this task so much in my script (searching
>> over the sample size and repeated for several different distribution
>> parameters) that it takes too many hours to run.
>>
>> I've searched for pre-existing code to do this in R and haven't found
>> anything. Perhaps I'm missing something. Is anyone aware of an R function
>> to
>> compute this probability?
>>
>> I've tried writing my own code using the fact that R's quantile() function
>> is a linear combination of 2 order statistics. Basically, I wrote down the
>> mathematical form of the joint pdf for the 2 order statistics (a function
>> of
>> the sample size and the distribution parameters) then performed a
>> pseudo-Monte Carlo integration (i.e. using Halton Draws rather than R's
>> random draws) over the region where the sample quantile exceeds the
>> threshold. In theory, this should work and it takes about 1000 times fewer
>> clock cycles to compute than the Brute Force approach. My problem is that
>> there is a significant discrepancy between the results using Brute Force
>> and
>> using this more efficient approach that I have coded up. I believe that
>> the
>> problem is numerical error but it could be some programming bug;
>> regardless,
>> I have been unable to locate the source of this problem and have spent
>> over
>> 20 hours trying to identify it this weekend. Please, somebody help!!!
>>
>> So, again, my question: is there code in R for quickly evaluating the CDF
>> of
>> a Sample Quantile given the sample size and the parameters governing the
>> distribution of each iid point in the sample?
>>
>
> I think the answer to your question is no, but I think it's the wrong
> question.  Suppose Xm:n is the mth sample quantile in a sample of size n,
> and you want to calculate P(Xm:n > x).  Let X be a draw from the underlying
> distribution, and suppose P(X > x) = p.  Then the event Xm:n > x
> is the same as the event that out of n independent draws of X, at least
> n-m+1 are bigger than x:  a binomial probability.  And R can calculate
> binomial probabilities using pbinom().
>
> Duncan Murdoch
>
>

[[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] How to change dataframe to tables

[Lao Meng]

The data is in the attachment.

What I wanna get is:
, , Sex = Male
   Eye
HairBrown Blue Hazel Green
  Black32   1110 3
  Brown53   502515
  Red  10   10 7 7
  Blond 3   30 5 8
, , Sex = Female
   Eye
HairBrown Blue Hazel Green
  Black369 5 2
  Brown66   342914
  Red  167 7 7
  Blond 4   64 5 8


Then,how to change the dataframe to tables?

Many thanks for your help.

Best
"hair"  "eye"   "sex"   "freq"
"Black" "Brown" "Male"  32
"Black" "Blue"  "Male"  11
"Black" "Hazel" "Male"  10
"Black" "Green" "Male"  3
"Brown" "Brown" "Male"  38
"Brown" "Blue"  "Male"  50
"Brown" "Hazel" "Male"  25
"Brown" "Green" "Male"  15
"Red"   "Brown" "Male"  10
"Red"   "Blue"  "Male"  10
"Red"   "Hazel" "Male"  7
"Red"   "Green" "Male"  7
"Blond" "Brown" "Male"  3
"Blond" "Blue"  "Male"  30
"Blond" "Hazel" "Male"  5
"Blond" "Green" "Male"  8
"Black" "Brown" "Female"36
"Black" "Blue"  "Female"9
"Black" "Hazel" "Female"5
"Black" "Green" "Female"2
"Brown" "Brown" "Female"81
"Brown" "Blue"  "Female"34
"Brown" "Hazel" "Female"29
"Brown" "Green" "Female"14
"Red"   "Brown" "Female"16
"Red"   "Blue"  "Female"7
"Red"   "Hazel" "Female"7
"Red"   "Green" "Female"7
"Blond" "Brown" "Female"4
"Blond" "Blue"  "Female"64
"Blond" "Hazel" "Female"5
"Blond" "Green" "Female"8
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] xyplot formula

[Geoff Russell]

On Fri, Feb 18, 2011 at 5:42 PM, Phil Spector wrote:

> Geoff -
>   I think this will get you closer to a solution:
>

Yes, much closer, many thanks.

Geoff.


>
> newdf = reshape(df,varying=names(df)[-c(1,2)],direction='long',
>times=2000:2003,idvar=c('country','food'),v.names='X',
>timevar='year')
>
> xyplot(X~year|country*food,data=newdf)
>- Phil Spector
> Statistical Computing Facility
> Department of Statistics
> UC Berkeley
> spec...@stat.berkeley.edu
>
>
>
>
> On Fri, 18 Feb 2011, Geoff Russell wrote:
>
>  df=data.frame(country=c("A","A","A","B","B","B"),
>>   food=rep(c("Apples","Pears","Bananas"),2),
>>   X2000=c(4,5,6,7,6,8),
>>   X2001=c(4,5,6,7,6,8),
>>   X2002=c(4,5,6,7,6,8),
>>   X2003=c(4,5,6,7,6,8));
>>
>> I have data in the above form trying to get a plot of each fruit over time
>> year conditioned on country and food.
>>
>> I tried,
>>
>> xyplot(X2000+X2001+X2002+X2003~2000:2003|country*food,df)
>>
>> But think I need to transform the data, just not sure
>> how.  Any help gratefully received.
>>
>> Cheers,
>> Geoff.
>>
>>[[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.
>>
>>


-- 
6 Fifth Ave,
St Morris, S.A. 5068
Australia
Ph: 041 8805 184 / 08 8332 5069
http://perfidy.com.au

[[alternative HTML version deleted]]

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

Re: [R] xyplot formula

[Phil Spector]

Geoff -
   I think this will get you closer to a solution:

newdf = reshape(df,varying=names(df)[-c(1,2)],direction='long',
times=2000:2003,idvar=c('country','food'),v.names='X',
timevar='year')

xyplot(X~year|country*food,data=newdf)
- Phil Spector
 Statistical Computing Facility
 Department of Statistics
 UC Berkeley
 spec...@stat.berkeley.edu



On Fri, 18 Feb 2011, Geoff Russell wrote:


df=data.frame(country=c("A","A","A","B","B","B"),
   food=rep(c("Apples","Pears","Bananas"),2),
   X2000=c(4,5,6,7,6,8),
   X2001=c(4,5,6,7,6,8),
   X2002=c(4,5,6,7,6,8),
   X2003=c(4,5,6,7,6,8));

I have data in the above form trying to get a plot of each fruit over time
year conditioned on country and food.

I tried,

xyplot(X2000+X2001+X2002+X2003~2000:2003|country*food,df)

But think I need to transform the data, just not sure
how.  Any help gratefully received.

Cheers,
Geoff.

[[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] problem with rbind when data frame contains an date-time variable "POSIXt" "POSIXlt"

[Andrew Anglemyer]

Thanks again for your help.  I was sniffing around this solution forever,
but that did it cleanly and without any problems.



On Thu, Feb 17, 2011 at 10:21 PM,  wrote:

> The solution is probably to make the data-time columns POSIXct:
>
> > x <- read.table(textConnection("
> + ID event.date.time
> + 1 '2009-07-23 00:20:00'
> + 2 '2009-08-18 16:25:00'
> + 3 '2009-08-13 08:30:00'
> + "), header = TRUE)
> > y <- read.table(textConnection("
> + ID event.date.time
> + 4 '2009-08-25 10:25:00'
> + 5 '2009-08-10 06:20:00'
> + 6 '2009-10-09 08:20:00'
> + "), header = TRUE)
> > closeAllConnections()
> >
> > x <- within(x,
> + event.date.time <- as.POSIXct(as.character(event.date.time),
> + format = "%Y-%m-%d %H:%M:%S"))
> > y <- within(y,
> + event.date.time <- as.POSIXct(as.character(event.date.time),
> + format = "%Y-%m-%d %H:%M:%S"))
> > z <- rbind(x, y)
> > z
>  ID event.date.time
> 1  1 2009-07-23 00:20:00
> 2  2 2009-08-18 16:25:00
> 3  3 2009-08-13 08:30:00
> 4  4 2009-08-25 10:25:00
> 5  5 2009-08-10 06:20:00
> 6  6 2009-10-09 08:20:00
> >
>
> No problems.
>
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> On Behalf Of Andrew Anglemyer
> Sent: Friday, 18 February 2011 3:54 PM
> To: r-help@r-project.org
> Subject: [R] problem with rbind when data frame contains an date-time
> variable "POSIXt" "POSIXlt"
>
> I'm trying to rbind two data frames, both with the same columns names.  One
> of the columns is a variable with date-time and this variable is causing
> the
> rbind to fail--giving the error
> "Error in names(value[[jj]])[ri] <- nm :  'names' attribute [7568] must be
> the same length as the vector [9]"
>
> Is there a way to stack or rbind these two data frames even with this
> extended date-time variable?  The class of event.date.time in each data
> frame is POSIXt POSIXlt.
> x
> ID event.date.time
> 1 2009-07-23 00:20:00
> 2 2009-08-18 16:25:00
> 3 2009-08-13 08:30:00
>
> y
> ID event.date.time
> 4 2009-08-25 10:25:00
> 5 2009-08-10 06:20:00
> 6 2009-10-09 08:20:00
>
> I would like to get
>
> z
> ID event.date.time
> 1 2009-07-23 00:20:00
> 2 2009-08-18 16:25:00
> 3 2009-08-13 08:30:00
> 4 2009-08-25 10:25:00
> 5 2009-08-10 06:20:00
> 6 2009-10-09 08:20:00
>
> I've looked at stripping the dates and times, but it would be really
> helpful
> for my purposes to keep the extended variable date-time variable (have to
> eventually get 24 hours from baseline.date.time).
>
> thanks for any and all help!
> Andy
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

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

[R] xyplot formula

[Geoff Russell]

df=data.frame(country=c("A","A","A","B","B","B"),
food=rep(c("Apples","Pears","Bananas"),2),
X2000=c(4,5,6,7,6,8),
X2001=c(4,5,6,7,6,8),
X2002=c(4,5,6,7,6,8),
X2003=c(4,5,6,7,6,8));

I have data in the above form trying to get a plot of each fruit over time
year conditioned on country and food.

I tried,

xyplot(X2000+X2001+X2002+X2003~2000:2003|country*food,df)

But think I need to transform the data, just not sure
how.  Any help gratefully received.

Cheers,
Geoff.

[[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] problem with rbind when data frame contains an date-time variable "POSIXt" "POSIXlt"

[Bill.Venables]

The solution is probably to make the data-time columns POSIXct:

> x <- read.table(textConnection("
+ ID event.date.time
+ 1 '2009-07-23 00:20:00'
+ 2 '2009-08-18 16:25:00'
+ 3 '2009-08-13 08:30:00'
+ "), header = TRUE)
> y <- read.table(textConnection("
+ ID event.date.time
+ 4 '2009-08-25 10:25:00'
+ 5 '2009-08-10 06:20:00'
+ 6 '2009-10-09 08:20:00'
+ "), header = TRUE)
> closeAllConnections()
> 
> x <- within(x, 
+ event.date.time <- as.POSIXct(as.character(event.date.time),
+ format = "%Y-%m-%d %H:%M:%S"))
> y <- within(y, 
+ event.date.time <- as.POSIXct(as.character(event.date.time),
+ format = "%Y-%m-%d %H:%M:%S"))
> z <- rbind(x, y)
> z
  ID event.date.time
1  1 2009-07-23 00:20:00
2  2 2009-08-18 16:25:00
3  3 2009-08-13 08:30:00
4  4 2009-08-25 10:25:00
5  5 2009-08-10 06:20:00
6  6 2009-10-09 08:20:00
> 

No problems. 

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of Andrew Anglemyer
Sent: Friday, 18 February 2011 3:54 PM
To: r-help@r-project.org
Subject: [R] problem with rbind when data frame contains an date-time variable 
"POSIXt" "POSIXlt"

I'm trying to rbind two data frames, both with the same columns names.  One
of the columns is a variable with date-time and this variable is causing the
rbind to fail--giving the error
"Error in names(value[[jj]])[ri] <- nm :  'names' attribute [7568] must be
the same length as the vector [9]"

Is there a way to stack or rbind these two data frames even with this
extended date-time variable?  The class of event.date.time in each data
frame is POSIXt POSIXlt.
x
ID event.date.time
1 2009-07-23 00:20:00
2 2009-08-18 16:25:00
3 2009-08-13 08:30:00

y
ID event.date.time
4 2009-08-25 10:25:00
5 2009-08-10 06:20:00
6 2009-10-09 08:20:00

I would like to get

z
ID event.date.time
1 2009-07-23 00:20:00
2 2009-08-18 16:25:00
3 2009-08-13 08:30:00
4 2009-08-25 10:25:00
5 2009-08-10 06:20:00
6 2009-10-09 08:20:00

I've looked at stripping the dates and times, but it would be really helpful
for my purposes to keep the extended variable date-time variable (have to
eventually get 24 hours from baseline.date.time).

thanks for any and all help!
Andy

[[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] problem with rbind when data frame contains an date-time variable "POSIXt" "POSIXlt"

[Andrew Anglemyer]

I'm trying to rbind two data frames, both with the same columns names.  One
of the columns is a variable with date-time and this variable is causing the
rbind to fail--giving the error
"Error in names(value[[jj]])[ri] <- nm :  'names' attribute [7568] must be
the same length as the vector [9]"

Is there a way to stack or rbind these two data frames even with this
extended date-time variable?  The class of event.date.time in each data
frame is POSIXt POSIXlt.
x
ID event.date.time
1 2009-07-23 00:20:00
2 2009-08-18 16:25:00
3 2009-08-13 08:30:00

y
ID event.date.time
4 2009-08-25 10:25:00
5 2009-08-10 06:20:00
6 2009-10-09 08:20:00

I would like to get

z
ID event.date.time
1 2009-07-23 00:20:00
2 2009-08-18 16:25:00
3 2009-08-13 08:30:00
4 2009-08-25 10:25:00
5 2009-08-10 06:20:00
6 2009-10-09 08:20:00

I've looked at stripping the dates and times, but it would be really helpful
for my purposes to keep the extended variable date-time variable (have to
eventually get 24 hours from baseline.date.time).

thanks for any and all help!
Andy

[[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] recoding a data in different way: please help

[Dennis Murphy]

Hi:

This is as far as I could get:

df <- read.table(textConnection("
 Individual  Parent1  Parent2 mark1   mark2
 10   0   12  11
 20   0   11  22
 30   0   13  22
 40   0   13  11
 51   2   11  12
 61   2   12  12
 73   4   11  12
 83   4   13  12
 91   4   11  12
 10   1   4   11  12"), header = TRUE)
df2 <- transform(df, Parent1 = replace(Parent1, Parent1 == 0, NA),
 Parent2 = replace(Parent2, Parent2 == 0, NA))
df2 <- transform(df2, imark1p1 = df2$mark1[df2$Parent1],   # Parent 1's
mark1
  imark1p2 = df2$mark1[df2$Parent2],
# Parent 2's mark1
  imark2p1 = df2$mark2[df2$Parent1],
# Parent 1's mark2
  imark2p2 = df2$mark2[df2$Parent2])
# Parent 2's mark2

I created df2 so as not to overwrite the original in case of a mistake. At
this point, you have several sets of vectors that you can compare; e.g.,
mark1 with imark1p1 and imark1p2. Like Josh, I couldn't make heads or tails
out of what these logical tests were meant to output, but perhaps this gives
you a broader template with which to work. At this point, you can probably
remove the rows corresponding to the parents. I believe ifelse() is your
friend here - it can perform logical tests in a vectorized fashion. As long
as the tests are consistent from one individual to the next, it's likely to
be an efficient route.

HTH,
Dennis

On Thu, Feb 17, 2011 at 6:21 PM, Umesh Rosyara  wrote:

> Dear R users
>
> The following question looks simple but I have spend alot of time to solve
> it. I would highly appeciate your help.
>
> I have following dataset from family dataset :
>
> Here we have individuals and their two parents and their marker scores
> (marker1, marker2,and so on). 0 means that their parent information not
> available.
>
>
> Individual  Parent1  Parent2 mark1   mark2
> 10   0   12  11
> 20   0   11  22
> 30   0   13  22
> 40   0   13  11
> 51   2   11  12
> 61   2   12  12
> 73   4   11  12
> 83   4   13  12
> 91   4   11  12
> 10   1   4   11  12
>
> I want to recode mark1 and other mark2.and so on column by looking
> indvidual parent (Parent1 and Parent2).
>
> For example
>
> Take case of Individual 5, who's Parent 1 is 1 (has mark1 score 12) and
> Parent 2 is 2 (has mark1 score 11). Individual 5 has mark1 score 11.
> Suppose
> I have following condition to recode Individual 5's mark1 score:
>
> For mark1 variable, If Parent1 score "11" and Parent2 score "22" and recode
> indvidual 5's score, "12"=1, else 0
>If Parent1 score "12" and Parent2 score
> "22" and recode individual 5's score, "22"=1, "12"= 0.5, else 0
>.more conditions
>
> Similarly the pointer should move from individual 5 to n individuals at the
> end of the file.
>
>  Thank you in advance
>
> Umesh 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.
>

[[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] recoding a data in different way: please help

[Joshua Wiley]

Dear Umesh,

I could not figure out exactly what your recoding scheme was, so I do
not have a specific solution for you.  That said, the following
functions may help you get started.

?ifelse # vectorized and different from using if () statements
?if #
?Logic ## logical operators for your tests
## if you install and load the "car" package by John Fox
?recode # a function for recoding in package "car"

I am sure it is possible to string together some massive series of if
statements and then use a for loop, but that is probably the messiest
and slowest possible way.  I suspect there will be faster, neater
options, but I cannot say for certain without having a better feel for
how all the conditions work.

Best regards,

Josh

On Thu, Feb 17, 2011 at 6:21 PM, Umesh Rosyara  wrote:
> Dear R users
>
> The following question looks simple but I have spend alot of time to solve
> it. I would highly appeciate your help.
>
> I have following dataset from family dataset :
>
> Here we have individuals and their two parents and their marker scores
> (marker1, marker2,and so on). 0 means that their parent information not
> available.
>
>
> Individual      Parent1  Parent2         mark1   mark2
> 1        0       0       12      11
> 2        0       0       11      22
> 3        0       0       13      22
> 4        0       0       13      11
> 5        1       2       11      12
> 6        1       2       12      12
> 7        3       4       11      12
> 8        3       4       13      12
> 9        1       4       11      12
> 10       1       4       11      12
>
> I want to recode mark1 and other mark2.and so on column by looking
> indvidual parent (Parent1 and Parent2).
>
> For example
>
> Take case of Individual 5, who's Parent 1 is 1 (has mark1 score 12) and
> Parent 2 is 2 (has mark1 score 11). Individual 5 has mark1 score 11. Suppose
> I have following condition to recode Individual 5's mark1 score:
>
> For mark1 variable, If Parent1 score "11" and Parent2 score "22" and recode
> indvidual 5's score, "12"=1, else 0
>                                    If Parent1 score "12" and Parent2 score
> "22" and recode individual 5's score, "22"=1, "12"= 0.5, else 0
>                                    .more conditions
>
> Similarly the pointer should move from individual 5 to n individuals at the
> end of the file.
>
>  Thank you in advance
>
> Umesh 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.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.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] Integrate with an indicator function

[li li]

Hi Dennis,
  Thank you for your kind reply. Yes, essentially, we take integration
twice. However, I still have a few questions:
  First, if we consider doing a double integration at the end, since the
first integration in within the indicator function,
it seems to be difficult.

  Second, m1star is a univariate function, say, of x, since I already take
the integration. Then T is an univariate
function of x and func2 is a function of  x and c. It seems to make sense
if I integrate func2 with respect to x.

  Can you give more hint in terms of doing a double integration? Thank you
again.
   Hannah
2011/2/17 Dennis Murphy 

> Hi Hannah:
>
> You have a few things going on, but the bottom line is that numer and denom
> are both double integrals.
>
>  On Thu, Feb 17, 2011 at 1:06 PM, li li  wrote:
>
>> Hi all,
>>   I have some some problem with regard to finding the integral of a
>> function containing an indicator function.
>> please see the code below:
>>
>>
>> func1 <- function(x, mu){
>>  (mu^2)*dnorm(x, mean = mu, sd = 1)*dgamma(x, shape=2)}
>>
>
> curve(func1(x, 10), 0, 20)
> curve(func1(x, 5), 0, 20)
>
> both yield reasonable plots, so this is a function of one variable when mu
> is supplied. If you wanted to plot this as a function of mu, you could get a
> single curve by fixing x or integrating over x, which is what m1star appears
> to be doing.
>
> m1star <- function(x){
>>   integrate(func1, lower = 0, upper = Inf,x)$val}
>>
>
> u <- seq(0, 20, 0.05)
> v <- sapply(u, m1star)
> plot(u, v, type = 'l')
>
> yields a plot of m1star, which appears to marginalize func1 to make it a
> function of mu, from what I can tell.
>
>  T <- function(x){
>>   0.3*dnorm(x)/(0.3*dnorm(x)+0.7*m1star(x))}
>>
>
> plot(u, sapply(u, T), type = 'l')
>
> yields a plot of T, but having to use sapply() indicates that T also
> marginalizes a 2D function.
>
>
>
>> func2 <- function(x,c){(T(x) <=c)*0.3*dnorm(x)}
>> func3 <- function(x,c){(T(x) <= c)*(0.3*dnorm(x)+0.7*m1star(x))}
>>
>
> func2 uses T, which in turn uses m1star, so func2 is a marginalization of a
> 2D function whose support is on T(x) <= c. To integrate func2, you have to
> do the integration in m1star first, so basically you have a double integral
> to evaluate in numer. The same problem occurs in func3, since m1star() is a
> part of both T() and the convex combination. Therefore, denom also requires
> double integration.
>
>
>> numer <- function(c){
>>   integrate(func2, -Inf, Inf, c)$val}
>>
>> denom <- function(c){
>>   integrate(func3, lower, Inf,c)$val}
>>
>> Look into cubature or Rcuba for two packages that are capable of
> performing numerical double integration. An alternative solution is to use
> approxfun() to create a function object from evaluations of m1star and T,
> and then use the approxfun()s as the functions to be integrated in numer and
> denom. If you decide to go the approxfun route, make sure to make enough
> evaluations to reasonably capture the curvature in both m1star and T.
>
> HTH,
> Dennis
>
>>
>>
>> The error message is as below :
>>
>> > numer(0.5)
>> Error in integrate(func1, lower = 0, upper = Inf, x) :
>>  the integral is probably divergent
>>
>>[[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>

[[alternative HTML version deleted]]

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

[R] recoding a data in different way: please help

[Umesh Rosyara]

Dear R users 
 
The following question looks simple but I have spend alot of time to solve
it. I would highly appeciate your help.  
 
I have following dataset from family dataset : 
 
Here we have individuals and their two parents and their marker scores
(marker1, marker2,and so on). 0 means that their parent information not
available. 
 
 
Individual  Parent1  Parent2 mark1   mark2  
10   0   12  11 
20   0   11  22 
30   0   13  22 
40   0   13  11 
51   2   11  12 
61   2   12  12 
73   4   11  12 
83   4   13  12 
91   4   11  12 
10   1   4   11  12 
 
I want to recode mark1 and other mark2.and so on column by looking
indvidual parent (Parent1 and Parent2). 
 
For example 
 
Take case of Individual 5, who's Parent 1 is 1 (has mark1 score 12) and
Parent 2 is 2 (has mark1 score 11). Individual 5 has mark1 score 11. Suppose
I have following condition to recode Individual 5's mark1 score:
 
For mark1 variable, If Parent1 score "11" and Parent2 score "22" and recode
indvidual 5's score, "12"=1, else 0
If Parent1 score "12" and Parent2 score
"22" and recode individual 5's score, "22"=1, "12"= 0.5, else 0
.more conditions
 
Similarly the pointer should move from individual 5 to n individuals at the
end of the file. 
 
 Thank you in advance
 
Umesh 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] Transforming relational data

[Matthew Dowle]

Mathijs,

To my eyes you seem to have repeated back what is already done.

More R and less English would help. In other words if it is not 2.5
you need, what is it? Please provide some input and state what the
output should be (and what you tried already).

Matthew

-- 
View this message in context: 
http://r.789695.n4.nabble.com/Re-Transforming-relational-data-tp3307449p3311954.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] Using Weights in R

[Phil Spector]

Well, if you really want to produce what SAS does:


ageval = read.table(textConnection('AgeCat FINWT

+ 1 98
+ 2 62
+ 1 75
+ 3 39
+ 4 28
+ 2 47
+ 2 66
+ 4 83
+ 1 19
+ 3 50
+ '),header=TRUE)

one = aggregate(ageval$FINWT,ageval['AgeCat'],sum)
two = prop.table(as.table(one$x)) * 100
three = cumsum(one$x)
four = cumsum(two)
answer = data.frame(AgeCat=one$AgeCat,

+ Frequency=one$x,
+ Percent=as.numeric(two),
+ "Cumulative.Frequency"=three,
+ "Cumulative.Percent"=four)

answer

  AgeCat Frequency  Percent Cumulative.Frequency Cumulative.Percent
A  1   192 33.86243  192   33.86243
B  2   175 30.86420  367   64.72663
C  389 15.69665  456   80.42328
D  4   111 19.57672  567  100.0

You can probably find the specific part you want from the above code.

- Phil Spector
 Statistical Computing Facility
 Department of Statistics
 UC Berkeley
 spec...@stat.berkeley.edu

On Thu, 17 Feb 2011, Krishnan Viswanathan wrote:


I am new to R. I have a data set like this (given below is a fictional
dataset):

AgeCat FINWT
1 98
2 62
1 75
3 39
4 28
2 47
2 66
4 83
1 19
3 50

I need to calculate the weighted distribution of the variable AgeCat.

In SAS i can do:

proc freq data=ageval;
  tables agecat;
weight finwt;
run;

What or is there an equivalent in R?

TIA,
Krishnan

--
Krishnan Viswanathan
1101 High Meadow Dr
Tallahassee FL 32311

[[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] logsumexp function in R

[Spencer Graves]

On 2/17/2011 3:49 PM, William Dunlap wrote:

-Original Message-
From: r-help-boun...@r-project.org
[mailto:r-help-boun...@r-project.org] On Behalf Of Petr Savicky
Sent: Wednesday, February 16, 2011 10:46 PM
To: r-help@r-project.org
Subject: Re: [R] logsumexp function in R

On Wed, Feb 16, 2011 at 04:14:38PM -0500, Brian Tsai wrote:

Hi,

I was wondering if anyone has implemented a numerically

stable function to

compute log(sum(exp(x))) ?

Try the following

log(sum(exp(x - max(x + max(x)

Sometimes the log1p(x) function, which gives log(1+x)
accurately for small abs(x), helps a little more.
Compare the following 3 functions, which I think give
the same thing for 'ordinary' values of x:

f0<- function(x) log(sum(exp(x)))
f1<- function(x) log(sum(exp(x - max(x + max(x)
f2<- function(x) { x<- sort(x) # mainly so max(x)==x[length(x)]
 n<- length(x)
 log1p(sum(exp(x[-n] - x[n]))) + x[n]
   }



  That's great.  The following seems to be equivalent but takes 
less time by replacing an expensive sort with "which.max":



f2.0 <- function(x){
  xmax <- which.max(x)
  log1p(sum(exp(x[-xmax]-x[xmax])))+x[xmax]
}


 set.seed(1)
> system.time(f2(rnorm(1e7)))
   user  system elapsed
   8.330.258.69
> system.time(f2.0(rnorm(1e7)))
   user  system elapsed
   3.410.564.09
> system.time(f2(rnorm(1e7)))
   user  system elapsed
   7.990.278.35
> system.time(f2.0(rnorm(1e7)))
   user  system elapsed
   3.550.334.07


  Spencer


But for the following x f2 gives a more accurate result
than f1, which in turn often gives a more accurate result
than f0:
   >  x<- c(0, -50)
   >  exp(x)
   [1] 1.000e+00 1.928749847963918e-22
   >  f0(x)
   [1] 0
   >  f1(x)
   [1] 0
   >  f2(x) # log(1+epsilon) =~ epsilon
   [1] 1.928749847963918e-22
I don't think f2 should ever be less accurate.

expm1(x) is the inverse of log1p(x): it gives exp(x)-1
accurately for small abs(x).

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


If this is not sufficient, consider using CRAN package Rmpfr with
extended arithmetic.

Petr Savicky.

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

[R] Using Weights in R

[Krishnan Viswanathan]

I am new to R. I have a data set like this (given below is a fictional
dataset):

AgeCat FINWT
1 98
2 62
1 75
3 39
4 28
2 47
2 66
4 83
1 19
3 50

I need to calculate the weighted distribution of the variable AgeCat.

In SAS i can do:

proc freq data=ageval;
   tables agecat;
weight finwt;
run;

What or is there an equivalent in R?

TIA,
Krishnan

-- 
Krishnan Viswanathan
1101 High Meadow Dr
Tallahassee FL 32311

[[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] RCurl HTTP Post ?

[Hasan Diwan]

According to [1] and [2], using RCurl to post a form with basic
authentication is done using the postForm method. I'm trying to post
generated interpolation data from R onto an HTTP form. The call I'm using is
page <- postForm('http://our.server.com/dbInt/new', opts =
curlOptions=(userpwd="test:test", verbose=T), profileid = "-1",
value="1.801", type="history"). The page instance shows the HTTP response
500 screen and I get a nullpointerexception in the server logs. The line it
points to is dealing with getting an integer out of "profileid". Help?
Many thanks in advance...
-- 
Hasan Diwan
Developer
Economic Risk Management LLC
1. http://cran.r-project.org/web/packages/RCurl/RCurl.pdf
2. http://www.omegahat.org/RCurl/FAQ.html

[[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] sort by column and row names

[Jim Moon]

Yes, that is it.  Thank you very much, Phil. 
 

-Original Message-
From: Phil Spector [mailto:spec...@stat.berkeley.edu] 
Sent: Thursday, February 17, 2011 3:36 PM
To: Jim Moon
Cc: r-help@r-project.org
Subject: Re: [R] sort by column and row names

Jim -
I believe

dat[order(rownames(dat)),order(colnames(dat))]

is what you are looking for.
- Phil Spector
 Statistical Computing Facility
 Department of Statistics
 UC Berkeley
 spec...@stat.berkeley.edu


On Thu, 17 Feb 2011, Jim Moon wrote:

> Hello, All,
>
> How can one sort on column and row names.  For example:
>
> How can this
> X1   X3   X2
> X1   1 0  0
> X3   0 1  0
> X2   0 0  1
>
> become this?
> X1   X2   X3
> X1   1 0  0
> X2  0 1  0
> X3   0 0  1
>
>
> Thank you for your time!
>
> Jim
>
>
>   [[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] logsumexp function in R

[William Dunlap]

> -Original Message-
> From: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-project.org] On Behalf Of Petr Savicky
> Sent: Wednesday, February 16, 2011 10:46 PM
> To: r-help@r-project.org
> Subject: Re: [R] logsumexp function in R
> 
> On Wed, Feb 16, 2011 at 04:14:38PM -0500, Brian Tsai wrote:
> > Hi,
> > 
> > I was wondering if anyone has implemented a numerically 
> stable function to
> > compute log(sum(exp(x))) ?
> 
> Try the following
> 
>log(sum(exp(x - max(x + max(x)

Sometimes the log1p(x) function, which gives log(1+x)
accurately for small abs(x), helps a little more.
Compare the following 3 functions, which I think give
the same thing for 'ordinary' values of x:

f0 <- function(x) log(sum(exp(x)))
f1 <- function(x) log(sum(exp(x - max(x + max(x)
f2 <- function(x) { x <- sort(x) # mainly so max(x)==x[length(x)]
n <- length(x)
log1p(sum(exp(x[-n] - x[n]))) + x[n]
  }
But for the following x f2 gives a more accurate result
than f1, which in turn often gives a more accurate result
than f0:
  > x <- c(0, -50)
  > exp(x)
  [1] 1.000e+00 1.928749847963918e-22
  > f0(x)
  [1] 0
  > f1(x)
  [1] 0
  > f2(x) # log(1+epsilon) =~ epsilon
  [1] 1.928749847963918e-22
I don't think f2 should ever be less accurate.

expm1(x) is the inverse of log1p(x): it gives exp(x)-1
accurately for small abs(x).

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com  

> 
> If this is not sufficient, consider using CRAN package Rmpfr with
> extended arithmetic.
> 
> Petr Savicky.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] sort by column and row names

[Phil Spector]

Jim -
   I believe

dat[order(rownames(dat)),order(colnames(dat))]

is what you are looking for.
- Phil Spector
 Statistical Computing Facility
 Department of Statistics
 UC Berkeley
 spec...@stat.berkeley.edu


On Thu, 17 Feb 2011, Jim Moon wrote:


Hello, All,

How can one sort on column and row names.  For example:

How can this
X1   X3   X2
X1   1 0  0
X3   0 1  0
X2   0 0  1

become this?
X1   X2   X3
X1   1 0  0
X2  0 1  0
X3   0 0  1


Thank you for your time!

Jim


[[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] sort by column and row names

[Jim Moon]

Hello, All,

How can one sort on column and row names.  For example:

How can this
X1   X3   X2
X1   1 0  0
X3   0 1  0
X2   0 0  1

become this?
X1   X2   X3
X1   1 0  0
X2  0 1  0
X3   0 0  1


Thank you for your time!

Jim


[[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] convert the sas file into csv in R

[Joshua Wiley]

Hi yf,

On Thu, Feb 17, 2011 at 2:35 PM, yf  wrote:
>
> i am trying to convert sas file into csv. I used write.csv(a,
> file="cool.csv") but nothing come out. i don't know why. Thanks.
>
>
> library(Hmisc)
> a<- sasxport.get("C:\\Users\.")

Let's back up a little bit.  What type of file are you trying to read
into R with this statement?  Did you make sure it was a SAS XPORT
format?  SAS files often come in the native .sas7bdat format, which R
cannot read (it is proprietary).

At this point, try taking a look at your object "a".  At the console
you can type:

class(a)
str(a)

and see what the output says.  Is it a data frame?  A list? or  ?
Even if the data is read into R properly, write.csv() may not be able
to properly write all classes of data.  You might also try:

print(a)

to make sure there is actually data to be written and it is not just
an empty dataset.

> write.csv(a, file="cool.csv")

Assuming "a" is a matrix or data frame or can readily be coerced to
such, this line should be fine.

Just as a note, if it is a .sas7bdat file, Revolution Analytics now
provides software that can read it into R.  There is normally a
charge, but it is free for academics, which, judging by your email
address, you are.
http://www.revolutionanalytics.com/news-events/news-room/2011/revolution-analytics-unlocks-sas-data-for-r.php

Cheers,

Josh

> --
> View this message in context: 
> http://r.789695.n4.nabble.com/convert-the-sas-file-into-csv-in-R-tp3311769p3311769.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.

-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.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] convert the sas file into csv in R

[Nordlund, Dan (DSHS/RDA)]

> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
> project.org] On Behalf Of yf
> Sent: Thursday, February 17, 2011 2:35 PM
> To: r-help@r-project.org
> Subject: [R] convert the sas file into csv in R
> 
> 
> i am trying to convert sas file into csv. I used write.csv(a,
> file="cool.csv") but nothing come out. i don't know why. Thanks.
> 
> 
> library(Hmisc)
> a<- sasxport.get("C:\\Users\.")
> write.csv(a, file="cool.csv")

What do you mean "nothing come out"?  Where did you look for the file?  It 
should have been written to your current working directory.  You might try 
specifying a complete path so you know where you expect to find the CSV file.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] convert the sas file into csv in R

[yf]

i am trying to convert sas file into csv. I used write.csv(a,
file="cool.csv") but nothing come out. i don't know why. Thanks.


library(Hmisc)
a<- sasxport.get("C:\\Users\.")
write.csv(a, file="cool.csv")
-- 
View this message in context: 
http://r.789695.n4.nabble.com/convert-the-sas-file-into-csv-in-R-tp3311769p3311769.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] Problem with writing a file in UTF-8

[tpklein]

Hello,

I am working with a data frame containg character strings with many special
symbols from various European languages.  When writing such character
strings to a file using the UTF-8 encoding, some of them are converted in a
strange way.  See the following example, run in R 2.12.1 on Windows 7:

out <- file( description="out.txt", open="w", encoding="UTF-8")
write( x="äöüßæűŁ", file=out )
close( con=out )

The last two symbols in the character string are converted to "uL" while all
other characters are not changed (which is what I want).  How to explain
this?  Does it have something to do with my locale?  And is there a way to
work around this problem? -- Any help would be greatly appreciated.

Thomas
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Problem-with-writing-a-file-in-UTF-8-tp3311721p3311721.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] cv.glmnet errors

[Brian Tsai]

Hi,

I am trying to do multinomial regression using the glmnet package, but the
following gives me an error (for no reason apparent to me):

library(glmnet)
cv.glmnet(x=matrix(c(1,2,3,4,5,6,1,2,3,4,5,6),
nrow=6),y=as.factor(c(1,2,1,2,3,3)),family='multinomial',alpha=0.5,
nfolds=2)

The error i get is:
Error in if (outlist$msg != "Unknown error") return(outlist) :
  argument is of length zero


If i change the number of folds to 1, i get a seg fault:
 *** caught segfault ***
address 0x0, cause 'memory not mapped'

Traceback:
 1: .Fortran("lognet", parm = alpha, nobs, nvars, nc, as.double(x), y,
offset, jd, vp, ne, nx, nlam, flmin, ulam, thresh, isd, maxit, kopt, lmu
= integer(1), a0 = double(nlam * nc), ca = double(nx * nlam * nc),
ia = integer(nx), nin = integer(nlam), nulldev = double(1), dev =
double(nlam), alm = double(nlam), nlp = integer(1), jerr = integer(1),
PACKAGE = "glmnet")
 2: lognet(x, is.sparse, ix, jx, y, weights, offset, type, alpha, nobs,
nvars, jd, vp, ne, nx, nlam, flmin, ulam, thresh, isd, vnames, maxit,
HessianExact, family)
 3: glmnet(x[!which, ], y_sub, lambda = lambda, offset = offset_sub,
weights = weights[!which], ...)
 4: cv.glmnet(x = matrix(c(1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6), nrow =
6), y = as.factor(c(1, 2, 1, 2, 3, 3)), family = "multinomial",
alpha = 0.5, nfolds = 1)

Possible actions:



any ideas?


Brian.

[[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] Integrate with an indicator function

[Dennis Murphy]

Hi Hannah:

You have a few things going on, but the bottom line is that numer and denom
are both double integrals.

On Thu, Feb 17, 2011 at 1:06 PM, li li  wrote:

> Hi all,
>   I have some some problem with regard to finding the integral of a
> function containing an indicator function.
> please see the code below:
>
>
> func1 <- function(x, mu){
>  (mu^2)*dnorm(x, mean = mu, sd = 1)*dgamma(x, shape=2)}
>

curve(func1(x, 10), 0, 20)
curve(func1(x, 5), 0, 20)

both yield reasonable plots, so this is a function of one variable when mu
is supplied. If you wanted to plot this as a function of mu, you could get a
single curve by fixing x or integrating over x, which is what m1star appears
to be doing.

m1star <- function(x){
>   integrate(func1, lower = 0, upper = Inf,x)$val}
>

u <- seq(0, 20, 0.05)
v <- sapply(u, m1star)
plot(u, v, type = 'l')

yields a plot of m1star, which appears to marginalize func1 to make it a
function of mu, from what I can tell.

T <- function(x){
>   0.3*dnorm(x)/(0.3*dnorm(x)+0.7*m1star(x))}
>

plot(u, sapply(u, T), type = 'l')

yields a plot of T, but having to use sapply() indicates that T also
marginalizes a 2D function.


> func2 <- function(x,c){(T(x) <=c)*0.3*dnorm(x)}
> func3 <- function(x,c){(T(x) <= c)*(0.3*dnorm(x)+0.7*m1star(x))}
>

func2 uses T, which in turn uses m1star, so func2 is a marginalization of a
2D function whose support is on T(x) <= c. To integrate func2, you have to
do the integration in m1star first, so basically you have a double integral
to evaluate in numer. The same problem occurs in func3, since m1star() is a
part of both T() and the convex combination. Therefore, denom also requires
double integration.


> numer <- function(c){
>   integrate(func2, -Inf, Inf, c)$val}
>
> denom <- function(c){
>   integrate(func3, lower, Inf,c)$val}
>
> Look into cubature or Rcuba for two packages that are capable of performing
numerical double integration. An alternative solution is to use approxfun()
to create a function object from evaluations of m1star and T, and then use
the approxfun()s as the functions to be integrated in numer and denom. If
you decide to go the approxfun route, make sure to make enough evaluations
to reasonably capture the curvature in both m1star and T.

HTH,
Dennis

>
>
> The error message is as below :
>
> > numer(0.5)
> Error in integrate(func1, lower = 0, upper = Inf, x) :
>  the integral is probably divergent
>
>[[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] Populate a list / recursively set list values to NA

[Gene Leynes]

Works perfectly, thank you!

On Thu, Feb 17, 2011 at 4:11 PM, Erik Iverson  wrote:

> Gene,
>
> ?rapply is a recursive version of ?lapply, and should work.
>
> rapply(masterlist, function(x) x*NA, how = "replace")
>
>
> --Erik
>
> Gene Leynes wrote:
>
>> Hello all,
>>
>> Maybe I'm being thick, but I was trying to figure out a simple way to
>> create
>> a list with the same dimension of another list, but populated with NA
>> values.
>>
>>masterlist = list(
>>aa=list(
>>a=matrix(rnorm(100),10,10),
>>b=matrix(rnorm(130),13,10),
>>c=matrix(rnorm(140),14,10)),
>>bb=list(
>>a=matrix(rnorm(100),10,10),
>>b=matrix(rnorm(110),11,10),
>>c=matrix(rnorm(120),12,10)),
>>cc=list(
>>a=matrix(rnorm(100),10,10),
>>b=matrix(rnorm(100),10,10),
>>c=matrix(rnorm(100),10,10)),
>>dd=matrix(rnorm(225),15,15))
>>str(masterlist)
>>## Yay!  A new list that's just like the master list!
>>newlist = masterlist
>>## BOO!  Can't just set all the values to NA!
>>newlist = masterlist * NA
>>
>>## This works... but requires knowledge of the structure
>>newlist = masterlist
>>newlist[['aa']] = lapply(masterlist[['aa']], function(x) x * NA)
>>newlist[['bb']] = lapply(masterlist[['bb']], function(x) x * NA)
>>newlist[['cc']] = lapply(masterlist[['cc']], function(x) x * NA)
>>newlist[['dd']] = masterlist[['dd']] * NA
>>str(newlist)
>>
>> At this point I thought a recursive function would be a good idea, but I
>> wasn't sure how to keep the nice list names.
>>
>> I feel sure that this must be something that others have already done.
>>
>> Thank you!!!
>>
>> PS:
>> I found an email in the archive about it, but that email un-helpfully
>> instructed the user to LEARN HOW TO USE EMAIL.  Not sure how to use email
>> to
>> initialize a list. so if someone could explain to me how to use emails to
>> initialize / populate a list, that would be cool also.
>> Also, I did type ?list and ?lapply, and I did read the help files.  In
>> fact,
>> I learned about the very helpful vector("list", length) command yesterday!
>> Pretty cool!
>> I'm not opposed to looking at the help, but please provide some hint about
>> what direction to take with the help.  Even a small hint will do!!
>>
>> If you want to yell at me about reading the help or whatever, that's fine.
>> But please do so offline so that the archives will be more valuable for
>> future users.
>>
>>[[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] tikzDevice compiling problem

[Yihui Xie]

I don't think tikz() can capture the plots in rgl. To export rgl
graphics, you need to use rgl.snapshot() or rgl.postscript().

Regards,
Yihui
--
Yihui Xie 
Phone: 515-294-2465 Web: http://yihui.name
Department of Statistics, Iowa State University
2215 Snedecor Hall, Ames, IA



On Wed, Feb 16, 2011 at 9:19 AM, cuass  wrote:
>
> Hi guys.
> The tex file which compile the following graphic has many problems. I don't
> know what's happening. Any input would be really appreciated.
>
> # Cargo el prgrama que produce el código en Latex
> require(tikzDevice)
>
> # Establezco directorio del programa
> setwd('/Users/fabiangarcia/Documents/R')
>
> # El siguiente programa produce el archivo tex
> tikz('CobbGRAF.tex', standAlone = TRUE, width=5, height=5)
>
> # La gráfica de la función de utilidad
> f = function(x, y) ((y)^1*(x)^1)
> x = seq(0,5,len=40)
> y = seq(0,5,len=40)
> z = outer(x, y, f)
>  showsurface = function(x, y, z)
>  persp3d(x,y,z, col="blue", alpha=0.3, axes=
>  F)+{
>  contours = contourLines(x,y,z)
>  for (i in 1:length(contours)) {
>  with(contours[[i]], lines3d(x, y, level, col="darkred"))
>  }
>  }
>  open3d()
>  showsurface(x,y,z)
>
>  # Cierro el device
> dev.off()
>
> # Compilo el archivo tex
> tools::texi2dvi('CobbGRAF.tex',pdf=F)
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/tikzDevice-compiling-problem-tp3309028p3309028.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

Re: [R] using the bootstrap method for latent class analysis

[David Duffy]

My main question is: can the bootstrap procedure sometimes be less 
precise than the non-bootstrap procedure when carrying out latent class 
analysis?


I am asking this question because I carried out some analyses and when I 
*did not* use the bootstrap procedure I found adequate fit for one of 
the models but when I did use the bootstrap none of the models indicated

adequate fit to the data.


It should not usually lead to this behaviour, though I am curious what 
model fit measure you are using in the bootstrap situation.  You need to 
tell us what packages you are using, and send us a small self-contained

example.

Cheers, David Duffy.

--
| David Duffy (MBBS PhD) ,-_|\
| email: dav...@qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  / *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Populate a list / recursively set list values to NA

[Erik Iverson]

Gene,

?rapply is a recursive version of ?lapply, and should work.

rapply(masterlist, function(x) x*NA, how = "replace")


--Erik

Gene Leynes wrote:

Hello all,

Maybe I'm being thick, but I was trying to figure out a simple way to create
a list with the same dimension of another list, but populated with NA
values.

masterlist = list(
aa=list(
a=matrix(rnorm(100),10,10),
b=matrix(rnorm(130),13,10),
c=matrix(rnorm(140),14,10)),
bb=list(
a=matrix(rnorm(100),10,10),
b=matrix(rnorm(110),11,10),
c=matrix(rnorm(120),12,10)),
cc=list(
a=matrix(rnorm(100),10,10),
b=matrix(rnorm(100),10,10),
c=matrix(rnorm(100),10,10)),
dd=matrix(rnorm(225),15,15))
str(masterlist)
## Yay!  A new list that's just like the master list!
newlist = masterlist
## BOO!  Can't just set all the values to NA!
newlist = masterlist * NA

## This works... but requires knowledge of the structure
newlist = masterlist
newlist[['aa']] = lapply(masterlist[['aa']], function(x) x * NA)
newlist[['bb']] = lapply(masterlist[['bb']], function(x) x * NA)
newlist[['cc']] = lapply(masterlist[['cc']], function(x) x * NA)
newlist[['dd']] = masterlist[['dd']] * NA
str(newlist)

At this point I thought a recursive function would be a good idea, but I
wasn't sure how to keep the nice list names.

I feel sure that this must be something that others have already done.

Thank you!!!

PS:
I found an email in the archive about it, but that email un-helpfully
instructed the user to LEARN HOW TO USE EMAIL.  Not sure how to use email to
initialize a list. so if someone could explain to me how to use emails to
initialize / populate a list, that would be cool also.
Also, I did type ?list and ?lapply, and I did read the help files.  In fact,
I learned about the very helpful vector("list", length) command yesterday!
Pretty cool!
I'm not opposed to looking at the help, but please provide some hint about
what direction to take with the help.  Even a small hint will do!!

If you want to yell at me about reading the help or whatever, that's fine.
But please do so offline so that the archives will be more valuable for
future users.

[[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] which functions are being debugged?

[Shi, Tao]

Hi Henrique,

Thanks for the reply!

I posted the question in a hurry yesterday.  After a bit research today, I 
found 
that back in 2009, Andrew actually asked the same question and there was a good 
discussion about it:

http://r.789695.n4.nabble.com/how-do-you-know-which-functions-are-being-debugged-td906109.html


Steve's approach of keeping a record of the debugged functions yourself is more 
proper, but your function works well for me for now.  


Thanks!

...Tao






From: Henrique Dallazuanna 

Cc: r-help@r-project.org
Sent: Wed, February 16, 2011 4:56:01 PM
Subject: Re: [R] which functions are being debugged?

Try this:

isDebug <- lapply(sapply(search(), lapply, ls), function(f)sapply(f, 
function(.f)tryCatch(isdebugged(.f), error = function(...)FALSE)))
which(unlist(isDebug))




Hi list,
>
>Is there a function that can let me know which functions are being debugged?  I
>know I'm probably not doing a very good job of keeping track of things, but it
>does get messier when you dig into different layers of a function.  I know
there
>is "isdebugged", but it only works on one function at a time.
>
>Thanks!
>
>...Tao
>
>
>
>
>   [[alternative HTML version deleted]]
>
>__
>R-help@r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>and provide commented, minimal, self-contained, reproducible code.
>


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



  
[[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] RGtk2 on Debian Testing

[Lorenzo Isella]

Dear All,
I am running Debian testing on my system for the amd64 architecture,
When trying to install the RGtk package I get this error


> install.packages('RGtk2')
Installing package(s) into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
trying URL 
'http://rm.mirror.garr.it/mirrors/CRAN/src/contrib/RGtk2_2.20.8.tar.gz'

Content type 'application/x-gzip' length 2637806 bytes (2.5 Mb)
opened URL
==
downloaded 2.5 Mb

* installing *source* package ‘RGtk2’ ...
checking for pkg-config... /usr/bin/pkg-config
checking pkg-config is at least version 0.9.0... yes
checking for INTROSPECTION... no
checking for GTK... no
configure: error: GTK version 2.8.0 required
ERROR: configuration failed for package ‘RGtk2’
* removing ‘/usr/local/lib/R/site-library/RGtk2’

The downloaded packages are in
‘/tmp/RtmpMTHLGF/downloaded_packages’
Warning message:
In install.packages("RGtk2") :
  installation of package 'RGtk2' had non-zero exit status

Does anyone know why there is a mismatch between my GTK and the one 
required by R?
Should I enable some particular R repositories (I know that the previous 
Debian testing was released a few days ago, but I do not know if this is 
relevant).

Any suggestion is welcome.
Cheers

Lorenzo

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Populate a list / recursively set list values to NA

[Gene Leynes]

Hello all,

Maybe I'm being thick, but I was trying to figure out a simple way to create
a list with the same dimension of another list, but populated with NA
values.

masterlist = list(
aa=list(
a=matrix(rnorm(100),10,10),
b=matrix(rnorm(130),13,10),
c=matrix(rnorm(140),14,10)),
bb=list(
a=matrix(rnorm(100),10,10),
b=matrix(rnorm(110),11,10),
c=matrix(rnorm(120),12,10)),
cc=list(
a=matrix(rnorm(100),10,10),
b=matrix(rnorm(100),10,10),
c=matrix(rnorm(100),10,10)),
dd=matrix(rnorm(225),15,15))
str(masterlist)
## Yay!  A new list that's just like the master list!
newlist = masterlist
## BOO!  Can't just set all the values to NA!
newlist = masterlist * NA

## This works... but requires knowledge of the structure
newlist = masterlist
newlist[['aa']] = lapply(masterlist[['aa']], function(x) x * NA)
newlist[['bb']] = lapply(masterlist[['bb']], function(x) x * NA)
newlist[['cc']] = lapply(masterlist[['cc']], function(x) x * NA)
newlist[['dd']] = masterlist[['dd']] * NA
str(newlist)

At this point I thought a recursive function would be a good idea, but I
wasn't sure how to keep the nice list names.

I feel sure that this must be something that others have already done.

Thank you!!!

PS:
I found an email in the archive about it, but that email un-helpfully
instructed the user to LEARN HOW TO USE EMAIL.  Not sure how to use email to
initialize a list. so if someone could explain to me how to use emails to
initialize / populate a list, that would be cool also.
Also, I did type ?list and ?lapply, and I did read the help files.  In fact,
I learned about the very helpful vector("list", length) command yesterday!
Pretty cool!
I'm not opposed to looking at the help, but please provide some hint about
what direction to take with the help.  Even a small hint will do!!

If you want to yell at me about reading the help or whatever, that's fine.
But please do so offline so that the archives will be more valuable for
future users.

[[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] barplot with errorbars

[Marc Schwartz]

Aldi,

Thanks for the reference to the R News article. That seems like a lifetime 
ago...  :-)

As an exercise in learning to use R's base graphics and some of the concepts of 
manipulating and annotating base R plots, it is a decent intro.

That being said, unless pressured (eg. weapon pointing at cranium), I would 
never actually use that particular format (the barplots) today in presenting 
data. The same goes for the barplot2() function in the gplots CRAN package, 
which I wrote, but have literally not used in years. That function automates 
barplots and confidence interval plotting, amongst other enhancements (eg. log 
scaled axes), to the base barplot() function.

If I may, allow me to offer some additional comments to Maria.

First, note that barplots are best used to present discrete counts or perhaps 
percentages, not continuous data. The y axis (in the case of a vertical 
barplot) is intended to start at 0 to allow the appropriate perception of 
count/percentage based differences within or between groups, as represented by 
the heights of the bars. 

The addition of confidence intervals to barplots, especially only the upper 
bound, which are also known as dynamite plots 
(http://biostat.mc.vanderbilt.edu/wiki/Main/DynamitePlots), are generally 
frowned upon these days, though there are communities that still prefer that 
format for count data.

For continuous data however, point plots with confidence intervals are much 
preferred. Other alternatives include boxplots and violin plots, if one is 
wanting to visually present more information on the distribution of a 
continuous variable between groups, beyond point estimates for the means and 
CIs.

Frank has a nice intro document on statistical graphics here:

  
http://biostat.mc.vanderbilt.edu/twiki/pub/Main/StatGraphCourse/graphscourse.pdf

that you might want to review.

I would urge Maria to consider alternatives to using a barplot to present this 
data. There are various custom functions to aid with that (searching the list 
archives using rseek.org should yield good results) and the ?arrows and 
?segments functions typically serve as the basis for drawing the CIs around the 
points representing the means.

HTH,

Marc Schwartz



On Feb 17, 2011, at 2:38 PM, Aldi Kraja wrote:

> Hi Toby and Maria,
> 
> I did a check on Toby's suggestion and it is not there:
> 
> R version 2.12.1 (2010-12-16)
>> ??barploterrbar
> No help files found with alias or concept or title matching
> 'barploterrbar' using fuzzy matching.
> 
> Also I went to the following location which does not exist.
> http://www.dkfz.de/abt0840/whuber
> 
> To Maria: here is a suggestion: to do the barplot with confidence 
> intervals / errorbars:
> 
> look at the following R news under www.r-project.org
> "An introduction to using R's base graphics." (see following citation)
> on page 4 you will find explanations and code how to build bars with se.
> 
> [62]  Marc Schwartz. R Help Desk: An introduction to using R's base 
> graphics. /R News/, 3(2):2-6, October 2003. [ bib | http 
>  | .pdf 
>  ]
> 
> 
> HTH,
> 
> Aldi
> 
> On 2/17/2011 10:08 AM, Toby Marthews wrote:
>> If you google "barplot with error bars" you immediately find 
>> http://svitsrv25.epfl.ch/R-doc/library/prada/html/barploterrbar.html .
>> Toby.
>> 
>> 
>> From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] On Behalf 
>> Of Lathouri, Maria [m.lathour...@imperial.ac.uk]
>> Sent: 17 February 2011 16:00
>> To: r-help@r-project.org
>> Subject: [R] barplot with errorbars
>> 
>> Dear all
>> 
>> I have six variables of the average metal concentrations
>> 
>> Var1 4.77
>> Var2 23.5
>> Var3 5.2
>> Var4 12.3
>> Var5 42.1
>> Var6 121.2
>> 
>> I want to plot them as a barplot with error bars. Could you help me?
>> 
>> Cheers
>> Maria

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] need tick marks by 10s or 20s up to 200 on x axis

[Joshua Wiley]

Hi Julie,

Here is a little example.  First I make a histogram (it looks odd
because there isn't much data), then add the x axis in manually
specifying at what points to place tick marks.

## histogram suppressing x axis
hist(c(0, 200), xaxt = "n")
## add axis with custom ticks
axis(1, at = seq(0, 200, 20))

## For documentation see
?par
?axis

## the documentation for par is long and is never fun to go through,
but it really explains a lot and R graphics will seem less mysterious
if you schedule an hour to sit down and read through it

Cheers,

Josh


On Thu, Feb 17, 2011 at 12:58 PM, Julie McWhorter
 wrote:
>
> I haven't been able to get any code to make more tick marks than the default 
> by 50.
> Thanks!  Code:
> hist(OMY$FL, main = "2010 Oncorhynchus mykiss Fork Length Frequencies at Buck 
> Creek Reach 1", include.lowest = TRUE, col = "blue", border = "white", breaks 
> = 140,xlab = "Fork Length  in  mm", ylab = "Frequency", xlim = range(0,200), 
> ylim = range (0,15), xbar=10)
>        [[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.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.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] sampling

[Dennis Murphy]

Hi:

A couple more approaches to consider:

# Utility function to extract two rows from a data frame
# Meant to be applied to each data subset
sampler <- function(d) if(nrow(d) > 2) d[sample(1:nrow(d), 2, replace =
FALSE), ] else d

library(plyr)
> ddply(x, 'id', sampler)
  id v1 V2
1  1  2 13
2  1  1 12
3  2  4 15
4  2  6 17
5  3  8 19
6  3 10 21
7  4 11 22
8  4 12 23

library(data.table)
dtx <- data.table(x, key = 'id')
> dtx[, sampler(.SD), by = 'id']
 id v1 V2
[1,]  1  1 12
[2,]  1  3 14
[3,]  2  5 16
[4,]  2  7 18
[5,]  3  9 20
[6,]  3 10 21
[7,]  4 11 22
[8,]  4 12 23


HTH,
Dennis


On Wed, Feb 16, 2011 at 8:35 PM, yf  wrote:

>
> I want to sample from the ID. For each ID, i want to have 2 set of data. I
> try the sample() function but it didn't work.
>
> > x<-data.frame(id=c(1,1,1,2,2,2,2,3,3,3,4,4), v1=c(1:12), V2=c(12:23))
> > x
>   id v1 V2
> 1   1  1 12
> 2   1  2 13
> 3   1  3 14
> 4   2  4 15
> 5   2  5 16
> 6   2  6 17
> 7   2  7 18
> 8   3  8 19
> 9   3  9 20
> 10  3 10 21
> 11  4 11 22
> 12  4 12 23
> --
> View this message in context:
> http://r.789695.n4.nabble.com/sampling-tp3310184p3310184.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] need tick marks by 10s or 20s up to 200 on x axis

[Julie McWhorter]

I haven't been able to get any code to make more tick marks than the default by 
50.
Thanks!  Code:
hist(OMY$FL, main = "2010 Oncorhynchus mykiss Fork Length Frequencies at Buck 
Creek Reach 1", include.lowest = TRUE, col = "blue", border = "white", breaks = 
140,xlab = "Fork Length  in  mm", ylab = "Frequency", xlim = range(0,200), ylim 
= range (0,15), xbar=10)
[[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] need tick marks by 10s or 20s up to 200 on x axis

[juliemcwh]

Sorry, I have read other posts on tick marks but nothing is working.

Here is my R code.  I was able to turn off the x axis at one point (I
removed that code), but never got the tick marks on by 20s.  It always goes
to the default of 50, 100, 150, 200.

Thanks.

R code: hist(OMY$FL, main = "2010 Oncorhynchus mykiss Fork Length
Frequencies at Buck Creek Reach 1", include.lowest = TRUE, col = "blue",
border = "white", breaks = 140,
xlab = "Fork Length  in  mm", ylab = "Frequency", xlim = range(0,200), ylim
= range (0,15), xbar=10)
-- 
View this message in context: 
http://r.789695.n4.nabble.com/need-tick-marks-by-10s-or-20s-up-to-200-on-x-axis-tp3311627p3311627.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] Credible intervals for Bayesian Model Averaging using the BAS package

[matthewnkim]

Dear all,

I am using the BAS package to make predictions using the predict.bma
function.  Does anyone know how to calculate credible intervals for
predictions from objects of class bma?  It seems like all of the
required information is available from the bma class objects and/or
output from the predict.bma function.

Sincerely,

Matt Bowser

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] summing 15 minute precip data to daily

[Gabor Grothendieck]

On Thu, Feb 17, 2011 at 3:39 PM, Stephen Sefick  wrote:
> Janet:
>
> The zoo package and aggregate.zoo should do the trick.  I have done this
> many times with these tools.

Here is some code. Suppose precip is the data.  We convert it to a zoo
object with date/times assuming the first one is at the start of
2000-01-01.  Then we aggregate them by day:

> library(zoo)
> library(chron)
> precip <- 1:200 # sample data
> z <- zooreg(precip, start = as.chron("2000-01-01 00:00:00"), freq = 96)
> z.ag <- aggregate(z, as.Date)
> z.ag
2000-01-01 2000-01-02 2000-01-03
  4656  13872   1572


-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

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

Re: [R] barplot with errorbars

[William Revelle]

As the FAQ 
(http://cran.r-project.org/doc/FAQ/R-FAQ.html#How-can-I-put-error-bars-or-confidence-bands-on-my-plot_003) 
says,

there are about 5 different packages that will draw error bars.



At 2:38 PM -0600 2/17/11, Aldi Kraja wrote:

Hi Toby and Maria,

I did a check on Toby's suggestion and it is not there:

R version 2.12.1 (2010-12-16)
 > ??barploterrbar
No help files found with alias or concept or title matching
'barploterrbar' using fuzzy matching.

Also I went to the following location which does not exist.
http://www.dkfz.de/abt0840/whuber

To Maria: here is a suggestion: to do the barplot with confidence
intervals / errorbars:

look at the following R news under www.r-project.org
"An introduction to using R's base graphics." (see following citation)
on page 4 you will find explanations and code how to build bars with se.

[62]Marc Schwartz. R Help Desk: An introduction to using R's base
graphics. /R News/, 3(2):2-6, October 2003. [ bib | http
 | .pdf
 ]


HTH,

Aldi

On 2/17/2011 10:08 AM, Toby Marthews wrote:
 If you google "barplot with error bars" you immediately find 
http://svitsrv25.epfl.ch/R-doc/library/prada/html/barploterrbar.html 
.

 Toby.

 
 From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] 
On Behalf Of Lathouri, Maria [m.lathour...@imperial.ac.uk]

 Sent: 17 February 2011 16:00
 To: r-help@r-project.org
 Subject: [R] barplot with errorbars

 Dear all

 I have six variables of the average metal concentrations

 Var1 4.77
 Var2 23.5
 Var3 5.2
 Var4 12.3
 Var5 42.1
 Var6 121.2

 I want to plot them as a barplot with error bars. Could you help me?

 Cheers
 Maria





--



[[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] Integration with an Indicator Function in R

[li li]

Hi all,
   I have some some problem with regard to finding the integral of a
function containing an indicator function.
please see the code below:

func1 <- function(x, mu){
  (mu^2)*dnorm(x, mean = mu, sd = 1)*dgamma(x, shape=2)}
m1star <- function(x){
   integrate(func1, lower = 0, upper = Inf,x)$val}
T <- function(x){
   0.3*dnorm(x)/(0.3*dnorm(x)+0.7*m1star(x))}

func2 <- function(x,c){(T(x) <=c)*0.3*dnorm(x)}
func3 <- function(x,c){(T(x) <= c)*(0.3*dnorm(x)+0.7*m1star(x))}
numer <- function(c){
   integrate(func2, -Inf, Inf, c)$val}

denom <- function(c){
   integrate(func3, lower, Inf,c)$val}



The error message is as below :

> numer(0.5)
Error in integrate(func1, lower = 0, upper = Inf, x) :
  the integral is probably divergent



Thank you very much!
   Hannah

[[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] Integrate with an indicator function

[li li]

Hi all,
   I have some some problem with regard to finding the integral of a
function containing an indicator function.
please see the code below:


func1 <- function(x, mu){
  (mu^2)*dnorm(x, mean = mu, sd = 1)*dgamma(x, shape=2)}
m1star <- function(x){
   integrate(func1, lower = 0, upper = Inf,x)$val}
T <- function(x){
   0.3*dnorm(x)/(0.3*dnorm(x)+0.7*m1star(x))}

func2 <- function(x,c){(T(x) <=c)*0.3*dnorm(x)}
func3 <- function(x,c){(T(x) <= c)*(0.3*dnorm(x)+0.7*m1star(x))}
numer <- function(c){
   integrate(func2, -Inf, Inf, c)$val}

denom <- function(c){
   integrate(func3, lower, Inf,c)$val}



The error message is as below :

> numer(0.5)
Error in integrate(func1, lower = 0, upper = Inf, x) :
  the integral is probably divergent

[[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] Data frame sampling

[Hosack, Michael]

R users,

I have been trying to write a program in R that will extract rows from a data 
frame and combine the rows 
into a new smaller data frame while meeting several criteria. I would greatly 
appreciate any advice 
that could help me get started down the right path. What I want to do is to 
extract two rows WEEK 
(26 weeks total) from the data frame by sampling without replacement from the 
variable STRATA_NUM (1:40) 
and using this output to index the data frame and extract rows, moving from one 
week to the next. I 
created STRATA_NUM to use as an index to identify unique sampling frame strata 
(e.g. unique DOW, SITE, 
and TOD permutations---40 total). This would result in the entire STRATA_NUM 
vector being sampled by 
the end of WEEK 20, at which point I would need the code to resample again 
without replacement from 
the full vector for the remaining 6 weeks (i.e. 12 more rows extracted). One 
final condition is that 
DOW's must not be duplicated for any given WEEK (e.g. I cannot sample two 
Tuesdays in any one week). 
The reasoning for my approach is that I want even coverage of all 40 unique 
spatiotemporal strata in 
my final data frame. 

Thank you,

Mike

### The following code extracts two rows per week randomly without replacement 
and ensures that no two 
### DOW are duplicated per week, but fails to incorporate even coverage of 
strata. 

set.seed(3000)
WEEKDAY_STRATA1 <- do.call("rbind", lapply(split(WEEKDAY_STRATA, 
WEEKDAY_STRATA$WEEK), function(x) x[sample(nrow(x), 10, replace=FALSE),])) 
l <- do.call("rbind", 
lapply(split(WEEKDAY_STRATA1,WEEKDAY_STRATA1$WEEK),function(y)!duplicated(y[,1])))
M <- as.matrix(l,byrow = FALSE)
WEEKDAY_STRATA2 <- WEEKDAY_STRATA1[as.vector(t(M)),]
WEEKDAY_STRATA3 <- do.call("rbind", lapply(split(WEEKDAY_STRATA2, 
WEEKDAY_STRATA2$WEEK), function(x) x[sample(nrow(x), 2, replace=FALSE),])) 
WEEKDAY_STRATA3


### I only included the first three weeks for obvious reasons

WEEKDAY_STRATA <-
structure(list(DOW = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L), .Label = c("Mon", "Tue", "Wed", "Thu", "Fri"), class = "factor"), 
SITE = c(101, 101, 102, 102, 103, 103, 104, 104, 101, 101, 
102, 102, 103, 103, 104, 104, 101, 101, 102, 102, 103, 103, 
104, 104, 101, 101, 102, 102, 103, 103, 104, 104, 101, 101, 
102, 102, 103, 103, 104, 104, 101, 101, 102, 102, 103, 103, 
104, 104, 101, 101, 102, 102, 103, 103, 104, 104, 101, 101, 
102, 102, 103, 103, 104, 104, 101, 101, 102, 102, 103, 103, 
104, 104, 101, 101, 102, 102, 103, 103, 104, 104, 101, 101, 
102, 102, 103, 103, 104, 104, 101, 101, 102, 102, 103, 103, 
104, 104, 101, 101, 102, 102, 103, 103, 104, 104, 101, 101, 
102, 102, 103, 103, 104, 104, 101, 101, 102, 102, 103, 103, 
104, 104), TOD = structure(c(2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L), .Label = c("Morn", "Aftn"
), class = "factor"), STRATA_NUM = c(2L, 1L, 12L, 11L, 22L, 
21L, 32L, 31L, 4L, 3L, 14L, 13L, 24L, 23L, 34L, 33L, 6L, 
5L, 16L, 15L, 26L, 25L, 36L, 35L, 8L, 7L, 18L, 17L, 28L, 
27L, 38L, 37L, 10L, 9L, 20L, 19L, 30L, 29L, 40L, 39L, 2L, 
1L, 12L, 11L, 22L, 21L, 32L, 31L, 4L, 3L, 14L, 13L, 24L, 
23L, 34L, 33L, 6L, 5L, 16L, 15L, 26L, 25L, 36L, 35L, 8L, 
7L, 18L, 17L, 28L, 27L, 38L, 37L, 10L, 9L, 20L, 19L, 30L, 
29L, 40L, 39L, 2L, 1L, 12L, 11L, 22L, 21L, 32L, 31L, 4L, 
3L, 14L, 13L, 24L, 23L, 34L, 33L, 6L, 5L, 16L, 15L, 26L, 
25L, 36L, 35L, 8L, 7L, 18L, 17L, 28L, 27L, 38L, 37L, 10L, 
9L, 20L, 19L, 30L, 29L, 40L, 39L), DATE = structure(c(14732, 
14732, 14732, 14732, 14732, 14732, 14732, 14732, 14733, 14733, 
14733, 14733, 14733, 14733, 14733, 14733, 14734, 14734, 14734, 
14734, 14734, 14734, 14734, 14734, 14735, 14735, 14735, 14735, 
14735, 14735, 14735, 14735, 14736, 14736, 14736, 14736, 14736, 
14736, 14736, 14736, 14739, 14739, 14739, 14739, 14739, 14739, 
14739, 14739, 14740, 14740, 14740, 14740, 14740, 14740, 14740, 
14740, 14741, 14741, 14741, 14741, 14741, 14741, 14741, 14741, 
14742, 14742, 14742, 14742, 14742, 14742

[R] [R-pkgs] New version of rms package on CRAN

[Frank Harrell]

A new version of rms is now available on CRAN for Linux/UNIX.  I expect  
Mac and Windows versions to be available in a day or so. This version  
works with and requires the newest version of Therneau's survival package.


More information is at http://biostat.mc.vanderbilt.edu/Rrms


Changes in version 3.2-0 (2011-02-14)
   * Changed to be compatible with survival 2.36-3 which is now required
   * Added logLik.rms and AIC.rms functions to be compatible with standard  
R

   * Fixed oos.loglik.Glm
   * Fixed bootcov related to nfit='Glm'
   * Fixed (probably) old bug in latexrms with strat predictors

Frank
--
Frank E Harrell Jr Professor and Chairman  School of Medicine
   Department of Biostatistics Vanderbilt University

___
R-packages mailing list
r-packa...@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-packages

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

Re: [R] A very basic line-plot question

[Phil Spector]

The reason that it's easy to produce dotcharts or barplots in R with
data like yours, and difficult to produce a line plot, is because 
dotcharts and barplots are appropriate for your data, whereas a line

plot is not.  Since the values on the x-axis represent names,
and those names are not measured on any kind of ordinal scale,
connecting the points on a plot would be at best meaningless and at 
worst deceptive.


Of course, one of R's charms is that it provides ample rope to hang
yourself it that's what you really want to do.  You could create a 
line plot with a custom axis labeled in what ever way you want, to 
produce a whatever sort of plot you wanted.  I'll leave it to you to 
figure out the details.

- Phil Spector
 Statistical Computing Facility
 Department of Statistics
 UC Berkeley
 spec...@stat.berkeley.edu


On Thu, 17 Feb 2011, world peace wrote:


Hi All

I have data like this

 tom   randy mike dan doug
height   150   152 155  134 141

I am trying to create a line plot, with names on X-axis and height measure
on Y. how can i get it through R.
I could get several versions which are close (dotchart, bargraph), but not
quite the same thing.

I am looking for something like
http://www.statsoft.com/Portals/0/blog/line_plot.jpg
with option to add more lines, such as weight measurement for above data.

Thanks,
A

[[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] barplot with errorbars

[Aldi Kraja]

Hi Toby and Maria,

I did a check on Toby's suggestion and it is not there:

R version 2.12.1 (2010-12-16)
 > ??barploterrbar
No help files found with alias or concept or title matching
'barploterrbar' using fuzzy matching.

Also I went to the following location which does not exist.
http://www.dkfz.de/abt0840/whuber

To Maria: here is a suggestion: to do the barplot with confidence 
intervals / errorbars:

look at the following R news under www.r-project.org
"An introduction to using R's base graphics." (see following citation)
on page 4 you will find explanations and code how to build bars with se.

[62]Marc Schwartz. R Help Desk: An introduction to using R's base 
graphics. /R News/, 3(2):2-6, October 2003. [ bib | http 
 | .pdf 
 ]


HTH,

Aldi

On 2/17/2011 10:08 AM, Toby Marthews wrote:
> If you google "barplot with error bars" you immediately find 
> http://svitsrv25.epfl.ch/R-doc/library/prada/html/barploterrbar.html .
> Toby.
>
> 
> From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] On Behalf 
> Of Lathouri, Maria [m.lathour...@imperial.ac.uk]
> Sent: 17 February 2011 16:00
> To: r-help@r-project.org
> Subject: [R] barplot with errorbars
>
> Dear all
>
> I have six variables of the average metal concentrations
>
> Var1 4.77
> Var2 23.5
> Var3 5.2
> Var4 12.3
> Var5 42.1
> Var6 121.2
>
> I want to plot them as a barplot with error bars. Could you help me?
>
> Cheers
> Maria
>
>
>

-- 



[[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] A very basic line-plot question

[Jinyan Huang]

dat<-c(150,152,155,134,141)
plot(dat,type="o",ylim=c(100,160),xlab="Names",ylab="Height")

On Thu, Feb 17, 2011 at 7:45 PM, world peace  wrote:
> Hi All
>
> I have data like this
>
>              tom   randy mike dan doug
> height       150   152     155  134 141
>
> I am trying to create a line plot, with names on X-axis and height measure
> on Y. how can i get it through R.
> I could get several versions which are close (dotchart, bargraph), but not
> quite the same thing.
>
> I am looking for something like
> http://www.statsoft.com/Portals/0/blog/line_plot.jpg
> with option to add more lines, such as weight measurement for above data.
>
> Thanks,
> A
>
>        [[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] sampling

[David Winsemius]

On Feb 17, 2011, at 1:33 PM, andrija djurovic wrote:


This is, maybe, not the best solution but I hope it will help you:

x<-data.frame(id=c(1,1,1,2,2,2,2,3,3,3,4,4), v1=c(1:12), V2=c(12:23))

do.call("rbind",by(x,x$id,function(x) x[c(sample(nrow(x),2)),]))

Andrija



Another way (and note that by is just a wrppare for tapply):

> tapply(1:nrow(x), x$id, sample, 2)
$`1`
[1] 2 3

$`2`
[1] 5 4

$`3`
[1] 10  8

$`4`
[1] 11 12

> x[unlist( tapply(1:nrow(x), x$id, sample, 2) ), ]
   id v1 V2
2   1  2 13
3   1  3 14
5   2  5 16
6   2  6 17
9   3  9 20
8   3  8 19
12  4 12 23
11  4 11 22



On Thu, Feb 17, 2011 at 6:39 PM, yf  wrote:



But i need for each id have two data.
Like...

x

 id v1 V2
1   1  1 12
2   1  2 13

4   2  4 15
5   2  5 16


8   3  8 19
9   3  9 20

11  4 11 22
12  4 12 23

So should write sample( if sample id >2  ,2). I don't know how to  
write

(if
sample id >2). Thanks.
--
View this message in context:
http://r.789695.n4.nabble.com/sampling-tp3310184p3311253.html
Sent from the R help mailing list archive at Nabble.com.

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



[[alternative HTML version deleted]]

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


David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] summing 15 minute precip data to daily

[Stephen Sefick]

Janet:

The zoo package and aggregate.zoo should do the trick.  I have done this
many times with these tools.
HTH,

Stephen

On Thu, 2011-02-17 at 11:56 -0800, Janet Choate wrote:
> Hi all,
> i'm sure there is an easy way to do this, but i'm stumped, so any help would
> be appreciated.
> 
> i have a single column of data for precipitation every 15 minutes over a
> year.  i want to sum the precip to daily data.
> so the first 96 records = the first day, the second 96 records = the second
> day, and so on
> is there a way to write a for loop that would sum the data to daily, and
> write each value to a second object so i end up with a file of daily precip?
> 
> thanx,
> Janet
> 
>   [[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] Urgent Request

[Ben Bolker]

muhammad mohsin  yahoo.com> writes:

> Hope you will be fine. I am student of Ph.D and doing '
> some work on distribution. 
> I developed a new distribution and having some problems in estimating their 
> parameters by MLE. I used R-program and  used "maxLik" function (maxLik: A 
> Package for Maximum Likelihood Estimation in R)
> But there is some problem, it is 
> not estimated the parameters properly. 
> I also write an e-mail to the author of 
> this paper but he could not solve my problem. 
> His function works well for simple 
> and known distribution but does not work for a new function.  
> Can anybody spare 
> some time for me?  I really need your help. 
> Please inform me so that I can send 
> you the material.
> Waiting for a quick reply 
> Best Regards,
> 
>  Muhammad Mohsin

  Dear Mr Mohsin,

  You are free to post a (self-contained/reproducible and preferably
small/minimal) example here and see if it interests anyone sufficiently
for them to volunteer time to see if they can find the problem.  With respect,
though, if you are a PhD student in statistics then this is part of your
training, and it should really fall to you, or to your supervisor or other
people at your institution, to work out how to solve it. You are (much)
more likely to get useful help from this group if you can narrow 
your problem down to a specific
point, and if you can indicate what steps you have tried to take to
solve your problem for yourself.  Maximum likelihood estimation is
in general a challenging computational problem -- just because a
general-purpose function or package exists doesn't mean it can solve
all problems easily. You may have to work harder to understand 
the particular structure of your optimization problem and what 
methods will work for it.

  For a start, you might try other optimization algorithms (see e.g.
the 'optimx' package, which may be on R-forge rather than CRAN
[I don't remember], as well as the Optimization task view on CRAN).

  good luck
Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] missing values in party::ctree

[Andrew Ziem]

After ctree builds a tree, how would I determine the direction missing values 
follow by examining the BinaryTree-class object?  For instance in the example 
below Bare.nuclei has 16 missing values and is used for the first split, but 
the missing values are not listed in either set of factors.   (I have the same 
question for missing values among numeric [non-factor] values, but I assume the 
answer is similar.)


> require(party)
> require(mlbench)
> data(BreastCancer)
> BreastCancer$Id <- NULL
> ct <- ctree(Class ~ . , data=BreastCancer, controls = ctree_control(maxdepth 
> = 1))
> ct

 Conditional inference tree with 2 terminal nodes

Response:  Class 
Inputs:  Cl.thickness, Cell.size, Cell.shape, Marg.adhesion, Epith.c.size, 
Bare.nuclei, Bl.cromatin, Normal.nucleoli, Mitoses 
Number of observations:  699 

1) Bare.nuclei == {1, 2}; criterion = 1, statistic = 488.294
  2)*  weights = 448 
1) Bare.nuclei == {3, 4, 5, 6, 7, 8, 9, 10}
  3)*  weights = 251 
> sum(is.na(BreastCancer$Bare.nuclei))
[1] 16
> nodes(ct, 1)[[1]]$psplit
Bare.nuclei == {1, 2}
> nodes(ct, 1)[[1]]$ssplit
list()



Based on below, the answer is node 2, but I don't see it in the object.

> sum(BreastCancer$Bare.nuclei %in% c(1,2,NA))
[1] 448
> sum(BreastCancer$Bare.nuclei %in% c(1,2))
[1] 432
> sum(BreastCancer$Bare.nuclei %in% c(3:10))
[1] 251


Andrew

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] summing 15 minute precip data to daily

[Joshua Wiley]

Hi Janet,

One relatively simple way would be to transofrm the data into a 96 x
Ndays matrix and use colSums().  Of course, lets say on one day, the
measurement tool had technical difficulties and missed two
observations, then you only have 94 observations for that day, you
will need a fancier solution that deals with time not number of
observations.  Below is an example.

Cheers,

Josh


## Imaginary precipitation data for 7 days
set.seed(10)
x <- rnorm(96*7, 1, .1)
## An alternate way you may have the data stored
xalt <- data.frame(precip = x)

## Assuming _no_ missing observations
colSums(matrix(x, nrow = 96))

## alternate version
colSums(matrix(xalt$precip, nrow = 96))


On Thu, Feb 17, 2011 at 11:56 AM, Janet Choate  wrote:
> Hi all,
> i'm sure there is an easy way to do this, but i'm stumped, so any help would
> be appreciated.
>
> i have a single column of data for precipitation every 15 minutes over a
> year.  i want to sum the precip to daily data.
> so the first 96 records = the first day, the second 96 records = the second
> day, and so on
> is there a way to write a for loop that would sum the data to daily, and
> write each value to a second object so i end up with a file of daily precip?
>
> thanx,
> Janet
>
>        [[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.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

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

[R] summing 15 minute precip data to daily

[Janet Choate]

Hi all,
i'm sure there is an easy way to do this, but i'm stumped, so any help would
be appreciated.

i have a single column of data for precipitation every 15 minutes over a
year.  i want to sum the precip to daily data.
so the first 96 records = the first day, the second 96 records = the second
day, and so on
is there a way to write a for loop that would sum the data to daily, and
write each value to a second object so i end up with a file of daily precip?

thanx,
Janet

[[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] A very basic line-plot question

[world peace]

Hi All

I have data like this

  tom   randy mike dan doug
height   150   152 155  134 141

I am trying to create a line plot, with names on X-axis and height measure
on Y. how can i get it through R.
I could get several versions which are close (dotchart, bargraph), but not
quite the same thing.

I am looking for something like
http://www.statsoft.com/Portals/0/blog/line_plot.jpg
with option to add more lines, such as weight measurement for above data.

Thanks,
A

[[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] Indentify polygons that are on the border of a shapefile

[Leonardo Monasterio]

Dear R users,

I would like to know how to indentify the polygons that are located on the
border of a map (i.e.shapefile).

Do you have any suggestion on how to do it?
Thank you very much,
Leo  Monasterio.

[[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] Multivariate BLUP

[Muhammad Yaseen]

*Dear All,*
*I'm trying to do Multivariate BLUP in R. I'd highly appreciate if someone
can share R code and data for Multivariate BLUP. Thanks*
*
*
*Regards!
*
-- 
*

Muhammad Yaseen
*

[[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] Predictive Analytics with R, PMML and ADAPA

[MZ]

This is a presenation from the R Users Group (Bay Area)

Covers building predictive analytic models in R, exporting in PMML and
using ADAPA for model deployment and execution.  Introduction to the
Predictive Model Markup Language (PMML) standard and how it helps to
overcome memory and speed limitations of R and makes models available
for operational deployment and integration via web services.

http://prezi.com/nyxnpa4kqbdo/predictive-modeling-with-r-pmml-and-adapa/


Covering the same topics, this is also available as a video:
http://adapasupport.zementis.com/2010/12/predictive-analytics-with-r-pmml-adapa.html

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Pre-allocation of matrices is LESS efficient?

[Duncan Murdoch]

On 17/02/2011 11:02 AM, Alex F. Bokov wrote:

Motivation: during each iteration, my code needs to collect tabular data (and 
use it only during that iteration), but the rows of data may vary. I thought I 
would speed it up by preinitializing the matrix that collects the data with 
zeros to what I know to be the maximum number of rows. I was surprised by what 
I found...

# set up (not the puzzling part)
x<-matrix(runif(20),nrow=4); y<-matrix(0,nrow=12,ncol=5); foo<-c();

# this is what surprises me... what the?
>  system.time(for(i in 1:10){n<-sample(1:4,1);y[1:n,]<-x[1:n,];});
user  system elapsed
   1.510   0.000   1.514
>  system.time(for(i in 1:10){n<-sample(1:4,1);foo<-x[1:n,];});
user  system elapsed
   1.090   0.000   1.085

These results are very repeatable. So, if I'm interpreting them correctly, 
dynamically allocating 'foo' each time to whatever the current output size is 
runs faster than writing to a subset of a preallocated 'y'? How is that 
possible?


The expression

y[1:n,]<-x[1:n,]


creates a new temporary variable to hold the result of the expression 
x[1:n,], then copies the elements of it to y[1:n,].


The expression

foo <- x[1:n,]

creates the same temporary, and then binds foo to it without doing any 
copying.  Much less work.



And, more generally, I'm sure other people have encountered this type of 
situation. Am I reinventing the wheel? Is there a best practice for storing 
temporary loop-specific data?


Storing the value of an expression in a new variable will always be 
faster than copying it into part of an existing variable.


Duncan Murdoch

P.S.  You might be aware of this, but there's one other thing that might 
be a surprise to you:  x[1:1,] will be a vector, while x[1:n,] will be
a matrix for n>1.  Use the "drop=FALSE" argument if you always want a 
matrix result.



Thanks.

PS:  By the way, though I cannot write to foo[,] because the size is different 
each time, I tried writing to foo[] and the runtime was worse than either of 
the above examples.

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

[andrija djurovic]

This is, maybe, not the best solution but I hope it will help you:

x<-data.frame(id=c(1,1,1,2,2,2,2,3,3,3,4,4), v1=c(1:12), V2=c(12:23))

do.call("rbind",by(x,x$id,function(x) x[c(sample(nrow(x),2)),]))

Andrija

On Thu, Feb 17, 2011 at 6:39 PM, yf  wrote:

>
> But i need for each id have two data.
> Like...
> > x
>   id v1 V2
> 1   1  1 12
> 2   1  2 13
>
> 4   2  4 15
> 5   2  5 16
>
>
> 8   3  8 19
> 9   3  9 20
>
> 11  4 11 22
> 12  4 12 23
>
> So should write sample( if sample id >2  ,2). I don't know how to write
>  (if
> sample id >2). Thanks.
> --
> View this message in context:
> http://r.789695.n4.nabble.com/sampling-tp3310184p3311253.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Pre-allocation of matrices is LESS efficient?

[Douglas Bates]

On Thu, Feb 17, 2011 at 10:02 AM, Alex F. Bokov
 wrote:
> Motivation: during each iteration, my code needs to collect tabular data (and 
> use it only during that iteration), but the rows of data may vary. I thought 
> I would speed it up by preinitializing the matrix that collects the data with 
> zeros to what I know to be the maximum number of rows. I was surprised by 
> what I found...
>
> # set up (not the puzzling part)
> x<-matrix(runif(20),nrow=4); y<-matrix(0,nrow=12,ncol=5); foo<-c();

There is no purpose in initializing foo here.  Your assignment in the
second version overwrites any assignment here.

> # this is what surprises me... what the?
>> system.time(for(i in 1:10){n<-sample(1:4,1);y[1:n,]<-x[1:n,];});
>   user  system elapsed
>  1.510   0.000   1.514

This version performs extraction from x and assignment into a
submatrix of y.  The second version performs only the extraction and
assignment to a name in the evaluation environment, which is a much
faster operation.

>> system.time(for(i in 1:10){n<-sample(1:4,1);foo<-x[1:n,];});
>   user  system elapsed
>  1.090   0.000   1.085
>
> These results are very repeatable. So, if I'm interpreting them correctly, 
> dynamically allocating 'foo' each time to whatever the current output size is 
> runs faster than writing to a subset of a preallocated 'y'? How is that 
> possible?
>
> And, more generally, I'm sure other people have encountered this type of 
> situation. Am I reinventing the wheel? Is there a best practice for storing 
> temporary loop-specific data?
>
> Thanks.
>
> PS:  By the way, though I cannot write to foo[,] because the size is 
> different each time, I tried writing to foo[] and the runtime was worse than 
> either of the above examples.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] censoring symbols on survfit plot

[threshold]

Hi, when ploting Kaplan-Meier estimate curves as below, the censoring symbols
(crosses) to not change thickness along the lines 
plot(survfit(surv ~ I(x>=cut.off) ),lty=c(1,2), lwd=2)

is there any strightforward way to make it happen? thanks
robert

-- 
View this message in context: 
http://r.789695.n4.nabble.com/censoring-symbols-on-survfit-plot-tp3311283p3311283.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] Urgent Request

[muhammad mohsin]

Dear Colleagues,
Hope you will be fine. I am student of Ph.D and doing some work on 
distribution. 
I developed a new distribution and having some problems in estimating their 
parameters by MLE. I used R-program and  used "maxLik" function (maxLik: A 
Package for Maximum Likelihood Estimation in R) But there is some problem, it 
is 
not estimated the parameters properly. I also write an e-mail to the author of 
this paper but he could not solve my problem. His function works well for 
simple 
and known distribution but does not work for a new function.  Can anybody spare 
some time for me?  I really need your help. Please inform me so that I can send 
you the material.
Waiting for a quick reply 
Best Regards,

 Muhammad Mohsin
PhD Research Fellow
University of Klagenfurt, Department of Statistics
University st. 65-67, 9020
Klagenfurt, Austria, Europe
University E-mail Address:mmoh...@edu.uni-klu.ac.at
Phone No Office 0043 (0)463 27003129
Mobile No: 0043 (0)676 7218836 


  
[[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] incomplete final line

[Joshua Wiley]

On Thu, Feb 17, 2011 at 12:35 AM, mipplor  wrote:
>
> thx for your suggestions , i have made it csv file,and it looks:
>
>
>> data<‐read.table('E:/my documents/r/1.csv',header=TRUE)

It looks like you could benefit from (re)reading my previous email.
You either need to specify the sep = argument in read.table() or use
read.csv(), which is a wrapper for read.table() with defaults for CSV
files.

Josh

> Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings,
> :
>  line 1 did not have 2 elements
>
>
> is there wrong with the data?
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/incomplete-final-line-tp3309405p3310357.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] exact logistic regression

[Łukasz Ręcławowicz]

I believe that this code will work (...for very small) samples, but let some
correct me if there is something wrong.

require(logistf);require(combinat)
permY<-permn(data$y)
ntimes<-length(permY)
results<-matrix(nrows=ntimes,ncols=number_of_coefficients)
for(i in 1:ntimes){
results[i,]<-logistf(unlist(permY[i])~factor_A+factor_B,data=data)}
observed<-logistf(y~factor_A+factor_B,data=data)$coefficients
# Exact p-values will be:
sum(results[,1]>=observed[1])/ntimes # One for intercept
sum(results[,2]>=observed[2])/ntimes
sum(results[,3]>=observed[3])/ntimes
-- 
Mi³ego dnia

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

[yf]

But i need for each id have two data. 
Like...
> x
   id v1 V2
1   1  1 12
2   1  2 13

4   2  4 15
5   2  5 16


8   3  8 19
9   3  9 20

11  4 11 22
12  4 12 23

So should write sample( if sample id >2  ,2). I don't know how to write  (if
sample id >2). Thanks.
-- 
View this message in context: 
http://r.789695.n4.nabble.com/sampling-tp3310184p3311253.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] How to speed up a for() loop

[Phil Spector]

Simona -
   I don't think preallocating your random variables would make
the code run any faster.
   A very simple change that would speed things up a little would
be to replace

 simPD.vec[i]=length(R.vec[R.vec

Dear all,

Does anyone have any idea on how to speed up the for() loop below.
Currently it takes approximately 2 minutes and 30 seconds.

Because of the size of Nsim and N, simulating a multivariate normal
(instead of simulating Nsim times a vector of N normal distributions)
would require too much memory space.

Many thanks for your kind help,

Simona


N=3000
PD=runif(N,0,1)
cutoff.=qnorm(PD)

rho.=0.1

Nsim=10
simPD.vec=0*(1:Nsim)

systemic = rnorm(Nsim,0,1);

for (i in 1:Nsim)
{
R.vec=sqrt(rho.)*systemic[i]+sqrt(1-rho.)*rnorm(N,0,1)
simPD.vec[i]=length(R.vec[R.vec

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 and replace all the elements in a data frame

[Gong-Yi Liao]

You may write as this:

for (i in 1:nrow(x)){
   for (j in 1:ncol(x)){
   if (!is.na(x[i, j])) {
  if(x[i, j] == 'A') {x2[i, j] <- 'A/A'} else{
  if(x[i, j] == 'T') {x2[i, j] <- 'T/T'} else{
   if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{
  if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{x2[i, j]
<- NA}
 }
  }
  }
   }
   }
}


On Thu, Feb 17, 2011 at 11:54 AM, Josh B  wrote:

> Hi all,
>
> I'm having a problem once again, trying to do something very simple.
> Consider
> the following data frame:
>
> x <- read.table(textConnection("locus1 locus2 locus3
> A T C
> A T NA
> T C C
> A T G"), header = TRUE)
> closeAllConnections()
>
> I am trying to make a new data frame, replacing "A" with "A/A", "T" with
> "T/T",
> "G" with "G/G", and "C" with "C/C." Note also the presence of an "NA"
> (missing
> data) in the data frame, which should be carried over to the new data
> frame.
>
> Here is what I am trying, which fails miserably:
>
> x2 <- data.frame(matrix(nrow = nrow(x), ncol = ncol(x)))
>
> for (i in 1:nrow(x)){
>for (j in 1:ncol(x)){
>if(x[i, j] == 'A') {x2[i, j] <- 'A/A'} else{
>if(x[i, j] == 'T') {x2[i, j] <- 'T/T'} else{
> if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{
>if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{x2[i, j] <-
> NA}
>}
>   }
>   }
>}
> }
>
> I get the following error message:
> Error in if (x[i, j] == "A") { : missing value where TRUE/FALSE needed
>
> So what am I doing wrong? If you can provide me with specific code that
> fixes
> the problem and gets the job done, that would be the most useful.
>
>
> Thanks very much in advance for your help!
>
> Sincerely,
> ---
> Josh Banta, Ph.D
> Center for Genomics and Systems Biology
> New York University
> 100 Washington Square East
> New York, NY 10003
> Tel: (212) 998-8465
> http://plantevolutionaryecology.org
>
>
>
>[[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.
>



-- 
Gong-Yi Liao

Department of Statistics
University of Connecticut
215 Glenbrook Road  U4120
Storrs, CT 06269-4120

860-486-9478

[[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] removing lower and upper quantiles from an arry

[Greg Snow]

In addition to the other answers that you received you can also do:

library(TeachingDemos)
i[ quantile(i,.25) %<% i %<% quantile(i,.75) ]

This may or may not be more readable than the others.  Also note that 
precomputing both quantiles in one step may be faster than calling quantile 
twice.

You could also do a partial sort of your data and just pull out the middle 
section (though you would probably lose the ordering in the data).

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


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
> project.org] On Behalf Of Maas James Dr (MED)
> Sent: Thursday, February 17, 2011 3:09 AM
> To: r-help@r-project.org
> Subject: [R] removing lower and upper quantiles from an arry
> 
> I'm trying to work out the simplest way to remove the upper and lower
> quantiles, in this case upper and lower 25% from an array.  I can do it
> in two steps but when I try it in one, it fails.  Is there something
> simple missing from my syntax or are there other simple elegant way to
> accomplish this?
> 
> Thanks
> 
> J
> 
> > i <-1:20
> > i2 <- i[i > i3 <- i[i>quantile(i,.25)]
> > i4 <- i[quantile(i,.25)< i > quantile(i,.75)]
> Error: unexpected '>' in "i4 <- i[quantile(i,.25)< i >"
> 
> ===
> Dr. Jim Maas
> University of East Anglia
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] color.scale error

[Alaios]

Dear all
when I call color.scale like this:


require('plotrix')
colcolor<-color.scale(c(range_sr,sr),extremes=c("red","blue"))


Error in if (min(reds) < 0 || max(reds) > 1) reds <- rescale(reds, c(0,  : 
  missing value where TRUE/FALSE needed


range_sr 
[1] -10.00  44.02977813958366

str(sr)
 num [1:100, 1:100] 2.54 2.71 2.89 2.95 3.03 ...


this was working on simpler data sets before. 

How can I understand more what was the problem in depth.? I have read the 
source code of color.scale but unfortunately I did not make any progress to 
understand this error message. Is it possible to raise an exception before the 
error happens and see what went wrong?

Regards
Alex

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 and replace all the elements in a data frame

[baptiste auguie]

Hi,

You could use car::recode to change the levels of the factors,

library(car)
transform(x, locus1 = recode(locus1, "'A' = 'A/A' ; else = 'T/T'"),
locus2 = recode(locus2, "'T'='T/T' ; 'C' = 'C/C'"),
locus3 = recode(locus3, "'C'='C/C' ; 'G' = 'G/G'"))


HTH,

baptiste


On 17 February 2011 17:54, Josh B  wrote:
> Hi all,
>
> I'm having a problem once again, trying to do something very simple. Consider
> the following data frame:
>
> x <- read.table(textConnection("locus1 locus2 locus3
> A T C
> A T NA
> T C C
> A T G"), header = TRUE)
> closeAllConnections()
>
> I am trying to make a new data frame, replacing "A" with "A/A", "T" with 
> "T/T",
> "G" with "G/G", and "C" with "C/C." Note also the presence of an "NA" (missing
> data) in the data frame, which should be carried over to the new data frame.
>
> Here is what I am trying, which fails miserably:
>
> x2 <- data.frame(matrix(nrow = nrow(x), ncol = ncol(x)))
>
> for (i in 1:nrow(x)){
>    for (j in 1:ncol(x)){
>        if(x[i, j] == 'A') {x2[i, j] <- 'A/A'} else{
>            if(x[i, j] == 'T') {x2[i, j] <- 'T/T'} else{
>                 if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{
>                    if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{x2[i, j] <- NA}
>                }
>           }
>       }
>    }
> }
>
> I get the following error message:
> Error in if (x[i, j] == "A") { : missing value where TRUE/FALSE needed
>
> So what am I doing wrong? If you can provide me with specific code that fixes
> the problem and gets the job done, that would be the most useful.
>
>
> Thanks very much in advance for your help!
>
> Sincerely,
> ---
> Josh Banta, Ph.D
> Center for Genomics and Systems Biology
> New York University
> 100 Washington Square East
> New York, NY 10003
> Tel: (212) 998-8465
> http://plantevolutionaryecology.org
>
>
>
>        [[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] Find and replace all the elements in a data frame

[Henrique Dallazuanna]

Try this:

xNew <- as.data.frame(mapply(paste, x, x, sep = "/"))
xNew[is.na(x)] <- NA
xNew

On Thu, Feb 17, 2011 at 2:54 PM, Josh B  wrote:

> Hi all,
>
> I'm having a problem once again, trying to do something very simple.
> Consider
> the following data frame:
>
> x <- read.table(textConnection("locus1 locus2 locus3
> A T C
> A T NA
> T C C
> A T G"), header = TRUE)
> closeAllConnections()
>
> I am trying to make a new data frame, replacing "A" with "A/A", "T" with
> "T/T",
> "G" with "G/G", and "C" with "C/C." Note also the presence of an "NA"
> (missing
> data) in the data frame, which should be carried over to the new data
> frame.
>
> Here is what I am trying, which fails miserably:
>
> x2 <- data.frame(matrix(nrow = nrow(x), ncol = ncol(x)))
>
> for (i in 1:nrow(x)){
>for (j in 1:ncol(x)){
>if(x[i, j] == 'A') {x2[i, j] <- 'A/A'} else{
>if(x[i, j] == 'T') {x2[i, j] <- 'T/T'} else{
> if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{
>if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{x2[i, j] <-
> NA}
>}
>   }
>   }
>}
> }
>
> I get the following error message:
> Error in if (x[i, j] == "A") { : missing value where TRUE/FALSE needed
>
> So what am I doing wrong? If you can provide me with specific code that
> fixes
> the problem and gets the job done, that would be the most useful.
>
>
> Thanks very much in advance for your help!
>
> Sincerely,
> ---
> Josh Banta, Ph.D
> Center for Genomics and Systems Biology
> New York University
> 100 Washington Square East
> New York, NY 10003
> Tel: (212) 998-8465
> http://plantevolutionaryecology.org
>
>
>
>[[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



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

[[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 and replace all the elements in a data frame

[Sarah Goslee]

Josh, you've made it far too complicated. Here's one simpler way (note
that I changed your read.table statement to make the values NOT factors,
since I wouldn't think you want that).

> x <- read.table(textConnection("locus1 locus2 locus3
+ A T C
+ A T NA
+ T C C
+ A T G"), header = TRUE, as.is=TRUE)
> closeAllConnections()
>
> x2 <- x
> x2[x2 == "A"] <- "A/A"
> x2[x2 == "T"] <- "T/T"
> x2[x2 == "G"] <- "G/G"
> x2[x2 == "C"] <- "C/C"
> x2
  locus1 locus2 locus3
1A/AT/TC/C
2A/AT/T   
3T/TC/CC/C
4A/AT/TG/G

If you do for some reason want a factor, you'll need to adjust the
levels for each
column before doing this.

Sarah

On Thu, Feb 17, 2011 at 11:54 AM, Josh B  wrote:
> Hi all,
>
> I'm having a problem once again, trying to do something very simple. Consider
> the following data frame:
>
> x <- read.table(textConnection("locus1 locus2 locus3
> A T C
> A T NA
> T C C
> A T G"), header = TRUE)
> closeAllConnections()
>
> I am trying to make a new data frame, replacing "A" with "A/A", "T" with 
> "T/T",
> "G" with "G/G", and "C" with "C/C." Note also the presence of an "NA" (missing
> data) in the data frame, which should be carried over to the new data frame.
>
> Here is what I am trying, which fails miserably:
>
> x2 <- data.frame(matrix(nrow = nrow(x), ncol = ncol(x)))
>
> for (i in 1:nrow(x)){
>    for (j in 1:ncol(x)){
>        if(x[i, j] == 'A') {x2[i, j] <- 'A/A'} else{
>            if(x[i, j] == 'T') {x2[i, j] <- 'T/T'} else{
>                 if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{
>                    if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{x2[i, j] <- NA}
>                }
>           }
>       }
>    }
> }
>
> I get the following error message:
> Error in if (x[i, j] == "A") { : missing value where TRUE/FALSE needed
>
> So what am I doing wrong? If you can provide me with specific code that fixes
> the problem and gets the job done, that would be the most useful.
>
>
> Thanks very much in advance for your help!
>
> Sincerely,
> ---
> Josh Banta, Ph.D
> Center for Genomics and Systems Biology
> New York University
> 100 Washington Square East
> New York, NY 10003
> Tel: (212) 998-8465
> http://plantevolutionaryecology.org
>
>
>
>        [[alternative HTML version deleted]]




-- 
Sarah Goslee
http://www.functionaldiversity.org

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 and replace all the elements in a data frame

[Josh B]

Hi all,

I'm having a problem once again, trying to do something very simple. Consider 
the following data frame:

x <- read.table(textConnection("locus1 locus2 locus3
A T C
A T NA
T C C
A T G"), header = TRUE)
closeAllConnections()

I am trying to make a new data frame, replacing "A" with "A/A", "T" with "T/T", 
"G" with "G/G", and "C" with "C/C." Note also the presence of an "NA" (missing 
data) in the data frame, which should be carried over to the new data frame.

Here is what I am trying, which fails miserably:

x2 <- data.frame(matrix(nrow = nrow(x), ncol = ncol(x)))

for (i in 1:nrow(x)){
for (j in 1:ncol(x)){
if(x[i, j] == 'A') {x2[i, j] <- 'A/A'} else{
if(x[i, j] == 'T') {x2[i, j] <- 'T/T'} else{
 if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{
if(x[i, j] == 'G') {x2[i, j] <- 'G/G'} else{x2[i, j] <- NA}
}
   }
   }
}
}

I get the following error message:
Error in if (x[i, j] == "A") { : missing value where TRUE/FALSE needed

So what am I doing wrong? If you can provide me with specific code that fixes 
the problem and gets the job done, that would be the most useful. 


Thanks very much in advance for your help!

Sincerely,
---
Josh Banta, Ph.D
Center for Genomics and Systems Biology
New York University
100 Washington Square East
New York, NY 10003
Tel: (212) 998-8465
http://plantevolutionaryecology.org


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

[rex.dwyer]

I hate to sound like David "Have You Read The Posting Guide?" Winsemius, but 
there's no way for anyone to know what you are trying to accomplish here 
without a lot more information.  You don't show us the output you expect and 
the output you got.  I would expected "relatedness" to be on a scale from 0 to 
1, but it's clear that you'll get values >1 in this program.
To use R effectively, you need to rephrase your computation as a matrix 
computation.  People generally use R at least partly to avoid debugging index 
computations in for-loops.  For-loops are also much slower than the 
corresponding matrix operations in R.  If you want to use for-loops, you can 
always put in some prints and trace what's going on, just like the old days!

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of Val
Sent: Wednesday, February 16, 2011 3:14 PM
To: r-h...@stat.math.ethz.ch
Subject: [R] covar

Hi all,

I want to construct relatedness among individuals and have a look at the
following script.

#
rm(list=ls())

N=5
id   = c(1:N)
dad = c(0,0,0,3,3)
mom  = c(0,0,2,1,1)
sex  = c(2,2,1,2,2) # 1= M and 2=F

   A=diag(nrow = N)
   for(i in 1:N){
  for(j in i:N)  {
 ss = dad[j]
 dd = mom[j]
 sx = sex[j]
  if( ss > 0 && dd > 0 )
{
  if(i == j)
   { A[i,i] = 1 + 0.5*A[ss,dd] }
 else
  { A[i,j] = A[i,ss] + 0.5*(A[i,dd])
A[j,i] = A[i,j] }
}

  } #inner for loop
 } # outer for loop
  A

If the sex is male ( sex=1)  then I want to set A[i,i]=0.5*A[ss,dd]
If it is female ( sex=2) then A[i,i] = 1 + 0.5*A[ss,dd]


How do I do it ?

I tried several cases but it did not work from me. Your assistance is
highly  appreciated  in advance

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.




message may contain confidential information. If you are not the designated 
recipient, please notify the sender immediately, and delete the original and 
any copies. Any use of the message by you is prohibited. 
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Transforming relational data

[mathijsdevaan]

Thanks for helping me out so generously. After reading the vignettes and the
other info I still have a question (sorry I am a R novice):

I am not so much trying to construct time series (although it comes very
close). Rather for each pair (Bi,Bj) in project (An) I am trying to sum up
the values of v for (Bi,Bj) where CHello. One (of many) solution might be:

>require(data.table)
>DT = data.table(read.table(textConnection("A  B  C
>1 1  a  1999
>2 1  b  1999
>3 1  c  1999
>4 1  d  1999
>5 2  c  2001
>6 2  d  2001"),head=TRUE,stringsAsFactors=FALSE))

>firststep = DT[,cbind(expand.grid(B,B),v=1/length(B)),by=C][Var1!=Var2]
>setkey(firststep,Var1,Var2)
>grp3 = c("a","b","d")
>firststep[J(expand.grid(grp3,grp3)),nomatch=0][,sum(v)]
># 2.5

>If I guess the bigger picture correctly, this can be extended
>to make a time series of prior familiarity by including
>the year in the key.

>If you decide to try this, please make sure to grab the latest
>(recent) version of data.table from CRAN (v1.5.3). Suggest that
>you run it first to confirm it does return 2.5, then break it
>down and run it step by step to see how each part works. You
>will need some time to read the vignettes and ?data.table
>(which has recently been improved) but I hope you think it is
>worth it. Support is available at maintainer("data.table").

>HTH
>Matthew


>>On Mon, 14 Feb 2011 09:22:12 -0800, mathijsdevaan wrote:
>> Hi,
>>
>> I have a large dataset with info on individuals (B) that have been
>> involved in projects (A) during multiple years (C). The dataset contains
>> three columns: A, B, C. Example:
>>
>>A  B  C
>> 1 1  a  1999
>> 2 1  b  1999
>>3 1  c  1999
>> 4 1  d  1999
>>5 2  c  2001
>>6 2  d  2001
>> 7 3  a  2004
>> 8 3  c  2004
>> 9 3  d  2004
>>
>> I am interested in how well all the individuals in a project know each
>> other. To calculate this team familiarity measure I want to sum the
>> familiarity between all individual pairs in a team. The familiarity
>> between each individual pair in a team is calculated as the summation of
>> each pair's prior co-appearance in a project divided by the total number
>> of team members. So the team familiarity in project 3 = (1/4+1/4) +
>> (1/4+1/4+1/2) + (1/4+1/4+1/2) = 2,5 or a has been in project 1 (of size
>> 4) with c and d > 1/4+1/4 and c has been in project 1 (of size 4) with 1
>> and d > 1/4+1/4 and c has been in project 2 (of size 2) with d > 1/2.
>>
>> I think that the best way to do it is to transform the data into an
>> edgelist (each pair in one row/two columns) and then creating two
>> additional columns for the strength of the familiarity and the year of
>> the project in which the pair was active. The problem is that I am stuck
>> already in the first step. So the question is: how do I go from the
>> current data structure to a list of projects and the familiarity of its
>> team members?
>>
>> Your help is very much appreciated. Thanks!

-- 
View this message in context: 
http://r.789695.n4.nabble.com/Re-Transforming-relational-data-tp3307449p3311101.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] barplot with errorbars

[Toby Marthews]

If you google "barplot with error bars" you immediately find 
http://svitsrv25.epfl.ch/R-doc/library/prada/html/barploterrbar.html .
Toby.


From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] On Behalf Of 
Lathouri, Maria [m.lathour...@imperial.ac.uk]
Sent: 17 February 2011 16:00
To: r-help@r-project.org
Subject: [R] barplot with errorbars

Dear all

I have six variables of the average metal concentrations

Var1 4.77
Var2 23.5
Var3 5.2
Var4 12.3
Var5 42.1
Var6 121.2

I want to plot them as a barplot with error bars. Could you help me?

Cheers
Maria

[[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] Dependency on R-Forge Package

[Prof Brian Ripley]

On Thu, 17 Feb 2011, Uwe Ligges wrote:




On 17.02.2011 15:29, Damian Betebenner wrote:
In building a package, is it possible to make the package depend upon 
another package that is only available on R-Forge (not CRAN). For example, 
by doing something in the DESCRIPTION file
I'd like to add a dependency to my package such that when the user install 
my package it will automatically install this other package from R-forge as 
well.



Not completely, the user has to choose R-forge as a repository (which is not 
the default but easily selectable).


And of course on R-forge being up (and its admins tell me they never 
intended it to be used to distribute packages other than to testers).


However, you can if you want to get your package to install this other 
r-forge package: a package which does that is HSAUR2.  That's 
reasonably sure to work if the r-forge package does not contain 
compiled code (just using install.packages(type="source")).  Otherwise 
you depend on R-forge having built the binary for the end user's 
platform.


For the record, some of us are pretty unhappy about CRAN packages 
having off-CRAN dependencies (and have talked about adding metadata 
to say where such packages should be located) and would be very 
unhappy about an R-forge one.




Uwe Ligges






Any help greatly appreciated.

Damian


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



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

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

Re: [R] VAR with HAC

[Marta Lachowska]

Thank you for your hint!

I see that there was a thread discussing implementation of what I
wanted to do (Newey-West standard errors in a VAR context), but that
there is a conflict due to how the type = "const" is defined in the VAR
command: https://stat.ethz.ch/pipermail/r-sig-finance/2009q2/004272.html
that appears not to be resolved. 
 
Best, 
 
Marta
 

>>> "Pfaff, Bernhard Dr."  2/17/2011
4:31 AM >>>
Hello Marta,

have you read ?coeftest and ? VAR carefully enough? The function does
expect a lm/glm object for x as argument. Hence, the following does
work:

library(vars)
data(Canada)
myvar <- VAR(Canada, p = 2, type = "const")
lapply(myvar$varresult, coeftest)

Best,
Bernhard 

> -Ursprüngliche Nachricht-
> Von: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-project.org] Im Auftrag von Marta Lachowska
> Gesendet: Mittwoch, 16. Februar 2011 16:50
> An: r-help@r-project.org
> Betreff: [R] VAR with HAC
> 
> 
> Hello,
> I would like to estimate a VAR model with HAC corrected 
> standard errors. I tried to do this by using the sandwich 
> package, for example: 
>  
> > library(vars)
> > data(Canada)
> > myvar = VAR(Canada, p = 2, type = "const") coeftest(myvar, vcov = 
> > vcovHAC)
> Error in umat - res : non-conformable arrays
>  
> Which suggests that this function is not compatible with the 
> VAR command. Has anyone tried to modify the code to get HAC 
> corrected standard errors with VAR? Any suggestions are welcome. 
>  
> Thank you. 
>  
> Marta
> 
> [[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.
> 
*
Confidentiality Note: The information contained in this ...{{dropped:14}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Multi-response MCMCglmm (gaussian and zapoisson)

[Susanne Lachmuth]

Dear MCMCglmm users,

I am currently struggling with the specification of a proper prior and model 
formula for a multi-response MCMCglmm with two of the three response variables 
being Gaussian and the third being za-poisson. The model includes several fixed 
effects and three nested random effects.
In general, I would prefer to fit a model with a fixed effect of trait and 
suppressed intercept for getting trait specific intercepts. Further, I wish to 
measure covariances between the response variables and consequently to specify 
completely parameterized covariance matrices. 

My current model looks like this:

prior<-list(R=list(V=diag(4),nu=0.004),G=list(G1=list(V=diag(4)/4,n=4),G2=list(V=diag(4)/4,n=4),G3=list(V=diag(4)/4,n=4)))

m1<-MCMCglmm(fixed=cbind(resp1,resp2,resp3) ~
trait:(expl1 + expl2+ expl3+ expl4+ expl5...)+trait-1,
random=~us(trait):rand1+us(trait): rand1: rand2+us(trait): rand1: rand2: rand3,
rcov=~us(trait):units,family=c("gaussian","gaussian","zapoisson"),nitt=2,burnin=5000,thin=10,prior=prior,data=dat)

However, for zero-altered models it is recommended in the MCMCglmm course notes 
to constrain over-dispersion to be the same for both processes (the 
zero-alteration and poisson process) by using a trait by unit interaction in 
the R-structure. Additionally, the intercept should not be suppressed for 
getting the differences between the zero-altered regression coefficients and 
the Poisson regression coefficients. This allows identification of 
zero-inflation or zero-deflation in response to the explanatory variables.
I therefore fit an additional model for the zapoisson response only, looking 
like this:

prior2<-list(R=list(V=diag(1)/4,nu=0.002),G=list(G1=list(V=diag(1)/4,nu=0.002),G2=list(V=diag(1)/4,nu=0.002),G3=list(V=diag(1)/4,nu=0.002)))

m2<-MCMCglmm(fixed=resp3 ~
trait:(expl1 + expl2+ expl3+ expl4+ expl5….)+trait,
 random=~trait:rand1+trait: rand1: rand2+ trait: rand1: rand2: rand3,
rcov=~trait:units, family="zapoisson",nitt=2,burnin=1000,thin=10,data=dat, 
prior=prior2)

The results show that there is “significant” zero-inflation and deflation in 
response to some variables.

My main questions are:
Are the two priors specified correctly?
Does it make sense to include the zapoisson response (resp3) in model m1 and is 
the model formula (in particular the R-structure) appropriate?
An alternative might be to analyze resp1 and resp2 (both size variables of 
plants) with model m1 and fit an extra model (like m2) for resp 3 (Number of 
flowers) with the size parameters (i.e. responses of m1) as covariates. We are 
interested in a trade-off between growth and reproduction.

Any help would be greatly appreciated.

Many thanks,

Susanne Lachmuth




Susanne Lachmuth
Department of Plant Ecolgy
Martin Luther University Halle-Wittenberg 
Germany

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

[R] How to speed up a for() loop

[simona.costanzo]

Dear all,

Does anyone have any idea on how to speed up the for() loop below.
Currently it takes approximately 2 minutes and 30 seconds.

Because of the size of Nsim and N, simulating a multivariate normal
(instead of simulating Nsim times a vector of N normal distributions)
would require too much memory space.

Many thanks for your kind help,

Simona


N=3000
PD=runif(N,0,1)
cutoff.=qnorm(PD)

rho.=0.1

Nsim=10
simPD.vec=0*(1:Nsim)

systemic = rnorm(Nsim,0,1);

for (i in 1:Nsim)
{
R.vec=sqrt(rho.)*systemic[i]+sqrt(1-rho.)*rnorm(N,0,1)
simPD.vec[i]=length(R.vec[R.vecVisit our website at http://www.ubs.com 

This message contains confidential information and is intended only 
for the individual named. If you are not the named addressee you 
should not disseminate, distribute or copy this e-mail. Please 
notify the sender immediately by e-mail if you have received this 
e-mail by mistake and delete this e-mail from your system. 

E-mails are not encrypted and cannot be guaranteed to be secure or 
error-free as information could be intercepted, corrupted, lost, 
destroyed, arrive late or incomplete, or contain viruses. The sender 
therefore does not accept liability for any errors or omissions in the 
contents of this message which arise as a result of e-mail transmission. 
If verification is required please request a hard-copy version. This 
message is provided for informational purposes and should not be 
construed as a solicitation or offer to buy or sell any securities 
or related financial instruments. 

UBS Limited is a company limited by shares incorporated in the United 
Kingdom registered in England and Wales with number 2035362. 
Registered office: 1 Finsbury Avenue, London EC2M 2PP.  UBS Limited 
is authorised and regulated by the Financial Services Authority. 

UBS AG is a public company incorporated with limited liability in 
Switzerland domiciled in the Canton of Basel-City and the Canton of 
Zurich respectively registered at the Commercial Registry offices in 
those Cantons with Identification No: CH-270.3.004.646-4 and having 
respective head offices at Aeschenvorstadt 1, 4051 Basel and 
Bahnhofstrasse 45, 8001 Zurich, Switzerland.  Registered in the 
United Kingdom as a foreign company with No: FC021146 and having a 
UK Establishment registered at Companies House, Cardiff, with No:  
BR 004507.  The principal office of UK Establishment: 1 Finsbury Avenue, 
London EC2M 2PP.  In the United Kingdom, UBS AG is authorised and 
regulated by the Financial Services Authority.

UBS reserves the right to retain all messages. Messages are protected 
and accessed only in legally justified cases. __
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] VGAM in R

[Weichao Bao]

I am using VGAM package to model a multinormial outcome with spline
function. But from the following result, I only can get p-value for variable
with spline function, such as s(As). However, there is no p-value available
for variable with linear function, such as SEX, Parity. Does anybody know
why?


> summary(fit1)
Call:
vgam(formula = mr~ s(As) + s(Pb) + SEX + parity,
family = "multinomial", data = paper)

Number of linear predictors:3

Names of linear predictors: log(mu[,1]/mu[,4]), log(mu[,2]/mu[,4]),
log(mu[,3]/mu[,4])

Dispersion Parameter for multinomial family:   1

Residual Deviance:  3292.768 on 6423.567 degrees of freedom

Log-likelihood: -1646.384 on 6423.567 degrees of freedom

Number of Iterations:  10

DF for Terms and Approximate Chi-squares for Nonparametric Effects

  Df Npar Df Npar Chisq  P(Chi)
(Intercept):1  1
(Intercept):2  1
(Intercept):3  1
s(As):11 2.2 1.5061 0.51558
s(As):21 2.8 3.7633 0.26187
s(As):31 3.8 1.8600 0.74002
s(Pb):11 2.4 7.7796 0.03172
s(Pb):21 2.8 8.5664 0.03091
s(Pb):31 3.3 5.6044 0.16036
SEX:1 1
SEX:2 1
SEX:3 1
parity:1   1
parity:2   1
parity:3   1


Thanks a lot.
Weichao

[[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] ggplot2, 'se' variable in geom_errorbar's limits?

[Eric Fail]

Dear R-list

I'm working with with geom_errorbar; specifically I'm trying to
reproduce the example Hadley Wickham have on
http://had.co.nz/ggplot2/geom_errorbar.html (all in the button of the
page) where he makes an nice plot with errorbars and then draw lines
between the points.

What confuses me is the 'limits' he defines for the errorbars from the
se variable.

First he creates a dataset,

df <- data.frame(
  trt = factor(c(1, 1, 2, 2)),
  resp = c(1, 5, 3, 4),
  group = factor(c(1, 2, 1, 2)),
  se = c(0.1, 0.3, 0.3, 0.2)
)

# library(ggplot2)

and then he creates some limits from the se variables.

limits <- aes(ymax = resp + se, ymin=resp - se)

[elements omitted]

# and then he creates the plot (I'm interested in).

p <- ggplot(df, aes(colour=group, y=resp, x=trt))
p + geom_line(aes(group=group)) + geom_errorbar(limits, width=0.2)

I can (of course) get Hadley's example to run, but I can't do it on my
data as I don't have a 'se' variable/don't know how to create it. I
have a group variable, a treatment variable, and a response variable,
but no se variable.

Could anyone out there explain how I create a 'se' variable in my data?

I'm sure my reasoning is the one that is off, and not ggplot2 (I'm a big fan).

Your help is appreciated!

Thanks,
Eric

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

Re: [R] calling pairs of variables into a function

[rex.dwyer]

Try putting d,e,f in a list:
Xxx = list(d,e,f)

For (I in 1:length(xxx))
For (j in 1:length(xxx))
If (i!=j) bigfunction(xxx[[i]],xxx[[j]])

(bad indentation, caps thanks to outlook)


-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of squamous
Sent: Wednesday, February 16, 2011 7:11 PM
To: r-help@r-project.org
Subject: [R] calling pairs of variables into a function


Hi,
I have imported three text files into R using read.table.  Their variables
are called d, e and f.

I want to run a function on all the possible combinations of these three
files.  The only way I know how to do that is like this:

bigfunction(d,e)
bigfunction(d,f)
bigfunction(e,d)
bigfunction(e,f)
bigfunction(f,e)
bigfunction(f,d)

Is there an easier way?  I will have five files later on, so it would be
useful to know!  I'd imagine I can use a loop somehow, and I have installed
a package (gregmisc) so that typing permutations(3,2) gives all the possible
pairs of three numbers, but I don't know how to combine these things to make
it work.
--
View this message in context: 
http://r.789695.n4.nabble.com/calling-pairs-of-variables-into-a-function-tp3309993p3309993.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.




message may contain confidential information. If you are not the designated 
recipient, please notify the sender immediately, and delete the original and 
any copies. Any use of the message by you is prohibited. 
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] [BioC] Make.cdf.package error

[elodiem]

Hi everybody, 

I tried to analyze a custom Affymetrix 3'-biased Array. So I wanted to make
a cdf package. (My CDF file size is 1.12Go). 

I tried several methods but the same error occured 

Method 1

> #Set the working directory 
> setwd("D:/Analyse R/Cel files") 
> #library to create cdf env 
> library("makecdfenv") 
>#Create cdf environment 
>pkgpath <- tempdir() 
>make.cdf.package("file.cdf", 
   cdf.path=getwd(), 
   compress=FALSE, species = "Homo_sapiens", 
   package.path = pkgpath, 
   verbose=TRUE) 
dir(pkgpath) 

Method 2
#Set the working directory 
setwd("D:/Analyse R/Cel files") 
#library to create cdf env 
library("makecdfenv") 
#Create cdf environment 
env <- make.cdf.env("file.cdf") 


>Reading CDF file. 
>Error in .Call("reaD file", as.character(file), as.integer(3),
as.integer(compress), :promise already under evaluation: recursive default
argument reference or earlier problems? 

The two methods take about one hour and use all the RAM of my computer
(8Go). 

Is it a limit of memory? 
Is it a non appropriated method for my array? 

Thanks in advance for your response, 
Elodie 
-- 
View this message in context: 
http://r.789695.n4.nabble.com/BioC-Make-cdf-package-error-tp3310927p3310927.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] deriving waveform values

[whizevans]

Im am curerntly trying to use R software to extract uncompressed waveform
values with I have acquired from IDLreadGLAS tool.
I am very lost with this & know I have to assess the GLA01 product data,
however, any advice on the matter atall would be brilliant, particularly on
define thresholds and connecting GLA01 and GLA14 data (or any data for that
matter)
-- 
View this message in context: 
http://r.789695.n4.nabble.com/deriving-waveform-values-tp3310999p3310999.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] Categorical Variables and Machine Learning

[Andrew Ziem]

Try the function ctree() in the package party or earth() in earth.   You can 
use factor variable as is, or you can transform the factor to binary variables 
(i.e., is_P is 0 or 1, is_D is 0 or 1).  In the second case, you can use any 
algorithm, and earth() automatically transforms factors to binary features.

However, you may find 120 variables is not much data.


Andrew



-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of Lorenzo Isella
Sent: Thursday, February 17, 2011 7:14 AM
To: r-help
Subject: [R] Categorical Variables and Machine Learning

Dear All,
Please consider a dataframe like the one below (I am showing only a few
rows).

>  role degree strength weight count disparity intermittency
>P 10   82  18017 2  2.317073  5.550314e-05
>P  7  529   434560  5.178466  6.904488e-03
>P  8  609   438210  6.204535  1.141031e-03
>D 42  230   691088  1.791153  6.367583e-03

You have a categorical variable (the role variable) which can assume
only a few values ("P","D","C","N","A") referring to different
individuals for whom you collect some extra properties (namely, degree,
strength, weight, disparity and intermittency, like in the table above).
My goal is to find the most suitable property (or combination of
properties) to guess the role of an individual. It looks like a typical
machine learning problem, but I have categorical variables to predict.
I am drowning in the wealth of R packages for machine learning, but I
really would like something simple and easy to use (consider that the
dataset covers only 120 individuals, so performance is not a problem).
Any suggestion is appreciated.
Cheers

Lorenzo

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Pre-allocation of matrices is LESS efficient?

[Alex F. Bokov]

Motivation: during each iteration, my code needs to collect tabular data (and 
use it only during that iteration), but the rows of data may vary. I thought I 
would speed it up by preinitializing the matrix that collects the data with 
zeros to what I know to be the maximum number of rows. I was surprised by what 
I found...

# set up (not the puzzling part)
x<-matrix(runif(20),nrow=4); y<-matrix(0,nrow=12,ncol=5); foo<-c();

# this is what surprises me... what the?
> system.time(for(i in 1:10){n<-sample(1:4,1);y[1:n,]<-x[1:n,];});
   user  system elapsed 
  1.510   0.000   1.514 
> system.time(for(i in 1:10){n<-sample(1:4,1);foo<-x[1:n,];});
   user  system elapsed 
  1.090   0.000   1.085

These results are very repeatable. So, if I'm interpreting them correctly, 
dynamically allocating 'foo' each time to whatever the current output size is 
runs faster than writing to a subset of a preallocated 'y'? How is that 
possible?

And, more generally, I'm sure other people have encountered this type of 
situation. Am I reinventing the wheel? Is there a best practice for storing 
temporary loop-specific data?

Thanks.

PS:  By the way, though I cannot write to foo[,] because the size is different 
each time, I tried writing to foo[] and the runtime was worse than either of 
the above examples.

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

[Pfaff, Bernhard Dr.]

Hello Marta,
 
arrg, sorry, I have not carefully enough read your message. Well, if you follow 
the cited thread further down, you will find: 
(and hereby directly quoting from: 
https://stat.ethz.ch/pipermail/r-sig-finance/2009q2/004274.html) 

<> 
On Tue, 9 Jun 2009, Matthieu Stigler wrote:

> Hi
>
> I wasn't aware of the fact that HAC is not designed for time series model 
> (thanks Achim!). But nevertheless I think that HC is still usable, well at > 
> least I saw it in couple of papers dealing with times series.
>
> So if you still want to use an HC, two solutions:
>
> A solve the problem (workaround):

No, there is no "problem" at least not from the "sandwich" point of view. 
If you want "sandwich" to cooperate with "varest" objects, you just need 
to provide the appropriate methods (essentially, bread() and estfun()) for 
"varest" objects. This is a clean and non-invasive solution and not very 
difficult to implement given that all this is OLS.
<>
 
and further down in this reply, Achim is referring you to the package's 
vignette:

<> 
vignette("sandwich-OOP", package = "sandwich")
<>

Best,
Bernhard
 
 
 
 




Von: Marta Lachowska [mailto:ma...@upjohn.org] 
Gesendet: Donnerstag, 17. Februar 2011 17:01
An: Pfaff, Bernhard Dr.; r-help@r-project.org
Betreff: Re: AW: [R] VAR with HAC


Thank you for your hint!

I see that there was a thread discussing implementation of what I 
wanted to do (Newey-West standard errors in a VAR context), but that there is a 
conflict due to how the type = "const" is defined in the VAR command: 
https://stat.ethz.ch/pipermail/r-sig-finance/2009q2/004272.html that appears 
not to be resolved. 
 
Best, 
 
Marta
 

>>> "Pfaff, Bernhard Dr."  2/17/2011 
4:31 AM >>>
Hello Marta,

have you read ?coeftest and ? VAR carefully enough? The function does 
expect a lm/glm object for x as argument. Hence, the following does work:

library(vars)
data(Canada)
myvar <- VAR(Canada, p = 2, type = "const")
lapply(myvar$varresult, coeftest)

Best,
Bernhard 

> -Ursprüngliche Nachricht-
> Von: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-project.org] Im Auftrag von Marta Lachowska
> Gesendet: Mittwoch, 16. Februar 2011 16:50
> An: r-help@r-project.org
> Betreff: [R] VAR with HAC
> 
> 
> Hello,
> I would like to estimate a VAR model with HAC corrected 
> standard errors. I tried to do this by using the sandwich 
> package, for example: 
>  
> > library(vars)
> > data(Canada)
> > myvar = VAR(Canada, p = 2, type = "const") coeftest(myvar, vcov = 
> > vcovHAC)
> Error in umat - res : non-conformable arrays
>  
> Which suggests that this function is not compatible with the 
> VAR command. Has anyone tried to modify the code to get HAC 
> corrected standard errors with VAR? Any suggestions are welcome. 
>  
> Thank you. 
>  
> Marta
> 
> [[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.
> 
*
Confidentiality Note: The information contained in this message,
and any attachments, may contain confidential and/or privileged
material. It is intended solely for the person(s) or entity to
which it is addressed. Any review, retransmission, dissemination,
or taking of any action in reliance upon this information by
persons or entities other than the intended recipient(s) is
prohibited. If you received this in error, please contact the
sender and delete the material from any computer.
*



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

[Mohamed Lajnef]

Hi Maria,

Look at  barplot function

? barplot  for help

Regards
M
Le 17/02/11 17:00, Lathouri, Maria a écrit :
> Dear all
>
> I have six variables of the average metal concentrations
>
> Var1 4.77
> Var2 23.5
> Var3 5.2
> Var4 12.3
> Var5 42.1
> Var6 121.2
>
> I want to plot them as a barplot with error bars. Could you help me?
>
> Cheers
> Maria
>
>   [[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.
>


-- 

Mohamed Lajnef,IE INSERM U955 eq 15#
Pôle de Psychiatrie#
Hôpital CHENEVIER  #
40, rue Mesly  #
94010 CRETEIL Cedex FRANCE #
mohamed.laj...@inserm.fr   #
tel : 01 49 81 31 31 (poste 18467) #
Sec : 01 49 81 32 90   #
fax : 01 49 81 30 99   #




[[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] barplot with errorbars

[Lathouri, Maria]

Dear all

I have six variables of the average metal concentrations

Var1 4.77
Var2 23.5
Var3 5.2
Var4 12.3
Var5 42.1
Var6 121.2

I want to plot them as a barplot with error bars. Could you help me?

Cheers
Maria

[[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] does range of values in array include a third value?

[David Winsemius]

On Feb 17, 2011, at 10:36 AM, Maas James Dr (MED) wrote:



I'm using the range command to get the minimum and maximum values of  
an array as in


x <- range(array_y)

which gives me two values such as

[1] -2 9

I need to be able to test if this range of values includes a third  
value.  For example I'd like to query


1) does the range of -2 to 9 include 3, answer TRUE
2) does the range of -2 to 9 include -6, answer FALSE?

All values could be negative or positive.  Is there a R function  
that will test this or do I need to programme it?  I have searched  
but not found one.


There is a function that handles intervals well, findInterval:

> findInterval(0, c(-2, 9))
[1] 1
> findInterval(-3, c(-2, 9))
[1] 0
> findInterval(10, c(-2, 9))
[1] 2

So:

isxInRange_y <- function(x, y) findInterval(x, range(y)) == 1

If you want to omit NA's which would otherwise poison the effort, you  
need to wrap the y argument in na.omit().


--
David.



Thanks

J

===
Dr. Jim Maas
University of East Anglia

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


David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] does range of values in array include a third value?

[Peter Ehlers]

On 2011-02-17 07:36, Maas James Dr (MED) wrote:


I'm using the range command to get the minimum and maximum values of an array 
as in

x<- range(array_y)

which gives me two values such as

[1] -2 9

I need to be able to test if this range of values includes a third value.  For 
example I'd like to query

1) does the range of -2 to 9 include 3, answer TRUE
2) does the range of -2 to 9 include -6, answer FALSE?

All values could be negative or positive.  Is there a R function that will test 
this or do I need to programme it?  I have searched but not found one.


Here too you can use findInterval():

  findInterval( -2, range(array_y) ) == 1

You may want to set the 'rightmost.closed' argument to TRUE.

Peter Ehlers



Thanks

J

===
Dr. Jim Maas
University of East Anglia

[[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] does range of values in array include a third value?

[Maas James Dr (MED)]

I'm using the range command to get the minimum and maximum values of an array 
as in

x <- range(array_y)

which gives me two values such as

[1] -2 9

I need to be able to test if this range of values includes a third value.  For 
example I'd like to query

1) does the range of -2 to 9 include 3, answer TRUE
2) does the range of -2 to 9 include -6, answer FALSE?

All values could be negative or positive.  Is there a R function that will test 
this or do I need to programme it?  I have searched but not found one.

Thanks

J

===
Dr. Jim Maas
University of East Anglia

[[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] Dependency on R-Forge Package

[Uwe Ligges]

On 17.02.2011 15:29, Damian Betebenner wrote:

In building a package, is it possible to make the package depend upon another 
package that is only available on R-Forge (not CRAN). For example, by doing 
something in the DESCRIPTION file
I'd like to add a dependency to my package such that when the user install my 
package it will automatically install this other package from R-forge as well.



Not completely, the user has to choose R-forge as a repository (which is 
not the default but easily selectable).


Uwe Ligges






Any help greatly appreciated.

Damian


[[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] Saturated model in binomial glm

[Giovanni Petris]

Dear Bill,

Thank you very much for your careful discussion of the issue.

It is not surprising that the deviance is the same whether you fit the
model using a factor response with weights or individual 0/1 responses.
I think this happens because the fitted probabilities in the saturated
models are zero or one (depending on the response being a failure or a
success) in both cases. The degrees of freedom are clearly different.

However, if you consider as response the two-column matrix of successes
and failures then you have a different deviance, with different degrees
of freedom:

> covariates <- subset(UCBAd, subset = Admit == "Admitted",
+  select = -c(Admit, Freq))
> AdmitReject <- as.matrix(cbind(subset(UCBAd, subset = Admit == "Rejected",
+   select = Freq),
+subset(UCBAd, subset = Admit == "Admitted",
+   select = Freq)))
> UCBAd_matrix <- glm(AdmitReject ~ Gender + Dept, family = binomial,
+ data = covariates) ## matrix as response
> UCBAd_matrix$deviance
[1] 20.20428
> UCBAd_matrix$df.residual
[1] 5
> coef(UCBAd_matrix)
 (Intercept) GenderFemaleDeptBDeptCDeptDDeptE 
 -0.58205140  -0.09987009   0.04339793   1.26259802   1.29460647   1.73930574 
   DeptF 
  3.30648006 

In this case the fitted probabilities in the saturated model are the
observed proportions. The fitted coefficients are clearly the same.

The deviance (and degrees of freedom) resulting from this way of fitting
the model _can_ be used as a goodness-of-fit statistic, based on
standard arguments, while the deviance (and degrees of freedom)
resulting from fitting the model as individual 0/1 responses _can't_.

What about deviance and degrees of freedom that I get using the "factor
with weights" approach? The number of parameters in the saturated model
does not increase with the number of observations, as it happens in the
individual 0/1 responses case, but the fitted probabilities in the
saturated model are fixed to zero and one. I guess it is not clear to me
that the standard likelihood ratio test asymptotic theory holds in this
case. Could anybody clarify?

Thank you in advance,
Giovanni


On Thu, 2011-02-17 at 10:50 +1100, bill.venab...@csiro.au wrote:
> This is a very good question.  You have spotted something that not
> many people see and it is important.
> 
> The bland assertion that the "deviance can be used as a test of fit"
> can be seriously misleading.
> 
> For this data the response is clearly binary, "Admitted" (success) or
> "Rejected" (failure) and the other two factors are explanatory
> variables.  
> 
> Any binomial model can be fitted by expressing the data as a binary
> response.  In this case there are
> 
> > sum(UCBAd$Freq)
> [1] 4526
> > 
> 
> 4526 trials, corresponding to the individual applicants for admission.
> We can expand the data frame right out to this level and fit the model
> with the data in that form, and in this case the weights will be the
> default, ie. all 1's.
> 
> We can also *group* the data into subsets which are homogeneous with
> respect to the explanatory variables.  
> 
> The most succinct grouping would be into 12 groups corresponding to
> the 2 x 6 distinct classes defined by the two explanatory variables.
> In this case you would specify the response either as a two-column
> matrix of successes/failures, or as a proportion with the totals for
> each of the 12 cases as weights.
> 
> Another succince grouping is into 2 x 2 x 6 classes as you do in your
> example.  In this case your response is the factor and the weights are
> the frequencies.
> 
> For all three cases a) the estimates will be the same, and so the
> predictions will be identical, b) the deviance will also be the same,
> but c) the degrees of freedom attributed to the deviance will be
> different.
> 
> The reason for c) is, as you have intuited, the saturated model is
> different.  Literally, the saturated model is a model with one mean
> parameter for each value taken as an observation when the model is
> fitted.  So the saturated model is *not* invariant with respect to
> grouping.
> 
> Let's look at two of these cases computationally:
> 
> 
> > UCB_Expanded <- UCBAd[rep(1:24, UCBAd$Freq), 1:3] ## expand the data frame
> > dim(UCB_Expanded)
> [1] 45263
> 
> >  ### now fit your model
> 
> > m1 <- glm(Admit ~ Gender + Dept, binomial, UCBAd, weights = Freq)
> 
> >  ### and the same model using the binary data form
> 
> > m2 <- glm(Admit ~ Gender + Dept, binomial, UCB_Expanded)
> 
> >  ### as predicted, the coefficients are identical (up to round off)
> 
> > coef(m1)
>  (Intercept) GenderFemaleDeptBDeptCDeptDDeptE 
>  -0.58205140  -0.09987009   0.04339793   1.26259802   1.29460647   1.73930574 
>DeptF 
>   3.30648006 
> > coef(m2)
>  (Intercept) GenderFemaleDeptBDeptCDeptDDept

[R] does range of values in array include a third value?

[Maas James Dr (MED)]

I'm using the range command to get the minimum and maximum values of an array 
as in

x <- range(array_y)

which gives me two values such as 

[1] -2 9

I need to be able to test if this range of values includes a third value.  For 
example I'd like to query 

1) does the range of -2 to 9 include 3, answer TRUE
2) does the range of -2 to 9 include -6, answer FALSE?

All values could be negative or positive.  Is there a R function that will test 
this or do I need to programme it?  I have searched but not found one.

Thanks

J

===
Dr. Jim Maas
University of East Anglia

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

[Uwe Ligges]

On 16.02.2011 22:38, poisontonic wrote:



Uwe Ligges-3 wrote:


If the available space got too fragmented, there is not single 3.8 block
of memory available any more



Is there anything I can do to prevent this?


If you did it after a fresh reboot: I don't see a way to prevent it.
Nevertheless, I doubt you really have that much memory free in that 
case. Have you inspected how much memory was already allocated by R?


uwe Ligges


I've restarted and rerun the

whole thing straight up, and still the error...?

Ben


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