Traduções desta página:

Ferramentas do usuário

Ferramentas do site


05_curso_antigo:r2017:alunos:trabalho_final:iaraottoni:rotina

Trabalho Final

compara_coefs <- function (data,plota)
  {
## Retira-se os possíveis NA do vetor e cria-se um objeto para os novos vetores
nadata <- na.omit(data)

x <- data$x ## Atribui-se os dados de x a um objeto
y <- data$y ## Atribui-se os dados de y a um objeto
expo <- lm(y~x)  ## Aplica-se a regressão nos termos desejados
expo ## Checa-se os dados obtidos

## Utiliza-se a função Segmented para identificar o ponto de inflexão
library(segmented) ## Ativa-se o pacote Segmented
seg <- segmented(expo, seg.Z = ~x) ##  Segmented é usado em função de x
seg ## Pode-se identificar o ponto de inflexão estimado

## Será necessário retirar o ponto de inflexão para separar os dados
seg ## Checa-se os dados
## O ponto de inflexão se encontra no dado de psi, pontanto:
seg$psi ## Identifica-se onde estão os pontos
coef <- seg$psi ## E atribui-se a um objeto para que os dados desejados sejam extraídos
psi <- coef [1,2] ## Extrai-se o ponto de inflexão
psi ## Ponto de inflexão extraído dos dados

## Agora será necessário separar os dados que estão antes e depois do psi
## Utilizando um indexação simples será possível separar os dados
data1 <- data[data$x<=psi,]  ## Separa-se os dados menores ou iguais ao ponto de inflexão
data1 ## Checa-se os dados
data2 <- data[data$x>psi,] ## Separa-se os dados maiores que ponto de inflexão
data2 ## Checa-se os dados

## Portanto, agora temos x e y para os dados anteriores ao ponto de inflexão e, x e y para os 
## posteriores ao ponto de inflexão
## data1: x, y
## data2: x, y
## É necessário então comparar as inclinações das retas dos conjuntos 1 e 2
## Para isso serão criadas duas regressões
lm1 <- lm (y~x, data1) ## Regressão para o primeiro conjunto de dados, já aplicando a um objeto
lm1 ## Checa-se os dados
lm2 <- lm (y~x, data2) ## Regressão para o segundo conjunto de dados, já aplicando a um objeto
lm2 ## Checa-se os dados

## Os coeficientes de inclinação das retas:
summary(lm1)$coefficients ## Checa-se os dados
summary(lm2)$coefficients  ## Checa-se os dados
s1 <- summary(lm1)$coefficients ## Atribui-se um objeto para que se possa extrair os coeficientes para o conjunto de dados 1
s2 <- summary(lm2)$coefficients ## Atribui-se um objeto para que se possa extrair os coeficientes para o conjunto de dados 2

## Para comparar as inclinações das retas será utilizada uma fórmula apresentada em: 
  ## Paternoster, R., Brame, R., Mazerolle, P., & Piquero, A. R. 1998. Using the Correct Statistical Test for the Equality of Regression Coefficients. Criminology, 36(4), 859–866
## Os componetes da fórmula foram separados para que diminua o erro.

## É necessário subtrair os coeficientes de inclinação e encontra a diferença entre eles
b1 <- s1[2,1] ## Coeficiente da reta 1 (<=psi)
b2 <- s2[2,1] ## Coeficiente da reta 2 (>psi)
difb <- (b2-b1) ## Diferença entre as inclinações
difb ## Checando

## Para a fórmula também e necessário extrair o Standard Error
erro1 <- s1[2,2] ## Standard Error da reta 1 (<=psi)
erro1 ## Checando
erro2 <- s2[2,2] ## Standard Error da reta 2 (>psi)
erro2 ## Checando

## Eleva-se o erro ao quadrado para as duas retas
std1 <- erro1^2 ## Standard Error elevado ao quadrado para a reta 1 (<=psi)
std1  ## Checando
std2 <- erro2^2 ## Standard Error elevado ao quadrado para a reta 2 (>psi)
std2 ## Checando

div <- (std1+std2)^1/2 ## Soma-se os erros ao quadrado e retira-se a raiz
div ## Checando

## A fórmula final para o cálculo da inclinação apresenta-se como 
z <- difb/div ## Diferença entre as inclinações encontrada e atribuída ao objeto z
z ## Checando

## Se plota=="RTS" os pontos das duas retas são plotados separadamente
if (plota=="RTS") ## Condição RTS para que sejam produzidos dois gráficos separados com as regressões
{
  par(mfrow=c(1,2)) ## Condição para se emparelhar os gráficos
  plot(y~x,data1) ## Gráfico para o primeiro conjunto de dados
  abline(lm1, col=2) ## Linha de tendência dos dados da data1
  plot(y~x,data2) ## Gráfico para o segundo conjunto de dados
  abline(lm2, col=4) ## Linha de tendência dos dados da data2
}

if (plota=="PSI")  ## Condição RTS para que seja produzido um gráfico com a linha de tendência considerando o ponto de inflexão
{
  ## Para plotar os dados completos com o ponto de inflexão:
  par(mfrow=c(1,1)) ## Informa que será plotado apenas um gráfico
  plot(x,y) ## Plota-se os dados para o conjunto de dados total
  slope(seg) ## Informa os valores das regressões separadamente
  plot(seg, add = T, col=2) ## Plota a linha de tendência considerando o ponto de inflexão
}

print(paste('Ponto de Inflexão = ',psi,sep = '')) ## Retorna o ponto de inflexão
print(paste('b1 = ',b1,sep = '')) ## Retorna a inclinação da reta 1
print(paste('b2 = ',b2,sep = '')) ## Retorna a inclinação da reta 2 
print(paste('z = ',z,sep = '')) ## Retorna z
}
05_curso_antigo/r2017/alunos/trabalho_final/iaraottoni/rotina.txt · Última modificação: 2020/08/12 06:04 (edição externa)