# This program 
# 1) determines the average birth rate for aphids that have previously reproduced
# 2) determines the average birth rate for aphids that have not previously reproduced
# 3) plots the average reproduction versus age with the 95% confidence interval
# 	and the fraction of aphids surviving up to the corresponding day

# The program requires one data file.

# predation.csv is a file that contains the predation data.
#	The first row contains column headers.
#	The first column contains the ID numbers for the cages.
#	Rows 2 and 3 are the before and after counts of aphids for the larva in the first cage.
#	After all larvae data, there is a blank row.
#	The following rows are before and after counts for adults.
#	No data should be included on the day the coccinellid died.

# DATA

data <- read.csv("predation.csv",header=TRUE)

# INITIALIZATION

blankrow <- which(is.na(data[,1]))

data <- 
rows <- 
days <- 

larvaecages <- 
adultcages <- 
larvae <- matrix(0,nrow=,ncol=)
adults <- matrix(0,nrow=,ncol=)

mu <- rep(0,days)					# average for larvae for each day
sigma <- rep(0,days)				# standard deviation for larvae for each day
allvalues <- rep(0,0)				# list of all adult predation values

# PREPARE DATA FILES

	
# COMPUTE MEANS AND 95% CONFIDENCE INTERVALS FROM DATA

for (day in 1:days)
	{
	thisday <- larvae[,day]					# selects the appropriate column
	values <- thisday[which(is.na(thisday)==FALSE)]	# removes the NA entries
	mu[day] <- mean(values)
	sigma[day] <- sd(values)
	}
conf95 <- 1.96*sigma/sqrt(cages)

# FUNCTIONS 

# errorbars is a low-level plotting command that adds error bars to a set of plotted points.
#	The first argument is the set of x values.
#	The second argument is the set of y values.
#	The third argument is the height of the bar above and below the point.
#	The standard optional arguments col and lwd are accepted.
#	The optional argument 'cap' is a multiplier for the cap length.

errorbars <- function(x,y,dy,col="black",lwd=1,cap=1)
	{
	n <- length(x)
	top <- y+dy
	btm <- y-dy
	d <- cap*0.01*(max(x)-min(x))
	for (i in 1:n)
		{
		lines( c(x[i]  ,x[i]  ), c(btm[i],top[i]), col=col, lwd=lwd )
		lines( c(x[i]-d,x[i]+d), c(btm[i],btm[i]), col=col, lwd=lwd )
		lines( c(x[i]-d,x[i]+d), c(top[i],top[i]), col=col, lwd=lwd )
		}
	}

# OUTPUT

par(mar=c(5.1,4.1,4.1,4.1)) 
plot(1:days, mu, pch=15, 
		xlim=(c(0,days+1)), ylim=(c(0,10)), xlab="days", ylab="fecundity")
errorbars(1:days, mu, conf95, cap=1.5)
points(1:days, 10*survival, pch=16, col="blue3")
lines(1:days, 10*survival, lty=3, col="blue3")
axis(4, at=2*0:5, labels=.2*0:5, col="blue3", col.axis="blue3")
mtext("survival", side=4, line=2.8, col="blue3")

m
m1
c95