Quick-R: a great R tutorial site

October 23, 2007 · Posted in R bloggers · Comments Off 

Quick-R, by Robert Kabacoff, is a wonderful R introduction site. It covers data management, basic and advanced statistics, and graphing in R, and it is aimed at an audience that has previous experience using other packages (such as SAS or Stata) that work differently than R.


Tags:

Approximate sunrise and sunset times

October 17, 2007 · Posted in R bloggers · Comments Off 
This function is not perfect, but it does a reasonable job estimating sunrise and sunset times for my field site. If more accurate data are required, try here.
suncalc<-function(d,Lat=48.1442,Long=-122.7551){
  ## d is the day of year
  ## Lat is latitude in decimal degrees
  ## Long is longitude in decimal degrees (negative == West)
  
  ##This method is copied from:
  ##Teets, D.A. 2003. Predicting sunrise and sunset times.
  ##  The College Mathematics Journal 34(4):317-321.
 
  ## At the default location the estimates of sunrise and sunset are within
  ## seven minutes of the correct times (http://aa.usno.navy.mil/data/docs/RS_OneYear.php)
  ## with a mean of 2.4 minutes error.

  ## Function to convert degrees to radians
  rad<-function(x)pi*x/180
  
  ##Radius of the earth (km)
  R=6378
  
  ##Radians between the xy-plane and the ecliptic plane
  epsilon=rad(23.45)

  ##Convert observer's latitude to radians
  L=rad(Lat)

  ## Calculate offset of sunrise based on longitude (min)
  ## If Long is negative, then the mod represents degrees West of
  ## a standard time meridian, so timing of sunrise and sunset should
  ## be made later.
  timezone = -4*(abs(Long)%%15)*sign(Long)

  ## The earth's mean distance from the sun (km)
  r = 149598000

  theta = 2*pi/365.25*(d-80)

  z.s = r*sin(theta)*sin(epsilon)
  r.p = sqrt(r^2-z.s^2)

  t0 = 1440/(2*pi)*acos((R-z.s*sin(L))/(r.p*cos(L)))
  
  ##a kludge adjustment for the radius of the sun
  that = t0+5 

  ## Adjust "noon" for the fact that the earth's orbit is not circular:
  n = 720-10*sin(4*pi*(d-80)/365.25)+8*sin(2*pi*d/365.25)

  ## now sunrise and sunset are:
  sunrise = (n-that+timezone)/60
  sunset = (n+that+timezone)/60

  return(list("sunrise" = sunrise,"sunset" = sunset))
}

Tags: , ,

Convert polar coordinates to Cartesian

October 14, 2007 · Posted in R bloggers · Comments Off 
When I want to calculate the coordinates of a location (e.g., a nest or burrow) based on distance and bearing from a grid point, this function helps me avoid writing down SOH-CAH-TOA every time. Just note that the bearing in this case is from the grid point (known location) to the unknown location.

polar2cart<-function(x,y,dist,bearing,as.deg=FALSE){
  ## Translate Polar coordinates into Cartesian coordinates
  ## based on starting location, distance, and bearing
  ## as.deg indicates if the bearing is in degrees (T) or radians (F)
  
  if(as.deg){
    ##if bearing is in degrees, convert to radians
    bearing=bearing*pi/180
  }
  
  newx<-x+dist*sin(bearing)  ##X
  newy<-y+dist*cos(bearing)  ##Y
  return(list("x"=newx,"y"=newy))
}


##Example

oldloc=c(0,5) 
bearing=200 #degrees
dist=5

newloc<-polar2cart(oldloc[1],oldloc[2],dist,bearing,TRUE)
plot(oldloc[1],oldloc[2],xlim=c(-10,10),ylim=c(-10,10))
points(newloc$x,newloc$y,col="red")

Tags: ,

If you don’t use R to fit multilevel models, there’s no point in reading this one

October 6, 2007 · Posted in R bloggers · Comments Off 
David Afshartous writes, What is the difference between sim() and mcsamp(), if any? If there is a difference, what should I expect in terms of difference in results? My reply:...

Tags: , ,

GillespieSSA 0.3-1 released

October 4, 2007 · Posted in R bloggers · Comments Off 

I recently rolled up the new version of the GillespieSSA package, version 0.3-1. The tar ball of the new version is posted on its web page (here). I also submitted it to CRAN so in (due time) it should appear on the official R package list.

The release consists of a number of bug fixes (among others typos in the package URL and demo models). The biggest change for the end user is, however, a now and cleaner way of passing model parameters to the main interface function ssa(). Rather than having to define model parameters in the global environment or hard coding their values into the propensity functions they can now be passed as a formal argument in the form of a named vector directly to ssa(). I owe a thanks to Ben Bolker for suggesting this deuglification solution. An added benefit of this revised feature is that model definitions are actually shorter now. The whole model + running it can now easily be a one liner.

