mensagem de erro em simulacao

#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
participantes (1)
-
Maurício Lordêlo