Antonio Olinto wrote:
> Hello,
>
> I'm analyzing some fish length-frequency data to access relative age and
> growth
> information and I would like to do something different from FiSAT / ELEFAN
> analysis.
>
> I found a package named MIX in http://www.math.mcmaster.ca/peter/mix/mix.html
> but it's compiled for R 2.0.0
>
> So I have two questions:
> Is there any problem to run this package in R 2.7.0? - probably yes …
> And, is there any other package in R that can be used to analyze and to
> separate
> mixture distributions?
>
> Thanks for any help. Best regards,
>
> Antonio Olinto
> Sao Paulo Fisheries Institute
> Brazil
This problem is not too difficult. Maybe you can write your own code.
Say, you enter the number of fish you measured, n, and a two-column
dataframe with the mid point of the length class in the first column
(call it l) and the observed relative frequency of each length class in
the second column (call it obs_prob). Then using the multinomial
likelihood for the number of fish in each length class as the random
variable (the approach in Peter Macdonald's Mix), minimize
log_likelihood=-n*sum(elementwise product(obs_prob,log(pred_prob)))
where
pred_prob=(p1*0.3989423/s1)*exp(-0.5*square((l-m1)/s1))
+(p2*0.3989423/s2)*exp(-0.5*square((l-m2)/s2))
+(p3*0.3989423/s3)*exp(-0.5*square((l-m3)/s3))
where s1, s2, s3, m1, m2, m3, p1 and p2 (p3=1-p1+p2) are the parameters.
Do you know your number of components (cohorts) in the mixture?
In the model above it is three. If you don't know it, try several models
with different number of components and the AIC may let you know which
one is the best working model. Don't use the likelihood ratio test as
one of the p parameters is on the boundary of parameter space if the
null is true.
Also, you should bound the parameters (m1<m2<m3, 0<p1<1, 0<p2<1, and
establish the condition p3=1-p1+p2.
I wrote ADMB code to do this if you want it. You can translate it to R.
Below you can find some example code to do the graphics.
R.
cont <-
c(4,15,32,44,62,69,61,99,114,119,106,108,89,88,95,99,91,84,137,190,224,297,368,484,566,629,446,349,377,405,438,367,215,115,36,10,4,1,0,0,1)
tot <- sum(cont)
rel.freq <- rep(cont,each=10)/tot/10
lect <- seq(9.1,50,0.1)
prop1 <- 0.0219271
prop2 <- 0.121317
prop3 <- 0.0328074
prop4 <- 0.315534
prop5 <- 0.203071
prop6 <- 0.3053439
sigma1 <- 1.50760
sigma2 <- 2.82352
sigma3 <- 1.40698
sigma4 <- 3.00063
sigma5 <- 1.41955
sigma6 <- 1.99940
media1 <- 13.4148
media2 <- 19.2206
media3 <- 24.5748
media4 <- 31.9498
media5 <- 34.6998
media6 <- 39.7016
prot1<-(1/10)*(prop1*(1/sqrt(2*pi))/sigma1)*exp(-0.5*((lect-(media1-0.5))/sigma1)^2)
prot2<-(1/10)*(prop2*(1/sqrt(2*pi))/sigma2)*exp(-0.5*((lect-(media2-0.5))/sigma2)^2)
prot3<-(1/10)*(prop3*(1/sqrt(2*pi))/sigma3)*exp(-0.5*((lect-(media3-0.5))/sigma3)^2)
prot4<-(1/10)*(prop4*(1/sqrt(2*pi))/sigma4)*exp(-0.5*((lect-(media4-0.5))/sigma4)^2)
prot5<-(1/10)*(prop5*(1/sqrt(2*pi))/sigma5)*exp(-0.5*((lect-(media5-0.5))/sigma5)^2)
prot6<-(1/10)*(prop6*(1/sqrt(2*pi))/sigma6)*exp(-0.5*((lect-(media6-0.5))/sigma6)^2)
prot.tot<-prot1+prot2+prot3+prot4+prot5+prot6
plot(lect,rel.freq,type="l",lwd=3,xlab="",ylab="",ylim=c(0,0.01))
lines(lect,prot1)
lines(lect,prot2)
lines(lect,prot3)
lines(lect,prot4)
lines(lect,prot5)
lines(lect,prot6)
lines(lect,prot.tot,lwd=2)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
