Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

Link para esta página de comparações

Ambos lados da revisão anteriorRevisão anterior
Próxima revisão
Revisão anterior
artigos:ernesto3:simpj [2008/10/08 21:58] paulojusartigos: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 66: Linha 67:
 Comentário adicional: para fechar o modelo seria bom que os p's já fosse gerados de forma a somar 1 sem a necessidade de "padronizar". Imagino que deve ser possível impor uma restrição nos coeficientes <latex>\beta</latex> para garantir isto.\\ Comentário adicional: para fechar o modelo seria bom que os p's já fosse gerados de forma a somar 1 sem a necessidade de "padronizar". Imagino que deve ser possível impor uma restrição nos coeficientes <latex>\beta</latex> para garantir isto.\\
 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%  -- 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.
 +
  
  
  

QR Code
QR Code artigos:ernesto3:simpj (generated for current page)