Modelo de Regressão Linear Simples

Nesta aula computacional, verificaremos empiricamente as propriedades estudadas sobre os estimadores vistos em sala de aula.

Modelo

Considere o modelo de regressão linear simples: \[Y_i = \beta_0 + \beta_1 x_i + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2), \quad i=1,\ldots,n.\] Suponha \(\beta_0\; , \; \beta_1\; , \; \sigma^2\; , \; \mu_1 = \beta_0 + \beta_1 x_1, \; \ldots\; , \; \mu_n = \beta_0 + \beta_1 x_n\) desconhecidos. Considere os seguintes estimadores não viesados para estes parâmetros: \[\hat{\beta}_0 = \bar{Y}-\bar{x}\hat{\beta}_1,\] \[\hat{\beta}_1 = \frac{\sum_{i=1}^n{Y_ix_i - n\bar{x}\bar{Y}}}{\sum_{i=1}^n{(x_i-\bar{x})^2}},\] \[\hat{\sigma}^2 = \frac{\sum_{i=1}^n{(Y_i-\hat{Y}_i)^2}}{n-2},\] \[\hat{\mu}_i = \hat{y}_i = \hat{\beta}_0 + x_i \hat{\beta}_1, \quad i=1, \ldots, n.\]

Estimadores não tendenciosos

Verifique empiricamente que \[\hat{\beta_0} \sim N\left(\beta_0 \quad , \quad \sigma^2\left(\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^n{(x_i-\bar{x})^2}}\right)\right),\]

\[\hat{\beta_1} \sim N\left(\beta_1 \quad , \quad \displaystyle\frac{\sigma^2}{\sum_{i=1}^n{(x_i-\bar{x})^2}}\right),\] \[E[\hat{\sigma}^2] = \sigma^2,\]

\[\hat{Y}_i \sim N\left( \beta_0 + \beta_1 x_i \quad , \quad \sigma^2\left(\frac{1}{n}+\frac{(x_i-\bar{x})^2}{\sum_{i=1}^n{(x_i-\bar{x})^2}}\right)\right).\]

Motivação

Para verificar empiricamente as distribuições dos estimadores no modelo de regressão linear, precisamos observar como esses estimadores se comportam quando aplicamos o mesmo procedimento a diferentes conjuntos de dados. Esse tipo de abordagem é chamado de simulação, e é uma ferramenta poderosa para compreender, na prática, propriedades teóricas como a distribuição amostral de um estimador.

A ideia é a seguinte: vamos gerar vários conjuntos de dados artificialmente, sempre seguindo o mesmo modelo probabilístico, e para cada conjunto vamos calcular os estimadores dos parâmetros desconhecidos: \(\hat{\beta}_0\), \(\hat{\beta}_1\), \(\hat{\sigma}^2\), \(\hat{\mu}_1\), \(\ldots\), \(\hat{\mu}_n\). Assim, ao final, teremos uma amostra de estimadores — por exemplo, 2000 valores distintos de \(\hat{\beta}_1\), cada um correspondente a uma simulação — e poderemos estudar, por exemplo, a média, a variância e a forma da distribuição empírica desses valores.

Para tornar isso mais concreto, considere o seguinte exemplo:

Uma pesquisadora deseja investigar se existe uma relação entre o nível de glicose e a quantidade de cafeína ingerida 1h antes do exame. Para isso, \(n\) pessoas participarão da pesquisa por \(n_{rep}\) dias. Cada pessoa \(i\) beberá sempre a mesma quantidade de cafeína 1h antes de realizar o exame, sendo essa quantidade denotada por \(x_i\). O nível de glicose observado para a \(i\)-ésima pessoa no \(j\)-ésimo dia, será denotado por \(Y_{ij}\). Dessa forma, o nosso estudo empírico será equivalente ao seguinte modelo: \[Y_{ij} = \beta_0 + \beta_1 x_i + e_{ij}, \quad e_{ij} \sim N(0, \sigma^2), \quad i=1,\ldots,n, \quad j=1, \ldots, n_{rep}.\]

Nesse modelo:

\(\beta_0\) representa o nível médio de glicose das pessoas quando não há consumo de cafeína (\(x_i = 0\)),
\(\beta_1\) representa o efeito médio da cafeína sobre a glicose,
\(\sigma^2\) é a variância da glicose após remover o efeito da cafeína, assumida constante.

Para investigar empiricamente as distribuições dos estimadores, simularemos \(n_{rep}=2000\) conjuntos de dados com \(n = 100\) pessoas cada, considerando os seguintes valores reais dos parâmetros: \[\beta_0 = 10, \quad \beta_1 = 2, \quad \sigma^2 = 2, \quad x_i \stackrel{iid}{\sim} N(5,1).\] Ou seja,

