[R] convert delimited strings with ranges to numeric

2013-08-14 Thread Chris Stubben
Is there an easy way to convert character strings with comma-separated 
numbers and ranges to a numeric vector?


x<-  "2,5-7,10,12-15"

[1]  2  5  6  7 10 12 13 14 15

Thanks,
Chris


--

Chris Stubben

Los Alamos National Lab
Bioscience Division
MS M888
Los Alamos, NM 87545

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] [XML packages] how to get the sub-node according to the sub-node's attribute?

2013-07-12 Thread Chris Stubben

My question is:
Is there a function that can get the sub-node according to the sub-node's 
attribute ?
like that,in the mtcar.xml ,there is sub-node as follow:
30.4   4  95.1 113 3.77 1.513 16.90  1  15
2
I want to get the specific sub-node according to the attribute "id="Lotus 
Europa".
Is there a simple function in the xml package?



Just use xpathSApply with brackets to specify the specific id...

doc <- xmlParse(system.file("exampleData", "mtcars.xml", package="XML"))

xpathSApply(doc, "//record[@id='Lotus Europa']", xmlValue)
[1] "30.4   4  95.1 113 3.77 1.513 16.90  1  152"

#OR?
xpathSApply(doc, "//record[@id='Lotus Europa']")
getNodeSet(doc, "//record[@id='Lotus Europa']")

# And to list all Ids...

xpathSApply(doc, "//record", xmlGetAttr, "id")
[1] "Mazda RX4"   "Mazda RX4 Wag"   "Datsun 710"   .   


--

Chris Stubben

Los Alamos National Lab
Bioscience Division
MS M888
Los Alamos, NM 87545

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

2013-05-08 Thread Chris Stubben
Sorry to answer my own question - I guess here's one way to read this 
table.  Other suggestions are still welcome.


Chris

--

x<-htmlParse("
abX
YZ
c
")

# split by rows
z <- getNodeSet(x, "//tr")

# create empty data.frame - probably not the best solution...
t1<- data.frame(matrix(NA, nrow = 3,  ncol = 2 ))

for (i in 1:3){
  rowspan <- as.numeric( xpathSApply(z[[i]], ".//td", xmlGetAttr, 
"rowspan", 1) )

  val <- xpathSApply(z[[i]], ".//td", xmlValue)

  # fill values into empty cells
  n <- which(is.na(t1[i,]))
  t1[ i ,n] <- val

  if( any(rowspan > 1) ){
 for(j in 1:length( rowspan ) ){
if(rowspan[j] > 1){
## repeat value down column
  t1[ (i+1):(i+ ( rowspan[j] -1) ) , n[j] ]   <- val[j]
}
 }
  }
}


t1
 X1 X2
1 ab  X
2 ab YZ
3  c YZ


If you are interested, I used this code in the pmcTable function at 
https://github.com/cstubben/pubmed .  To get  Table 1, this now works...


doc<-pmc("PMC3544749")  # downloads XML from OAI service
t1 <- pmcTable(doc,1)  # parse table... also saves caption and footnotes 
to attributes

t1[1:4,1:4]
  Category Gen Name Rv 
number  Description
1 Lipids and Fatty Acid Metabolism kasBRv2246 
3-oxoacyl-[acyl-carrier protein] synthase 2 kasb
2   Mycolic acid synthesismmaA4   Rv0642c  
Methoxy mycolic acid synthase 4
3   Mycolic acid synthesis pcaA   Rv0470cMycolic acid 
synthase (cyclopropane synthase)
4   Mycolic acid synthesis pcaA   Rv0470cMycolic acid 
synthase (cyclopropane synthase)





--

Chris Stubben

Los Alamos National Lab
Bioscience Division
MS M888
Los Alamos, NM 87545

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] create unique ID for each group

2013-05-07 Thread Chris Stubben

Yes, I tried, but the order of the IDs in dat1 and dat2 is not exactly the
same, I simplify the data here. So in dat2, it may have records for ID=0002
first then ID=0001, also I have more than two categories under ID col


I should have looked at the question more closely, sorry.   Unique ids in raw 
datasets
are pretty important, especially if observations are split into different files and 
you are trying to join them later.  How do you know for ID 0001 and obs 1 that height 
is  3.2 and not 2.6, especially if order in the two files are "not exactly the same".


Chris



--

Chris Stubben

Los Alamos National Lab
Bioscience Division
MS M888
Los Alamos, NM 87545

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] recode categorial vars into binary data

2013-05-07 Thread Chris Stubben


First off, stop using cbind() when it is not needed. You will not see the reason 
when the columns are all numeric but you will start experiencing pain and puzzlement 
when the arguments are of mixed classes. The data.frame function will do what you want. 
(Where do people pick up this practice anyway?)


Maybe from help( data.frame)?

It's in most of the  examples and is not needed ...

L3 <- LETTERS[1:3]
(d <- data.frame(cbind(x=1, y=1:10), fac=sample(L3, 10, replace=TRUE)))

## The same with automatic column names:

data.frame(cbind(  1,   1:10), sample(L3, 10, replace=TRUE))



Chris




--

Chris Stubben

Los Alamos National Lab
Bioscience Division
MS M888
Los Alamos, NM 87545

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] create unique ID for each group

2013-05-07 Thread Chris Stubben

> I want to merge dat1 and dat2 based on "ID" in order

Have you tried merge(dat1, dat2) ?

If ID is the common column (and no others), then that should be all you 
need to join (see ?merge).  And then order if needed.


Chris



--

Chris Stubben

Los Alamos National Lab
Bioscience Division
MS M888
Los Alamos, NM 87545

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

2013-05-06 Thread Chris Stubben
I'm trying to read html tables with lots of rowspan attributes, for 
example...


