Update to the R code found in the appendix of:

Brunner, J.L. K. LoGiudice, and R.S. Ostfeld. 2008. Estimating reservoir competence of Borrelia burgdorferi hosts: prevalence and infectivity, sensitivity and specificity. Journal of Medical Entomology 45(1): 139-147. [PDF]

The anneal() function (found in the Likelihood package here) was changed recently and now requires that the parameter names be specified in a separate list called "var". If these names are instead specified in the "par" list, as the original code had it, you will likely get an error such as:

Error in anneal(model = EstPrevInf, par = pars, source_data = ExampleData, : anneal: All values in par must be numeric.

The parts of the code that have been changed are in bold.

EstPrevInf= function(P,I,Pos,Tested) { LL=log(P* dbinom(Pos, Tested, prob=I) +(1-P) * dbinom(Pos, Tested, prob=0)) #to prevent the function from returning -inf, which #causes errors in the optimization routine, those rows #with -inf are replaced with a very low value LL[is.infinite(LL)] = -99999 return(LL) }
PDF.wrapper=function(predicted) { #just to keep anneal() happy return(predicted) }
pars= list(P=0.5, I=0.5) vars = list(Pos= "Pos", Tested= "Tested", predicted="predicted") par_lo= list(P=0, I=0) par_hi= list(P=1, I=1)
ExampleData = data.frame(cbind( Tested = c(3, 10, 25, 16, 15, 12, 7, 7, 3, 3, 2), Pos =c(3, 0, 13, 13, 8, 10, 7 , 5, 0, 0, 0) ))
Results = anneal(model=EstPrevInf,par=pars, var = vars, source_data=ExampleData, pdf=PDF.wrapper, par_lo, par_hi, dep_var="Pos", initial_temp=1) #lower temp converges better

And simililarly, in the case with multiple host species:

EstPrevInfSpec = function(T1, P1, T2, P2, T3, P3, Spec, Pos, Tested, Spp) { LL = numeric(length(Pos)) #empty vector of the right length #First calculate the likelihoods of the subset of #individuals that are species "One" LL[which(Spp == "One")] = log( #Positive tests from vectors fed on infected animals P1*dbinom(Pos[which(Spp == "One")], Tested[which(Spp == "One")], prob = (T1 + (1 - T1)*(1 - Spec)) ) #Positives tests from vectors fed on uninfected animals + (1-P1)*dbinom(Pos[which(Spp == "One")], Tested[which(Spp == "One")], prob = (1-Spec) ) ) #Then for those that are species "Two" LL[which(Spp == "Two")] = log( P2*dbinom(Pos[which(Spp == "Two")], Tested[which(Spp == "Two")], prob = (T2 + (1 - T2)*(1 - Spec)) ) + (1-P2)*dbinom(Pos[which(Spp == "Two")], Tested[which(Spp == "Two")], prob = (1-Spec) ) ) #Then finally those that are species "Three" LL[which(Spp == "Three")] = log( P3*dbinom(Pos[which(Spp == "Three")], Tested[which(Spp == "Three")], prob = (T3 + (1 - T3)*(1 - Spec)) ) + (1-P3)*dbinom(Pos[which(Spp == "Three")], Tested[which(Spp == "Three")], prob = (1-Spec) ) ) #to prevent the function from returning -inf, #which causes errors in the optimization routine, those rows #with -inf are replaced with a very low value LL[is.na(LL)] = -99999 LL[is.infinite(LL)] = -99999 return(LL) }
pars = list(T1 = 0.5, P1 = 0.5, T2 = 0.5, P2 = 0.5, T3 = 0.5, P3 = 0.5, Spec = 0.5) vars = list(Pos = "Pos", Tested = "Tested", Spp = "Spp", predicted = "predicted") par_lo = list(T1 = 0, P1 = 0, T2 = 0, P2 = 0, T3 = 0, P3 = 0, Spec = 0) par_hi = list(T1 = 1, P1 = 1, T2 = 1, P2 = 1, T3 = 1, P3 = 1, Spec = 1)
MultiSppData = data.frame( cbind( Tested = c(3, 10, 25, 16, 15, 12, 7, 7, 3, 3, 2, 15, 16, 22, 25, 18, 18, 16, 21, 23, 5, 9, 10, 12, 13, 12, 8, 9, 17, 13, 6, 5, 16, 3, 1, 1), Pos = c(3, 0, 13, 13, 8, 10, 7, 5, 0, 0, 1, 5, 1, 6, 0, 2, 4, 3, 0, 0, 0, 1, 2, 1, 1, 7, 6, 6, 7, 4, 3, 4, 7, 2, 0, 1) ) )
MultiSppData$Spp = c(rep("One", 11), rep("Two", 12), rep("Three", 13) )
MultiSppResults = anneal(model = EstPrevInfSpec, par = pars, var = vars, source_data = MultiSppData, pdf = PDF.wrapper, par_lo, par_hi, dep_var = "Pos", initial_temp = 1)


Last updated February 2008