Aula de Modelos Lineares I

Modelo Normal com média e variância constantes

1 Contextualização

Uma fábrica de peças automotivas busca garantir a padronização do diâmetro de um componente essencial para o funcionamento dos motores. O diâmetro ideal desse componente é de aproximadamente 40 mm, porém, pequenas variações são inevitáveis no processo de fabricação. Estima-se que um desvio padrão de 3 mm seja aceitável. Para monitorar a qualidade da produção, foram coletadas medidas de 1000 peças selecionadas aleatoriamente, com o objetivo de verificar se o processo está operando dentro dos padrões esperados.

Para modelar essas medições, assumimos que o diâmetro das peças segue a distribuição normal com média \(\theta = 40\)mm e variância \(\sigma^2 = 3^2 = 9\). Dessa forma, cada medida pode ser expressa como:

\[ Y_i = \theta + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2), \quad i = 1, \dots, n=1000, \]

onde

• \(Y_i\) representa o diâmetro da \(i\)-ésima peça.
• \(\theta = 40\)mm é o valor ideal do diâmetro.
• \(e_i\) é um erro aleatório, assumido como normalmente distribuído com média zero e variância \(\sigma^2 = 9\) mm\(^2\).

2 Dados Simulados

Simule um conjunto de dados para o problema descrito acima.

#Simulando dados
n       = 1000  #tamanho da amostra
theta   = 40    #media populacional
sigma2  = 9 #variancia
e       = rnorm(n, 0, sqrt(sigma2)) #efeito aleatorio
y       = theta + e #variavel resposta

3 Analisando os dados

Realize uma análise exploratória da variável resposta \(Y\) e dos erros aleatórios \(e\).

3.1 Medidas resumo

Calcule medidas descritivas como média, mediana, variância e desvio padrão.

# Medidas descritivas
mean(y)  # Média

## [1] 39.99922

median(y)  # Mediana

## [1] 39.9236

var(y)  # Variância

## [1] 8.9281

sd(y)  # Desvio padrão

## [1] 2.987993

3.2 Histogramas

Visualize a distribuição da variável resposta e do erro aleatório por meio de histogramas e curvas de densidade ajustadas.

#analisando as densidades de Y e de e
par(mfrow=c(1,2))
hist(y, freq=F, main="", cex.lab=1.4, cex.axis=1.4, ylab="densidade", xlab="Y")
curve(dnorm(x, mean(y), sd(y)), col="red", add=T, lwd=2)
#analisando a densidade do erro aleatorio
hist(e, freq=F, main="", cex.lab=1.4, cex.axis=1.4, ylab="densidade", xlab="erro aleatorio")
curve(dnorm(x, mean(e), sd(e)), col="red", add=T, lwd=2)
abline(v=0, lwd=2, lty=3, col="red")

3.3 Normalidade (QQ Plots)

Verifique a normalidade da variável resposta e do erro aleatório utilizando QQ Plots.

#outra forma de verificar normalidade
par(mfrow=c(1,2))
qqnorm(y, bty="n", xlab="quantis teóricos", ylab="quantis amostrais",
       cex.lab=1.4, cex.axis=1.4, main="QQ Plot da variável resposta")
qqline(y, col="red", lwd=2)
qqnorm(e, bty="n", xlab="quantis teóricos", ylab="quantis amostrais",
       cex.lab=1.4, cex.axis=1.4, main="QQ Plot dos erros aleatórios")
qqline(e, col="red", lwd=2)

3.4 Normalidade (Shapiro-Wilk e Kolmogorov-Smirnov)

Teste a normalidade da variável resposta e do erro aleatório usando os testes de Shapiro-Wilk e Kolmogorov-Smirnov.

#testando normalidade 
shapiro.test(y)

## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.9975, p-value = 0.1306

ks.test(y, "pnorm", mean(y), sd(y))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.021047, p-value = 0.7675
## alternative hypothesis: two-sided

shapiro.test(e)

## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.9975, p-value = 0.1306

ks.test(e, "pnorm", mean(e), sd(e))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  e
## D = 0.021047, p-value = 0.7675
## alternative hypothesis: two-sided

3.5 Independência e dispersão

Faça um gráfico de dispersão para a variável resposta e outro para os erros aleatórios. Analise os gráficos.

#analisando a dispersao
par(mfrow=c(1,2))
plot(y, cex.lab=1.4, cex.axis=1.4, xlab="unidades", ylab="Y")
plot(e, cex.lab=1.4, cex.axis=1.4, xlab="unidades", ylab="e")

3.6 Simetria

Avalie a simetria com boxplots.