x<-htmlParse("
  abX
  YZ
  c
")

readHTMLTable(x, which=1)
 V1   V2
1 abX
2 YZ 
3  c 

Does anyone know how to use the rowspan attributes and repeat cell 
values  to format a table like this?


 V1   V2
1 abX
2 ab   YZ
3  c   YZ

Also, the actual tables I'm using are large, for example, this one has 
206 rows and rowspan attributes ranging from 2-14 scattered in all 8 
columns, so the shifted rows in t1 are not very useful right now.


t1 <- readHTMLTable( 
"http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3544749/table/T1";, which=1)


Thanks,
Chris












t1<-readHTMLTable( 
"http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3544749/table/T1";, which=1)



--

Chris Stubben

Los Alamos National Lab
Bioscience Division
MS M888
Los Alamos, NM 87545

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

2012-04-13 Thread Chris Stubben
Sorry if I was not clear.  I wanted to remove the superscripts using xpath
queries if possible.  For example this will get p nodes with superscripts,
but how do I remove the superscripts if there are many matching nodes and
different superscripts?

xpathSApply(doc, "//p[sup]", xmlValue) 
[1] "Cata"


Chris

--
View this message in context: 
http://r.789695.n4.nabble.com/Remove-superscripts-from-HTML-objects-tp4550738p4555370.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] Remove superscripts from HTML objects

2012-04-11 Thread Chris Stubben
Is there some way to remove superscripts from objects returned by
html/xmlParse (XML package)?

h <- "CataDog"
doc <- htmlParse(h)
 xpathSApply(doc, "//p", xmlValue)
[1] "Cata" "Dog"

I could probably remove the   tags from the "h" object above, but I'd
rather just work with the results from htmlParse if possible (and not use
readLines to load raw HTML first).

Thanks,
Chris Stubben
 


--
View this message in context: 
http://r.789695.n4.nabble.com/Remove-superscripts-from-HTML-objects-tp4550738p4550738.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] Getting codebook data into R

2012-02-13 Thread Chris Stubben
Just to follow up on Dan's code - once you have a data.frame listing column
positions, then it's just a couple steps to download the file...

x <- data.frame(name=c('caseid', 'nbrnaliv', 'babysex',
'birthwgt_lb','birthwgt_oz','prglength',
'outcome', 'birthord',  'agepreg',  'finalwgt'),
begin = c(1, 22, 56, 57, 59, 275, 277, 278, 284, 423),
end =  c(12, 22, 56, 58, 60, 276, 277, 279, 287, 440)
)


x$width <- x$end - x$begin + 1
x$skip <-  (-c(x$begin[-1]-x$end[-nrow(x)]-1,0))

widths <- c(t(x[,4:5]))
widths <- widths[widths!=0]

ftp<-
"ftp://ftp.cdc.gov/pub/Health_Statistics/NCHS/Datasets/NSFG/2002FemPreg.dat";
# drop the n=10 option to get all lines
y<- read.fwf(ftp, widths, n=10)
names(y) <- x$name
y
   caseid nbrnaliv babysex birthwgt_lb birthwgt_oz prglength outcome
birthord agepreg  finalwgt
1   11   1   8  1339   1   
13316  6448.271
2   11   2   7  1439   1   
23925  6448.271
3   23   1   9   239   1   
11433 12999.542
4   21   2   7   039   1   
21783 12999.542
5   21   2   6   339   1   
31833 12999.542
6   61   1   8   938   1   
12700  8874.441
7   61   2   9   940   1   
22883  8874.441
8   61   2   8   642   1   
33016  8874.441
9   71   1   7   9    39   1   
12808  6911.880
10  71   2   6  1035   1   
23233  6911.880


Chris Stubben




Daniel Nordlund-4 wrote
> 
>> -Original Message-
> 
>> I've been trying to get some data from the National Survey for Family
>> Growth
>> into R - however, the data is in a .dat file and the data I need doesn't
>> have any spaces or commas separating fields - rather you have to look
>> into
>> the codebook and what number of digits along the line the data you need
>> is.
>> The data I want are the following, where 1,12,int means that the data I'm
>> interested starts in column 1 and finishes in column 12 and is an
>> integer.
>> 
>> ('caseid', 1, 12, int),
>>  ('nbrnaliv', 22, 22, int),
>> ('babysex', 56, 56, int),
>> ('birthwgt_lb', 57, 58, int),
>> ('birthwgt_oz', 59, 60, int),
>> ('prglength', 275, 276, int),
>> ('outcome', 277, 277, int),
>> ('birthord', 278, 279, int),
>> ('agepreg', 284, 287, int),
>> ('finalwgt', 423, 440, float)
>> 
>> How can I do this using R? I've written a python programme which
>> basically
>> does it but it'd be nicer if I could skip the Python bit and just do it
>> using R. Cheers for any help.
>> 
> 
> 
> Dan
> 
> 
> 


--
View this message in context: 
http://r.789695.n4.nabble.com/Getting-codebook-data-into-R-tp4374331p4386135.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] Markov transition matrices , missing transitions for certain years

2011-04-19 Thread Chris Stubben
To keep the table dimensions the same, try changing the columns to factors...


boxes$y97<-factor(boxes$y97, 1:4)
boxes$y98<-factor(boxes$y98, 1:4)
boxes$y99<-factor(boxes$y99, 1:4)
...

> table(boxes$y98, boxes$y97)
   
1 2 3 4
  1 1 0 1 0
  2 0 0 0 0
  3 0 0 0 0
  4 1 0 0 1
> table(boxes$y99, boxes$y98)
   
1 2 3 4
  1 1 0 0 1
  2 0 0 0 0
  3 0 0 0 0
  4 1 0 0 1

Chris


Abby_UNR wrote:
> 
> My problem is that there are 16 possible transitions, but not all possible
> transitions occur at each time step. Therefore, don't think I could do
> something easy like create a table for each time step and add them
> together, for example:
> 
>> t1.boxes <- table(boxes$y98, boxes$y97)
>> t1.boxes
>
> 1 3 4
>   1 1 1 0
>   4 1 0 1
>> t2.boxes <- table(boxes$y99, boxes$y98)
>> t2.boxes
>
> 1 4
>   1 1 1
>   4 1 1
>> 
> 
> 


--
View this message in context: 
http://r.789695.n4.nabble.com/Markov-transition-matrices-missing-transitions-for-certain-years-tp3459072p3460651.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] rotate column names in large matrix

2010-11-16 Thread Chris Stubben

Just increase the margins on the left side and add the rownames

x <- cor(matrix(rnorm(600), 60, 100))
rownames(x)<-paste("row", 1:100)
op<-par(mar=c(1,5,1,1), xpd=TRUE)
image(t(x[nrow(x):1,]), axes=FALSE)
 text(-0.01, seq(0,1,length=nrow(x) ), rownames(x), pos = 2,  offset = 0,
cex = .7)

Another option is to search the R graphics manual at
http://rgm2.lab.nig.ac.jp/RGM2/images.php and maybe you'll find a function
in all those packages that plots what you need.   Try searching for "symnum"
or "image key" or "matrix plot" (in quotes).  For "matrix plot", there's a
function in the plotrix package that will add a key or even show values in
the "cells"..

library(plotrix)
color2D.matplot(x, show.legend=TRUE, red=1, blue=0)

Chris




Lara Poplarski wrote:
> 
> Chris, could you please suggest how to modify what you sent to also show
> the
> same labels as (horizontal) row names? I have not yet mastered the details
> of R graphics...
> 

-- 
View this message in context: 
http://r.789695.n4.nabble.com/rotate-column-names-in-large-matrix-tp3043493p3045366.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] rotate column names in large matrix

2010-11-15 Thread Chris Stubben

You could display the matrix as an image then add column names (rotated used
srt) if you really want them. 

x <- cor(matrix(rnorm(600), 60, 100))  

# set margins with extra space at top and xpd=TRUE to write outside plot
region
op<-par(mar=c(1,1,5,1), xpd=TRUE)

# display image without the 90 degree counter clockwise rotation
image(t(x[nrow(x):1,]), axes=FALSE)

## add 100 column names
y<-paste("column", 1:100)
text( seq(0,1,length=100) , 1.01, y,  pos = 2, srt = 270, offset=0, cex=.7)


Chris Stubben





Lara Poplarski wrote:
> 
> I have a large (1600*1600) matrix generated with symnum, that I am using
> to
> eyeball the structure of a dataset.
> 
> I have abbreviated the column names with the abbr.colnames option. One way
> to get an even more compact view of the matrix would be to display the
> column names rotated by 90 degrees.
> 
> Any pointers on how to do this would be most useful. Any other tips for
> displaying the matrix in compact form are of course also welcome.
> 
> 
> 

-- 
View this message in context: 
http://r.789695.n4.nabble.com/rotate-column-names-in-large-matrix-tp3043493p3043927.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] R eat my data

2010-05-25 Thread Chris Stubben

Gene names often have single quotes like

5'-methylthioadenosine phosphorylase  
ATP synthase B' chain
ppGpp 3'-pyrophosphohydrolase

so maybe  try adding quote="" to the read table options.

Chris Stubben
-- 
View this message in context: 
http://r.789695.n4.nabble.com/R-eat-my-data-tp2230217p2230303.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] Using substitute in a print method

2010-03-19 Thread Chris Stubben

I've defined a new class for some really large tables, and I'd like to use
substitute() in the print function.  For example, if I define a new class
below and create a print method for it, I can't figure out how to return the
object name by just typing "d". 

d<-diag(2)
class(d)<-c("diag")
print.diag<-function(x,...)
  {
# do something then
print(paste( "A diagonal matrix",  substitute(x)) )
  }

> d
[1] "A diagonal matrix 1" "A diagonal matrix 0" "A diagonal matrix 0" "A
diagonal matrix 1"

##  This prints what I'd want, but I'm not sure why its different.  
print.diag(d)
[1] "A diagonal matrix d"

Thanks for any help,

Chris Stubben

-- 
View this message in context: 
http://n4.nabble.com/Using-substitute-in-a-print-method-tp1606635p1606635.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] popbio and stochastic lambda calculation

2010-02-18 Thread Chris Stubben


Shawn Morrison-2 wrote:
> 
> # The paper reports a 95% CI of 0.79 - 1.10
> # "My" reproduced result for the CIs is much larger, especially on the 
> upper end. Why would this be?
> # The authors report using the 'delta' method (Caswell, 2001) to 
> calculate the CI in which the
> 
> 

Shawn,

I probably can't help much with the vitalsim example, but I would check Box
8.10 in Morris and Doak (2002).

I do have a few ideas about the delta method below.

# List the vital rates

vr<-list(cub= 0.64, yly = 0.67, sub=0.765, adt=0.835, mx=0.467)

# and the matrix using an expression in R

el<- expression(
 0,  0,  0,  0,  adt*mx, 
 cub,0,  0,  0,  0, 
 0,  yly,0,  0,  0, 
 0,  0,  sub,0,  0, 
 0,  0,  0,  sub,adt)

# this should get the projection matrix

A<-matrix( sapply( el, eval, vr), nrow=5, byrow=TRUE)

lambda(A)
[1] 0.9534346

# use the vitalsens function to get the vital rate sensitivites and
  save the second column

vitalsens(el, vr)
estimate sensitivity elasticity
cub0.640   0.1236186 0.08298088
yly0.670   0.1180835 0.08298088
sub0.765   0.2068390 0.16596176
adt0.835   0.7628261 0.66807647
mx 0.467   0.1694131 0.08298088

sens<-vitalsens(el, vr)[,2]


# I'm not sure about the covariance matrix next.  In Step 7 in Slakski et al
2007 ("Calculating the variance of the finite rate of population change from
a matrix model in Mathematica") they just use the square of the standard
errors, so I'll do the same...

se<-list(cub= 0.107, yly = 0.142, sub=0.149, adult=0.106, mx=0.09)
cov<-diag(unlist(se)^2)

## and then the variance of lambda from step 8
var<-t(sens) %*% ( cov%*%sens)
[,1]
[1,] 0.008176676

# and the confidence intervals

lambda(A) - 1.96*sqrt(var)
lambda(A) + 1.96*sqrt(var)

CI of 0.78 and 1.13 is close to paper

Hope that helps,

Chris




-- 
View this message in context: 
http://n4.nabble.com/popbio-and-stochastic-lambda-calculation-tp1557647p1560745.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] splitting scientific names into genus, species, and subspecies

2009-11-04 Thread Chris Stubben


Mark W. Miller wrote:
> 
> I have a list of scientific names in a data set.  I would like to split
> the names into genus, species and subspecies.  Not all names include a
> subspecies.  Could someone show me how to do this?
> 

strsplit should work for your example...

data.frame( 
  genus=sapply(strsplit(aa, " "), "[", 1),
species=sapply(strsplit(aa, " "), "[", 2),
subspecies=sapply(strsplit(aa, " "), "[", 3)   ## will be NA for missing
subsp 
 ) 

However, scientific names are often pretty messy - I often have datasets
like this...
x
 [1] "Aquilegia caerulea James var. caerulea"   
  
 [2] "Aquilegia caerulea James var. ochroleuca Hook."   
  
 [3] "Aquilegia caerulea James var. pinetorum (Tidestrom) Payson ex Kearney
& Peebles"
 [4] "Aquilegia caerulea James" 
  
 [5] "Aquilegia chaplinei Standl."  
  
 [6] "Aquilegia chaplinei Standley ex Payson"   
  
 [7] "Aquilegia chrysantha Gray var. chrysantha"
  
 [8] "Aquilegia chrysantha Gray"   

So I first strip out author names using strsplit and use grep to find
subspecies/variety abbreviations 

noauthor<-function(x){
  ## split name into vector of separate words
  y<-strsplit(x, " ")
  sapply(y, function(x){  
n<-grep( "^var\\.$|^ssp\\.$|^var$|^f\\.$",x)
# apply a function to paste together the first and second elements
# plus element after matching var., spp., f. (or and others) 
# use sort in case the name includes both var and spp -sometimes happens
paste( x[sort(c(1:2, n,n+1))], collapse=" ")  })}


noauthor(x[1:8])
[1] "Aquilegia caerulea var. caerulea"
[2] "Aquilegia caerulea var. ochroleuca"  
[3] "Aquilegia caerulea var. pinetorum"   
[4] "Aquilegia caerulea"  
[5] "Aquilegia chaplinei" 
[6] "Aquilegia chaplinei" 
[7] "Aquilegia chrysantha var. chrysantha"
[8] "Aquilegia chrysantha"


Chris






-- 
View this message in context: 
http://old.nabble.com/splitting-scientific-names-into-genus%2C-species%2C-and-subspecies-tp26204666p26205654.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] Matrix multiplication and random numbers

2009-09-09 Thread Chris Stubben


RFish wrote:
> 
> I new to using R and am struggling with some matrix multiplication. 
> 

I'm not sure what you're trying to print, but you could place this vector in
an expression

mat3<-expression(c(0,rnorm(1,0.6021,0.0987),0,0,0,0,0,0,0,rnorm(1,0.6021,0.0987),0,0,0,0,1.9,0,0,rnorm(1,0.6021,0.0987),0,0,0,4.8,0,0,0,rnorm(1,0.6021,0.0987),0,0,9.7,0,0,0,0,rnorm(1,0.6021,0.0987),0,18,0,0,0,0,0,rnorm(1,0.6021,0.0987),32.6,0,0,0,0,0,0))

# and then evaluate to get a new matrix each time
matrix(eval(mat3), nrow=7)

#I think this may be easier to follow. First create a matrix of zeros, stick
in fertilities and then add random survival probabilities each time


mat3<-diag(0,7)
#fertilities
mat3[1,3:7]<-c(1.9, 4.8, 9.7, 18, 32.6)
# random survival on sub-diagonal
mat3[row(mat3)==col(mat3)+1]<-rnorm(6,0.6021,0.0987)


# and if you want to project the population over 10 time steps in a loop ?

n<-matrix(c(500,0,0,0,0,0,0))

popsize <- matrix(numeric(7 * 10), nrow = 7)
 for (i in 1:10) {
popsize[, i] <- n
 mat3[row(mat3)==col(mat3)+1]<-rnorm(6,0.6021,0.0987)
 n <- mat3 %*% n
}


 [,1] [,2] [,3] [,4]  [,5]  [,6]   [,7] 
[,8]
[1,]  500   0.   0. 531.6256 709.89940 940.19337 1697.52862
3403.6610
[2,]0 352.5116   0.   0. 298.97874 424.71160  561.32525
1027.1605
[3,]0   0. 279.8029   0.   0.0 231.45988  316.83352 
424.8883
[4,]0   0.   0. 147.8957   0.0   0.0  136.36804 
220.7370
[5,]0   0.   0.   0.  96.92715   0.00.0 
108.6551
[6,]0   0.   0.   0.   0.0  69.875270.0   
0.
[7,]0   0.   0.   0.0000   0.0   0.0   65.86229   
0.


Chris Stubben

-- 
View this message in context: 
http://www.nabble.com/Matrix-multiplication-and-random-numbers-tp25365184p25369495.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] Stochastic (transition) matrices: how to determine distributions and variance?

2009-08-31 Thread Chris Stubben


Jonathan Greenberg-2 wrote:
> 
> I am trying to determine the distribution and variance for a classic 
> stochastic (transition) matrix 
> 

I have a few suggestions that I've tried with matrix population models
(so probabilities are survival rates plus dead fates).

# here's a test example

A<-matrix(c(.1,.3,.6,.2,.4,.4, 0,0,1), nrow=3)
n<-c(100, 300, 200)

A %*% n
 [,1]
[1,]   70
[2,]  150
[3,]  380
 
# columns should add to 1, correct?
colSums(A)


#  two possible solutions

# option 1.   sample A from a multinomial distribution 


apply(A, 2, function(x) rmultinom(1, size = 1000, prob=x) /1000)
  [,1]  [,2] [,3]
[1,] 0.090 0.1950
[2,] 0.322 0.4010
[3,] 0.588 0.4041


n1<-matrix(numeric(1000*3), ncol=3)
for(i in 1:1000)
{
A2<-apply(A, 2, function(x) rmultinom(1, size = 1000, prob=x) /1000)
n1[i,] <-A2  %*% n
}

apply(n1, 2, quantile, c(0.025,0.975))
  [,1]  [,2] [,3]
2.5%  62.3 140.5 370.4975
97.5% 78.3 159.6 389.8000


# option 2.  Since I usually don't know anything about the distribution of
values in A, I
would resample transitions using a bootstrap


# Get number of each transition

a2<-data.frame( round ( as.table(t(A)*n)))
  Var1 Var2 Freq
1AA   10
2BA   60
3CA0
4AB   30
5BB  120
6CB0
7AC   60
8BC  120
9CC  200

## expand into 600 rows (so A->A repeats 10 times, B->A repeats 60 times,
etc)

r <- rep(row.names(a2), a2$Freq)
y <- a2[r, 1:2 ]
rownames(y)<-1:nrow(y)

##  this gets the original matrix
prop.table( table(y[,2], y[,1]), 2)
   
# so resample rows in y 

z <- sample(nrow(y), replace = TRUE)
y1<-y[z,]
prop.table( table(y1[,2], y1[,1]), 2)  %*% n

## and repeat 1000 times

n2<-matrix(numeric(1000*3), ncol=3)

for(i in 1:1000)
{
z <- sample(nrow(y), replace = TRUE)
y1<-y[z,]
n2[i,] <-prop.table( table(y1[,2], y1[,1]), 2)  %*% n
}

# the CIs are much wider here

apply(n2, 2, quantile, c(0.025,0.975))
  [,1] [,2] [,3]
2.5%  55.89959 132.6868 361.2091
97.5% 85.11097 169.3798 399.7032


Chris Stubben
-- 
View this message in context: 
http://www.nabble.com/Stochastic-%28transition%29-matrices%3A-how-to-determine-distributions-and-variance--tp25209297p25227303.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] Long character not truncated by str() after adding comment

2009-07-22 Thread Chris Stubben

If I add a comment to a long character string, I can't get the truncated
string to display in the str() output.  Anyone know how to get str to
compactly display a truncated string and the comment? 

longch <- paste(rep(letters,100), collapse="")
str(longch)  # truncated
comment(longch)<-"100 ABCs"
str(longch) # not truncated


Thanks,

Chris Stubben
-- 
View this message in context: 
http://www.nabble.com/Long-character-not-truncated-by-str%28%29-after-adding-comment-tp24611767p24611767.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] sampling quadrats of increasing size in a list of matrices

2009-06-30 Thread Chris Stubben


Jens Oldeland wrote:
> 
> I am looking for a solution on how to sample increasing sizes of 
> "quadrats" in a list of matrices. Each matrix is a binary raster map and 
> I want to know if there is a 1 or only 0 s in the sampling unit 
> (quadrat, e.g. 4x4 cells). 
> 


You could just use a few loops to index the submatrices/quadrats, maybe
something like this.  

set.seed(7)
x<-matrix(sample(c(0,1), 256, replace=TRUE, prob=c(.8,.2)),  nrow=16,
byrow=TRUE)
# matrix size (16)
n<-nrow(x)
#quadrat sizes (1  2  4  8 16)
q<-  2^(0:sqrt(n))
# store results in list
quads<-vector("list", length(q) )
names(quads)<-q
for (i in 1:length(q) ) 
{
   z<-c()
   for(j in 1:(n-q[i]+1) )
   {
  for(k in 1:(n-q[i]+1))
  {
 y<-x[j:(j+q[i]-1), k:(k+q[i]-1)]
 #do something with submatrix
 z<-c(z, ifelse( sum(y)>0, 1,0))
 #z<-c(z, sum(y) ) 

  }
   }
   quads[[i]]<-z
}

#empty quads
sapply(quads, function(x) sum(x==0))
  1   2   4   8  16 
206  91   3   0   0 
# total number of quadrats sampled
sapply(quads, length)
#or
(n-q+1)^2


For random sampling, try replacing 1:(n-q[i]+1)  with  sample( 1:(n-q[i]+1),
replace=TRUE) ?

Also, if you have a list of matrices, then maybe combine the code above into
a function and then apply it to your list using lapply.


Chris Stubben

-- 
View this message in context: 
http://www.nabble.com/sampling-quadrats-of-increasing-size-in-a-list-of-matrices-tp24270434p24276950.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 calculate means of matrix elements

2009-05-18 Thread Chris Stubben


dxc13 wrote:
> 
>   If these three matrices are given in separate text files, how can I
> write code that will get this result I need?
> 


If you have matrices in separate text files like mat1.txt, mat2.txt,
mat3.txt, you could load them into a list using a loop

x<- vector('list', 3)

for ( i in 1:3)
{
   ## you may to change some default options for read.table
   x[[i]]<-as.matrix(read.table(paste( 'mat', i, '.txt', sep='')))
}

Then write a function to calculate the mean of a list of matrices

mean(x)


# Paste this function into R before running mean(x) – its also included in
the popbio package

mean.list<-function (x, ...) 
{
if (!all(sapply(x, is.matrix))) 
stop("'x' must be a list containing matrices")
dims <- sapply(x, dim)
n <- dims[1, 1]
p <- dims[2, 1]
if (!all(n == dims[1, ]) || !all(p == dims[2, ])) 
stop("the matrices must have the same dimensions")
mat <- matrix(unlist(x), n * p, length(x))
    mm <- matrix(rowMeans(mat, ...), n, p)
dimnames(mm) <- dimnames(x[[1]])
mm
}


Chris Stubben

-- 
View this message in context: 
http://www.nabble.com/how-to-calculate-means-of-matrix-elements-tp23607694p23609472.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] Looking for a quick way to combine rows in a matrix

2009-05-13 Thread Chris Stubben

You can automate this step
>key[key == "AT"] <- "TA"

## create a function to reverse a string -- see strsplit help page for this
strReverse function
reverse <- function(x) sapply(lapply(strsplit(x, NULL), rev), paste,
collapse="")

key <- rownames(a)
# combine rownames with reverse (rownames)
n<-cbind(key, rev=reverse(key))
 key  rev 
[1,] "AA" "AA"
[2,] "AT" "TA"
[3,] "TA" "AT"
[4,] "TT" "TT"

# Now just sort the values in the rows   (apply returns column vectors so I
also use t() ) and then run do.call on first column
 n<-t(apply(n,1, sort))

