I am looking to optimize the fit of a model that describes the amount of litter collected in a network of .5m^2 'litter traps' in a plot of mapped trees of known diameter and species. The model of choice has two factors, allometric scaling of litter production, and exponential decay in litter travel distance.
tree1.litter = alpha*gamma^2 * DBH^Beta/(2*pi) * exp(-gamma*z-delta*DBH)
However, our trap data contains input from multiple trees (this is the "missing level" referred to in title):
Obs.Litter = tree1.litter + tree2.litter + ... + treej.litter + error
So far had very mixed results on even simulated data. It seems like with enough combinations of diameters and distances the functions should be somewhat well constrained. This analysis has been performed in an article I'm copy-catting. I've also tried the analysis on the log(Obs.Litter), which I think is the way to go. But I am not sure that the way I've coded the log version would have resulted in something that you would expect to perform any better.
At this point I suppose I'm just looking for any sort of advice (code based or conceptual) from someone more experienced with fitting nonlinear regressions or model fitting problems with this type of "hidden process". Code for data simulation and the various likelihoods are included below. I've had a bit more success with estimating these parameters with a Bayesian hierarchical model in OpenBUGS, with informative priors only.
library(plyr)
########################
##Generate Data#########
########################
alpha = 5
Beta = 2
gamma = .2
delta = .02
n = 600 #Number of trees
N.trap = 45 #Number of litter traps
D = rlnorm(n, 2)+5 #generate diameters
Z = runif(n, 0, 25) #generate distances
trap.id = sort(sample(1:N.trap, size = n, replace = T)) #assign trees to traps
tree.lit = (2*pi)^-1*alpha*gamma^2*D^Beta * exp(-gamma*Z-delta*D) #generate litter
log.tree.lit = -(2*pi) + log(alpha) + 2*log(gamma) + Beta*DBH -gamma*Z - delta*D
litter = data.frame(D=D, Z = Z, trap.id = trap.id, tree.lit = tree.lit)
data = ddply(litter, .(trap.id), summarize, trap.lit = sum(tree.lit), n.trees=length(trap.id) )
trap.lit = data[,2]
#####################################
##### Maximum Likelihood Optimization
#####################################
library(bbmle)
log.Litter.Func<-function(alpha, Beta, gamma, delta, sigma, D, Z, N.trap, trap.id, Obs.Litter){
log.Expected.Litter.tree = -log(2*pi) + log(alpha) + 2*log(gamma) + Beta*log(D) -gamma*Z - delta*D
log.Expected.Litter.trap = rep(0, N.trap)
for(i in 1:N.trap){
log.Expected.Litter.trap[i] <- sum(exp(log.Expected.Litter.tree[trap.id==i]))
}
-sum(dlnorm(log(Obs.Litter), log.Expected.Litter.trap, sd=sigma, log=T))
}
Litter.Func<-function(alpha, Beta, gamma, delta, sigma, D, Z, N.trap, trap.id, Obs.Litter){
Expected.Litter.tree = 1/(2*pi) * alpha * gamma^2 * D^Beta *exp(-gamma*Z - delta*D)
Expected.Litter.trap = rep(0, N.trap)
for(i in 1:N.trap){
Expected.Litter.trap[i] <- sum(Expected.Litter.tree[trap.id==i])
}
-sum(dnorm(Obs.Litter, Expected.Litter.trap, sd=sigma, log=T))
}
log.fit<-mle2(log.Litter.Func,
start = list(alpha = 5,gamma = .2,Beta = 2,delta = .02, sigma = 1),
#upper = list(alpha = 20,gamma = 1,Beta = 4,delta = .2,sigma = 20),
#lower = list(alpha = .002,gamma = .002,Beta = .0002,delta = .000000002,sigma = .020),
#method="L-BFGS-B",
data=list(D=D, Z=Z, N.trap=N.trap, trap.id=litter$trap.id, Obs.Litter=trap.lit)
)
fit<-mle2(Litter.Func,
start = list(alpha = 5,gamma = .2,Beta = 2,delta = .02, sigma = 1),
#upper = list(alpha = 20,gamma = 1,Beta = 4,delta = .2,sigma = 20),
#lower = list(alpha = .002,gamma = .002,Beta = .0002,delta = .000000002,sigma = .020),
#method="L-BFGS-B",
data=list(D = D,Z = Z,N.trap=N.trap, trap.id=litter$trap.id,Obs.Litter = trap.lit)
)
