Montando função de verossimilhança: um exemplo para dados censurados
## Verossimilhança para o parâmetro de média da normal
## (assumindo variância conhecida igual a 16)
## i) Assumindo só ter observações completas
dat <- c(23, 24, 27, 20, 32, 26, 28)
## um primeira forma de escrever a função
f1 <- function(mu) {
sum(dnorm(dat, mean=mu, sd = 4, log=TRUE))
}
f1(25)
mu.vals <- seq(20, 32, length=100)
plot(mu.vals, sapply(mu.vals, f1), type="l")
## algumas opções no gráfico
plot(mu.vals, sapply(mu.vals, f1), type="l", xlab=expression(mu),
ylab=expression(l(mu)), main=expression(paste("log-verossimilhança de ", mu)))
## agora mais geral colocando os dados como um argumento da função
## para permitir receber outros conjuntos de dados
f1 <- function(mu, dados) {
sum(dnorm(dados, mean=mu, sd = 4, log=TRUE))
}
f1(25, dados=dat)
plot(mu.vals, sapply(mu.vals, f1, dados=dat), type="l", xlab=expression(mu),
ylab=expression(l(mu)), main=expression(paste("log-verossimilhança de ", mu)))
## agora vetorizando a função para permitir receber vários valores de mu
f1 <- function(par, dados) {
lVero <- function(mu) sum(dnorm(dados, mean=mu, sd = 4, log=TRUE))
sapply(par, lVero)
}
f1(25, dados = dat)
plot(mu.vals, f1(mu.vals, dados=dat), type="l", xlab=expression(mu),
ylab=expression(l(mu)), main=expression(paste("log-verossimilhança de ", mu)))
## Estimador de MV de mu
## a) solução analitica
mean(dat)
abline(v=mean(dat), col=3, lwd=4) # ver o gráfico
## embora nao necessario até aqui, a solucao numerica pode ser obtida com:
mu.est1 <- optimize(f1,c(10,50), dados = dat, maximum=T)
mu.est1
abline(v=mu.est1$max) ## ver o gráfico
## Introduzindo agora as informações censuradas
f2 <- function(par, dados) {
sapply(par,
function(mu){
## obs sem censura
semC <- sum(dnorm(dados, mean=mu, sd = 4, log=TRUE))
## obs com censura a direita
Cdir <- sum(pnorm(c(30, 30), mean=mu, sd=4, low=F, log=TRUE))
logL <- semC + Cdir
return(logL)
}
)
}
## veja como mudou a curva
lines(mu.vals, f2(mu.vals, dados=dat), col=2)
mu.vals <- seq(20, 35, length=100)
## fazendo num novo gráfico
plot(mu.vals, f2(mu.vals, dados=dat), ty="l")
## agora escrevendo a função de forma mais geral com um argumento
## para receber os pontos de censura
f2 <- function(par, dados, cenD) {
sapply(par,
function(mu){
## obs sem censura
semC <- sum(dnorm(dados, mean=mu, sd = 4, log=TRUE))
## obs com censura a direita
Cdir <- sum(pnorm(cenD, mean=mu, sd=4, low=F, log=TRUE))
logL <- semC + Cdir
return(logL)
}
)
}
plot(mu.vals, f2(mu.vals, dados=dat, cenD=rep(30, 2)), ty="l", xlab=expression(mu),
ylab=expression(l(mu)), main=expression(paste("log-verossimilhança de ", mu)))
## e a estimativa de MV é agora obtida somente por maximização numérica
mu.est2 <- optimize(f2,c(10,50), dados=dat, cenD=rep(30, 2), maximum=T)
mu.est2
abline(v=mu.est2$max)
## Agora definindo uma função geral que aceita diferentes tipos de censura
## (a direita, esquerda ou intervalar)
fcens <- function(par, dados=NULL, cenD=NULL, cenE=NULL, cenI=NULL) {
sapply(par,
function(mu){
## obs sem censura
if(is.null(dados)) lVero1 <- 0
else lVero1 <- sum(dnorm(dados, mean=mu, sd = 4, log=TRUE))
## obs com censura a direita
if(is.null(cenD)) lVero2 <- 0
else lVero2 <- sum(pnorm(cenD, mean=mu, sd=4, lower=F, log=TRUE))
## obs com censura a esquerda
if(is.null(cenE)) lVero3 <- 0
else lVero3 <- sum(pnorm(cenE, mean=mu, sd=4, log=TRUE))
## obs com censura intervalar
if(is.null(cenI)) lVero4 <- 0
else lVero4 <- sum(apply(cenI, 1,
function(int)
log(diff(pnorm(int, mean=mu, sd=4)))))
logL <- lVero1 + lVero2 + lVero3 + lVero4
return(ifelse(logL==0, NULL, logL))
}
)
}
plot(mu.vals, fcens(mu.vals, dados=dat, cenD=rep(30, 2),
cenI= cbind(rep(28, 3), rep(30,3))), ty="l", xlab=expression(mu),
ylab=expression(l(mu)), main=expression(paste("log-verossimilhança de ", mu)))
## escrevendo a função com um único argumento para censura:
## censura deve estar em uma matrix com
## [-Inf, x] para censura a esquerda
## [x, Inf] para censura a direita
## [x1, x2] para censura intervalar
fcens <- function(par, dados=NULL, censurados=NULL) {
sapply(par,
function(mu){
## obs sem censura
if(is.null(dados)) lVero1 <- 0
else lVero0 <- sum(dnorm(dados, mean=mu, sd = 4, log=TRUE))
## obs com censura a direita
if(is.null(censurados)) lVero1 <- 0
else lVero1 <- sum(apply(censurados, 1,
function(xc)
log(diff(pnorm(xc, mean=mu, sd=4)))))
logL <- lVero0 + lVero1
return(ifelse(logL==0, NULL, logL))
}
)
}
## note que o gráfico não se altera
lines(mu.vals, fcens(mu.vals, dados=dat, cen = censMat), col=2)
plot(mu.vals, fcens(mu.vals, dados=dat, cen = censMat), ty="l", xlab=expression(mu),
ylab=expression(l(mu)), main=expression(paste("log-verossimilhança de ", mu)))
## No exemplo
censMat <- rbind(c(30, Inf), c(28, 30))
censMat
## aora replicando as censuras (duas do primeiro tipo de 3 do segundo)
censMat <- censMat[c(rep(1, 2), rep(2, 3)),]
censMat
## e a estimativa de MV por maximização numérica pode ser obtida como antes
mu.est3 <- optimize(fcens,c(10,50), dados=dat, censurados=censMat, maximum=T)
mu.est3
abline(v=mu.est3$max)
## as 3 estimativas obtidas foram portanto
c(mu.est1$max, mu.est2$max, mu.est3$max)
Exemplo de simulação: teste aleatorizado para comparação de médias de duas amostras
## Exemplo do teste-t
mandibulas <- data.frame(
femeas = c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112),
machos = c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111)
)
with(mandibulas, t.test(femeas, machos))
## uma alternativa ao teste t: teste aleatorizado
## jutando todos os dados (sob H0)
todos <- unlist(mandibulas)
## inicialmente mostrando uma forma não eficiente de fazer (for() deve ser evitado)
Nsim <- 999
dm <- numeric(Nsim)
for(i in 1:Nsim){
st <- sample(todos)
dm[i] <- mean(st[1:10]) - mean(st[11:20])
}
## outra forma
dm <- replicate(Nsim, function(dados){
st <- sample(dados)
return(mean(st[1:10]) - mean(st[11:20]))})
hist(dm, prob=T)
lines(density(dm))
rug(dm)
## acrescentando a diferença real ao gráfico das simuladas
dreal <- with(mandibulas, mean(machos) - mean(femeas))
dreal
abline(v=dreal, col=2)
## p-valor estimado (pelas simulações)
pvalor <- sum(dm < dreal)/(Nsim+1) ## p/ teste unilateral
pvalor
2* pvalor # p/ teste bilateral