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
}