#verificando a simetria
par(mfrow=c(1,2))
boxplot(y, main="Boxplot dos diâmetros", ylab="Diâmetro (mm)")
boxplot(e, main="Boxplot dos erros aleatórios", ylab="e (mm)")

4 Ajustando os dados assumindo \(\sigma^2\) conhecido.

Considere \(\sigma^2\) conhecido.

4.1 Estimativa pontual

Forneça uma estimativa pontual para o parâmetro \(\theta\).

#Estimativa pontual para theta
mean(y)

## [1] 39.99922

4.2 Estimativa intervalar

Construa um intervalo de confiança para \(\theta\) assumindo \(\sigma^2\) conhecido e um nível de confiança de 97%.

#Estimativa intervalar para theta
gama      = 0.97
alfa      = 1-gama
z_critico = qnorm(1-alfa/2)
#z_critico = qnorm(alfa/2, lower.tail=F)
erro_padrao_theta = sqrt(sigma2/ n)
round(c(mean(y) - z_critico*erro_padrao_theta,
            mean(y) + z_critico*erro_padrao_theta),3)

## [1] 39.793 40.205

4.3 Função de log-verossimilhança

Represente graficamente a função de log-verossimilhança de \(\theta\).

#função de log-verossimilhança
funcao_log_vero = function(y, theta, sigma2){
  sum(dnorm(y, theta, sqrt(sigma2), log=T))
}

#aplicando a funcao
theta_val   = seq(10,70, l=1000)
log_vero    = NULL
for(j in 1:length(theta_val)){
  log_vero[j] = funcao_log_vero(y, theta_val[j], sigma2)
}

plot(theta_val, log_vero, type="l", bty="n", xlab=expression(theta), ylab="log verossimilhança", lwd=2)
abline(v=theta, lwd=2, col="blue")
abline(v=mean(y), lwd=2, col="red", lty=3)
legend("topright", legend=c(expression("Verdadeiro "*theta), "Estimativa da média amostral"), 
       col=c("blue", "red"), lty=c(1, 3), lwd=2, bty="n")

4.4 Erros aleatórios

Encontre as estimativas pontuais dos erros aleatórios e compare com os verdadeiros.

#Estimando os erros aleatorios
echapeu = y - mean(y)
plot(e, echapeu, lwd=2, bty="n")
lines(e,e, col="red", lwd=2)

5 Ajustando os dados assumindo \(\sigma^2\) desconhecido.

Considere \(\sigma^2\) desconhecido.

5.1 Estimativa pontual

Forneça uma estimativa pontual para o parâmetro \(\theta\) e uma para \(\sigma^2\).

# Estimativas pontuais de theta e sigma^2
theta_chapeu = mean(y)  # Estimador pontual de theta
c(theta_chapeu, theta)

## [1] 39.99922 40.00000

sigma2_chapeu = var(y)  # Estimador pontual de sigma^2
c(sigma2_chapeu, sigma2)

## [1] 8.9281 9.0000

5.2 Estimativa intervalar

Construa um intervalo de confiança para \(\theta\) e outro para \(\sigma^2\) considerando um nível de confiança de 97%.

# Intervalos de confiança (97%)
gama = 0.97
alfa = 1-gama  # Nível de significância

# IC para theta (usando distribuição t de Student)
t_critico = qt(1 - alfa/2, df = n - 1)
erro_padrao_theta = sqrt(sigma2_chapeu / n)
IC_theta = c(theta_chapeu - t_critico * erro_padrao_theta, 
              theta_chapeu + t_critico * erro_padrao_theta)
IC_theta

## [1] 39.79388 40.20456

# IC para sigma^2 (usando distribuição qui-quadrado)
chi_inf = qchisq(alfa/2, df=n-1)
chi_sup = qchisq(1 - alfa/2, df=n-1)
IC_sigma2 = c((n-1) * sigma2_chapeu / chi_sup, (n-1) * sigma2_chapeu / chi_inf)
IC_sigma2

## [1] 8.119698 9.861061

5.3 Função de log-verossimilhança perfilada

Represente graficamente a função de log-verossimilhança de \(\theta\) considerando a estimativa pontual de \(\sigma^2\). E depois represente graficamente a função de log-verossimilhança de \(\sigma^2\) considerando a estimativa pontual de \(\theta\).

#função de log-verossimilhança
funcao_log_vero = function(y, theta, sigma2){
  sum(dnorm(y, theta, sqrt(sigma2), log=T))
}

