jmsigner / amt

37 stars 13 forks source link

Understanding log_RSS #52

Closed anjankatna closed 2 years ago

anjankatna commented 2 years ago

Hi, I have a query on correctly interpreting the RSS values, as I am getting different results when using the naive estimate and the model vs the results from the log_RSS function. For simplicity, lets only consider the landcover class variable in my model, with the reference as agriculture. Model coefficient for fallow = -1.533 [ exp(coeff) = 0.215]. == Thus, more likely to select agri.

Using naive estimates, I weight it with the available locations for fallow and agri, using the following formula: 1/(exp(coef(m1.1$model)["landuseCFallow"])*a.flw/a.agri) which gives me a value of 5.09 (i.e. 5 times more likely to use agri???)

However, when I use the log_RSS function to do the same, with x1 as fallow and x2 as the reference class agriculture, log_RSS = 0.1417 exp(log_RSS) = 1.15 (so this says that the animal is 1.15 times more likely to use fallow?)

Or am using the function incorrectly. Please help!

bsmity13 commented 2 years ago

Hello,

Your explanations sound correct. I'd have to see your model structure and your code to figure out the discrepancy. If you have no interactions between land cover and any other terms, you should get the same result with the function vs. manually with the coefficients.

Can you share more details? At least the model formula and your x1 and x2 objects for the log_rss() function would be helpful.

anjankatna commented 2 years ago

Hi, Thank you for your reply. Please see the model below (this is a pilot analysis on a subset of my data which I am doing to get my ssf and amt understanding right, my final model may have some differences) - jungle cat data for 1 month, the reference landcover class is agriculture

m1.1 <- ssf_dat.jc.jan18 %>% fitclogit(case ~ landuseC + sl_ + costa + logsl + dist.poultry + dist.roads + dist.settle + landuseC:(sl_ + logsl + costa) + cluster(id) + strata(stepid), method = 'approximate', model = TRUE)

the coefficient for fallow is -1.53382

Here are the x1 and x2

-data.frame for s1 s1 <- data.frame( landuseC = factor("Fallow", levels = levels(ssfdat.jc.jan18$landuseC)), sl = 100, logsl = log(100), costa = 1, dist.poultry = 1, dist.roads = 1, dist.settle = 1)

-data.frame for s2 #reference location s2 <- data.frame( landuseC = factor("Agriculture", levels = levels(ssfdat.jc.jan18$landuseC)), sl = 100, logsl = log(100), costa = 1, dist.poultry = 1, dist.roads = 1, dist.settle = 1)

lr1 <- log_rss(m1.1, x1 = s1, x2 = s2, ci = "se")

The manual calculation values: Agri == avail = 13963, used = 2203 Fallow == avail = 12716, used = 2811

bsmity13 commented 2 years ago

Thanks for sharing! Yeah, the discrepancy is caused by the interaction between landuseC and the movement parameters. Even though sl_, log_sl_, and cos_ta_ are the same between s1 and s2, and thus the main effects of those parameters will cancel out, the interaction terms are not the same. I.e., there is a term for sl_:landuseCAgriculture and sl_:landuseCFallow (or perhaps agriculture is the reference level, but the idea holds), but those terms being non-zero is what is altering your results (which is correct, by the way).

If you want to calculate by hand and check, it should be something along the lines of:

1/(exp(coef(m1.1$model)["landuseCFallow"] + coef(m1.1$model)["landuseCFallow:sl_"] * 100 + coef(m1.1$model)["landuseCFallow:log_sl_"] * log(100) + coef(m1.1$model)["landuseCFallow:cos_ta_"] * 1) * a.flw/a.agri)

(I'm assuming agriculture is the reference level. Your names may be slightly different, but hopefully that gives you the idea.)

I hope that helps!

anjankatna commented 2 years ago

yes, thanku very much. This was very helpful. Just a minor confirmation - you meant the results of the log_RSS function are correct.

bsmity13 commented 2 years ago

Oh sorry -- yes, I would trust the results of the log_rss() function.

anjankatna commented 2 years ago

great. that solves my RSS related doubts. Thanks a tonne. I had a few other doubts with my ssf model, might ask them in a separate thread.