[R-br] Problema numérico no gnm

[Davi]

Boa noite a todos!

Estou tentando simular de um modelo não-linear generalizado e depois ajustar
esse modelo aos dados gerados na simulação. A variável resposta pode seguir
Poisson ou Binomial, não tem problema.

 

A média é modelada por: n*(1-exp(-(b0+b1*x))*(exp(-gx) ), em que x é uma
variável quantitativa, n pode ser uma variável offset, no caso da Poisson,
ou o número de ensaios Bernoulli, no caso da binomial.

 

Aí vai o exemplo que vem dando erro, que eu suponho ter a ver com o optim,
devido aos valores muito pequenos dos parâmetros b0, b1 e g e do valor da
offset “n” ser muito grande.

 

Agradeço a todos pela atenção! Abraços.

 

 

# no R

 

require(gnm)

 

x  <- sort( rep(0:5,3) )

b0 <- 1.5*1e-6

b1 <- 3*1e-6

g  <- 3*1e-1

n  <- rep(1e+8,length(x))

p  <- (1 - exp(-(b0+b1*x)) ) * exp(-g*x) ; p

mu <- n*p

 

## Funções não lineares

# Parametrização original

margo.betas <- function(x){

  list(predictors=list(beta0=1,beta1=1,gamma=1),

       variables=list(substitute(x)),

       term=function(predictors,variables){

 
pred<-paste("(1-exp(",variables[1],"+",variables[2],"*",predictors[1],"))*(e
xp(",

                     variables[3],"*",predictors[1],"))",sep="")

       })

}

class(margo.betas) <- "nonlin"

 

# Reparametrização, com b0 = 1/(ln k0), b1=1/(ln k1) e g=1/(ln k2)

margo.ks <- function(x){

  list(predictors=list(k0=1,k1=1,k2=1),

       variables=list(substitute(x)),

       term=function(predictors,variables){

 
pred<-paste("(1-",variables[1],"*(",variables[2],"^",predictors[1],"))*(",

                     variables[3],"^",predictors[1],")",sep="")

       })

}

class(margo.ks) <- "nonlin"

 

 

######## ABORDAGEM POISSON

set.seed(0119)

y <-rpois(length(x),lambda=mu)

 

eta.margo1   <- formula(y ~ -1 + Mult(offset(n), margo.ks(x) ) )

eta.margo2   <- formula(y ~ -1 + Mult(offset(n), margo.betas(x) ) )

 

gnm(eta.margo1,family=poisson,tolerance=1e-50)

gnm(eta.margo2,family=poisson,tolerance=1e-50)

 

 

######## ABORDAGEM BINOMIAL

set.seed(0119)

y            <-rbinom(length(x),size=n,prob=p)

 

eta.margo1   <- formula(cbind(y,n-y) ~ -1+ margo.ks(x) )

eta.margo2   <- formula(cbind(y,n-y) ~ -1+ margo.betas(x) )

 

 

gnm(eta.margo1,family=binomial)

gnm(eta.margo2,family=binomial)

 

 

___________________________________________________

Davi Butturi-Gomes

Doutorando em Estatística e Experimentação Agronômica - ESALQ /USP -
Piracicaba

Mestre em Biometria - IB/UNESP - Botucatu

Ecólogo - IB/UNESP - Rio Claro

Skype: davibg87

(19)8206.6565

 

-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20130508/3af95fb4/attachment.html>