[R] apply() and sample()

Patrick Presi ppresi at gmail.com
Mon Apr 16 14:51:22 CEST 2012


Hi All,

I am developing an epidemiological model at animal level.

> head(all.herds)
  ID herd age status status2 dry dry.per pasture id
1  1    1   3      s       s lac       4       1  1
2  2    1   8      s       s lac       4       1  2
3  3    1  30      s       s lac       4       1  3
4  4    1  10      s       s lac       4       1  4
5  5    1  20      s       s lac       4       1  5
6  6    1  57    gtb     gtb lac       4       1  6

I am trying to find a way to speed up this loop (9 sec so far per
iteration for 20000 animals)

all.herds[all.herds[,p]=="s",f]	<-
	apply(all.herds[all.herds[,p]=="s",],1,function(x)
sample(c("s","eco","gto","gtb","gtc","sub","sdy"),1,prob=p.inf[ifelse((x[6]=="lac"
| (x[6]=="dry" &
x[7]==4)),as.numeric(x[8]),max.pasture+as.numeric(x[8])),2:8]))

the idea is to define the status of each animal based on the herd it
belongs to and the lactation phase. For each herd there are different
probabilities associated to each status.

legend to the code:
p and f refer to present situation or future (at each iteration it
changes using this:

	p     <- ((leng+1)%%2)+4   		# situation at t
	f     <- (leng%%2)+4       		# situation at t+1

where leng is my counter.

c("s","eco","gto","gtb","gtc","sub","sdy") are the different status an
animal can get based on probability thar are valid at a herd level and
that are different for cow in lactating or dry period (+ dry.per=4 =>
behave as lact)

> head(p.inf)
  status         s   eco   gto gtb gtc        sub        sdy
1    lac 0.9980000 0.001 0.001   0   0 0.00000000 0.00000000
2    lac 0.9849608 0.001 0.001   0   0 0.00000000 0.01303924
3    lac 0.9980000 0.001 0.001   0   0 0.00000000 0.00000000
4    lac 0.9980000 0.001 0.001   0   0 0.00000000 0.00000000
5    lac 0.9735968 0.001 0.001   0   0 0.02440319 0.00000000
6    lac 0.9980000 0.001 0.001   0   0 0.00000000 0.00000000

> dim(p.inf)
[1] 2000    8  => two lines per herd (1 for lact one for dry)

max.pasture = 1000

(herd and pasture are exchangeable, I use this in order to simulate
animal movement to summer pasture)


Thank you in advance for your help

Patrick



More information about the R-help mailing list