Como criar um conjunto de dados com probabilidade condicional?

8

Suponha que uma determinada doença ( ) tenha uma prevalência de . Suponha também que um determinado sintoma ( ) tenha uma prevalência (na população geral = pessoas com essa doença D e pessoas sem essa doença [provavelmente com outra doença, mas não é importante]) de . Em uma pesquisa anterior, descobriu-se que a probabilidade condicional (a probabilidade de apresentar o sintoma , dada a doença é de ).D31000S51000P(S|D)=30%SD30%

Primeira pergunta : poderia ser interpretado como equivalente à prevalência do sintoma no grupo de pessoas com a doença ?P(S|D)SD

Segunda pergunta : quero criar em R um conjunto de dados, o que mostra que:

P(D|S)=P(S|D)P(D)P(S)
Com meus dados fictícios, podemos calcular , que é interpretado dessa maneira : dado um paciente com o sintoma , a probabilidade de ele ter a doença é .P(D|S)=0.18SD18%

Como fazer isso? Se eu usar simplesmente a samplefunção, meu conjunto de dados não possui as informações que :P(S|D)=30%

symptom <- sample(c("yes","no"), 1000, prob=c(0.005, 0.995), rep=T)
disease <- sample(c("yes","no"), 1000, prob=c(0.002, 0.998), rep=T)

Então, minha pergunta é: como criar um bom conjunto de dados, incluindo a probabilidade condicional que desejo?

EDIT : Também postei a mesma pergunta no stackoverflow.com ( /programming/7291935/how-to-create-a-dataset-with-conditional-probability ), porque, na minha opinião, minha pergunta é herdado do programa de linguagem R, mas também da teoria estatística.

Tommaso
fonte
3
A cortesia comum é indicar que você postou em outro site da SE. stackoverflow.com/questions/7291935/…
Brandon Bertelsen
1
Sinalizei sua pergunta no SO para migração. Por favor, não faça postagens cruzadas!
chl

Respostas:

11

Você conhece as seguintes probabilidades marginais

                Symptom        Total
                Yes     No
Disease Yes      a       b     0.003
        No       c       d     0.997  
Total           0.005   0.995  1.000

e que a/(a+b) = 0.3assim isso se torna

                Symptom        Total
                Yes     No
Disease Yes     0.0009  0.0021 0.003
        No      0.0041  0.9929 0.997  
Total           0.005   0.995  1.000

e de fato a/(a+c) = 0.18como você declarou.

Então, no R você poderia codificar algo como

diseaserate <- 3/1000
symptomrate <- 5/1000
symptomgivendisease <- 0.3

status  <- sample(c("SYDY", "SNDY", "SYDN", "SNDN"), 1000, 
            prob=c(diseaserate * symptomgivendisease,
                   diseaserate * (1-symptomgivendisease),
                   symptomrate - diseaserate * symptomgivendisease,
                   1 - symptomrate - diseaserate * (1-symptomgivendisease)),
            rep=TRUE)
symptom <- status %in% c("SYDY","SYDN")
disease <- status %in% c("SYDY","SNDY")

embora você deva observar que 1000 é uma amostra pequena quando um dos eventos tem uma probabilidade de 0,0009 de acontecer.

Henry
fonte
Solução incrível, funciona muito bem! Agora eu posso criar um conjunto de dados mostrando o que a fórmula de Bayes pode calcular. Muito obrigado!
Tommaso
Lhe disse que alguém viesse com algo mais elegante;)
fomite
@henry, ficaria muito feliz se você pudesse dar uma olhada na minha nova pergunta aqui: stats.stackexchange.com/questions/15202/… . É uma generalização desta questão, com 2 sintomas.
Tommaso
3

A tablefunção retorna um objeto semelhante a uma matriz:

> symptom <- sample(c("yes","no"), 100, prob=c(0.2, 0.8), rep=TRUE)
> disease <- sample(c("yes","no"), 100, prob=c(0.2, 0.8), rep=TRUE)
> dataset <- data.frame(symptom, disease)
> dst_S_D <-with(dataset, table(symptom, disease))
> dst_S_D
       disease
symptom no yes
    no  65  13
    yes 17   5

Então o Pr (D | S = "sim") =

> probD_Sy <- dst_S_D[2, 2]/sum(dst_S_D[2, ] )
> probD_Sy
[1] 0.2272727

Alterei o problema porque na primeira vez em que o executei com seus parâmetros, obtive:

> dst_S_D <-with(dataset, table(symptom, disease)); dst_S_D
       disease
symptom   no  yes
    no  9954   22
    yes   24    0

E eu pensei que um Pr (D | S = "yes") de 0 era bastante chato. Se você deseja executar isso muitas vezes, construa uma função e use essa função com a replicatefunção

Aqui está um método para construir um conjunto de dados que aplica uma probabilidade diferente de doença no grupo sintomático que é 3 vezes maior do que o usado no grupo assintomático:

