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