2016-07-08 5 views
1

J'essaie de faire une bootstrap à partir d'un modèle ajusté (pente maximale de la courbe). Bien que le code que je pense puisse être plus efficace, je pense que cela fonctionne généralement (des suggestions pour le rendre plus efficace seraient également les bienvenues). Le problème que j'ai est à la fin quand j'ai mes valeurs d'échantillon finales où il y a Inf dans l'une des cellules qui m'empêche d'obtenir un intervalle de confiance. Je ne suis pas sûr que ce soit un problème de bootstrap ou simple de la façon dont je l'ai codé. Exemple de code:Intervalle de confiance bootstrap avec "Inf" dans les estimations finales - paquet boot/dplyr

library(boot) 
library(dplyr) 
df <- data.frame (id=rep(1:10,each=10), 
         time1=rep(1:10,10), 
         ff=runif(100, 100, 150), 
         gg=runif(100, 120, 170)) 



    set.seed(10) 

#function 
    maxx<-function(formula, data,varr,indices) { 
     data <- data[indices,] # allows boot to select sample 
     fit <- lm(formula, data=data) 
     data$fit<-fit$fitted.values 
     data<-filter(data,time1>5) 
     data1<-data %>% group_by_(varr) %>% mutate(derivative = c(NA,diff(fit)/diff(time1))) %>% 
     group_by(id) %>% 
     slice(which.max(derivative)) %>% 
     ungroup() %>% 
     as.data.frame() 
     return(mean(data1$derivative,na.rm = TRUE)) 
    } 

res <- boot(data=df, statistic=maxx, 
        R=10, formula=ff~gg,varr="id") 

cela fonctionne, mais ne peut pas obtenir S'en raison d'un Inf dans l'une des cellules

R>res 

ORDINARY NONPARAMETRIC BOOTSTRAP 


Call: 
boot(data = df, statistic = maxx, R = 10, formula = ff ~ gg, 
    varr = "id") 


Bootstrap Statistics : 
    original bias std. error 
t1* 0.1726803  Inf   NaN 
     res 
     res$t 

On peut le voir en explorant davantage:

R>res$t 
      [,1] 
[1,] 0.50399242 
[2,] 0.52171509 
[3,] 0.04568459 
[4,] 1.41317481 
[5,] 0.39741115 
[6,] 0.10703703 
[7,] 0.52206909 
[8,] 0.91624253 
[9,]  Inf 
[10,] 0.05076168 


R>boot.ci(res, type="normal") 
Error in ci.out[[4L]] : subscript out of bounds 

Peut-être que je manquant quelque chose, mais je pensais return(mean(data1$derivative,na.rm = TRUE)) résoudrait tous les problèmes avec NA. Est-ce que quelqu'un a des suggestions s'il vous plaît? Je suppose que c'est une petite solution. C'est la première fois que je fais un bootstrap alors excuses si je fais quelque chose de naïf. Si quelqu'un savait comment rendre le code global plus efficace, ce serait génial car je suis en train d'utiliser un grand modèle à effets aléatoires plutôt que le simple modèle présenté. Merci

Répondre

1

Parce que boot rééchantillons avec le remplacement, vous pouvez obtenir des valeurs répétées de time1 pour une donnée id dans un échantillon donné. Lorsque cela se produit, votre calcul de dérivée diff(fit)/diff(time) évalue à 0/0, ce qui devrait retourner NaN. Les valeurs NaN ne devraient pas être intrinsèquement problématiques pour le reste de votre fonction maxx (bien que je ne sois pas tout à fait clair sur ce que vous essayez de faire), mais je suppose qu'en raison d'une inexactitude en virgule flottante diff(fit) renvoie parfois des valeurs non nulles qui sont divisées par zéro pour produire Inf ou -Inf. Les fonctions which.max et mean n'ignorent pas Inf, et donc votre fonction maxx renvoie occasionnellement Inf.