do.call(rbind, by(a, n[,1], colSums)) 
   V1 V2 V3 V4
AA  1  5  9 13
AT  5 13 21 29
TT  4  8 12 16


I often need to combine reverse complement DNA strings, so you could do that
too 

# DNA complement
comp <-  function(x) chartr("ACGT", "TGCA", x)

n<-cbind(key, rev=reverse(comp(key)))  
 n<-t(apply(n,1, sort))
do.call(rbind, by(a, n[,1], colSums)) 
   V1 V2 V3 V4
AA  5 13 21 29   
AT  2  6 10 14
TA  3  7 11 15


Chris Stubben


jholtman wrote:
> 
> Try this:
> 
>> key <- rownames(a)
>> key[key == "AT"] <- "TA"
>> do.call(rbind, by(a, key, colSums))
>V2 V3 V4 V5
> AA  1  5  9 13
> TA  5 13 21 29
> TT  4  8 12 16
> 
> 
> On Mon, May 11, 2009 at 4:53 PM, Crosby, Jacy R
> wrote:
> 
>> I'm working with genotype data in a frequency table:
>>
>> > a=matrix(1:16, nrow=4)
>> > rownames(a)=c("AA","AT","TA","TT")
>> > a
>>   [,1] [,2] [,3] [,4]
>> AA159   13
>> AT26   10   14
>> TA37   11   15
>> TT48   12   16
>>
>> 'AT' and 'TA' are essentially the same, and I'd like to combine (add) the
>> rows to reflect this. The final matrix should be:
>>
>>   [,1] [,2] [,3] [,4]
>> AA159   13
>> AT513   21   29
>> TT48   12   16
>>
>> Is there a fast way to do this?
>>
>> Thanks in advance!
>>
>> Jacy Crosby
>> jacy.r.cro...@uth.tmc.edu
>>
> 

-- 
View this message in context: 
http://www.nabble.com/Looking-for-a-quick-way-to-combine-rows-in-a-matrix-tp23491348p23525634.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] Download daily weather data

2009-02-27 Thread Chris Stubben


tlevine wrote:
> 
> The NOAA has very promising tabular forecasts
> (http://forecast.weather.gov/MapClick.php?CityName=Ithaca&state=NY&site=BGM&textField1=42.4422&textField2=-76.5002&e=0&FcstType=digital),
> but I can't figure out how to import them.
> 

Sometimes you can just use gsub to get html into R.   Try this.  

## read html 
wf<-readLines("http://forecast.weather.gov/MapClick.php?CityName=Ithaca&state=NY&site=BGM&textField1=42.4422&textField2=-76.5002&e=0&FcstType=digital";,
warn=FALSE) 

## you do need to find the rows with data (rows 18 and 19), grep may help
here
## mark end of lines
x<-gsub("", "\n", wf[18:19])
## remove white space...
x<-gsub(" ", "_", x)
## remove degree symbol or add comment.char="" to read.table below
x<-gsub("°", "", x)
## remove html tags
x<-gsub("<[^>]*>", " ", x)

y<-read.table(con<-textConnection(x), fill=TRUE)
close(con)

#now just reshape this mess 
# join rows 1,16 and 17-32 (skip 1st two rows with day and hour)

z<-data.frame(rbind( t(y[3:16,-1]), t(y[19:32,-1]) ), row.names=NULL,
stringsAsFactors=FALSE)
names(z)<-substr(y[3:16,1], 1,4)

head(z)
  Temp Dewp Wind Wind.1 Wind.2 Gust Sky_ Pcpn Rel. Thun Rain Snow Free Slee
1   49   43   43 18  S99   40   80   --  Chc   --   --   --
2   50   43   44 18SSW99   40   77   --  Chc   --   --   --
3   50   43   44 17 SW99   40   77   --  Chc   --   --   --
4   48   43   42 16 SW99   40   83   --  Chc   --   --   --
5   45   43   38 15WSW   100   40   93   --  Chc   --   --   --
6   42   42   35 14  W   100   35  100   --  Chc   --   --   --

# plot
plot(z["Temp"])


Chris Stubben
-- 
View this message in context: 
http://www.nabble.com/Download-daily-weather-data-tp22233373p22252517.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] how to save an updated dataset

2008-06-11 Thread Chris Stubben

I wrote a package which includes a number of genome sequencing project
statistics on the web like http://www.ncbi.nlm.nih.gov/genomes/lproks.cgi. I
included some generic functions to summarize, plot, and update the tables
with the most recent version

data(lproks)   
update(lproks)   
 [1] "lproks successfully updated, 7 new genomes added"

I usually save the dataset back to my package data directory...

save(lproks, file="/path/to/genomes/data/lproks.rda")

but I may one day put the package on CRAN (or BioConductor), so now I need
to know where the package/data directory is located, if the user has
permission to save to that directory, and probably some other complications
I'm missing.  Any suggestions?  

Thanks,
Chris Stubben

-- 
View this message in context: 
http://www.nabble.com/how-to-save-an-updated-dataset-tp17782584p17782584.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] SQL INSERT using RMySQL

2008-04-11 Thread Chris Stubben

Greg,

If you have a MySQL table with an auto_increment field, you could just
insert a NULL value into that column and the database will increment the key
(it may not work in SQL STRICT mode, I'm not sure).  I don't think there's
any way to specify which columns you want to load data into using
dbWriteTable yet, but that would be a nice feature since LOAD data now
allows that (and SET syntax and other options).  

Try this code below - I used cbind(NA, x) to insert a null into the first
column.

Chris

> dbSendQuery(con, "create table tmp (id int not null auto_increment primary
> key, a char(1), b char(1))")
 
> x<-data.frame( a=letters[1:3], b=letters[4:6])
> x
  a b
1 a d
2 b e
3 c f
> dbWriteTable(con, "tmp", cbind(NA, x), row.name=FALSE, append=TRUE)
[1] TRUE
> dbWriteTable(con, "tmp", cbind(NA, x), row.name=FALSE, append=TRUE)
[1] TRUE
> dbReadTable(con, "tmp")
  id a b
1  1 a d
2  2 b e
3  3 c f
4  4 a d
5  5 b e
6  6 c f





Gregory. R. Warnes wrote:
> 
> Hi All,
> 
> I've finally gotten around to database access using R.  I'm happily  
> extracting rows from a MySQL database using RMySQL, but am having  
> problems appending rows to an existing table.
> 
> What I *want* to do is to append rows to the table, allowing the  
> database to automatically generate primary key values.  I've only  
> managed to add rows by using
> 
>   dbWriteTable( con, "past_purchases", newRecords, overwrite=FALSE,  
> append=TRUE, ...)
> 
> And this only appears to properly append rows (as opposed to  
> overwriting them) IFF
> 1) the row names for newRecords are new unique primary key values,
> 2) the argument row.names is TRUE.
> 
> If row.names is FALSE, the records will not be appended, even if  
> newRecords contains a column (named 'id') of unique values that  
> corresponding to the primary key (named 'id').
> 
> It appears that in this case, the row names on the data frame are  
> still being used for the primary key, and since overwrite is FALSE,  
> the new records are being silently dropped.
> 
> 
> I did manage to get things working by doing the following:
> 
> ## get the last used id value (primary key)
> maxId <- dbGetQuery(con, "SELECT MAX(id) FROM past_purchases")[1,1]
> maxId
> if(is.na(maxId)) maxId <- -1
> 
> ## add the new unique primary keys as row names
> rownames(fulldata) <- maxId + 1:nrow(fulldata)
> 
> ## now write out the data
> dbWriteTable(con, "past_purchases", value=fulldata, overwrite=FALSE,  
> append=TRUE, row.names=TRUE)
> 
> 
> Is there a better way to accomplish this task?  (Session info is below)
> 
> Thanks!,
> 
> -Greg
> 
> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/SQL-INSERT-using-RMySQL-tp16640280p16644954.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] projection.matrix() {popbio} for more than one matrix

2008-02-24 Thread Chris Stubben

Michelle,

I would probably run a loop as well and save all matrices to a single list
(see hudsonia and calathea on working with lists of matrices).

First, run the example(test.census) to get the stage-fate data frame "trans"
and then run this code to save the matrices into a list "all".

years<-2001:2002   #or unique(trans$year)
all<-vector("list", length(years) )
names(all)<-years

## loop through years
for (s.year in years)
{
   ## Add individual fertilities 
   test.trans <- subset(trans, year==s.year)
   seedlings<-nrow(subset(test.census, year==s.year+1 & stage=="seedling"))
   test.trans$seedling<-test.trans$fruits/sum(test.trans$fruits) * seedlings
## save test.trans to another list if needed for bootstrapping,
calculating pooled matrix etc.
## test.trans  
 
   ## either one will work to get list index 
   all[[as.character(s.year)]] <-projection.matrix(test.trans)
   ##all[[ s.year- (years[1]-1) ]]  <-projection.matrix(test.trans)

}

all


$`2001`
  
   seedling vegetative reproductive
  seedling  0.00.0   1.6667
  vegetative0.50.5   0.
  reproductive  0.00.5   0.6667

$`2002`
  
   seedling vegetative reproductive
  seedling  0.00.0   0.6667
  vegetative0.20.5   0.
  reproductive  0.00.0   0.