#função de log-verossimilhança perfilada I
theta_seq   = seq(30,50, l=250)
log_vero    = NULL
for(j in 1:length(theta_seq)){
  log_vero[j] = funcao_log_vero(y, theta_seq[j], sigma2_chapeu) 
}
plot(theta_seq, log_vero, type="l", bty="n", xlab=expression(theta), ylab="log verossimilhança", lwd=2)
abline(v=theta, lwd=2, col="blue")
abline(v=theta_chapeu, lwd=2, col="red", lty=3)
legend("topright", legend=c(expression("Verdadeiro "*theta), "Estimativa da média amostral"), 
       col=c("blue", "red"), lty=c(1, 3), lwd=2, bty="n")

#função de log-verossimilhança perfilada II
sigma2_seq   = seq(7,12, l=250)
log_vero    = NULL
for(j in 1:length(sigma2_seq)){
  log_vero[j] = funcao_log_vero(y, theta_chapeu, sigma2_seq[j]) 
}
plot(sigma2_seq, log_vero, type="l", bty="n", xlab=expression(sigma^2), ylab="log verossimilhança", lwd=2)
abline(v=sigma2, lwd=2, col="blue")
abline(v=sigma2_chapeu, lwd=2, col="red", lty=3)
legend("topright", legend=c(expression("Verdadeiro "*sigma), "Estimativa da média amostral"), 
       col=c("blue", "red"), lty=c(1, 3), lwd=2, bty="n")

5.4 Função de log-verossimilhança

Represente graficamente a função de log-verossimilhança de \(\theta\) e \(\sigma^2\).

#criando valores para os parametros desconhecidos
#theta_seq = seq(20, 60, length = 100)
#sigma2_seq = seq(5,25, length = 100)

#Combina todos os valores de theta_seq com os de sigma2_seq
#Exemplo: se theta_seq tem os valores 1 e 2 e sigma2_seq tem os valores 3 e 4
#os pares de pontos são (1,3), (1,4), (2,3), (2,4) 
#e depois calcula a log-verossimilhanca para cada par de pontos
#log_vero_matrix = outer(theta_seq, sigma2_seq, Vectorize(funcao_log_vero))

log_vero_matrix = matrix(NA, length(theta_seq),length(sigma2_seq))
for(ind1 in 1:length(theta_seq)){
  for(ind2 in 1:length(sigma2_seq)){
    log_vero_matrix[ind1,ind2] = funcao_log_vero(y, theta_seq[ind1], sigma2_seq[ind2])
  }
}

#procurando pelo ponto de máximo
ind = which(log_vero_matrix==max(log_vero_matrix), arr.ind = T)
c(theta_seq[ind[1]], sigma2_seq[ind[2]])

## [1] 39.959839  8.927711

c(theta_chapeu, sigma2_chapeu)

## [1] 39.99922  8.92810

c(theta, sigma2)

## [1] 40  9

#Criando o gráfico 3D
persp_plot = persp(theta_seq, sigma2_seq, log_vero_matrix, theta=30, phi=30, 
                    xlab=expression(theta), ylab=expression(sigma^2), 
                    zlab="Log-Verossimilhança", col="lightblue", border="blue")

#Encontrando o ponto de maximo (espera-se que a funcao de log-verossimilhanca atinja seu valor maximo quando theta e sigma2 assumirem os seus respectivos valores verdadeiros)
max_log_vero = funcao_log_vero(y, theta, sigma2)
true_point = trans3d(theta, sigma2, max_log_vero, persp_plot)
points(true_point$x, true_point$y, pch=19, col="red", cex=1.5)

# Criando o gráfico interativo de superfície
# install.packages("plotly") # Instalar e carregar o pacote plotly, se necessário
library(plotly)

fig <- plot_ly(
  x = theta_seq, 
  y = sigma2_seq, 
  z = t(log_vero_matrix), 
  type = "surface"
)

# Adicionando o ponto verdadeiro (red point)
fig <- fig %>% add_trace(
  x = theta,
  y = sigma2,
  z = max_log_vero,
  type = "scatter3d",
  mode = "markers",
  marker = list(color = "red", size = 5)
)

# Ajustando o layout para os eixos
fig <- fig %>% layout(
  title = "Função de Log-Verossimilhança",
  scene = list(
    xaxis = list(title = "theta"),  
    yaxis = list(title = "sigma2"),
    zaxis = list(title = "Log-Verossimilhança")
  )
)


# Exibindo o gráfico interativo
fig

5.5 Erros aleatórios

Encontre as estimativas pontuais dos erros aleatórios e compare com os verdadeiros.

#Estimando os erros aleatorios
echapeu = y - mean(y)
plot(e, echapeu, lwd=2, bty="n")
lines(e,e, col="red", lwd=2)