Fixaremos os valores dos parâmetros \(n\), \(n_{rep}\), \(\beta_0\), \(\beta_1\), \(\sigma^2\) e geraremos as covariáveis \(x_i\) para \(i = 1, \ldots, n\) conforme descrito acima; (Note que estas medidas não variam com a replicação e por isso podem ser realizadas uma única vez!)
Para cada replicação \(j\):
1. geraremos o erro \(e_{ij} \stackrel{iid}{\sim} N(0,\sigma^2)\) para \(i = 1, \ldots, n\);
2. obteremos a variável resposta \(Y_{ij} = \beta_0 + \beta_1 x_i + e_{ij}\) para \(i = 1, \ldots, n\);
3. calcularemos os estimadores \(\hat{\beta}_0\), \(\hat{\beta}_1\), \(\hat{\sigma}^2\), \(\hat{Y}_1, \ldots, \hat{Y}_n\) e guardaremos estes valores.

Após repetir a 2ª etapa \(n_{rep}\) vezes, teremos amostras empíricas desses estimadores e poderemos verificar se:

eles são centrados em torno dos valores verdadeiros (verificando o não viés),
se as variâncias empíricas ficam próximas das variâncias teóricas,
suas distribuições empíricas se aproximam das distribuições teóricas fornecidas.

Solução

# Definindo parâmetros reais
n_rep   = 2000    # Número de replicacoes
n       = 100     # Tamanho da amostra
beta0   = 10      # Intercepto
beta1   = 2       # Coeficiente de regressão
sigma2  = 2       # variancia

# Vetores para armazenar as estimativas dos parametros
beta0chapeu = NULL
beta1chapeu = NULL
sigma2chapeu = NULL
sigma2chapeu_viesado = NULL
ychapeu       = matrix(NA, n, n_rep)

# Gerando a covariavel e a media mu
x       = rnorm(n, mean = 5, sd = 1) 
mu      = beta0 + beta1 * x

# Calculando algumas constantes
xbarra  = mean(x)
Sxx     = sum((x-xbarra)^2)

# Simulação de dados e ajuste do modelo
for (ite in 1:n_rep) {
    # Gerando erro aleatório
    e       = rnorm(n, mean = 0, sd = sqrt(sigma2)) 
    # Gerando a variavel resposta
    y       = mu + e
  
  # Ajustando o modelo de regressão linear simples
  ajuste = lm(y ~ x)
  # Armazenando as estimativas dos parametros
  beta0chapeu[ite]  = ajuste$coefficients[1]  # Estimativa de beta0
  beta1chapeu[ite]  = ajuste$coefficients[2]  # Estimativa de beta1
  sigma2chapeu[ite] = summary(ajuste)$sigma^2 # Estimativa de sigma2
  sigma2chapeu_viesado[ite] = sum((ajuste$residuals)^2)/n # Estimativa de sigma2

  ychapeu[,ite]     = ajuste$fitted.values

}

#####################################
# Verificando a media dos estimadores
#####################################
media_beta0chapeu  = mean(beta0chapeu)
media_beta1chapeu  = mean(beta1chapeu)
media_sigma2chapeu = mean(sigma2chapeu)
media_sigma2chapeu_viesado = mean(sigma2chapeu_viesado)

# Resultados
cat("Valor médio da amostra de beta0chapeu:", round(media_beta0chapeu,4), "Valor verdadeiro de beta0:", beta0, "\n")

## Valor médio da amostra de beta0chapeu: 10.0244 Valor verdadeiro de beta0: 10

cat("Valor médio da amostra de beta1chapeu:", round(media_beta1chapeu,4), "Valor verdadeiro de beta1:", beta1, "\n")

## Valor médio da amostra de beta1chapeu: 1.9964 Valor verdadeiro de beta1: 2

cat("Valor médio da amostra de sigma2chapeu:", round(media_sigma2chapeu,4), "Valor verdadeiro de sigma2:", sigma2, "\n")

## Valor médio da amostra de sigma2chapeu: 2.0104 Valor verdadeiro de sigma2: 2

cat("Valor médio da amostra de sigma2chapeu_viesado:", round(media_sigma2chapeu_viesado,4), "Valor verdadeiro de sigma2:", sigma2, "\n")

## Valor médio da amostra de sigma2chapeu_viesado: 1.9702 Valor verdadeiro de sigma2: 2

############################################
# Verificando a variancia dos estimadores
############################################

#obtendo a variancia empirica
var_beta0chapeu  = var(beta0chapeu)
var_beta1chapeu  = var(beta1chapeu)

#calculando a variancia teorica
var_t_beta0chapeu   = sigma2*(1/n+mean(x)^2/Sxx)
var_t_beta1chapeu   = sigma2/Sxx

# Resultados
cat("Variância empírica de beta0chapeu:", round(var_beta0chapeu,4), "Variância teórica de beta0chapeu:", round(var_t_beta0chapeu,4), "\n")

## Variância empírica de beta0chapeu: 0.5077 Variância teórica de beta0chapeu: 0.5183

cat("Variância empírica de beta1chapeu:", round(var_beta1chapeu,4), "Variância teórica de beta1chapeu:", round(var_t_beta1chapeu,4), "\n")

## Variância empírica de beta1chapeu: 0.019 Variância teórica de beta1chapeu: 0.0191

