#read in data set for presence/absence of two plants
presabs.plants=read.table("http://www.public.iastate.edu/~pcaragea/S40608/Compute/plants.txt",header=TRUE)
attach(presabs.plants)


#start by defining neighbors
library(spatstat)#need this to compute distances between all points in a data

dist.neighbors=function(data,distance)
#distance is set by user
{
#data has the x and y on the first two columns. Could have more than two columns, or exactly two columns, it doesn't matter.
#produces a "list" object with a vector of the neighbors for each location
#distance defines the range within which we define neighbors
nd=length(data[,1])
nbmat=list(array)
dist=pairdist(data[,1:2])
  for(i in 1:nd){
  nbmat[[i]]=0
 	for(j in 1:nd) if((dist[i,j]<distance)&(i!=j)) nbmat[[i]]=c(nbmat[[i]],j)
 	nbmat[[i]]=nbmat[[i]][nbmat[[i]]!=0]}
return(nbmat)		
}


#example of running the function:
nb.plants=dist.neighbors(presabs.plants,13)

####################################################
#define canonical parameter for a constant mean model
canon.const = function(kappa, eta, ys, nbrs){
	Ais = NULL
	cnt = 0
	repeat{
		cnt = cnt+1
		index = nbrs[[cnt]]
		tnbrs = ys[index]
		Ais[cnt] = log(kappa)-log(1-kappa)+eta*sum(tnbrs-kappa)
		if(cnt==length(ys)){break}
	}
	return(Ais)
}

#compute (-1)*(log pseudo likelihood) from cannonical parameters
logplik = function(pars, ys, nbrs){
  kappa = pars[1]
  eta = pars[2]

  Ai = canon.const(kappa, eta, ys, nbrs)
  Bi = log(1+exp(Ai))
  logplikelihood = sum(ys*Ai-Bi)
  return(-1*logplikelihood)
}

#example of using the two functions above:
#say you wanted to compute the log pseudolikelihood function for some arbitrary values of kappa(=0.5) and eta(=0.8)
logplik(pars=c(0.5,0.8),ys=S68,nbrs=nb.plants)

#since you want to estimate parameters by minimizing the negative pseudo likelihood function, you should use one of the many numerical optimization routines available in R. I will show you how to use "optim"
#under "par" you need to give the starting values for the two parameters kappa and eta. Make sure these starting values are within the domain (i.e. kappa is between 0 and 1, there are no restrictions on eta, really.)
#fn is the function to be minimized
#ys and nbrs are the arguments needed for logplik
#hessian=TRUE returns the observed information matrix, which allows computation of standard errors assuming independence
Model1=optim(par=c(0.3,0.5), fn=logplik, ys=S68, nbrs=nb.plants,hessian=TRUE)
#to get estimated parameters:
Model1$par
#to get their std.errors:
sqrt(diag(solve(Model1$hessian)))


#finally, we can compute fitted values:
const.fitted = function(pars, ys, nbrs){ 
#produce predictions and mspe
  kappa = pars[1]
  eta = pars[2]
  Ai = canon.const(kappa, eta, ys, nbrs)
  ps = exp(Ai)/(1+exp(Ai))
  mspe = mean((ys-ps)^2)
  list(mspe=mspe,preds=ps)
}

#here is an example of usage for this function:
Fits=const.fitted(Model1$par,S68,nb.plants)
Fits$preds #returns the actual predictions at each location
Fits$mspe  #returns the mean squared prediction errors



################### Model with covariates. 
#Here, only one covariate, called "covs"

canon.covar = function(pars, ys,covs, nbrs){
	beta0=pars[1]
	beta1=pars[2]
	eta=pars[3]
	Ais = NULL
	cnt = 0
	repeat{
		cnt = cnt+1
		index = nbrs[[cnt]]
		tnbrs = ys[index]
		kappajs = beta0+beta1*covs[index]
        Ais[cnt] = beta0+beta1*covs[cnt]+ eta*sum(tnbrs-kappajs)
		if(cnt==length(ys)){break}
	}
	return(Ais)
}

#and the log-likelihood function:
logplikcovar = function(pars, ys, covs, nbrs){

   Ais=canon.covar(pars,ys,covs,nbrs)
   Bis = log(1+exp(Ais))

  logpliklihood = sum(ys*Ais-Bis)
  return(-1*logpliklihood)
}


#esimate model parameters:
Model2=optim(par=c(0.5,-1,0.5),fn=logplikcovar,ys=S68,covs=H68,nbrs=nb.plants,hessian=TRUE)
#to get estimated parameters:
Model2$par
#to get their std.errors:
sqrt(diag(solve(Model2$hessian)))



#finally, we can compute fitted values:
covar.fitted = function(pars, ys, covs, nbrs){ 
#produce predictions and mspe
  Ai = concanon.covar(pars, ys, covs, nbrs)
  ps = exp(Ai)/(1+exp(Ai))
  mspe = mean((ys-ps)^2)
  list(mspe=mspe,preds=ps)
}

Fits.covar=covar.fitted(Model2$par,S68,H68,nb.plants)
Fits.covar$preds #returns the actual predictions at each location
Fits.covar$mspe  #returns the mean squared prediction errors


####################################
#directional models: NS and EW neighbors
#allows for EW and NS directional dependence. Constant mean
canon.dir = function(pars, ys, nbrs){
	kappa=pars[1]
	eta1=pars[2]
	eta2=pars[3]
	
	Ais = NULL
	cnt = 0
	repeat{
		cnt = cnt+1
		yis=ys[cnt]
		index = nbrs[[cnt]]
		tnbrs = ys[index]
  		   tnbs1<-tnbrs[abs(tnbrs-yis)<=1]#EW
		   tnbs2<-tnbrs[abs(tnbrs-yis)>1]#NS
   		   ts1<-sum(ys[tnbs1])
           ts2<-sum(ys[tnbs2])

      		Ais<-log(kappa)-log(1-kappa)+eta1*(ts1-kappa)+eta2*(ts2-kappa)
		if(cnt==length(ys)){break}
	}
	return(Ais)
}

#and the log-likelihood function:
logplikdir = function(pars, ys, nbrs){

   Ais=canon.dir(pars,ys,nbrs)
   Bis = log(1+exp(Ais))

  logpliklihood = sum(ys*Ais-Bis)
  return(-1*logpliklihood)
}

Model3=optim(par=c(0.5,0.5,0.5),fn=logplikdir,ys=S68,nbrs=nb.plants,hessian=TRUE)
#to get estimated parameters:
Model3$par


#finally, we can compute fitted values:
dir.fitted = function(pars, ys, nbrs){ 
#produce predictions and mspe
  Ai = canon.dir(pars, ys, nbrs)
  ps = exp(Ai)/(1+exp(Ai))
  mspe = mean((ys-ps)^2)
  list(mspe=mspe,preds=ps)
}

Fits.dir=dir.fitted(Model3$par,S68,nb.plants)
Fits.dir$preds #returns the actual predictions at each location
Fits.dir$mspe  #returns the mean squared prediction errors