1

Estou tentado estimar os limites inferiores e superiores usando o intervalo de confiança de 95%, mas o R retorna a seguinte mensagem:

Warning message:
In rnorm(nrow(df), media/-0.832, media/-0.399) : NAs produzidos

Data frame que estou utilizando:

y<-c(-1.0, -1.0, -1.0, -1.0, 0.769, 0.623, -1.0, 0.327, -1.0, -0.638,
   -1.0, -1.0, -0.618, -1.0, -0.670, -1.0, 0.028, -1.0, -1.0, -1.0,
   0.235, -0.286, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.148,
   0.857, -0.918)
df<-data.frame(y)

n<-length(df$y)
media<-mean(df$y)
var<-var(df$y)

ic<-media+qt(c(0.025, 0.975), df=n-1)*sqrt(var/length(df$y))

media<-abs(mean(df$y))
df$lim_inf<-df$y-rnorm(nrow(df), media/-0.832, media/-0.399)
df$lim_sup<-df$y+abs(rnorm(nrow(df), media/-0.832, media/-0.399))

Como faço para o R não retorna NAs?

2
  • 2
    Eu não consigo entender o que está acontecendo aqui. Por exemplo, porque calcular media<-abs(mean(df$y))? Existe alguma restrição para a média ser não-negativa? Porque, sem calcular o valor absoluto dela, seu valor é negativo. 29/06/2020 às 15:41
  • 2
    rnorm(nrow(df), media/-0.832, media/-0.399): como é que o desvio padrão pode ser negativo? 29/06/2020 às 15:46

2 Respostas 2

2

Não sei se o seguinte é o que a pergunta pede.

intconf <- function(x, nivel = 0.95, normal = FALSE){
  qq <- c((1 - nivel)/2, 1 - (1 - nivel)/2)
  n <- length(x)
  xbar <- mean(x)
  s2 <- var(x)
  ic <- if(normal && (n >= 30)){
    xbar + qnorm(qq)*sqrt(s2/n)
  } else {
    xbar + qt(qq, df = n - 1)*sqrt(s2/n)
  }
  setNames(ic, c("lim_inf", "lim_sup"))
}

Esta função calcula intervalos de confiança para a média. Se a amostra x for pequena, n < 30, utiliza a distribuição t de Student, caso contrário e se o argumento normal = TRUE utiliza a distribuição normal. O valor 30 é muitas vezes considerado para decidir se a amostra é pequena (não vem no artigo da Wikipedia sobre este tema).

Com os dados da pergunta,

intconf(df$y)
#   lim_inf    lim_sup 
#-0.8479307 -0.4110068

Como a amostra tem 32 elementos pode-se usar a distribuição normal.

nrow(df)
#[1] 32

intconf(df$y, normal = TRUE)
#   lim_inf    lim_sup 
#-0.8394098 -0.4195277

Para ver estes resultados, pode-se traçar um histograma com as linhas verticais correspondentes aos intervalos de confiança.

icdf <- rbind(
  intconf(df$y),
  intconf(df$y, normal = TRUE)
)

brks <- seq(min(df$y), max(df$y), length.out = 7)
hist(df$y, breaks = brks, freq = FALSE)
abline(v = icdf, col = c("red", "blue"))
legend("top", legend = c("t", "normal"), lty = 1, col = c("red", "blue"), horiz = TRUE)

inserir a descrição da imagem aqui

1

Rui Barradas já deu uma boa resposta sobre como calcular o intervalo de confiança, esta é para entender o que estava fazendo de errado. Porque não é apenas um erro de programação, mas de entendimento do IC.

A mensagem de erro

O erro que o R está retornando vem de usar um valor negativo para o desvio padrão na função rnorm:

> rnorm(5, -.4, .1)
[1] -0.4070113 -0.3860859 -0.3772370 -0.4746763 -0.3817434

> rnorm(5, -.4, -.1)
[1] NaN NaN NaN NaN NaN
Warning message:
In rnorm(5, -0.4, -0.1) : NAs produced

Isso porque desvio padrão não assume valores negativos. Mas por que está usando rnorm para calcular os limites? Por que usar o valor absoluto da média?

Cálculo do IC

Existe um motivo para os pacotes básicos do R não terem uma função para calcular IC. Ele é uma medida do quão confiável é a estimativa de um parâmetro, o modo com que é calculado depende do parâmetro e do modelo atribuído aos dados.

Sua fórmula está correta para mu estimado de uma amostra pequena (ou com grande variância) seguindo distribuição próxima da normal. Compare com o resultado de um teste t:

ic <- media + qt(c(0.025, 0.975), df = n-1) * sqrt(var/n)

> ic
[1] -0.8479307 -0.4110068

> t.test(dados$y)$conf.int
[1] -0.8479307 -0.4110068
attr(,"conf.level")
[1] 0.95

Mas esses são bons pressupostos para seus dados? Contraste com o resultado da função do @rui-barradas. E compare com o resultado usando diferentes métodos com a média estimada por bootstrap:

library(boot)

media.b <- boot(dados$y, function(x,i) mean(x[i]), 10000)

> boot.ci(media.b, conf = .95, type = c("norm","basic", "perc", "bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL :
boot.ci(boot.out = media.b, conf = 0.95, type = c("norm", "basic",
    "perc", "bca"))

Intervals :
Level      Normal              Basic
95%   (-0.8366, -0.4206 )   (-0.8527, -0.4365 )

Level     Percentile            BCa
95%   (-0.8224, -0.4063 )   (-0.8041, -0.3735 )
Calculations and Intervals on Original Scale

Adicionando limites aos dados

Uma vez que calculou o intervalo, bastaria adicionar/subtrair os valores de y:

df$lim_inf <- df$y - ic[1]
df$lim_sup <- df$y + ic[2]

Mas isso é errado. O intervalo de confiança aplica-se ao parâmetro estimado (a média, no caso), não aos valores individuais. Uma opção mais adequada é usar o erro padrão:

ep <- sqrt(var/n)

df$lim_inf <- df$y - ep
df$lim_sup <- df$y + ep
6
  • Carlos, estou tentando calcular o "ep" para os valores individuais: "ep<-sqrt(var/-1.0)", entretanto o R retorna: "non-numeric argument to binary operator". 4/07/2020 às 13:20
  • @MiltondePaula Porquê dividir por um valor negativo? A raíz quadrada vai dar... 4/07/2020 às 17:22
  • @Rui Barradas, perdão pelo equivoco. Meus valores individuais são resultados de um índice de preferência (Jacob's Index), e muitos são negativos, assim, preciso calcular o erro padrão para cada valor. 4/07/2020 às 17:37
  • @MiltondePaula Sejam quais forem os valores, a sua variância é sempre positiva (é a média de uma soma de quadrados). Dividindo por um valor negativo dá negativo e a raíz quadrada dá NaN. Se lhe está a dar esse erro, passa-se outra coisa, talvez var esteja a ser visto como uma closure, uma função, e não como um número real. 4/07/2020 às 17:43
  • @MiltondePaula sqrt(var/-1.0)Error in var/-1 : non-numeric argument to binary operator. A primeira parte da mensagem de erro diz em que operação, na divisão var/-1. 4/07/2020 às 17:45

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .