indsp.func <- function(sp, dfvec = c(4, 7, 10, 15, 20, 33))
{
# indsp.func:
# Calculates the data matrices of indices for different df.
# Arguments are the species data frame sp,  
# with column headers marked "site", "year", and "count", 
# and a list of required degrees of freedom (could be just one) for
# a selection of GAM fits. 
#
# All rows with NA's (missing values) must be removed from the data before
# applying this function.
# Result is a matrix with columns giving indices for each of the Nyears 
# years, with df as marked at the top of the column.
# 
# example: species data frame is called cb (for common bird):
# at command line:   indcb <- indsp.func(cb, c(4, 7, 10, 15, 20, 33))
#
	if(length(sp$count[is.na(sp$count)]) > 0) stop(
			"no missing data allowed")
	assign("sp", sp, frame = 1)
	sitemin <- min(sp$site)	# (any old site) 
	Nyears <- length(unique(sp$year))
	sp.new <- data.frame(year = seq(1, Nyears), site = rep(sitemin, Nyears)
		)	#
# sp.new is a new data frame with as few observations as possible to save
# time with extracting the indices.
	assign("sp.new", sp.new, frame = 1)	#
# fit.func fits the GAM and extracts the indices for a given value of df:
	fit.func <- function(dfval)
	{
		formula.name <- substitute(count ~ as.factor(site) + s(year, df
			 = dfval1) - 1, list(dfval1 = dfval))	#
# this waltz with formula.name is necessary so that both function gam()
# and predict.gam() understand how to find the value of dfval.
		gam.df <- gam(formula.name, family = poisson(link = log), data
			 = sp)	#
		pred.df <- predict.gam(gam.df, type = "terms", terms = labels(
			gam.df)[2], newdata = sp.new)
		pred.df <- as.vector(pred.df)
		ind.df <- exp(pred.df)/exp(pred.df[1])	#
# indices are scaled around the base year: chosen here as year 1. 
		cat("df ", dfval, " complete \n")
		ind.df
	}
	ind.all <- sapply(dfvec, fit.func)
	dimnames(ind.all) <- list(character(0), paste("df", as.character(dfvec),
		sep = ""))
	ind.all
}