Amakihi in Distance 1.0.3 vs DistWin 7.3

[erex]

Fitting 3 models to Amakihi in both R and Distwin:
Models hr(OBs+MAS), hn(OBs), hn(OBs+scale(MAS)) note that curiously, hn(OBs+MAS) fails

There are some degrees of freedom discrepancies, but note also that point estimates of density have best agreement for hn(OBs+scale(MAS) and least agreement for hn(OBs), intermediate for hr(OBs+MAS)

library(Distance)
data("amakihi")
amakihi$Area <- 1
amakihi$OBs <- relevel(amakihi$OBs, ref="TKP")
trunc <- 82.5
conv <- convert_units("meter", NULL, "hectare")
aukhr <- ds(data=amakihi, transect = "point", convert.units = conv,
 key="hr", formula=~OBs+MAS, truncation = trunc)
aukhn1 <- ds(data=amakihi, transect = "point", convert.units = conv,
 key="hn", formula=~OBs, truncation = trunc)
aukhn2 <- ds(data=amakihi, transect = "point", convert.units = conv,
 key="hn", formula=~OBs+scale(MAS), truncation = trunc)
# should fail
aukhn2 <- ds(data=amakihi, transect = "point", convert.units = conv,
 key="hn", formula=~OBs+MAS, truncation = trunc)
#

distwinhr <- read.table(text=
"
July 92 7.7366 9.58 173.17 6.4060 9.3436 
Dec 92 6.0996 12.11 110.67 4.8020 7.7479 
Apr 93 6.9409 7.85 117.23 5.9432 8.1060 
July 93 8.9561 9.06 157.23 7.4913 10.707 
Jan 94 4.9471 9.49 100.15 4.1000 5.9693 
Apr 94 7.5291 8.22 130.38 6.4013 8.8557 
Apr 95 7.9978 11.25 79.71 6.3984 9.9970
", stringsAsFactors=FALSE)
names(distwinhr) <- c("month", "year", "density", "cv", "df", "lcl", "ucl")
distwinhr$survey <- paste0(distwinhr$month, distwinhr$year)
rpack103hr <- aukhr$dht$individuals$D
rpack103hr <- rpack103hr[-8,] # remove total row
rpack103hr <- rpack103hr[c(6,1,3,7,2,4,5), ] # chronological ordering
plot(rpack103hr$Estimate,
 ylim=range(c(distwinhr$density, rpack103hr$Estimate)), xaxt="n",
 xlab="Survey", main="Amakihi hr(OBs+MAS)\nDistwin vs Distance 1.0.3")
axis(1, at=1:7, labels=distwinhr$survey)
points(distwinhr$density, pch=20)
legend("bottomright", legend = c("R", "DistWin"), pch=c(1,20))
pcthr <- data.frame(survey=distwinhr$survey,
 percenthr=round(rpack103hr$Estimate / distwinhr$density, 4))
pcthr

comphr <- data.frame(dfR=rpack103hr$df, dfD=distwinhr$df,
 lclR=rpack103hr$lcl, lclD=distwinhr$lcl,
 cvR=rpack103hr$cv, cvD=distwinhr$cv/100)
comphr


#================================= Only OBs covar, MAS fails
distwinhn <- read.table(text=
"
July 92 7.8171 9.27 168.17 6.5118 9.3840 
Dec 92 6.6367 11.64 99.29 5.2726 8.3537 
Apr 93 9.3903 8.09 128.31 8.0039 11.017 
July 93 9.4015 8.60 138.88 7.9338 11.141 
Jan 94 6.9919 9.81 110.57 5.7593 8.4884 
Apr 94 9.7917 8.52 136.87 8.2761 11.585 
Apr 95 9.4284 11.07 75.41 7.5674 11.747 
", stringsAsFactors=FALSE)
names(distwinhn) <- c("month", "year", "density", "cv", "df", "lcl", "ucl")
distwinhn$survey <- paste0(distwinhn$month, distwinhn$year)
rpack103hn <- aukhn1$dht$individuals$D
rpack103hn <- rpack103hn[-8,] # remove total row
rpack103hn <- rpack103hn[c(6,1,3,7,2,4,5), ] # chronological ordering
plot(rpack103hn$Estimate, 
 ylim=range(c(distwinhn$density, rpack103hn$Estimate)), xaxt="n",
 xlab="Survey", main="Amakihi hn(OBs)\nDistwin vs Distance 1.0.3")
axis(1, at=1:7, labels=distwinhn$survey)
points(distwinhn$density, pch=20)
legend("bottomright", legend = c("R", "DistWin"), pch=c(1,20))
pcthn <- data.frame(survey=distwinhn$survey,
 percenthr=round(rpack103hn$Estimate / distwinhn$density, 4))
pcthn

comphn <- data.frame(dfR=rpack103hn$df, dfD=distwinhn$df,
 lclR=rpack103hn$lcl, lclD=distwinhn$lcl,
 cvR=rpack103hn$cv, cvD=distwinhn$cv/100)
comphn


#================================= OBs+scale(MAS) covar
distwinhn <- read.table(text=
"
July 92 6.0353 9.09 164.77 5.0452 7.2197 
Dec 92 4.9281 11.48 95.51 3.9268 6.1847 
Apr 93 7.1205 7.88 118.53 6.0938 8.3203 
July 93 7.5425 8.50 134.89 6.3768 8.9213 
Jan 94 4.9985 9.50 100.53 4.1416 6.0326 
Apr 94 7.8371 8.28 131.74 6.6552 9.2290 
Apr 95 7.2941 10.93 72.03 5.8696 9.0643 
", stringsAsFactors=FALSE)
names(distwinhn) <- c("month", "year", "density", "cv", "df", "lcl", "ucl")
distwinhn$survey <- paste0(distwinhn$month, distwinhn$year)
rpack103hn <- aukhn2$dht$individuals$D
rpack103hn <- rpack103hn[-8,] # remove total row
rpack103hn <- rpack103hn[c(6,1,3,7,2,4,5), ] # chronological ordering
plot(rpack103hn$Estimate, 
 ylim=range(c(distwinhn$density, rpack103hn$Estimate)), xaxt="n",
 xlab="Survey", main="Amakihi hn(OBs+scale(MAS))\nDistwin vs Distance 1.0.3")
axis(1, at=1:7, labels=distwinhn$survey)
points(distwinhn$density, pch=20)
legend("bottomright", legend = c("R", "DistWin"), pch=c(1,20))
pcthn <- data.frame(survey=distwinhn$survey,
 percenthr=round(rpack103hn$Estimate / distwinhn$density, 4))
pcthn

comphn <- data.frame(dfR=rpack103hn$df, dfD=distwinhn$df,
 lclR=rpack103hn$lcl, lclD=distwinhn$lcl,
 cvR=rpack103hn$cv, cvD=distwinhn$cv/100)
comphn

[dill]

@erex could you post or e-mail me the Distance project file that you used to obtain the above results?

[erex]

This was the bog standard amakihi file that shipped with DistWin in Sample projects folder, I believe. But here's the version I just zipped up using DistWin.
forDave.zip

I "think" the labelling of the models described in this issue are:

R code name	DistWin label
auk	g7
aukhn1	l12
aukhn2	j10

[dill]

Thanks. I will look into this.

[dill]

In my copy of the default project (shipped with Distance 7.1 I think), analyses g7, l12 and j10 all use a model definition where the detection function estimation is stratified (i.e., one detection function fitted per stratum) so the results above are not comparable as in Distance we are only fitting one model.

This is not clear from the names of analyses, only once you look into the definitions (other analyses do include "by strat" in thier name, thus making one think that this is a naming convention -- it is not).

Looking at comparable analyses (letting one detection function be fitted then post-stratifying by year) there is still a discrepency between the results. I will look into this.

[dill]

Secondary discrepancy between the two analyses is that the Distance for Windows analysis includes 3 adjustment terms in it's analysis (along with the covariates) for l12. Removing these from the analysis gives an average p of 0.3462 in Distance for Windows and 0.3462969 in R (rounding as in presented results).

[erex]

Dave--I apologise for not comparing apples with apples.

I've revisited the amakihi and constructed the same detection function model both in DistWin and R: hazard rate with OBs factor covariate, with no DistWin adjustment terms and pooled detection function post-stratified to produce survey-specific estimates. Here is the comparison of the two pieces of software:

library(Distance)
trunc <- 82.5
data(amakihi)
cu <- convert_units("meter", NULL, "hectare")
amakihi$OBs <- relevel(amakihi$OBs, ref="TKP")
ERHROBs <- ds(data=amakihi, key="hr", truncation = trunc, transect = "point",
 formula=~OBs, convert.units=cu)
summary(ERHROBs)

WERHROBs <- read.table(text=
" 
July 92 5.5227 6.07 60.85 4.8918 6.2349 
Dec 92 5.4281 9.30 46.16 4.5031 6.5432 
Apr 93 8.7380 6.41 58.03 7.6867 9.9331 
July 93 7.7933 6.50 55.94 6.8424 8.8763 
Jan 94 6.5062 7.82 50.97 5.5619 7.6109 
Apr 94 9.1115 5.32 40.50 8.1833 10.145 
Apr 95 8.2198 9.69 45.53 6.7665 9.9852 
", stringsAsFactors=FALSE)
names(WERHROBs) <- c("month", "year", "density", "cv", "df", "lcl", "ucl")
WERHROBs$survey <- paste0(WERHROBs$month, WERHROBs$year)

den.est <- ERHROBs$dht$individuals$D
den.est <- den.est[-8,] # remove total row
den.est <- den.est[c(6,1,3,7,2,4,5), ] # chronological ordering

plot(den.est$Estimate, 
 ylim=range(c(WERHROBs$density, den.est$Estimate)), xaxt="n",
 xlab="Survey", main="Amakihi hr(OBs, no adjustments pooled detfn)\nDistwin vs Distance 1.0.3")
axis(1, at=1:7, labels=WERHROBs$survey)
points(WERHROBs$density, pch=20)
legend("bottomright", legend = c("R", "DistWin"), pch=c(1,20))

pcthn <- data.frame(survey=WERHROBs$survey,
 percenthr=round(den.est$Estimate / WERHROBs$density, 4))
print(pcthn) # do estimates differ by a constant fraction?

comphn <- data.frame(dfR=den.est$df, dfD=WERHROBs$df,
 lclR=den.est$lcl, lclD=WERHROBs$lcl,
 cvR=round(den.est$cv, 4), cvD=WERHROBs$cv/100)
print(comphn) # is there agreement between software in degrees of freedom or precision

> print(pcthn) # do estimates differ by a constant fraction?
 survey percenthr
1 July92 1.3547
2 Dec92 1.1200
3 Apr93 0.7734
4 July93 1.0818
5 Jan94 0.7734
6 Apr94 0.7734
7 Apr95 0.9429
> print(comphn) # is there agreement between software in degrees of freedom or precision
 dfR dfD lclR lclD cvR cvD
1 561.6757 60.85 5.569051 4.8918 0.1512 0.0607
2 164.9887 46.16 4.771005 4.5031 0.1232 0.0930
3 167.9178 58.03 5.722930 7.6867 0.0844 0.0641
4 167.2638 55.94 6.710579 6.8424 0.1160 0.0650
5 111.7762 50.97 4.165943 5.5619 0.0955 0.0782
6 163.4112 40.50 6.061277 8.1833 0.0764 0.0532
7 140.4718 45.53 6.415264 6.7665 0.0959 0.0969

[dill]

Thanks for this update.

Here there are differences in the estimated detection function parameters and resultant estimated average detectability between Distance for Windows and Distance (unlike in the half-normal case), but here we get an AIC of 10779.447 in Distance for Windows, whereas the above analysis gives an AIC of 10778.448 for Distance -- it seems that R gives us a "better" result than Distance for Windows.

[dill]

In my previous look at this it seemed that even when the estimated model matched (as in the half-normal case), there was a discrepency between estimated density. This might be because of a difference in the underlying data between the Distance for Windows project and the bundled R data.frame in flatfile format (though I didn't see a difference in number of transects visited per year between the two datasets).

[erex]

Dave Look also at the precision estimates and degrees of freedom associated with the survey-specific estimates. I don't think differences in the detection function parameters will have knock-on consequences for precision calculations.

[erex]

Dave--I double checked the data in the DistWin project and in the R package. There are the same number of points per stratum (year) as well as the same number of detections per stratum (year) when truncation is set to 82.5m.

Consequently, the stratum-specific encounter rates (and their CVs) match between software. I'm willing to assume the data being analysed are the same.