For example defining and running the stochastic version of the classical logistic growth model \frac{{\displaystyle dN}}{{\displaystyle dt}}=rN \left( 1-\frac{{\displaystyle N}}{\displaystyle K} \right) is now as simple as

library(GillespieSSA)
out <- ssa(x0=c(N=500),a=c(”b*{N}”, ”(d+(b-d)*{N}/K)*{N}”),nu=matrix(c(+1,-1),ncol=2),parms=c(b=2, d=1, K=1000),tf=25)

Here the per capita birth rate is b=2, death rate d=1, intrinsic growth rate r=b-d=1, carrying capacity K=1000. The Monte Carlo simulation is run for ten time units (t=10) using Gillespie’s Direct method (the default method in GillespieSSA) and starting with a population size of N=500.

One way of quickly visualizing the resulting time series is by

ssa.plot(out)

As expected the results look rather stochastically familiar,
ssaplot.png


Tags: , ,

Special Volume of JSS on Ecology and Ecological Modelling in R

October 3, 2007 · Posted in R bloggers · Comments Off 

The latest issue Journal of Statistical Software is a special volume on Ecology and Ecological Modelling in R. This has been in the woodwork for quite some time now and it is nice to finally to see this out “in print”. I am looking forward to many joyful moment thumbing through the papers in this issue. Here’s a run down of the contributions,

  • Introduction to the Special Volume on “Ecology and Ecological Modelling in R”. Thomas Kneib, Thomas Petzoldt
  • Analogue Methods in Palaeoecology: Using the analogue Package. Gavin L. Simpson
  • Maximum Likelihood Method for Predicting Environmental Conditions from Assemblage Composition: The R Package bio.infer. Lester L. Yuan
  • The ade4 Package: Implementing the Duality Diagram for Ecologists. Stephane Dray, Anne-Beatrice Dufour
  • Interactive Multivariate Data Analysis in R with the ade4 and ade4TkGUI Packages. Jean Thioulouse, Stephane Dray
  • Exploring Habitat Selection by Wildlife with adehabitat. Clement Calenge
  • The ecodist Package for Dissimilarity-based Analysis of Ecological Data. Sarah C. Goslee, Dean L. Urban
  • Statistical Methods for the Qualitative Assessment of Dynamic Models with Time Delay (R Package qualV). Stefani Jachner, K. Gerald van den Boogaart
  • simecol: An Object-Oriented Framework for Ecological Modeling in R. Thomas Petzoldt, Karsten Rinke
  • demogR: A Package for the Construction and Analysis of Age-structured Demographic Models in R. James Holland Jones
  • Estimating and Analyzing Demographic Models Using the popbio Package in R. Chris Stubben, Brook Milligan
  • Introducing untb, an R Package For Simulating Ecological Drift Under the Unified Neutral Theory of Biodiversity. Robin Hankin

Tags: , ,

Reorder factor levels

October 3, 2007 · Posted in R bloggers · Comments Off 
Very often, especially when plotting data, I need to reorder the levels of a factor because the default order is alphabetical. There must be many ways of reordering the levels; however, I always forget which package to look in. A direct way of reordering, using standard syntax is as follows:
## generate data
x = factor(sample(letters[1:5],100, replace=TRUE))

print(levels(x))  ## This will show the levels of x are "Levels: a b c d e"

## To reorder the levels:
## note, if x is not a factor use levels(factor(x))
x = factor(x,levels(x)[c(4,5,1:3)])

print(levels(x))  ## Now "Levels: d e a b c"
In this case, the level order could be set in the first line; however, if there are many levels (and you don't want to type out all of the levels explicitly), the above method is preferred. Again, if there is a better way to do this, please let me know in the comments.

Tags: , ,

R video tutorial number 2

October 2, 2007 · Posted in R bloggers · Comments Off 

READ TEXT FILES, RUN MODELS

rtut2.gif

The Decision Science News R video tutorials continue with number 2. (If you missed that last one, you will want to watch R Video Tutorial Number 1 first.) The Goldstein pedometer dataset can be downloaded from http://www.dangoldstein.com/flash/Rtutorial2/pedometer.csv

Topics covered this week include:

  • Tricking R into starting in your working directory
  • Reading in text files
  • Accessing columns in data frames
  • Creating histograms
  • Side effects and optional parameters of function calls
  • Fitting simple linear models
  • Adding least-squares and loess lines to plots

The commands in the tutorial are:

myData=read.table("pedometer.csv", header=TRUE, sep=",")
x=hist(myData$Steps,col="lightblue")
x=hist(myData$Steps,breaks=20,col="lightblue")
plot(myData$Steps ~ myData$Observation,col="blue")
myModel=lm(myData$Steps ~ myData$Observation)
summary(myModel)
lines(fitted(myModel))
lines(fitted(loess(myModel)),col="red")

Can’t view flash? Download movie If you see no image under Windows, download the TSSC Codec.

Tags: , , 

Diag| Memory: Current usage: 35967 KB
Diag| Memory: Peak usage: 36238 KB