Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
| Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior | ||
|
artigos:ernesto3:simpj [2008/10/08 22:14] paulojus |
artigos:ernesto3:simpj [2008/10/16 00:50] (atual) paulojus |
||
|---|---|---|---|
| Linha 1: | Linha 1: | ||
| ====== A joint model proposal???? ====== | ====== A joint model proposal???? ====== | ||
| + | - {{artigos:ernesto3:joint.rnw|Arquivo Sweave com o modelo proposto}} | ||
| ===== Ideias for simulation and model ===== | ===== Ideias for simulation and model ===== | ||
| Linha 67: | Linha 68: | ||
| Uma outra possibilidade é dividir as probs p1 e p2 por 2 e com isto garantir qq beta ser válido -- ver os efeitos disto pois coloca um limite de 0.5 em p1 e p2 o que pode nao ser uma boa ideia... | Uma outra possibilidade é dividir as probs p1 e p2 por 2 e com isto garantir qq beta ser válido -- ver os efeitos disto pois coloca um limite de 0.5 em p1 e p2 o que pode nao ser uma boa ideia... | ||
| - | Em uma ou outra situação ainda nao está claro como gerar uma situação onde uma das classes possa ter quase 100% | + | Em uma ou outra situação ainda nao está claro como gerar uma situação onde uma das classes possa ter quase 100% -- talvez pensar em algo como multiplicar as proporcoes ou na escala LRT |
| + | |||
| + | //**pensando melhor...**// | ||
| + | pensando melhor no que escrevi antes acho que me empolguei meio rápido demais: | ||
| + | - p3 precisa ser calculado associado aos betas tb caso contrário se for feito por diferenca como propus inicialmente composicoes (p1 e p2) que ficam restritas a serem menores que 0.5. Com isto acho que o modelo válido | ||
| + | - me parece ainda que uma outra idéia para gerar associação estava bem na nossa frente e nao percebemos: usar o modelo multinomial com a abundância como covárivel!!! - no paper isto pode ter certa vantagem por unificar com o diagnóstico que utilizamos. Seria bom fazer umas simu;ação [para ver as situações geradas. | ||
| + | |||
| + | ===== Alternativa usando CDA ===== | ||
| + | |||
| + | <code R> | ||
| + | require(lattice) | ||
| + | require(MASS) | ||
| + | require(geoR) | ||
| + | |||
| + | gs <- expand.grid((0:10)/10, (0:10)/10) | ||
| + | |||
| + | # basic gaussians to be used posteriorly | ||
| + | s2 <- 0.5 | ||
| + | Sig <- diag(c(1,1,1,1,1)) | ||
| + | Sig[1,4] <- 0.9 | ||
| + | Sig[4,1] <- 0.9 | ||
| + | Sig <- Sig*s2 | ||
| + | Sig | ||
| + | m0 <- mvrnorm(nrow(gs), c(2,0,0,0,0), Sig) | ||
| + | plot(as.data.frame(m0)) | ||
| + | cor(m0) | ||
| + | |||
| + | # Generate a spatial Gaussian process Y | ||
| + | phi <- 0.25 | ||
| + | sigmasq <- 1 | ||
| + | Y <- grf(grid=gs, cov.pars=c(sigmasq, phi)) | ||
| + | Y$data <- exp(m0[,1]+Y$data) | ||
| + | var(log(Y$data)) | ||
| + | |||
| + | ## option 1: age compositions independent from Y | ||
| + | # Generate age compositions independent from Y | ||
| + | comp.1 <- t(apply(cbind(exp(m0[,c(2,3)]),1),1,function(x) x/sum(x))) | ||
| + | dim(comp.1) | ||
| + | apply(comp.1, 2, range) | ||
| + | cor(cbind(comp.1, log(Y$data))) | ||
| + | plot(as.data.frame(cbind(comp.1, log(Y$data)))) | ||
| + | |||
| + | lr.1 <- data.frame(lr13 = log(comp.1[,1]/comp.1[,3]), lr23 = log(comp.1[,2]/comp.1[,3]), Y=log(Y$data)) | ||
| + | cor(lr.1) | ||
| + | cor(lr.1, met="spea") | ||
| + | plot(lr.1) | ||
| + | |||
| + | |||
| + | # Build abundance at age | ||
| + | Yi.1 <- comp.1*Y$data | ||
| + | dim(Yi.1) | ||
| + | cor(cbind(Yi.1, log(Y$data))) | ||
| + | plot(as.data.frame(cbind(Yi.1, log(Y$data)))) | ||
| + | |||
| + | ## option 2: age compositions dependent from Y | ||
| + | # Generate age compositions dependent from Y | ||
| + | comp.2 <- t(apply(cbind(exp(m0[,c(3,4)]),1),1,function(x) x/sum(x))) | ||
| + | apply(comp.1, 2, range) | ||
| + | cor(cbind(comp.2, log(Y$data))) | ||
| + | cor(cbind(comp.2, log(Y$data)), met="spea") | ||
| + | plot(as.data.frame(cbind(comp.2, log(Y$data)))) | ||
| + | |||
| + | lr.2 <- data.frame(lr13 = log(comp.2[,1]/comp.2[,3]), lr23 = log(comp.2[,2]/comp.2[,3]), Y=log(Y$data)) | ||
| + | cor(lr.2) | ||
| + | cor(lr.2, met="spea") | ||
| + | plot(lr.2) | ||
| + | |||
| + | lr.2a <- data.frame(lr12 = log(comp.2[,1]/comp.2[,2]), lr32 = log(comp.2[,3]/comp.2[,2]), Y=log(Y$data)) | ||
| + | cor(lr.2a) | ||
| + | cor(lr.2a, met="spea") | ||
| + | plot(lr.2a) | ||
| + | |||
| + | # Build abundance at age | ||
| + | Yi.2 <- comp.2*Y$data | ||
| + | dim(Yi.2) | ||
| + | cor(cbind(Yi.2, log(Y$data))) | ||
| + | plot(as.data.frame(cbind(Yi.2, log(Y$data)))) | ||
| + | |||
| + | df0 <- data.frame(abund=c(Yi.1,Yi.2), comp=c(comp.1,comp.2),Y=Y$data, opt=rep(c("indep","dep"), | ||
| + | rep(length(comp.1),2)), age=rep(rep(c(1,2,3),2),rep(length(comp.1)/3,6))) | ||
| + | |||
| + | # plots | ||
| + | xyplot(comp~log(Y)|age*opt,data=df0) | ||
| + | xyplot(abund~Y|age*opt,data=df0) | ||
| + | plot(as.geodata(Y)) | ||
| + | </code> | ||
| + | |||
| + | O problema que temos nestes casos em que simulamos Y e depois obtemos Yi por Y*pi é que a estrutura espacial dos Yi vai ser exactamente a mesma. Se acrescentarmos variabilidade com estrutura espacial temos que recalcular os pi ... | ||
| + | |||
| + | PJ: De fato mas para fazer diferente teria que simular de outra forma pois aqui só tem um Y geoestatistico gerado. O modelo multivariado conjunto para os LR seria a alternativa. Por outro lado eu nao vejo muito problema nesta fato pois as correlacoes entra as composicoes de certa forma tratam isto. | ||
| + | |||
| + | |||
