[R-br] mensagem de erro em simulacao

[Maurício Lordêlo]

#Caros,
#Estou fazendo uso da função "bild" pertencente ao pacote com mesmo nome.
#Nos dados que são usados como exemplo do pacote (locust) temos:
install.packages("bild")
require(bild)
locust2 <- bild(move ~ time+feed, data = locust,start = NULL, aggregate =
feed, dependence = "MC1R")
summary(locust2)
#O programa retorna os resultados sem problemas, como mostrado abaixo:
Call:
bild(formula = move ~ time + feed, data = locust, aggregate = feed,start =
NULL, dependence = "MC1R")
Number of profiles in the dataset:  24
Number of profiles used in the fit:  24
Log likelihood:  -1537.923
AIC:  3085.845
Coefficients:
            Label     Value Std. Error t value  p-value
(Intercept)     1 -1.048305 0.30730336  -3.411 0.000647
time            2  1.572124 0.14600334  10.768 0.000000
feed1           3 -3.161945 0.43072021  -7.341 0.000000
log.psi1        4  1.070874 0.09833972  10.890 0.000000
Random effect (omega):
      Value  Std. Error
-0.03532051  0.34195846
Message:  0

#Porém, se faço
warnings()
#Ele me retorna
Mensagens de aviso:
1: In sqrt(prob * (1 - prob)) : NaNs produzidos
2: In sqrt(prob * (1 - prob)) : NaNs produzidos
3: In sqrt(prob * (1 - prob)) : NaNs produzidos
4: In sqrt(prob * (1 - prob)) : NaNs produzidos
#Mesmo com estas mensagens consigo obter os resultados.
#O problema surge quando peço para simular muitos ajustes deste modelo.
#Por exemplo,
n.time=4   #numero de medidas no tempo
n.trat=3    #numero de tratamentos
n.rep=10  #numero de repeticoes
ns=15 #numero de simulacoes
n.ind <- n.trat * n.rep  #numero de individuos
n <- n.trat * n.rep * n.time   #numero total de registros
id <- rep(1:n.ind, each = n.time)   #identificacao do individuo
time <- rep(1:n.time, times = n.ind)   #tempo
trat <- rep(factor(LETTERS[1:n.trat]), each = n.rep*n.time) #tratamentos
dados=data.frame(id,time,trat) #estrutura do banco de dados
#funcao para gerar resposta binaria
z <- function()
{
 y <- integer(n)
 for(i in 1:n)
 y[i] <- rbinom(1,1,.5)
 return(y)
}
repet <- replicate(ns, z(), simplify = FALSE)

#funcao para obter as estimativas
coef_bild <- function(y)
{
dat <- data.frame(dados, y=z())
coefs=bild(y ~ trat+time, data = dat,
start=NULL,dependence="MC1R")@coefficients
coefs_bild=coefs[1:4]
}
#lista dos coeficientes ajustados
coefs_ajust <- do.call(rbind, lapply(repet, coef_bild))
#obtencao das medias dos coeficientes estimados
media_coefs <- apply(coefs_ajust, 2, mean)


#Com apenas ns=15 ele me dá as mensagens de aviso sem eu pedir (as mesmas
do exemplo para o conjunto de dados "locust"), mas me retorna o que preciso.
#Quando coloco ns=2000 (que é a quantidade que preciso)  demora bastante
tempo, não retorna o que preciso e emite a mensagem:
Error in solve.default(Info) :
  system is computationally singular: reciprocal condition number =
1.69179e-16
In addition: There were 50 or more warnings (use warnings() to see the
first 50)
#As 50 ou mais mensagens são as mesmas que aparecem quando rodo com uma
quantidade menor . Alguém tem alguma dica ou forma de conseguir estes
resultados para o número de simulações que desejo?
# Caso nao encontre uma solução, verifiquei que com ns=100 ele faz. Tentei
fazer um "for de 1 a 20, com ns=100" das funções que defini como "z" e
"coef_bild" e depois concatenar, mas não consegui. Em último caso, peço
ajuda a quem saiba fazer para este caso específico.
#Agradeço antecipadamente,
#Maurício
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20131019/01b9bb0f/attachment.html>