symptom <- sample(c("yes","no"), 10000, prob=c(0.02, 0.98), rep=TRUE)
dataset <- data.frame(symptom, disease=NA)
dataset$disease[dataset$symptom == "yes"] <- 
       sample(c("yes","no"), sum(dataset$symptom == "yes"), prob=c(0.15, 1-0.15), rep=TRUE)
dataset$disease[dataset$symptom == "no"] <- 
        sample(c("yes","no"), sum(dataset$symptom == "no"), prob=c(0.05, 1-0.05), rep=TRUE)
 dst_S_D <-with(dataset, table(symptom, disease)); dst_S_D
#       disease
symptom   no  yes
    no  9284  509
    yes  176   31
DWin
fonte
Truque perfeito, agradável e elegante! Eu adicionei algumas informações novas na minha resposta, para formalizar melhor o que estou procurando.
Tommaso
2

Eu diria que sua pergunta não é realmente muito dependente da linguagem R e é mais apropriada aqui, porque - para ser franco - a geração de dados como esse é principalmente uma tarefa estatística, e não de programação.

Primeira pergunta: p (S | D) é o risco de ter o sintoma S em uma população com doença D. Pode ser diretamente comparável à prevalência de certas advertências, como o sintoma que não afeta a duração da doença. Considere o seguinte exemplo: Um dos sintomas do SuperEbola é Morte instantânea, com p (Morte | Super Ebola) = 0,99. Aqui, sua prevalência do sintoma seria realmente extremamente baixa (de fato, 0,00), pois ninguém que você pode experimentar a doença tem o sintoma.

Segunda pergunta: eu voltaria a isso de uma maneira um pouco gradual. Primeiro, calcule o risco de linha de base do sintoma que você precisará para obter 0,15 em toda a população, levando em consideração que 0,03% da sua população terá uma taxa mais alta. Em seguida, gere essencialmente duas probabilidades:

  • Risco de doença = 0,003
  • Risco de sintoma = risco de base calculado + aumento relativo devido à doença * indicador binário do status da doença

Em seguida, gere dois números aleatórios uniformes. Se o primeiro for menor que 0,003, eles têm a doença. Isso é inserido no cálculo de risco do segundo e, se o número aleatório de cada indivíduo for menor que o risco, eles terão o sintoma.

Essa é uma maneira penosa e deselegante de fazer as coisas, e é provável que alguém venha com uma abordagem muito mais eficiente. Mas, nos estudos de simulação, encontro cada etapa do código, e é útil mantê-lo o mais próximo possível de ver um conjunto de dados no mundo real.

Fomite
fonte
Obrigado pela resposta; o exemplo do SuperEbola é realmente educativo e útil! O restante de sua resposta permanece bastante incerto, para mim, especialmente quando você diz "calcule o risco de linha de base do sintoma e precisará obter 0,15 em toda a população, levando em consideração que 0,03% da sua população estará em uma taxa mais alta " Como calcular esse risco de linha de base?
Tommaso
Honestamente, é uma dor de fazer. Se eu fosse você, mudaria um pouco o meu exemplo - em vez de afirmar que o risco geral na população é de 0,15, eu diria que o risco de referência nos não-doentes é, digamos, 0,15 ou 0,10 e, em seguida, determinarei o aumento quero o risco do doente e deixo o risco geral cair onde pode, em vez de tentar defini-lo. É consideravelmente mais fácil codificar, embora você possivelmente não tenha números tão limpos no final.
Fomite 03/09/11
0

Primeira pergunta:

Sim, claro que é quase a definição, embora você tenha algum erro associado ao tamanho da amostra. isto é, exatamente correto em um tamanho infinito de amostra.

Segunda questão:

Isso se chama Teorema de Bayes , mas presumo que você já saiba disso. Agora, dadas as informações que você forneceu, obtenho a probabilidade de P (D | S) como 0,18 ou 18%:

P(S|D)P(D)
----------
   P(S)

  0.3*(3/1000)
= ------------
    (5/1000)

= 0.18

Agora, infelizmente, eu não estou muito familiarizado com o R, então não posso realmente ajudá-lo com um programa exato. Mas certamente as quantidades de pessoas que se enquadram em cada grupo são bastante fáceis de calcular:

Para o seu conjunto de amostras 10000, você precisa:

  1. 50 pessoas com sintomas (população * P (S))
  2. 9 pessoas devem ter sintomas e a doença (50 * P (D | S))
  3. 21 pessoas com a doença e sem sintomas (população * P (D) = 30 e já temos 9)

O que deve tornar a geração de uma população adequada bastante trivial.


fonte
Sim, o valor verdadeiro é 0,18. Desculpe pela digitação incorreta. A segunda parte da sua resposta está correta, mas o problema é criar um conjunto de dados (em R) que realmente tenha 9 pessoas com doença e sintoma. A função "amostra" cria corretamente 50 e 30 "sim" para, respectivamente, sintoma e doença; mas não garante que 9 pessoas (em 30) também estejam no grupo de "sim-doença".
Tommaso
Novamente, com medo de que você precise de alguém mais familiarizado com R do que eu para ajudá-lo no uso dessa função de amostra. No entanto, você sempre pode gerar uma população muito maior e coletar aleatoriamente 10000 amostras.