There's also a loop in the demo(fillmore) example that uses the aq.census
data to create matrices for seven years.


Chris





Michelle DePrenger-Levin wrote:
> 
> Hello,
> 
>  
> 
> I am trying to use the projection.matrix( ) function and am following the
> example given. I have my data formatted very similar to the test.census
> example. 
> 
>  
> 
>> str(AsMi05mat)
> 
> `data.frame':   1854 obs. of  6 variables:
> 
>  $ Tag  : num  501 502 503 504 505 506 507 508 509 510 ...
> 
>  $ Year : int  1995 1995 1995 1995 1995 1995 1995 1995 1995 1995 ...
> 
>  $ Length   : num  34 37 11 24 7 44 4 7 12 20 ...
> 
>  $ Flowering: int  1 1 0 1 0 1 0 0 0 1 ...
> 
>  $ Fruits   : int  22 22 0 89 0 15 0 0 0 0 ...
> 
>  $ Stage: Factor w/ 5 levels "","Dead","Dormant",..: 4 4 5 4 5 4 5 5 5
> 4
> ...
> 
>  
> 
> The example data includes three years but only shows how to create a
> matrix
> from year 1 to 2. I have 13 years of data and would like to automate
> creating
> all matrices for each pair of years. 
> 
>  
> 
> I tried a for( ) loop but don't know how to assign new names for each
> matrix
> so end up with only the final comparison (2005 to 2006). I assume an
> apply( )
> function is the way to go but can't see how to do it. 
> 
>  
> 
> Thanks for any help!
> 
>  
> 
> Michelle
> 
>  
> 
> 
>   [[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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/projection.matrix%28%29-%7Bpopbio%7D-for-more-than-one-matrix-tp15646115p15673408.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] Population model question

2008-01-15 Thread Chris Stubben


-Halcyon- wrote:
> 
> My question is : How do I arbitrarily assign the animals being "FALSE"
> ***according to the uniformal chance generation*** to the array
> Dead.Animals, there being "TRUE"?
> 

Maybe I'm missing something, but aren't dead animals just the ones that
aren't alive?  Why do you even need the Dead.Animals array if you filled
Live.Animals correctly?

## 4 animals in cohort and 3 years 
Live.Animals<-matrix(c(T,T,F, T,T,T, T,T,F, F,F,F), nrow=4, byrow=TRUE)
 Live.Animals
  [,1]  [,2]  [,3]
[1,]  TRUE  TRUE FALSE
[2,]  TRUE  TRUE  TRUE
[3,]  TRUE  TRUE FALSE
[4,] FALSE FALSE FALSE

# Dead animals are NOT Live.Animals?

!Live.Animals
  [,1]  [,2]  [,3]
[1,] FALSE FALSE  TRUE
[2,] FALSE FALSE FALSE
[3,] FALSE FALSE  TRUE
[4,]  TRUE  TRUE  TRUE

# or maybe use rowSums to find dead animals 
dead<- which(rowSums(Live.Animals) < 3)
[1] 1 3 4

Chris
-- 
View this message in context: 
http://www.nabble.com/Population-model-question-tp14833982p14848821.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] stochastic growth rate (package popbio)

2008-01-14 Thread Chris Stubben


privalan wrote:
> 
> However, I cannot figure out the way to compute the asymptotic stochastic
> population growth rate using “popbio”. I would like to perform a
> stochastic model in which my demographic rates are sampled from beta
> distribution with known mean and SD. For example, rather than 0.50 (4th
> term of my matrix), I want to sample this parameter from a beta function
> at each time-step {e.g., betaval(0.5, 0.05)}. 
> 

I guess you can't do this right now, but it shouldn't be too difficult to
add an option to accept an expression as input to stoch.growth.rate or
stoch.projection in the next version of popbio.

Maybe this will get you started.

A<-expression(c(0.70, 0.70,0.35, betaval(0.5, 0.05)))

## stoch.growth.rate(A)  ## this doesn't work yet, so try...

matrix(eval(A), nrow=2,byrow=TRUE)
 
## code from simulation section in stoch.growth.rate 
##  just replace A with matrix(eval(A), nrow=2,byrow=TRUE) 
maxt<-1
r<-numeric(maxt)
n0<-c(1,1)

 for (t in 1:maxt)
   {
  n0<- matrix(eval(A), nrow=2,byrow=TRUE)  %*% n0
  N<-sum(n0)
  r[t]<-log(N)
  n0<-n0/N
   }
   loglsim<-mean(r)
   dse<-1.96*sqrt(var(r)/maxt)
   CI<-c(loglsim-dse, loglsim+dse)

 loglsim
#[1] 0.09982282
 CI
#[1] 0.0994700 0.1001757


# check
 eigen.analysis( matrix(c(0.70, 0.70,0.35,0.50), nrow=2,byrow=TRUE))$lambda1
 exp(loglsim)


Thanks,

Chris





-- 
View this message in context: 
http://www.nabble.com/stochastic-growth-rate-%28package-biopop%29-tp14800951p14806622.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] Diagonal matrix with off diagonal elements

2007-12-21 Thread Chris Stubben

Also try the odiag function in the demogR package

odiag( 1:5, -1)
 [,1] [,2] [,3] [,4] [,5] [,6]
[1,]000000
[2,]100000
[3,]020000
[4,]003000
[5,]000400
[6,]000050

Chris





Jonas Malmros wrote:
> 
> Hi, everyone
> 
> I wonder if there is a function in R with which I can create a square
> matrix with elements off main diagonal (for example one diagonal below
> the main diagonal).
> 
> Thanks in advance!
> 
> -- 
> Jonas Malmros
> Stockholm University
> Stockholm, Sweden
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Diagonal-matrix-with-off-diagonal-elements-tp14462310p14463831.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] Inserting a subsequence between values of a vector

2007-12-04 Thread Chris Stubben

You could use a combination of rle, cumsum and append.

> x <- c(1,1,1,2,2,3,3,3,3,3,4) 
> y<-rle(x)$lengths
> y
[1] 3 2 5 1
> z<-cumsum(y)[y>1]
> z
[1]  3  5 10
>
> for(i in rev(z)) x <- append(x, c(0,0,0), after = i)
> x
 [1] 1 1 1 0 0 0 2 2 0 0 0 3 3 3 3 3 0 0 0 4


Chris


Serguei Kaniovski-3 wrote:
> 
> 
> Hallo,
> 
> suppose I have a vector:
> 
> x <- c(1,1,1,2,2,3,3,3,3,3,4)
> 
> How can I generate a vector/sequence in which a fixed number of zeroes
> (say
> 3) is inserted between the consecutive values, so I get
> 
> 1,1,1,0,0,0,2,2,0,0,0,3,3,3,3,3,0,0,0,4
> 
> thanks a lot,
> Serguei
>   [[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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Inserting-a-subsequence-between-values-of-a-vector-tf4943930.html#a14154029
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] Generating these matrices going backwards

2007-11-15 Thread Chris Stubben

Sorry, I wasn't sure what you meant.  This way will return more than one
answer, right?

N<-c(1,2,1,3)
R<-c(1.75,3.5,1.75,1.3125)

## get all 126 combinations of five 0's and four 1's for matrix
cbn<-as.matrix(expand.grid( rep( list(0:1), 9)))
cbn<- cbn[rowSums(cbn)==4,] 

ans<-list()
ctr<-0
## loop through each combination 
for (i in 1:126){
x<-cbn[i,]
## replace 1's with N
x[which(x==1)]<-N
## create matrix
dim(x)<-c(3,3)
## calculate y and new R
y <-sum(x) * x / (rowSums(x)%o%colSums(x)) 
R1<-y[y[1:3,]>0] 
# check if equal to original R
if(identical(R, R1)) ans[[ctr<-ctr+1]]<-x
}
ans


[[1]]
 [,1] [,2] [,3]
[1,]001
[2,]103
[3,]020

[[2]]
 [,1] [,2] [,3]
[1,]001
[2,]020
[3,]103

[[3]]
 [,1] [,2] [,3]
[1,]020
[2,]001
[3,]103


Chris


francogrex wrote:
> 
> Hi, thanks but the way you are doing it is to assign the values of N in
> the x matrix, knowing from the example I have given where they are
> supposed to be. While the assumption is, you ONLY have values of N and R
> and do NOT know where they would be placed in the x and y matrix a-priori,
> but their position has to be derived from only the (N and R) dataframe you
> have.
> 
> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Generating-these-matrices-going-backwards-tf4807447.html#a13782676
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] Generating these matrices going backwards

2007-11-14 Thread Chris Stubben

Maybe something like this?

N<-c(1,2,1,3)
## create empty matrix
x<-diag(0,3)
## fill off diagonal
 x[row(x)==col(x)+1]<-N[1:2]
# fill 3rd column
 x[1:2,3]<-N[3:4]

Or create a function and return both x and y

mat3<-function(N)
{
 x<-diag(0,3)
# fill each element separately
x[2,1]<-N[1]
x[3,2]<-N[2]
x[1,3]<-N[3]
x[2,3]<-N[4]
dimnames(x)<-list(c("D1", "D2", "D3"), c("E1", "E2", "E3"))
y =sum(x) * x / (rowSums(x)%o%colSums(x)) 
list(x=x,y=y)
}
mat3(N)

$x
   E1 E2 E3
D1  0  0  1
D2  1  0  3
D3  0  2  0

$y
 E1  E2 E3
D1 0.00 0.0 1.7500
D2 1.75 0.0 1.3125
D3 0.00 3.5 0.




francogrex wrote:
> 
> I have generated the following:
> 
> x= 
> E1  E2  E3 
> D1  0   0   1 
> D2  1   0   3 
> D3  0   2   0 
> 
> y= 
> E1  E2  E3 
> D1  0   0   1.75 
> D2  1.750   1.3125 
> D3  0   3.5 0 
> 
> Where x and y are linked by: 
> y =sum(x) * x / (rowSums(x)%o%colSums(x)) 
> 
> N=x[x[1:3,]>0] 
> R=y[y[1:3,]>0] 
> 
> Now suppose I ONLY have N and R linked in this way below where each N 
> corresponds to an R 
> 
> N   R 
> 1   1.7500 
> 2   3.5000 
> 1   1.7500 
> 3   1.3125 
> 
> Is there a way to generate matrix "x" and matrix "y" having only the N 
> and the R as above? 
> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Generating-these-matrices-going-backwards-tf4807447.html#a13756103
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] Elasticity in Leslie Matrix

2007-10-29 Thread Chris Stubben

Sorry, I didn't think carefully about your question on my first reply, so
please ignore it.  

Chapter 9 in Morris and Doak (2002) cover vital rate elasticities and
sensitivities (see Box 9.1 for Matlab code
http://www.sinauer.com/PVA/vitalsens.m).   Simon, thanks for posting your
code, I will add a function called vitalsens to the next release of the
popbio package and output both sensitivities and elasiticities.

Here's the example for the emperor goose vital rates in chapter 9.

vl<-list(  Ss0=0.1357,  Ss1=0.8926,  Sf2=0.6388,  Sf3= 0.8943)

el<-expression(
0,  0,  Sf2*Ss1,Sf3*Ss1,
Ss0,0,  0,  0,
0,  Ss1,0,  0,
0,  0,  Ss1,Ss1)

## and the matching elasticities in table 9.1

elas.var(el, vl)

$Ss0
[1] 0.0797

$Ss1
[1] 0.92

$Sf2
[1] 0.00569

$Sf3
[1] 0.074


Also, the expression el should be formatted so this command returns the
projection matrix.

matrix(sapply(el, eval, vl), nrow=4, byrow=TRUE)


Thanks,

Chris



Simon Blomberg-4 wrote:
> 
> After a short exchange with the original questioner, I wrote the following
> function to calculate the elasticities of lower level variables in
> population transition matrices (Leslie matrices etc.) Perhaps it will be
> of use to others. There is no error-checking, so use with care. Users
> should consult Caswell (2001) for reference.
> 
> Cheers,
> 
> Simon.
> 
> # example values to construct leslie matrix
> 
> vl <- list(f1=1, f2=4, s0=.6, s=.4, v=.9, sigma=.5)
> 
> # Expressions for each matrix element
> F1 <- expression(sigma*s0*f1)
> F2 <- expression(sigma*s0*f2)
> S <- expression(s)
> V <- expression(v)
> el <- c(F1, F2, S, V)
> 
> elas.var <- function (elements, varlist) {
>   # elements should be a vector of expressions corresponding to the
> elements
>   # of the leslie matrix, in terms of the variables in varlist
> 
>   require(demogR)
>   res <- vector("list", length(varlist))
>   deriv.funcs <- sapply(elements, deriv, namevec=names(varlist),
> function.arg=TRUE)
>   devs <- lapply(deriv.funcs, function (x) do.call(x, varlist))
>   leslie.mat <-  matrix(as.numeric(devs), nrow=sqrt(length(elements)),
> byrow=TRUE)
>   eig <- eigen.analysis(leslie.mat)
> 
>   for (i in 1:length(varlist)) {
> derivs <- matrix( as.numeric(lapply(devs, function (x)
> [EMAIL PROTECTED])), nrow=sqrt(length(elements)), byrow=TRUE)
> res[[i]] <- varlist[[i]]/eig$lambda1*sum(derivs*eig$sensitivities)
> names(res)[i] <- names(varlist)[i]
>   }
>   res
> }
> 
> # example output
> 
>> elas.var(el, vl)
> $f1
> [1] 0.06671376
> 
> $f2
> [1] 0.2346064
> 
> $s0
> [1] 0.3013201
> 
> $s
> [1] 0.2346064
> 
> $v
> [1] 0.4640735
> 
> $sigma
> [1] 0.3013201
> 
> 
> Simon Blomberg, BSc (Hons), PhD, MAppStat. 
> Lecturer and Consultant Statistician 
> Faculty of Biological and Chemical Sciences 
> The University of Queensland 
> St. Lucia Queensland 4072 
> Australia 
> T: +61 7 3365 2506 
> email: S.Blomberg1_at_uq.edu.au
> 
> Policies:
> 1.  I will NOT analyse your data for you.
> 2.  Your deadline is your problem.
> 
> The combination of some data and an aching desire for 
> an answer does not ensure that a reasonable answer can 
> be extracted from a given body of data. - John Tukey.
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
>   [[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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Elasticity-in-Leslie-Matrix-tf4670533.html#a13470907
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] Elasticity in Leslie Matrix

2007-10-22 Thread Chris Stubben
> I would like to calculate elasticities and sensitivities of each parameters
> involved in the following transition matrix:
> 
> A <- matrix(c(
> sigma*s0*f1,  sigma*s0*f2,
>s, v
> ), nrow=2, byrow=TRUE,dimnames=list(stage,stage))
> 
> The command "eigen.analysis" avaliable in package "popbio" provides
> sensibility matrix and elasticity matrix (same dimension than A). I would
> like to know if there is a way to calculate separetely the elasticity of
> sigma, s0, f1, f2, s and v ?
> 

If you want the eigenvalue second derivatives, check the secder and fullsecder
functions in the demogR package (also elassens for partial derivatives of the
eigenvalue elasticities).  A description of that package is also found in the
recent special issue of Ecology in R in Journal of Statistical Software.  

Chris

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

2007-10-12 Thread Chris Stubben
> On 10/12/07, Ben Bolker  ufl.edu> wrote:
> >
> >   Trying to find a quick/slick/easily interpretable way to
> > collapse a data set.
> >

Another alternative for SQL fans is the sqldf package.  I used the MySQL driver
here since SQLite does not support standard deviation. 


sqldf("select BROOD, avg(TICKS) as 'TICKS.mean', stddev_samp(TICKS) as
'TICKS.sd', HEIGHT, YEAR, LOCATION
 from h
group by BROOD, HEIGHT, YEAR, LOCATION",  drv="MySQL")

  BROOD TICKS.mean TICKS.sd HEIGHT YEAR LOCATION
1   501  0 0.00465   95   32
2   502  0   NA472   95   36
3   503  1 1.73475   95   37



Chris

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] RMySQL NA/NULL value storage error

2007-09-28 Thread Chris Stubben


Adam Wilson-4 wrote:
> 
> I am running R 2.5.1, RMySQL 0.6 , and DBI 0.2-3 on Windows XP
> 
> Like others, I am having trouble with NA/Null value conversions between R
> and a MySQL database via DBI, but I could not find my exact problem in the
> archives.  Most of the time NA values in R get transferred correctly to
> the
> database and back again, unless (apparently) if they are in the last
> column
> of the table being saved.
> 

I have the R 2.5.0, RMySQL 0.6 , and DBI 0.2-3 on Linux and don't get this
problem.

What do you get running 

>  str(x)

>  dbGetQuery(test,"describe x")
   Field   Type Null Key Default Extra
 1X1 double  YES
 2X2 double  YES
 3X3 double  YES
 4X4 double  YES

>  dbGetInfo(test)$serverVersion
 [1] "5.0.41"


Chris


-- 
View this message in context: 
http://www.nabble.com/RMySQL-NA-NULL-value-storage-error-tf4531309.html#a12943167
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] 3d barplot in rgl

2007-09-26 Thread Chris Stubben
>Transition matrices are Markov transition matrices among different
> life stages of organisms -- in the simplest case (Leslie matrices,


Ben,

Thanks for your clear explanation and plot examples.  I like the dotplots alot
and added a few modifications below. Since I often compare rows and columns as
well as groups of elements representing growth (lower left trianlge), stasis
(diagonal), fertility etc., I like to preserve the matrix structure in the plot
if possible.


> 
> library(popbio)
> example(teasel)   ## shows 2 plots -- heatmap and stair-step
> A <- tea$sensitivities ##
> z <- as.data.frame.table(A)


#switch labels here- stage at time t in columns and fate at time t+1 in rows
names(z) <- c("To","From","Sensitivity")

# reverse order on rows to match matrix
z$To <- ordered(z$To, levels = rev(levels(z$To)))

dotplot(To~Sensitivity|From, data=z,
scales=list(alternating=c(1,0), tck=c(1,0)), layout=c(6,1), aspect=2)

# or log-scale 
dotplot(To~Sensitivity|From,data=z,scales=list(x=list(log=TRUE)), as.table=TRUE)


Thanks again,

Chris

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

2007-09-25 Thread Chris Stubben
Duncan Murdoch  stats.uwo.ca> writes:


> That demo gives you the basics of the code, so it shouldn't be too hard 
> to put your own together:  just strip out the counting part.
> 


Thanks, I did download the source to check the hist3d demo, but honestly it
didn't look very easy to simplify.  I switched to R from matlab, and 3d barplots
are the only thing I still use my student edition of matlab for anymore.  

I  would like to see a barplot3d added to your rgl package in the future, but
was hoping someone else had tried already.


Chris

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

2007-09-25 Thread Chris Stubben
hadley wickham  gmail.com> writes:

> 
> Why do you want a 3d barchart?  They are generally a bad way to
> present information as tall bars can obscure short bars, and it is
> hard to accurately read off the height of a bar.  While adding
> rotation can reduce some of these problems, why not create a graphic
> that your viewers can take in with a glance?
> 

3d barplots are a common way to display sensitvity/elasticity matrices in
stage-structured demography.  Here's a few other options from Caswell's Matrix
population models book (2001) - I definitely  prefer 3d barcharts to these
alternatives.


heatmap(log10(A[3:1,]), Rowv = NA,  Colv = NA, scale="none")

plot(log10(c(A)), type="s")


Thanks,

Chris

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

2007-09-17 Thread Chris Stubben

You could use a loop and fill small polygons with colors.

x<-seq(-3, 3, .01)
y<-eval(expression(x^3-3*x))
plot(x,y, type="n", las=1)   
n<-length(x)

# vertical bars
for(i in 1:n)
{
polygon(c(x[i], x[i], x[i+1], x[i+1]), c(min(y), y[i], y[i+1], min(y)),
border=0, col = rainbow(n)[i]) 
}

## or split each vertical bars into 100 smaller blocks
plot(x,y, type="n", las=1)
for(i in 1:n)
{
  y1<-seq(y[i], min(y), length.out=101)
  y2<-seq(y[i+1], min(y), length.out=101)
  for(j in 1:100)
  { 
polygon(c(x[i], x[i], x[i+1], x[i+1]), c(y1[j+1], y1[j], y2[j],
y2[j+1]), border=0, col = rev(heat.colors(100))[j]) 
  }
}

Chris Stubben




Van Dongen Stefan wrote:
> 
> Hi All,
>  
> I would like to fill the area under a curve with a gradient of colors. Are
> there any packages or trick I could use
>  
> Thanks
>  
> Stefan
>  
>  
> Stefan Van Dongen
> Antwerp
> 
> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/graphs-with-gradients-of-colors-tf4466780.html#a12743725
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.