[R] Data manipulation in columns (with apply?)
R Users, I have written a small simulation model in R which outputs a datafile consisting of ending population sizes for each simulation run (year). The data (see short data example below) is labeled by NUM (simulation run), sim (year) and N (yearly count). After searching the help files and coming up empty (probably because I used the wrong terms) I am appealing for some help for working with the output dataset. What I want to do is for each of the i simulation runs (NUM) I want to 1) take N(t+1)/N(t)=lambda(t) for each year (where in the below example t=1,...,10--total years of the simulation) 2) Sum lambda(t) and divide by t (e.g., output both the mean/se of lambda for each simulation run) 3) Take the mean of the mean(lambda's) (and associated stddev, min, max) over all NUM I think I have to write a function for use within an apply statement, but I am not quite there yet on the learning curve so most of my recent attempts in R have been useful learning experiences of what not to do... Any suggestions/direction is greatly appreciated. Bret Collier TX AM NUM sim N 1 1 466 1 2 450 1 3 473 1 4 531 1 5 515 1 6 502 1 7 471 1 8 460 1 9 458 1 10 434 2 1 289 2 2 356 2 3 387 2 4 440 2 5 457 2 6 466 2 7 467 2 8 449 2 9 387 2 10 394 3 1 367 3 2 400 3 3 476 3 4 508 3 5 478 3 6 501 3 7 513 3 8 505 3 9 492 3 10 465 platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 2 minor 3.0 year 2006 month 04 day24 svn rev37909 language R version.string Version 2.3.0 (2006-04-24) (yeah, I need to update) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Error with 'var' in JGR
R-Users, I had a strange occurrence this AM I could use some help with. 'var' seems to have ceased to work for me when I used JGR 1.3 interface. So, I ran a quick test (following pg 57 of Dalgaard) in both the basic RGui and JGR interface (see below). I have re-installed R 2.2.1 and JGR 1.3 but had no luck so far figuring this out? Anyone have a suggestion? Or did I just make a mistake somewhere (which is more than likely). Thanks, Bret TX AM Output from RGui: x-rnorm(50) x [1] -1.016219752 0.940801964 0.168647573 0.328696253 -0.600957761 [6] -0.420394338 -1.140567123 -0.755041336 0.907831188 0.198247025 [11] -0.065116449 0.732841048 -0.400194138 -1.490167845 -1.488611328 [16] 0.311665408 -1.534002497 -1.094357124 -1.717282302 0.217445299 [21] -1.878605987 -0.214200092 -0.096099830 -0.357325121 -0.118191356 [26] -0.214543361 -0.399768989 0.235104678 -0.363194200 -0.275338828 [31] -0.336677033 -0.009178995 1.294139880 1.442454681 1.192023689 [36] -0.521847342 -1.070860356 -0.756584142 0.747503883 0.939960584 [41] -1.444890590 1.218499328 0.343071542 1.303250666 -0.603718755 [46] 0.979527031 -3.709878278 -0.051090815 -0.452191654 -0.206564681 mean(x) [1] -0.226039 sd(x) [1] 0.9917835 var(x) [1] 0.9836345 Output from JGR: x-rnorm(50) x [1] 0.79961648 -0.72223557 -0.10589876 [4] -0.25367834 -0.53518039 0.18296671 [7] 0.33981503 0.29208966 1.15002574 [10] 0.84356083 -0.87841050 0.31345819 [13] 0.61348608 1.13257372 -0.14366714 [16] -1.28563266 1.39519683 -0.85124979 [19] -1.65043342 0.93956978 1.36692691 [22] 0.81240164 -0.44326482 0.20741846 [25] 0.38536713 -0.57994109 1.10457148 [28] -0.99307813 -2.31580515 -1.61582072 [31] 1.63273805 -1.49289425 -1.33675463 [34] -1.17789478 -0.42209334 -0.21208394 [37] -2.21572873 0.35724725 -0.24758106 [40] 0.60470892 1.71646244 1.02161560 [43] 2.93773901 -0.92356017 -0.38588476 [46] -0.08622336 -0.09622387 -0.91419002 [49] 0.93453732 -0.25280526 mean(x) [1] -0.02108243 sd(x) [1] 1.077132 var(x) Error in get(x, envir, mode, inherits) : variable var was not found version platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Estimating Daily Survival
R Users, I was wondering if someone might point me in the right direction. I am using a Cox model (survival package) to evaluate survival of pen-reared birds (time to event data collected daily) and I have been trying to determine how I can estimate a 'daily' survival rate and std error from the results of a Cox model? Using survfit (see input data below) I computed the predicted survivor function for a Cox model, but I have been unable to figure out how to estimate a 'daily' survival rate for an average individual? Should I be looking at survexp? Any suggestions would be appreciated. Thanks, Bret brettest1-coxph(Surv(Entry, Exit, Fate)~Sex , data=mydat) brettest1 coef exp(coef) se(coef) zp Sex -0.503 0.6050.524 -0.96 0.34 Likelihood ratio test=0.92 on 1 df, p=0.337 n= 19 survfit(brettest1)$surv [1] 0.9502716 0.9010356 0.8508737 0.7996929 0.7473818 0.6957615 0.5872297 0.5324263 0.4756545 [10] 0.4201535 0.3547775 0.2916481 0.2189030 0.1374015 mydat ID Year Dayrelease Agerelease Survivorship Entry Exit Fate Sex 1 16240 1996205 95 164 205 3691 0 2 16319 1996205 88 140 205 3451 0 3 16378 1996248108 100 248 3481 0 4 20383 1996241 98 204 241 4451 0 5 16324 1996219 90 227 219 4461 0 6 16327 1996219 90 497 219 7161 0 7 20373 1996241114 413 241 6541 0 8 20374 1996241111 211 241 4521 0 9 16241 1996205 95 234 205 4391 1 10 16321 1996219 90 118 219 3371 1 11 16323 1996219 90 180 219 3991 1 12 20375 1996241103 268 241 5091 1 13 20384 1996241 98 299 241 5401 1 14 20390 1996241 93 204 241 4451 1 15 20393 1996241 88 208 241 4491 1 16 16313 1996248122 512 248 7600 1 17 20378 1996241103 236 241 4770 1 18 20381 1996241101 329 241 5700 1 19 16328 1996219 90 224 219 4430 1 platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Survival Plots by Strata
All, I am struggling to create a survival plot using LTRC data for each year of a 10 year period. I have a set of individuals (birds) where 'entry' is the day of the year (1-365) they are released (let out of pens) into the wild (2 year data snip below). 'Entry' (e.g., day of year the first bird is released for each year) is highly variable, ranging from 48 to 250. When I create a plot using the below code I would like to remove the 'solid' line which represent S_hat=1 out to the LC point for each year and instead have the survival curves formatted in more of a 'hanging' style, where the LC day (e.g., min(Entry for each year)) is where each curve begins at S_hat=1 for each year (and does not extend back to the y-axis)? I could not find anything on this in the archives or MASS or Survival Analysis using S? Anyone have a suggestion on where to look? TIA, Bret #Code snip for R email. apc.coxfit1-coxph(Surv(Entry, Exit, Fate)~Sex + Agerelease + Dayrelease + strata(Year), data=mydat) coxfit.apc-survfit(apc.coxfit1) coxfit.apc plot(survfit(apc.coxfit1), conf.int=F, log=T, lty=c(1:2), col=c(1:2), xlim=c(205, 800)) #not run--first entry for this example is day 205 for 1996, 259 for 1997 mydat ID YearDayrelease Agerelease SurvivorshipEntry ExitFateSex 16240 1996205 95 164 205 369 1 0 16319 1996205 88 140 205 345 1 0 16378 1996248 108 100 248 348 1 0 20383 1996241 98 204 241 445 1 0 16324 1996219 90 227 219 446 1 0 16327 1996219 90 497 219 716 1 0 20373 1996241 114 413 241 654 1 0 20374 1996241 111 211 241 452 1 0 16241 1996205 95 234 205 439 1 1 16321 1996219 90 118 219 337 1 1 16323 1996219 90 180 219 399 1 1 20375 1996241 103 268 241 509 1 1 20384 1996241 98 299 241 540 1 1 20390 1996241 93 204 241 445 1 1 20393 1996241 88 208 241 449 1 1 16313 1996248 122 512 248 760 0 1 20378 1996241 103 236 241 477 0 1 20381 1996241 101 329 241 570 0 1 16328 1996219 90 224 219 443 0 1 16827 1997259 127 52 259 311 1 0 16828 1997259 127 216 259 475 1 0 16831 1997303 171 19 303 322 1 0 20466 1997289 149 31 289 320 1 0 20469 1997289 149 199 289 488 1 0 20483 1997289 134 18 289 307 1 0 16807 1997259 137 223 259 482 1 0 16809 1997259 137 1 259 260 1 0 16819 1997259 131 237 259 496 1 0 16829 1997303 171 7 303 310 1 0 20440 1997289 161 7 289 296 1 0 20470 1997289 148 257 289 546 1 0 20478 1997289 143 12 289 301 1 0 16817 1997259 130 85 259 344 1 0 20454 1997289 154 4 289 293 1 1 20459 1997289 153 335 289 624 1 1 20460 1997289 153 118 289 407 1 1 20465 1997289 150 31 289 320 1 1 20473 1997289 147 65 289 354 1 1 20484 1997289 133 58 289 347 1 1 20489 1997289 130 56 289 345 1 1 16808 1997259 137 137 259 396 0 1 16810 1997303 181 64 303 367 0 1 16816 1997259 130 1 259 260 0 1 20471 1997289 147 334 289 623 0 1 16826 1997259 127 20 259 279 0 1 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Survival Plots by Strata
Thomas, Thank you for the response. I will post what I come up with. Bret Thomas Lumley [EMAIL PROTECTED] 03/08/06 9:16 AM I don't see any easy way to do this. I think you may have to do the plotting yourself, based on the code in plot.survfit. -thomas On Wed, 8 Mar 2006, Bret Collier wrote: All, I am struggling to create a survival plot using LTRC data for each year of a 10 year period. I have a set of individuals (birds) where 'entry' is the day of the year (1-365) they are released (let out of pens) into the wild (2 year data snip below). 'Entry' (e.g., day of year the first bird is released for each year) is highly variable, ranging from 48 to 250. When I create a plot using the below code I would like to remove the 'solid' line which represent S_hat=1 out to the LC point for each year and instead have the survival curves formatted in more of a 'hanging' style, where the LC day (e.g., min(Entry for each year)) is where each curve begins at S_hat=1 for each year (and does not extend back to the y-axis)? I could not find anything on this in the archives or MASS or Survival Analysis using S? Anyone have a suggestion on where to look? TIA, Bret #Code snip for R email. apc.coxfit1-coxph(Surv(Entry, Exit, Fate)~Sex + Agerelease + Dayrelease + strata(Year), data=mydat) coxfit.apc-survfit(apc.coxfit1) coxfit.apc plot(survfit(apc.coxfit1), conf.int=F, log=T, lty=c(1:2), col=c(1:2), xlim=c(205, 800)) #not run--first entry for this example is day 205 for 1996, 259 for 1997 mydat IDYearDayrelease AgereleaseSurvivorshipEntry ExitFateSex 16240 1996205 95 164 205 369 1 0 16319 1996205 88 140 205 345 1 0 16378 1996248 108 100 248 348 1 0 20383 1996241 98 204 241 445 1 0 16324 1996219 90 227 219 446 1 0 16327 1996219 90 497 219 716 1 0 20373 1996241 114 413 241 654 1 0 20374 1996241 111 211 241 452 1 0 16241 1996205 95 234 205 439 1 1 16321 1996219 90 118 219 337 1 1 16323 1996219 90 180 219 399 1 1 20375 1996241 103 268 241 509 1 1 20384 1996241 98 299 241 540 1 1 20390 1996241 93 204 241 445 1 1 20393 1996241 88 208 241 449 1 1 16313 1996248 122 512 248 760 0 1 20378 1996241 103 236 241 477 0 1 20381 1996241 101 329 241 570 0 1 16328 1996219 90 224 219 443 0 1 16827 1997259 127 52 259 311 1 0 16828 1997259 127 216 259 475 1 0 16831 1997303 171 19 303 322 1 0 20466 1997289 149 31 289 320 1 0 20469 1997289 149 199 289 488 1 0 20483 1997289 134 18 289 307 1 0 16807 1997259 137 223 259 482 1 0 16809 1997259 137 1 259 260 1 0 16819 1997259 131 237 259 496 1 0 16829 1997303 171 7 303 310 1 0 20440 1997289 161 7 289 296 1 0 20470 1997289 148 257 289 546 1 0 20478 1997289 143 12 289 301 1 0 16817 1997259 130 85 259 344 1 0 20454 1997289 154 4 289 293 1 1 20459 1997289 153 335 289 624 1 1 20460 1997289 153 118 289 407 1 1 20465 1997289 150 31 289 320 1 1 20473 1997289 147 65 289 354 1 1 20484 1997289 133 58 289 347 1 1 20489 1997289 130 56 289 345 1 1 16808 1997259 137 137 259 396 0 1 16810 1997303 181 64 303 367 0 1 16816 1997259 130 1 259 260 0 1 20471 1997289 147 334 289 623 0 1 16826 1997259 127 20 259 279 0 1 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle
[R] JRG Console Output
All, I had a question about the JGR console and whether or not I can manipulate the location where line wrapping occurs. I have searched 'JGR' in the R listserve archives and attempted to find console manipulation on the JGR website to no avail and could use some direction. TIA, Bret As an example, the below output wraps every 4th value, leaving about 2/3 of the console empty. rnorm(20, 50, 17) [1] 43.42240 39.94807 8.94276 15.19369 [5] 82.56500 17.96678 80.58936 48.61693 [9] 75.92249 71.86615 53.39025 24.08080 [13] 23.92690 35.96344 61.29200 63.37290 [17] 77.71882 56.54847 70.16172 51.61530 where the below increases the 'after decimal' range 1 unit and cuts back to 3 values per line. rnorm(1000, 50, 17) [1] 64.805808 54.770722 46.552925 [4] 34.371987 61.971183 37.327260 [7] 60.403990 47.783683 51.744272 [10] 32.671062 31.395127 69.085938 [13] 77.554024 48.579639 48.111326 [16] 44.994786 65.241722 65.852035 [19] 53.254482 43.217719 30.255150 Reading in some data, the 'mydat' line extends across the console, but the output is wrapped around again? mydat-data.frame(ID, Year = factor(Year), Dayrelease, Agerelease, Survivorship, Entry, Exit, Fate, Gender) mydat ID Year Dayrelease Agerelease 1 16240 1996205 95 2 16319 1996205 88 3 16378 1996248108 . . . 576 3094 2005251136 577 2667 2005264861 578 3035 2005264497 Survivorship Entry Exit Fate Gender 1164 205 3691 0 2140 205 3451 0 3100 248 3481 0 4204 241 4451 0 5227 219 4461 0 version platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] JRG Console Output
The obvious answer was options(width=). Bret Bret Collier [EMAIL PROTECTED] 2/13/2006 12:24:22 PM All, I had a question about the JGR console and whether or not I can manipulate the location where line wrapping occurs. I have searched 'JGR' in the R listserve archives and attempted to find console manipulation on the JGR website to no avail and could use some direction. TIA, Bret As an example, the below output wraps every 4th value, leaving about 2/3 of the console empty. rnorm(20, 50, 17) [1] 43.42240 39.94807 8.94276 15.19369 [5] 82.56500 17.96678 80.58936 48.61693 [9] 75.92249 71.86615 53.39025 24.08080 [13] 23.92690 35.96344 61.29200 63.37290 [17] 77.71882 56.54847 70.16172 51.61530 where the below increases the 'after decimal' range 1 unit and cuts back to 3 values per line. rnorm(1000, 50, 17) [1] 64.805808 54.770722 46.552925 [4] 34.371987 61.971183 37.327260 [7] 60.403990 47.783683 51.744272 [10] 32.671062 31.395127 69.085938 [13] 77.554024 48.579639 48.111326 [16] 44.994786 65.241722 65.852035 [19] 53.254482 43.217719 30.255150 Reading in some data, the 'mydat' line extends across the console, but the output is wrapped around again? mydat-data.frame(ID, Year = factor(Year), Dayrelease, Agerelease, Survivorship, Entry, Exit, Fate, Gender) mydat ID Year Dayrelease Agerelease 1 16240 1996205 95 2 16319 1996205 88 3 16378 1996248108 . . . 576 3094 2005251136 577 2667 2005264861 578 3035 2005264497 Survivorship Entry Exit Fate Gender 1164 205 3691 0 2140 205 3451 0 3100 248 3481 0 4204 241 4451 0 5227 219 4461 0 version platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Can't get sample function from An Introduction to R to work
David,
If below is exactly what you typed, check your code again, I think you
are missing a '}' after the last 2 parentheses.
HTH,
Bret
David Groos [EMAIL PROTECTED] 7/15/2005 3:40:01 PM
I'm trying to figure out R, a piece at a time, hours at a time... I
was trying to copy the sample function in, An Introduction to R (for
version 2.1.0) by W. N. Venables, D. M. Smith, page 42. Section 10.1
Simple examples provides a sample function which I tried to duplicate
(I'm using Mac OS X 10.3.9, and R for Mac OS X Aqua GUI v1.11). The
following is what I typed and the last line is R's response when I hit
the return key after the penultimate line. I've re-checked and
re-typed the code many times to no avail. I wasn't able to find this
issue using search options, either. Any help is GREATLY appreciated!
twosam-function(y1, y2) {
+ n1-length(y1);n2 -length(y2)
+ yb1-mean(y1); yb2-mean(y2)
+ s1-var(y1);s2-var(y2)
+ s-((n1-1)*s1 + (n2-1)*s2)/(n1+n2-2)
+ tst-(yb1-yb2)/sqrt(s*(1/n1+1/n2))
Error: syntax error
David
[[alternative text/enriched version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Overlying a Normal Dist in a Barplot
R-Users, Hopefully someone can shed some light on these questions as I had little luck searching the archives (although I probably missed something in my search due to the search phrase). I estimated multinomial probabilities for some count data (number successful offspring) ranging from 0 to 8 (9 possible response categories). I constructed a barplot (using barplot2) and I want to overlay a normal distribution on the figure (using rnorm (1000, mean, sd)). My intent is to show that using a mean(and associated sd) estimated from discrete count data may not be a valid representation of the distribution of successful offspring. Obviously the x and y axes (as structured in barplot2) will not be equivalent for these 2 sets of information and this shows up in my example below. 1) Is it possible to somehow reconcile the underlying x-axis to the same scale as would be needed to overly the normal distribution (e.g. where 2.5 would fall on the normal density, I could relate it to 2.5 on the barplot)? Then, using axis (side=4) I assume I could insert a y-axis for the normal distribution. 2) Is lines(density(x)) the appropriate way to insert a normal distribution into this type of figure? Should I use 'curve'? If someone could point me in the right direction, I would appreciate it. TIA, Bret Example: testdata 00.196454948 10.063515510 20.149187592 30.237813885 40.282127031 50.066469719 60.001477105 70.001477105 80.001477105 x-rnorm(1000, 2.84, 1.57) barplot2(testdata, xlab=Fledgling Number, ylab=Probability, ylim=c(0, 1), col=black, border=black, axis.lty=1) lines(density(x)) --Version-- platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor0.1 year 2004 month11 day 15 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Plotting Ranges as Vertical Lines
R Users, I have not been able to find anything close to what I want searching R-help and I am hoping someone could point me in the right direction. The data consists of differences in length of individual salamanders collected at time of initial capture and last recapture (excerpt below). What I am trying to do is plot a vertical line for each individual (ID) representing the change in size between initial capture (CAPTURESVL) and last recapture (RECAPSVL) (the range of sizes would be on the y-axis), hence vertical line length would be equal to GROWTH (which can be + or -). I would like to make this plot over the time between recaptures (MONTHELAPSED) which would be on the x-axis. ID CAPMONTH CAPTURESVL CAPTUREMASSGR RECAPMONTH RECAPSVL RECAPMASSGR GROWTH MASSGR MONTHELAPSED 1SJC83 21 0.1 11 22 0.2 10.18 2SJC93 33 0.7 13 35 0.9 20.2 10 3 SJC103 19 0.5 9 20 0.2 1 -0.36 4 SJC114 36 1.1 13 41 1.5 50.49 5 SJC124 19 0.1 15 28 0.4 90.3 11 6 SJC176 24 0.3 11 26 0.4 20.15 7 SJC186 43 1.1 24 43 1.8 00.7 18 8 SJC237 16 0.1 15 20 0.1 40.08 9 SJC247 22 0.5 22 40 1.9 181.4 15 10 SJC258 23 0.4 15 27 0.6 40.27 11 SJC329 37 1.1 16 38 1.5 10.47 12 SJC349 31 0.8 12 33 1.1 20.33 13 SJC359 19 0.2 14 23 0.3 40.15 14 SJC53 11 41 1.2 12 40 1.4 -10.21 15 SJC55 11 41 1.3 12 40 1.7 -10.41 16 SJC60 11 39 1.5 22 41 1.8 20.3 11 17 SJC65 12 41 1.5 23 44 1.8 30.3 11 I am stuck on how to 1) create the vertical segments and 2) link them to x-axis values for MONTHELAPSED? I would have provided some example code, except I have not been able to figure out how to approach this yet. Since I am probably not looking in the right place, could someone point me in the right direction or towards some example figure code/help files that I must have missed in my searching? Thanks, Bret TX AM platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor0.1 year 2004 month11 day 15 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Cumulative Points and Confidence Interval Manipulation in barplot2
R-Users, I am working with gplots (in gregmisc bundle) plotting some posterior probabilities (using barplot2) of harvest bag limits for discrete data (x-axis from 0 to 12, data is counts) and I ran into a couple of questions whose solutions have evaded me. 1) When I create and include the confidence intervals, the lower bound of the confidence intervals for several of the posterior probabilities is below zero, and in those specific cases I only want to show the upper limit for those CI's so they do not extend below the x-axis (as harvest can not be 0). Also, comments on a better technique for CI construction when the data is bounded to be =0 would be appreciated. 2) I would also like to show the cumulative probability (as say a point or line) across the range of the x-axis on the same figure at the top, but I have been unable to figure out how to overlay a set of cumulative points over the barplot across the same range as the x-axis. Below is some example code showing the test data I am working on (xzero): xzero - table(factor(WWNEW[HUNTTYPE==DOVEONLY], levels=0:12)) xzero 0 1 2 3 4 5 6 7 8 9 10 11 12 179 20 9 2 2 0 1 0 0 0 0 0 0 n - sum(xzero) k - sum(table(xzero)) meantheta1 -((2*xzero + 1)/(2*n + k)) vartheta1 -((2*(((2*n)+k)-((2*xzero)+1)))*((2*xzero)+1))/2*n)+k)^2)*(((2*n)+k)+2)) stderr - sqrt(vartheta1) cl.l - meantheta1-(stderr*2)#Fake CI: Test cl.u - meantheta1+(stderr*2)#Fake CI: Test barplot2(meantheta1, xlab=WWD HARVEST DOVE ONLY 2001, ylab=Probability, ylim=c(0, 1),xpd=F, col=blue, border=black, axis.lty=1,plot.ci=TRUE, ci.u = cl.u, ci.l = cl.l) title(main=WHITE WING DOVE HARVEST PROBABILITIES: DOVE HUNT ONLY) I would greatly appreciate any direction or assistance, Thanks, Bret platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor0.1 year 2004 month11 day 15 language R *Note: I am working in Exmacs __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Measuring Variational Distance
R-Users, Does anyone know if there is a package or code somewhere in R that provides a measure of the maximum difference between 2 distributions defined on a common space? I think it is called variational distance? I have constructed several marginal distribution plots (proportion/sample where I defined bin-widths) showing the difference between 2 treatments and I have been trying to get a measure of the difference in area between the curves. I tried searches using variational distance on the R website and CRAN with no luck. TIA, Bret Collier Univ. Arkansas __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Histograms, density, and relative frequencies
R-users, I have been using R for about 1 year, and I have run across a couple of graphics problem that I am not quite sure how to address. I have read up on the email threads regarding the differences between density and relative frequencies (count/sum(count) on the R list, and I am hoping that someone could provide me with some advice/comments concerning my approach. I will admit that some of the underlying mathematics of the density discussion are beyond my current understanding, but I am looking into it. I have a data set (600,000 obs) used to parameterize a probabilistic causal model where each obs is a population response for one of 2 classes (either regs1 and regs2). I have been attempting to create 1 marginal probability plot with 2 lines (one for each class). Using my rather rough code, I created a plot that seems to adhere to the commonly used (although from what I can understand wrong) relative frequency histogram approach. My rough code looks like this: bk - c(0, .05, .1, .15, .2, .25,.3, .35, 1) par(mfrow=c(1, 1)) fawn1 - hist(MFAWNRESID[regs1], plot=F, breaks=bk) fawn2 - hist(MFAWNRESID[regs2], plot=F, breaks=bk) count1 - fawn1$counts/sum(fawn1$counts) count2 - fawn2$counts/sum(fawn2$counts) b - c(0, .05, .1, .15, .2, .25, .3, .35) plot(count1~b,xaxt=n, xlim=c(0, .5), ylim=c(0, .40), pch=., bty=l) lines(spline(count1~b), lty=c(1), lwd=c(2), col=black) lines(spline(count2~b), lty=c(2), lwd=c(2), col=black) axis(side=1, at=c(0, .05, .1, .15, .2, .25, .3, .35)) Using the above, I get frequency values for regs1 that look like this (which is the same as output for my probabilistic model): count1 [1] 1.213378e-01 3.454324e-01 3.365343e-01 1.580839e-01 3.342101e-02 [6] 4.698426e-03 4.488942e-04 4.322685e-05 First, count1 is the frequency of occurrence within range 0-0.05, but when plotted is the value at b=0 and does not really represent the range? Are there any suggestions on a technique to approach this? Next: Using the above code, the x-axis values end at 0.35, but the axis continues (because bk ends at 1)? While there is the chance of occurrence out past .35, it is low and I want to extend the lines to about .35 and clip the x-axis. But, I have been unable to figure out how to clip Could someone point me in the correct direction? TIA, Bret A. Collier Arkansas Cooperative Fish and Wildlife Research Unit Department of Biological Sciences University of Arkansas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Setting max values for rpois
R-users, I am simulating a birth process for 4 classes of individuals with l[i] being the average No. fetuses per individual. However, I need to bound the resulting values for each generated rpois to be =3 (no individual can have 3 offspring). I have not been able to figure out how to incorporate this into the below example. Any suggestions on integrating would be appreciated. recruit.f - c(12, 12, 25, 51) #No. females in each age class l - c(.05, 1.22, 1.6, 1.8) #mean No. fetuses for each age class x - sapply(lapply(1:4, function(i) rpois(recruit.f[i], l[i])), sum) TIA, Bret A. Collier Arkansas Cooperative Fish and Wildlife Research Unit University of Arkansas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Weird Error
R-Users,
I hope this is not a uniformed question, but I am a little lost.
I ran into a problem this morning and I was wondering if anyone had
seen it before. I was trying to summarize each column of a data set
(150,000 rows, ~50mb, so it was a relatively big file) imported from a text
file using the below code;
data.summary - read.csv(c:/summary.txt, sep=)
data.summary - as.matrix(data.summary)
my.summary - function(x){
return(c(min=min(x),max=max(x), mean=mean(x)))}
apply(data.summary, 2, my.summary)
And I got this weird error that I can not find out anything about?
Process R unknown signal at Wed Apr 14 08:17:22 2004
Have you seen anything like this before? Do you think it is the size of
the dataset that is causing the problem, since the same code works for
25000 rows (~17mb) and gives the correct results (I cross-checked in SAS
and EXCEL).
I was using R 1.8.0 in Xemacs.
TIA,
Bret Collier
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Matrix decomposition
Gideon, You might look into Caswell, H. 2001. Matrix population models 2nd edition. Sinauer Associates. Caswell's book covers a majority of matrix population modeling, including some good introductory information and has general matrix manipulation as an appendix. I think most of the work in this book was done in Matlab. HTH, Bret GIDEON WASSERBERG I am looking for a manual(s) or any kind of documentation, at the introductory level, regarding matrix algebra (specifically, matrix population models). Any help will be highly appreciated Gideon Gideon Wasserberg (Ph.D.) Wildlife research unit, Department of wildlife ecology, University of Wisconsin 218 Russell labs, 1630 Linden dr., Madison, Wisconsin 53706, USA. Tel.:608 265 2130, Fax: 608 262 6099 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Bret __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Evaluating AIC
R Users, I was just wondering if anyone has written a program (or if there is a package) out there that calculates the different derivations of AIC (e.g. AIC, AICc, QAIC, etc.) along with AIC differences (delta's), model likelihoods, Akaike weights and evidence ratio's (from Burnham and Anderson 2002). Just in a for instance if someone had the -2LL, sample sizes, parameter counts, and estimates of c_hat output from a program, is there a function out there that calculates the above information. I did not see anything on the help pages (or in packages, but I could have missed it) and I didn't want to re-invent the wheel. TIA, Bret A. Collier Arkansas Cooperative Fish and Wildlife Research Unit Department of Biological Sciences SCEN 513 University of Arkansas Fayetteville, AR 72701 [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Novice problems with write()
R-Users,
As a relatively new user of R, I have a quick (and probably
simple) question about using write(). I have a population simulation that
I am running and I want to output a set of variables for each run of the
simulation into a text file for use in another program. However, whenever
I attempt to use write(), the only output that I am able to get is the
final numbers from the simulation.
for example:
x - 5
for (i in 1:10){
z - x+i
print(z)
write(z, c:/test.txt)
}
In this simple case, with print(z) I can see that z has what I am looking
for, but all that is output for the write statement is 15; While this is
simplified, it shows my problem.
I searched the help files, and on the R website, but I could not find
anything addressing this. I suspect that it is my lack of knowledge and I
am missing something obvious (or should be using write.table). If anyone
could point me in the right direction I would appreciate it.
Thanks,
Bret Collier
Univ. Arkansas
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] New User Question regarding simulations
R-Users, Everyone, I am a new user of R and I was wondering if anyone could point me to a reference (book or article) that discusses writing population simulations in R (or S). Thanks in advance, Bret A. Collier Arkansas Cooperative Fish and Wildlife Research Unit Department of Biological Sciences SCEN 513 University of Arkansas Fayetteville, AR 72701 (479) 575-4720 [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