#Curiosidade:
#Mostrando como obtenho a variancia de beta0chapeu usando a funcao lm do R para o ajuste realizado na ultima replicacao
resultado = summary(ajuste)
coef(resultado)[1, "Std. Error"]^2

## [1] 0.452161

resultado$coefficients[1, "Std. Error"]^2

## [1] 0.452161

sigma2chapeu[n_rep]*(1/n+mean(x)^2/sum((x-mean(x))^2))

## [1] 0.452161

#################################################################
# Verificando a media e a variancia de ychapeu
#################################################################

media_ychapeu = NULL  #vetor de medias empiricas
var_ychapeu = NULL    #vetor de variancias empiricas
var_t_ychapeu = NULL  #vetor de variancias teoricas
for(i in 1:n){  
  media_ychapeu[i]  = mean(ychapeu[i,])
  var_ychapeu[i]    = var(ychapeu[i,]) 
  var_t_ychapeu[i]  = sigma2*(1/n+(x[i]-xbarra)^2/Sxx)
}

#comparando as estimativas das medias mu com as medias mu
plot(mu,media_ychapeu, bty="n")
lines(mu,mu,col="red", lwd=2, lty=3)

#comparando as variancias empiricas com as teoricas de ychapeu
plot(var_ychapeu, var_t_ychapeu, bty="n")
lines(var_ychapeu,var_ychapeu,col="red", lwd=2, lty=3)

#################################################################
#analisando as distribuicoes empiricas
#################################################################
hist(beta0chapeu, freq=F, main="", ylab="densidade", xlab=expression(hat(beta)[0]))
var_beta0chapeu_normal = sigma2*(1/n + xbarra^2/Sxx)
curve(dnorm(x, beta0, sqrt(var_beta0chapeu_normal)), add=T, col="red", lwd=2)

hist(beta1chapeu, freq=F, main="", ylab="densidade", xlab=expression(hat(beta)[1]))
var_beta1chapeu_normal = sigma2/Sxx
curve(dnorm(x, beta1, sqrt(var_beta1chapeu_normal)), add=T, col="red", lwd=2)

par(mfrow=c(3,3))
i = sort(sample(1:n, 9, replace=F))
for(inds in 1:9){
  hist(ychapeu[i[inds],], freq=F, main="", ylab="densidade", xlab=bquote(mu[.(i[inds])]))
  curve(dnorm(x, mu[i[inds]], sqrt(var_t_ychapeu[i[inds]])), add=T, col="red", lwd=2)
}

Demais propriedades

Verifique empiricamente que

\(\sum_{i=1}^n{\hat{e}_i} = 0, \quad \hat{e}_i = Y_i - \hat{Y}_i\)
\(\sum_{i=1}^n{x_i \hat{e}_i} = 0\)
\(\sum_{i=1}^n{\hat{Y}_i \hat{e}_i} = 0\)
\((\bar{x}, \bar{Y})\) é um ponto na regra ajustada (para isso plote de vermelho esse ponto no gráfico de dispersão com a reta ajustada em uma das simulações feitas)

Solução

# Vetores para armazenar as estimativas dos parametros
echapeu       = matrix(NA, n, n_rep)
ychapeu       = matrix(NA, n, n_rep)

# para armazenar as somas
soma_e = soma_xe = soma_ye = NULL

# Simulação de dados e ajuste do modelo
for (ite in 1:n_rep) {
    # Gerando erro aleatório
    e       = rnorm(n, mean = 0, sd = sqrt(sigma2)) 
    # Gerando a variavel resposta
    y       = mu + e
  
  # Ajustando o modelo de regressão linear simples
  ajuste = lm(y ~ x)


  ychapeu[,ite]     = ajuste$fitted.values
  echapeu[,ite]     = y-ychapeu[,ite]
  soma_e[ite]       = sum(echapeu[,ite])
  soma_xe[ite]      = sum(x*echapeu[,ite])
  soma_ye[ite]      = sum(ychapeu[,ite]*echapeu[,ite])
  
}

# Resultados 
cat("Soma de echapeu:", round(summary(soma_e),10), "\n")

## Soma de echapeu: 0 0 0 0 0 0

cat("Soma de xchapeu*echapeu:", round(summary(soma_xe),10), "\n")

## Soma de xchapeu*echapeu: 0 0 0 0 0 0

cat("Soma de ychapeu*echapeu:", round(summary(soma_ye),10), "\n")

## Soma de ychapeu*echapeu: 0 0 0 0 0 0

plot(x,y,main="Gráfico de dispersão", cex.lab=1.4, cex.axis=1.4)
lines(sort(x), beta0chapeu[n_rep]+beta1chapeu[n_rep]*sort(x), col="purple", lwd=2)
points(mean(x), mean(y), col="red", pch=19, cex=1.5)

Aula de Modelos Lineares I

Patrícia Lusié Velozo da Costa

2025-04-11

Modelo

Estimadores não tendenciosos

Motivação

Solução

Demais propriedades

Solução