[R-sig-eco] Help with distribution fitting and AIC in R

[Peter Solymos]

On Sat, Aug 6, 2011 at 4:06 AM, Lene Jung <ljkjar at hotmail.com> wrote:
> HI,
>
> I’m having several problems trying to fit distributions to data that I have
> sorted into a data frame, so the each ID has its own step length and turn
> angle.
>
>
>
> I can fit a Weibull distribution to step lengths with following code:
>
>
>
> ###create data frame###
>
> ID1 <- na.omit(ID)
>
> weib.test = data.frame(ID1, step1)
>
> set.seed(144)
>
>
>
> weib.dist<-rweibull(10000,shape=3,scale=8)
>
> names(weib.test)<-c('ID','steplength')
>
>
>
> g <- function(step1) {
>
>     require(MASS)
>
>     est <- fitdistr(step1, 'weibull')$estimate
>
>                AIC(est)
>
>     list(shape = est[1], scale = est[2])
>
>
>
>  }
>
> library(data.table)
>
> weib.test.dt <- data.table(weib.test, key = 'ID')
>
> weib.test.dt[, g(steplength), by = ID1]
>
>
>
> This gives me the following output:
>
>
>
>  ID1     shape    scale
>
> [1,]   1 0.8089580 131.1514
>
> [2,]   2 0.8393755 106.2641
>
> [3,]   3 0.8421661 196.7409
>
> [4,]   4 0.9096443 192.8869
>
> [5,]   5 0.8170656 172.3526
>
> [6,]   6 0.9815219 136.5133
>
> [7,]   7 0.9485611 133.9226
>
> [8,]   8 0.9279664 158.5760
>
> [9,]   9 0.7980447 250.8815
>
> [10,]  10 0.7882116 132.3324
>
> [11,]  11 0.6684602 176.5392
>
> [12,]  12 0.7636310 105.1567
>
> [13,]  13 0.7179191 126.0132
>
> [14,]  14 0.6833156 183.6390
>
> [15,]  15 0.9356103 171.3479
>
> [16,]  16 0.7865950 129.0388
>
> [17,]  17 0.8408600 157.9560
>
> [18,]  18 0.8023891 127.5567
>
>
>
> Which is exactly what I need – separate shape and scale parameters for each
> ID. However I also want to get AIC values for each ID, and I have tried
> several things such as
>
> AIC(weib.test.dt, by = ID1)
>
> or
>
> extractAIC(fit, scale, k = 2, ...)AIC(weib.test.dt)
>
> stopifnot(all.equal(AIC(weib.test.dt),
>
>                    AIC(logLik(weib.test.dt))))
>
>
>
>
>
> But I keep getting an eror:  Error in UseMethod("logLik") :
>
>  no applicable method for 'logLik' applied to an object of class
> "c('data.table', 'data.frame')"
>
>
>
> I’m assuming this had to do with the multiple ID’s? Should I just fit the
> distributions separately for each ID instead and then run AIC?
>

Lene,

You need to return the fitted object 'est', the error message tells
you that data.frame object class does not come with a logLik method,
unlike a fitdistr class returned by the function fitdistr(). You can
then use AIC, which assumes there is a logLik method with
log-likelihood and number of parameters defined. Alternatively, you
can return c(est$estimate, AIC(est)) that gives you a data frame with
AIC in 3rd column.

Careful reading of the help page for fitdistr would tell you the
structure of the return object.

>
>
> Furthermore, I would like to do the same thing as above with turn angles –
> but by fitting a wrapped Cauchy distribution to the several IDs. I have
> tried following code:
>
>
>
> ###wrapped cauchy fitting, simple mean and rho###
>
> mu0<- circ.mean(turn1)
>
> rho0 <- est.rho(turn1)
>
>
>
> wcauchy.test = data.frame(ID2, turn1)
>
> set.seed(144)
>
>
>
> wcauchy.dist<-rweibull(10000,mu0,rho0)
>
> names(wcauchy.test)<-c('ID','turnangle')
>
>
>
> z <- function(turn1) {
>
>     require(MASS)
>
>     est <-wrpcauchy.ml(turn1, mu0, rho0, acc=1e-015)$estimate
>
>     list(mu0 = est[1], rho0 = est[2])
>
>  }
>
>
>
> library(data.table)
>
> wcauchy.test.dt <- data.table(wcauchy.test, key = 'ID')
>
>
>
> wcauchy.test.dt[, z(turnangle), by = ID2]
>
>
>
> But I keep getting the error: Error in `[.data.table`(wcauchy.test.dt, ,
> z(turnangle), by = ID2) :
>
>  Type 0 in j column
>
>
>
> And of course I would like to calculate the AIC as well for distribution
> fit.
>
>
>
> Can anyone help me out with this?
>

There is no wrpcauchy.ml function in MASS, but in CircStats. It
returns only the estimates, and no logLik method is defined. From the
code or from the references, it should be possible to figure out PDF
for the wrapped Cauchy distribution, and from that it is easy to
calculate AIC.

Again, this is obvious from the help page.

Cheers,

Peter

>
>
> Thanks,
>
> Lene Jung Kjaer
>
>
>
> ___________________________
>
> Lene Jung Kjær
>
> Studsdalvej 20, Taulov
>
> 7000 Fredericia, Danmark
>
> Tlf: 29 86 96 14
>
> email: ljkjar at hotmail.com
>
>
>
>
>        [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>