TEXT S1: PROGRAMMING SCRIPT

Provided below are a set of S-Plus functions, custom-written by JM McPherson, to compute 
autologistic distribution models.  The main function is alogis; it in turn calls the 
other functions provided (neighbour, block, add.leftcol and add.rightcol).  By default, 
alogis uses an ordinary logistic regression as prior, but other priors can be supplied.  
It then iterates through a Gibbs sampler, updating blocks of independent cells in turn.  
These blocks are created using the function block.  The blocking pattern will depend on 
the neighbourhood structure used.  As set up currently, the neighbourhood structure 
corresponds to the Queen's case (eight neighbouring cells), and neighbours are identified 
by function neighbour.  Invalid neighbours (that were not sampled, contain no data, or 
fall outside the study region) are ignored in calculating the autocovariate.  To ensure 
that this is the case, the data matrix must be buffered on all sides. This happens 
automatically for the northern and southern edge of the study area; for West and East, 
add.leftcol and add.rightcol ensure a buffer.

alogis<-function(response,covariate,colnum,rownum,sim=20,burn=10,qual=NULL,prior=NULL){
	###function to compute autologistic regression
	###developed by Jana M. McPherson 
	###(jana.mcpherson@merton.oxon.org), 
	###primarily based on information in Augustin, Mugglestone & 
	###Buckland 1996 (Journal of Applied Ecology) and 1998 
	###(Environmetrics). The Gibbs' sampler borrows elements from 
	###Jennifer Hoeting's algorithm "autologit" (Hoeting, 
	###Leecaster & Bowden 2000, Journal of Agricultural Biological 
	###and Environmental Statistics).

        ###INPUT:
        #response:	a matrix with two columns specifying the location
	#		of observed responses submitted for model 
	#		training, as well as the response value (binary)
	#		observed; column 1 should give the location index
	#		and be named Record; column 2 should give the
	#		observed response and be named CODE
	#covariate:	a matrix with q=p+1 columns, where p is the number
        #		of candidate predictor variables. There will be as
        #		many rows as there are sites that you want 
        #		predictions for, including sites submitted for
        #		training in the variable response, above, and
        #		unsurveyed sites; column 1 should contain the
        #		location index, column 2-q should contain observed
        #		values of the candidate predictors at these sites;
        #		the matrix should have column names, with "Record"
        #		for the index variable, and appropriate names for
        #		each candidate predictor.
	#colnum:	number of columns in the matrix that the location
	#		indices are based on
	#rownum:	number of rows in the matrix that the location
        #		indices are based on
	#sim:		total number of iterations the Gibbs sampler
        #		should go through (default=20)
	#burn:	the number of Gibbs sampler iterations to discard,
        #		before outputs are stored for an average predicted
        #		response (default=10); you can set burn to NULL 
        #		and instead specify qual.
	#qual:	instead of fixing a burn-in period with burn, you 
        #		can specify a level of convergence ('quality') the 
        #		Gibbs sampler should achieve before outputs are 
        #		used to calculate an average response; indicate 
        #		the maximum percentage of lattice cells (among 
        #		those in covariate) you accept changes of more 
        #		than 0.05 in from one iteration to the next. E.g.
        #		if you set qual=10, Gibbs sampler outputs are 
        #		stored for calculation of an average response as 
        #		long as only 10% or less of the cells changed by 
	#		0.05 from the previous iteration. Default=NULL
	#prior:	this gives an opportunity to identify previously 
        #		selected predictors (including intercept in 1st 
        #		row) and corresponding parameter values; should 
        #		consist of a 2 column matrix, with predictor names 
        #		(equal to the relevant column names in covariate, 
        #		except for Intercept) in column 1, and estimated
        #		coefficients in column 2. If prior is NULL (the 
        #		default), alogis will compute priors by stepwise 
        #		glm, submitting all candidate predictors in 
        #		covariate			
	#############################################################

	###Important variables:
	cov.rec<-covariate[,"Record"]
	res.rec<-response[,"Record"]
	res.code<-response[,"CODE"]
	n<-nrow(covariate)	
        ###the number of sites to be predicted

	###Sites of interest (all those in covariate) should not touch 
        ###either of the lateral edges of the rectangular lattice 
        ###defined by colnum and rownum, as this is could lead to 
        ###wrapping (from left to right edge, or vice versa) when 
        ###calculating the autocovariate term later. So first we 
        ###should buffer sites of interest on the left or right 
        ###lattice edge, if necessary. 

	print("checking lateral edge")
	###Are any of the sites of interest on the lateral edges? 
        ###Check in cov.rec
	right.flag<-0
	left.flag<-0
	for(k in 1:n){
		if(cov.rec[k]%%colnum==0){right.flag<-1}
		if((cov.rec[k]-1)%%colnum==0){left.flag<-1}
	}
	###If there are sites on the left edge, shift indices 
        ###appropriately
	if(left.flag==1){
		tmp<-add.leftcol(cov.rec, colnum, rownum, addnum=1)
		cov.rec<-tmp$newpos
		tmp<-add.leftcol(res.rec, colnum, rownum, addnum=1)
		res.rec<-tmp$newpos
		colnum<-tmp$newcol
		print("left edge buffered")
	}
	###If there are sites on the right edge, shift indices 
        ###appropriately
	if(right.flag==1){
		tmp<-add.rightcol(cov.rec, colnum, rownum, addnum=1)
		cov.rec<-tmp$newpos
		tmp<-add.rightcol(res.rec, colnum, rownum, addnum=1)
		res.rec<-tmp$newpos
		colnum<-tmp$newcol
		print("right edge buffered")
	}

	###Next we need to make sure we have prior parameters to 
        ###initialise the model in the absence of an autocovariate. If 
        ###prior is NULL (as by default), derive prior parameters now.

	if(is.null(prior)){
		print("initialising prior")
		###put together a training dataset
		rec<-match(res.rec, cov.rec)
		trn<-covariate[rec,2:ncol(covariate)]	
        #exclude the "Record" column
		cov.names<-name.cols(trn)	
        #names of covariates
		trn<-cbind(res.code, trn)	
        #add the observed responses

        ###make sure trn is a data frame, otherwise glm doesn't 
        ###recognise data=trn
		if(!is.data.frame(trn)){trn<-as.data.frame(trn)}	
		dimnames(trn)[[2]]<-c("CODE",cov.names)

		###identify candidate predictors
		prds<-name.cols(trn)
		prds<-prds[2: length(prds)]		
        #this is to remove "CODE" from among predictors
		up.f<-paste(prds, collapse = "+")
		up.f<-paste("CODE~", up.f, collapse = "")

		###run stepwise glm
		basic.glm<-glm(CODE~1, family = binomial, data = trn, 
        maxit = 50)
		assign("trn", trn, frame=1)
		stepwise.glm<-step.glm(basic.glm, scope = list(upper= 
		up.f, lower= ~1), direction = "forward")
		prior.sum<-summary.glm(stepwise.glm)
		prior.vars<-name.rows(prior.sum$coefficients)
		prior.coef<-prior.sum$coefficients[,1]
		prior<-cbind(prior.vars, prior.coef)
		dimnames(prior)<-NULL
		print("priors initialised")
	}	

	###Now define some variables needed later
	tot<-colnum%*%rownum	
        ###the total number of sites in the (adjusted) lattice matrix

	###initialise g, which indexes irrelevant, unsurveyed and 
        ###surveyed sites in g 2=outside study area, 1=NOT surveyed 
        ###but to be predicted, 0=surveyed
	g<-rep(2, times=tot)	
	g[cov.rec]<-1
        ###initially put all sites of interest as NOT surveyed (1)
	g[res.rec]<-0
        ###then update to indicate training sites as surveyed (0)

	pars<-prior[,1]
        ###names of selected predictor variables
	pars.2<-pars[2:length(pars)]	
        ###this sheds name of the intercept
	pvals<-prior[,2]		
        ###parameter values corresponding to pars
	pvals<-c(pvals,0)	
        ###add zero as initial parameter for the autocovariate
	covs<-cbind(rep(1, n),covariate[pars.2], numeric(n))
        ###matrix of chosen covariates
	###with first column for intercept and last column for the 
        ###spatial autocovariate
	covs<-as.data.frame(covs)
	dimnames(covs)[[2]]<-c("Intercept", pars.2, "SPATIAL")
	last<-ncol(covs)
	print("variables defined")

	###	initialise prob.save, y, and obs.p
	prob.save<-rep(NA, times=tot)	
        ###this will store predicted probabilities
	prob.save[g==1]<-0	
        ###initially set all unsurveyed sites to probability zero
	prob.save[g==0]<-res.code	
        ###in surveyed sites use observed response
	y<-prob.save		
        ###this keeps the original set-up, to use in ifelse statements
	obs.prob<-prob.save[cov.rec]	
        ###this stores original observations for output
	not<-is.na(match(cov.rec, res.rec))	
        ###checks which sites were not surveyed
	obs.prob[not]<-NA	
        ###and sets them to NA
	print("prob.save initiated")

	###We now need to create a response surface for all sites of 
        ###interest, so that we can initiate the autologistic model.
        
	mcovs<-as.matrix(covs)	
        ###matrix version of covs for multiplication
	###update prob.save based on prior, but only if the site is 
        ###unsurveyed and a prediction is sought:
	prob.save<-ifelse(g==1,exp(mcovs%*%pvals)/
        (1+exp(mcovs%*%pvals)),ifelse(y==1,1,ifelse(y==0,0,NA)))
	###Unsurveyed sites can end up with NA values, because exp(a)
        ###returns 'inf' if a>709, and inf/(1+inf) returns NA. 
        ###Instead, a large number divided by itself plus 1 should
	###equal something very close to 1, so build in a correction 
        ###for that:
	prob.save[cov.rec]<-ifelse(is.na(prob.save[cov.rec]), 
        0.999999, prob.save[cov.rec])

	###we now have a response surface composed of observed and 
        ###simulated response values plus a value of NA outside the 
        ###study region where predictions are not sought
	print("preliminary updating of prob.save done")

	###next calculate the autocovariate.
	spatial<-neighbour(cov.rec, prob.save, colnum, rownum)
	covs[,last]<-spatial
	print("preliminary autocovariate calculated")

	###We now estimate new 'prior' parameter values by 
        ###incorporating the autocovariate as a predictor variable in 
        ###ordinary logistic regression. (This is equivalent to
        ###deriving parameter estimates by maximising the pseudo-
        ###likelihood.)

	###put together a training dataset
	CODE<-prob.save[res.rec]
	rec<-match(res.rec, cov.rec)
	trn<-cbind(CODE,covs[rec,2:last])	
        ###don't include the intercept column here
	trn<-as.data.frame(trn)	
        ###otherwise the data=trn argument in glm doesn't work
	cov.name<-name.cols(covs)
	dimnames(trn)[[2]]<-c("CODE", cov.name[2:last])

	###construct the formula
	prds<-name.cols(trn)
	prds<-prds[2: length(prds)]		
        ###this is to remove "CODE" from among predictors
	formul<-paste(prds, collapse = "+")
	formul<-paste("CODE~", formul, collapse = "")

	###run glm
	auto.glm<-glm(formula=formul, family=binomial, data=trn, 
        maxit = 50)
	new.prior<-summary.glm(auto.glm)
	print("autologistic model initiated")

	###update pvals
	pvals<-new.prior$coefficients[,1]
	print("new prior, constructed under inclusion of the 
        autocovariate")

	###Use the new coefficients to update prob.save before 
        ###initiating the Gibbs sampler

	mcovs<-as.matrix(covs)
        ###matrix version of covs, needed for multiplication
	###update prob.save only if the site is unsurveyed and a 
        ###prediction is sought:
	prob.save<-ifelse(g==1,exp(mcovs%*%pvals)/
        (1+exp(mcovs%*%pvals)),ifelse(y==1,1,ifelse(y==0,0,NA)))
	###Unsurveyed sites can end up with NA values, because exp(a) 
        ###returns 'inf' if a>709, and inf/(1+inf) returns NA. 
        ###Instead, a large number divided by itself plus 1 should
	###equal something very close to 1, so build in a correction 
        ###for that:
	prob.save[cov.rec]<-ifelse(is.na(prob.save[cov.rec]), 
        0.999999, prob.save[cov.rec])
		print("prob.save updated a second time")

	###Now set things up for the Gibbs sampler

	###initiate variables needed
	track.pvals<-matrix(0,sim+1,last)	
        ###saves parameter values at each iteration
	dimnames(track.pvals)<-list(NULL,name.cols(covs))
	track.pvals[1,]<-pvals
        ###store coefficients of the preliminary autologistic model
	probp<-numeric(n)
	###probability of presence, used later to calculate average
	prob.last<-prob.save[cov.rec]
	###used to monitor change in predictions
	cindex<-numeric(sim)
	###summarises change in predictions
	ccount<-0
	###counts how often specification 'qual' is met

	###identify the 4 'independent blocks' in the data
	blocks<-block(colnum, rownum)

	###filter relevant sites for block 1 (i.e. those we actually 
        ###want predictions for)
	b1<-blocks$b1
	aa<-match(b1,cov.rec)
	bb<-!is.na(aa)		
        ###true for useful elements of b1
	b1<-b1[bb]
        ###identifies appropriate index in prob.save, y and g
	s1<-aa[!is.na(aa)]	
        ###identifies appropriate index in covs
        
	###and block 2
	b2<-blocks$b2
	aa<-match(b2,cov.rec)
	bb<-!is.na(aa)		
        ###true for useful elements of b2
	b2<-b2[bb]
        ###identifies appropriate index in prob.save, y and g
	s2<-aa[!is.na(aa)]
        ###identifies appropriate index in covs

        ###and block 3
	b3<-blocks$b3
	aa<-match(b3,cov.rec)
	bb<-!is.na(aa)
        ###true for useful elements of b3
	b3<-b3[bb]
        ###identifies appropriate index in prob.save, y and g
	s3<-aa[!is.na(aa)]	
        ###identifies appropriate index in covs

        ###and block 4
	b4<-blocks$b4
	aa<-match(b4,cov.rec)
	bb<-!is.na(aa)
        ###true for useful elements of b4
	b4<-b4[bb]
        ###identifies appropriate index in prob.save, y and g
	s4<-aa[!is.na(aa)]
        ###identifies appropriate index in covs
	print("variables for Gibbs sampler initiated")

	###start the Gibbs sampler SIMULATION LOOP
	for(gg in 1:sim){
		print(paste("Gibbs iteration",gg))

	###UPDATE SPATIAL AUTOCOVARIATE for each independent 
        ###block
		###BLOCK1~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		###calculate the autocovariate
		spatial<-neighbour(b1, prob.save, colnum, rownum)
		covs[s1,last]<-spatial
		mcovs<-as.matrix(covs)
		###then update it based on current parameter values
		prob.save[b1]<-ifelse(g[b1]==1,exp(mcovs[s1,]%*%pvals)/
        	(1+exp(mcovs[s1,]%*%pvals)),ifelse(y[b1]==1,1,
        	ifelse(y[b1]==0,0,NA)))
		###Unsurveyed sites can end up with NA values, because 
        	###exp(a) returns 'inf' if a>709, and inf/(1+inf) 
        	###returns NA. Instead, a large number divided by itself 
        	###plus 1 should equal something very close to 1, so 
        	###build in a correction for that:
		prob.save[cov.rec]<-ifelse(is.na(prob.save[cov.rec]), 
        	0.999999, prob.save[cov.rec])

		###BLOCK 2~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		###calculate the autocovariate
		spatial<-neighbour(b2, prob.save, colnum, rownum)
		covs[s2,last]<-spatial
		mcovs<-as.matrix(covs)
		###then update it based on current parameter values
		prob.save[b2]<-ifelse(g[b2]==1,exp(mcovs[s2,]%*%pvals)/
        	(1+exp(mcovs[s2,]%*%pvals)),ifelse(y[b2]==1,1,
        	ifelse(y[b2]==0,0,NA)))
		###Unsurveyed sites can end up with NA values, because 
        	###exp(a) returns 'inf' if a>709, and inf/(1+inf) 
        	###returns NA. Instead, a large number divided by itself 
        	###plus 1 should equal something very close to 1, so 
        	###build in a correction for that:
        	prob.save[cov.rec]<-ifelse(is.na(prob.save[cov.rec]), 
        	0.999999, prob.save[cov.rec])

		###BLOCK 3~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		###calculate the autocovariate
		spatial<-neighbour(b3, prob.save, colnum, rownum)
		covs[s3,last]<-spatial
		mcovs<-as.matrix(covs)
		###then update it based on current parameter values
		prob.save[b3]<-ifelse(g[b3]==1,exp(mcovs[s3,]%*%pvals)/
        	(1+exp(mcovs[s3,]%*%pvals)),ifelse(y[b3]==1,1,
        	ifelse(y[b3]==0,0,NA)))
		###Unsurveyed sites can end up with NA values, because 
        	###exp(a) returns 'inf' if a>709, and inf/(1+inf) 
        	###returns NA. Instead, a large number divided by itself 
        	###plus 1 should equal something very close to 1, so 
        	###build in a correction for that:
		prob.save[cov.rec]<-ifelse(is.na(prob.save[cov.rec]), 
        	0.999999, prob.save[cov.rec])

		###BLOCK 4~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		###calculate the autocovariate
		spatial<-neighbour(b4, prob.save, colnum, rownum)
		covs[s4,last]<-spatial
		mcovs<-as.matrix(covs)
		###then update it based on current parameter values
		prob.save[b4]<-ifelse(g[b4]==1,exp(mcovs[s4,]%*%pvals)/
        	(1+exp(mcovs[s4,]%*%pvals)),ifelse(y[b4]==1,1,
        	ifelse(y[b4]==0,0,NA)))
		###Unsurveyed sites can end up with NA values, because 
        	###exp(a) returns 'inf' if a>709, and inf/(1+inf) 
        	###returns NA. Instead, a large number divided by itself 
        	###plus 1 should equal something very close to 1, so 
        	###build in a correction for that:
		prob.save[cov.rec]<-ifelse(is.na(prob.save[cov.rec]), 
        	0.999999, prob.save[cov.rec])
		###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        	print("autocovariates updated in blocks")

		###RE-CALIBRATE MODEL PARAMETERS
		###update the training dataset
		CODE<-prob.save[res.rec]
		rec<-match(res.rec, cov.rec)
		trn<-cbind(CODE,covs[rec,2:last])
        	###don't include the intercept column here
		trn<-as.data.frame(trn)	
        	###otherwise the data=trn argument in glm doesn't work
		cov.name<-name.cols(covs)
		dimnames(trn)[[2]]<-c("CODE", cov.name[2:last])
						
		###run glm
		auto.glm<-glm(formula=formul, family=binomial, data=trn, 
        	maxit = 50)
		new.prior<-summary.glm(auto.glm)
		###update pvals
		pvals<-new.prior$coefficients[,1]
		print("parameters updated")
        	###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		###MONITOR AND TRACK Gibbs sampler outputs

		###C-index: percentage of pixels that change by more 
        	###than 0.05 per iteration
		cindex[gg]<-sum(ifelse(abs(prob.last-
        	prob.save[cov.rec])>0.05,1,0))
		cindex[gg]<-100*(cindex[gg]/n)

		###save current predictions for next computation of 
        	###cindex
		prob.last<-prob.save[cov.rec]

		###store current parameter values
		track.pvals[gg+1,]<-pvals

		###add current predictions for calculation of average, 
        	###if appropriate
		if(!is.null(burn)){
			if(gg>burn){probp<-prob.save[cov.rec]+probp}
		}
		if(!is.null(qual)){
			if(cindex[gg]<=qual){
				probp<-prob.save[cov.rec]+probp
				ccount<-ccount+1
			}
		}

		###save model statistics for output if this is the final 
        	###iteration
		if(gg==sim){
			stats<-c(new.prior$df[2], new.prior$null.deviance, 
        		new.prior$deviance)
			names(stats)<-c("residual.df", "null.deviance", 
        		"residual.deviance")
		}
		print("monitored")
	}	###END OF GIBBS SAMPLER


	###OUTPUT for alogis

	###average predicted probability of occurrence
	if(!is.null(burn)){
		probp<-probp/(sim-burn)
		av.method<-"burn"
	}
	if(!is.null(qual)){
		if(ccount!=0) {probp<-(probp/ccount)}
		av.method<-"qual"
	}
	avg<-matrix(NA, nrow=n, ncol=2)
	avg[,1]<-cov.rec
	avg[,2]<-probp
	dimnames(avg)[[2]]<-c("Record", "Mean.Prob")

	###final predicted probability of occurrence, with location 
        ###index
	final<-matrix(NA, nrow=n, ncol=2)
	final[,1]<-cov.rec
	final[,2]<-prob.save[cov.rec]
	dimnames(final)[[2]]<-c("Record", "End.Prob")

	###difference between observed and predicted response
	###NB: this will be 0 for all surveyed and NA for unsurveyed 
        ###sites, unless you allow surveyed sites--or perhaps 
        ###observation of absence but not presence÷to be updated in 
        ###the preliminary model and/or Gibbs sampler.
	diffs<-matrix(NA, nrow=n, ncol=2)
	diffs[,1]<-cov.rec
	diffs[,2]<-final[,2]-obs.prob
	dimnames(diffs)[[2]]<-c("Record", "Difference")

	###observed probability with location index
	obsv<-matrix(NA, nrow=n, ncol=2)
	obsv[,1]<-cov.rec
	obsv[,2]<-obs.prob
	dimnames(obsv)[[2]]<-c("Record", "Obs.Prob")

	###chosen predictors
	cov.name<-name.cols(covs)

	list(final.p=final, mean.p=avg, obs.p=obsv, diff.p=diffs, 
        p.names=cov.name, p.coef=pvals, all.coefs=track.pvals, 
        meth=av.method, change=cindex, sims=sim, burns=burn, 
        quals=ccount, mstats=stats)		
}	###THE END OF ALOGIS





add.leftcol<-function(positn,colnums,rownums,addnum=1){
	###Developed by Jana M. McPherson
	###Called by alogis
	###This shifts location indices in position as if a column is 
        ###added to the left side of the matrix.
        
	###INPUT
	#positn:	should be a vector of location indices relevant to 
        #		a rectangular lattice defined by colnum and 
        #		rownum, where numbering of locations goes across 
        #		then down (i.e. row by row).
	#colnums:	should give the number of columns in the original 
        #		lattice.
	#rownums:	should give the number of rows in the original 
        #		lattice.
	#addnum:	indicates how many columns to add on the left; 
        #		default=1.

	new<-NULL
	for(k in 1:length(positn)){
		for(r in 1:rownums){
		###print(c(positn[k],colnums*r,colnums*(r-1)))
			if(positn[k]<=colnums*r && 
        positn[k]>colnums*(r-1)){
					new[k]<-positn[k]+(r*addnum)
			}
		}
	}
	colnums<-colnums+addnum
	list(newpos=new, newcol=colnums)
} ###END OF ADD.LEFTCOL



add.rightcol<-function(positn,colnums,rownums,addnum=1){
	###Developed by Jana M. McPherson
	###Called by alogis
	###This shifts location indices in position as if a column is 
        ###added to the right side of the matrix.
        
	###INPUT
	#positn:	should be a vector of location indices relevant to 
        #		a rectangular lattice defined by colnum and 
        #		rownum, where numbering of locations goes across 
        #		then down (i.e. row by row).
	#colnums:	should give the number of columns in the original 
        #		lattice.
	#rownums:	should give the number of rows in the original 
        #		lattice.
	#addnum:	indicates how many columns to add on the right; 
        #		default=1.
	
	new<-NULL
	for(k in 1:length(positn)){
		for(r in 1:rownums){
			#print(c(positn[k],colnums*r,colnums*(r-1)))
			if(positn[k]<=colnums*r && 
        positn[k]>colnums*(r-1)){
					new[k]<-positn[k]+((r-1)*addnum)
			}
		}
	}
	colnums<-colnums+addnum
	list(newpos=new, newcol=colnums)
} ###END OF ADD.RIGHTCOL




block<-function(colnums, rownums){
	#Developed by Jana M. McPherson
	#Called by alogis
	#This function divides image or matrix indices into 4 blocks.
	#Assuming a second-order neighbourhood structure, elements 
        #within each blocks are independent of each other, and thus 
        #can be updated simultaneously in a Gibbs sampling routine.
	#INPUTS:
	#colnums = number of columns in the raster image or matrix
	#rownums = number of rows in the raster image or matrix
	
	cols.1<-colnums/2
	cols.2<-ifelse(colnums%%2==0, colnums/2, 1+colnums/2)
	#cols.2 must be adjusted if there is an uneven number of 
        #columns
	
	rows.1<-rownums/2
	rows.2<-ifelse(rownums%%2==0, rownums/2, 1+rownums/2)
	#rows.2 must be adjusted if there is an uneven number of rows
	
	#initialise block.1
	strt<-NULL
	set<-matrix(data=NA, nrow=cols.2, ncol=rows.2)
	for(i in 1:(rows.2)){
		strt[i]<-ifelse(i==1, 1, strt[i-1]+(2*colnums))
		for(j in 1:(cols.2)){
			set[j,i]<-strt[i]+(2*(j-1))
		}
	}
	block.1<-as.vector(set)
	
	#initialise block.2
	strt<-NULL
	set<-matrix(data=NA, nrow=cols.1, ncol=rows.2)
	for(i in 1:(rows.2)){
		strt[i]<-ifelse(i==1, 2, strt[i-1]+(2*colnums))
		for(j in 1:(cols.1)){
			set[j,i]<-strt[i]+(2*(j-1))
		}
	}
	block.2<-as.vector(set)
	
	#initialise block.3
	strt<-NULL
	set<-matrix(data=NA, nrow=cols.2, ncol=rows.1) 
	for(i in 1:(rows.1)){
		strt[i]<-ifelse(i==1, colnums+1, strt[i-1]+(2*colnums))
		for(j in 1:(cols.2)){
			set[j,i]<-strt[i]+(2*(j-1))
		}
	}
	block.3<-as.vector(set)
	
	#initialise block.4
	strt<-NULL
	set<-matrix(data=NA, nrow=cols.1, ncol=rows.1) 
	for(i in 1:(rows.1)){
		strt[i]<-ifelse(i==1, colnums+2, strt[i-1]+(2*colnums))
		for(j in 1:(cols.1)){
			set[j,i]<-strt[i]+(2*(j-1))
		}
	}
	block.4<-as.vector(set)
	
	#save all output
	list(b1=block.1, b2=block.2, b3=block.3, b4=block.4)
}  ###END OF BLOCK




neighbour<-function(sites, probs, colnums, rownums){
	###Developed by Jana M. McPherson
	###Called by alogis
	###This calculates a spatial autocovariate with second order 
        ###neighbourhood. It sums occurrence (or probability thereof) 
        ###in 8 neighbouring cells, but ignores non-existing 
        ###neighbours and neighbours with value NA for probs. 
	###These neighbours are considered irrelevant, and so not 
        ###taken into account when the sum is divided by the total 
        ###number of neighbours. If a cell has no neighbours at all, 
        ###its autocovariate is set to 0, to avoid NA values in the 
        ###output (these would arise from dividing zero by zero).

        ###INPUT:
	###sites:	should list the location index (within a 
        #		rectangular lattice) of the cells for which the 
        #		autocovariate is to be calculated
	###probs:	should contain the variable (e.g. probability of 
        #		occurrence) based on which the autocovariate is to 
        #		be calculated for all sites in the lattice defined 
        #		by colnums and rownums.
	###colnum:	should state the number of columns in the 
        #		rectangular lattice. This assumes that site 
        #		indices in sites were numbered across then down, 
        #		i.e. row by row rather than column by column.
	###rownum:	should state the number of rows in the rectangular 
        #		lattice.

	###initialise variables
	max.site<-colnums*rownums
	spa<-numeric(length(sites))
		
	###check each site in turn, making sure that only valid 
        ###neighbours are considered
	for(ii in 1:length(sites)){
		s<-sites[ii]
		nbs<-NULL	
        ###empty this vector before each round

		###if statements below ensure that neighbours which 
        ###don't exist (i.e. fall outside the rectangular lattice)
        ###are not included. Not that this buffers sites at the 
        ###top or bottom image of the lattice, but could lead to 
        ###a wrap-around (from left to right edge and vice 
        ###versa) at lateral edges. So make sure that sites of 
        ###interest are buffered at these edges, with extra 
		###columns of NA if necessary. The alogis program does 
        ##this for you.	
		nb<-s-colnums-1	#top left neighbour
		if(nb<=max.site && nb>=1){nbs<-c(nbs, nb)}
		nb<-s-colnums		#top neighbour
		if(nb<=max.site && nb>=1){nbs<-c(nbs, nb)}
		nb<-s-colnums+1	#top right neighbour
		if(nb<=max.site && nb>=1){nbs<-c(nbs, nb)}
		nb<-s-1			#left neighbour
		if(nb<=max.site && nb>=1){nbs<-c(nbs, nb)}
		nb<-s+1			#right neighbour
		if(nb<=max.site && nb>=1){nbs<-c(nbs, nb)}
		nb<-s+colnums-1	#bottom left neighbour
		if(nb<=max.site && nb>=1){nbs<-c(nbs, nb)}
		nb<-s+colnums		#bottom neighbour
		if(nb<=max.site && nb>=1){nbs<-c(nbs, nb)}
		nb<-s+colnums+1	#bottom right neighbour
		if(nb<=max.site && nb>=1){nbs<-c(nbs, nb)}
		#print(nbs)
		
		valid.nbs<-sum(!is.na(probs[nbs]))
        ###number of valid neighbours (included and not NA)
		#print(valid.nbs)
		sum.nbs<-sum(probs[nbs], na.rm=T)
		#print(sum.nbs)
		if(sum.nbs==0){
			###if this pixel is completely isolated and has no 
        ###valid neighbours:
			spa[ii]<-0
		}
		if(sum.nbs>0){
			###if the pixel has any neighbours:
			spa[ii]<-sum.nbs/valid.nbs
		}
	}
	auto<-spa
}  ###END OF NEIGHBOUR