5

Gurus do R,

Tenho o seguinte data frame (Df) que estabelece a relação entre as variáveis X e Y:

     X     Y
 1  25  2457524
 2  25  2391693
 3  25  2450828
 4  25  2391252
 5  25  2444638
 6  25  2360293
 7  50  4693194
 8  50  4844527
 9  50  4835596
10  50  4878092
11  50  4809226
12  50  4722253
13  75  7142763
14  75  7182769
15  75  7135550
16  75  7173920
17  75  7216871
18  75  7076359
19  100 9496553
20  100 9537788
21  100 9405825
22  100 9439201
23  100 9609870
24  100 9707734
25  125 12031958
26  125 12027037
27  125 11935594
28  125 11930086
29  125 12154132
30  125 12096462
31  150 14298064
32  150 14396607
33  150 13964716
34  150 14221039
35  150 14283992
36  150 14042220

(Note que temos 7 níveis para a variável X com 6 pontos em cada nível)

Se ajustamos um modelo polinomial de 2º grau para estes dados obteremos o seguinte modelo:

Model<-lm(formula = Y ~ X + I(X^2))
print(Model)

Call:
lm(formula = Y ~ X + I(X^2))

Coefficients:
   (Intercept)     X       I(X^2)  
     -26588.12  97310.61   -14.02 

A representação gráfica deste modelo, que mais parece uma reta, é a seguinte:

Relação polinomial ente X e Y

Se quisermos utilizar o modelo para fazer a predição dos valores de "Y" a partir dos valores da variável "X" basta executar esta linha de código:

>predicted.intervals <- predict(Model,data.frame(x=X),interval='confidence',
+ level=0.95)

>predicted.intervals
        fit      lwr      upr
   1   2397413  2315346  2479481
   2   2397413  2315346  2479481
   3   2397413  2315346  2479481
   4   2397413  2315346  2479481
   5   2397413  2315346  2479481
   6   2397413  2315346  2479481
   7   4803887  4753705  4854070
   8   4803887  4753705  4854070
   9   4803887  4753705  4854070
   10  4803887  4753705  4854070
   11  4803887  4753705  4854070
   12  4803887  4753705  4854070
   13  7192834  7137649  7248019
   14  7192834  7137649  7248019
   15  7192834  7137649  7248019
   16  7192834  7137649  7248019
   17  7192834  7137649  7248019
   18  7192834  7137649  7248019
   19  9564252  9509067  9619438
   20  9564252  9509067  9619438
   21  9564252  9509067  9619438
   22  9564252  9509067  9619438
   23  9564252  9509067  9619438
   24  9564252  9509067  9619438
   25 11918144 11867961 11968326
   26 11918144 11867961 11968326
   27 11918144 11867961 11968326
   28 11918144 11867961 11968326
   29 11918144 11867961 11968326
   30 11918144 11867961 11968326
   31 14254507 14172440 14336574
   32 14254507 14172440 14336574
   33 14254507 14172440 14336574
   34 14254507 14172440 14336574
   35 14254507 14172440 14336574
   36 14254507 14172440 14336574

A pergunta que não quer calar:

Qual seria o linha(s) de código para fazermos a predição inversa, ou seja, neste modelo, prever "X" a partir dos dados da variável "Y"? Pesquisando no google já tentei vários pacotes e funções específicas mas infelizmente não obtive sucesso (talvez por falta de familiaridade com as funções). Será que algum de vocês poderiam me ajudar a desvendar este mistério? Grande abraço à todos.

6
  • 2
    Existe algum motivo especial para vc querer fazer isso? Pq vc não ajusta um modelo lm(X ~Y + I(Y^2)) e prevê X a partir de y diretamente? Estatisticamente, fazer isso que você quer é estranho pois X é considerada uma variável observada, não teria pq você querer prevê-la... Commented 17/11/2016 às 20:11
  • 1
    Em adição ao que o Daniel falou, eu ainda rodaria summary(Model), que é uma resposta mais interessante estatisticamente do que print(Model). Com ela tu vai testar as hipóteses dos três coeficientes do teu modelo serem iguais a zero. Dado o gráfico que estou vendo, eu aposto que o termo quadrático não é significante (i.e., p-valor > 0,05). Ou seja, tu tem dados que seguem um modelo linear, não quadrático. Commented 17/11/2016 às 20:25
  • Sim existe !!! Em suma, estou trabalhado com um classe especial dos modelos de regressão denominada de "calibração". Estes modelos são pouco conhecidos no meio acadêmico. Para estes modelos a predição é inversa, ou seja, primeiro você constrói o modelo e depois determina o X. Commented 17/11/2016 às 20:35
  • Olá Marcus!!! É sempre um prazer ver você por aqui...Meu problema não se resume em ver quais variáveis são importantes para o modelo. Isto eu já fiz no inicio da análise!!! Agora preciso prever o valor de "X". Ja verifiquei que o R tem uma função que faz isso denominada de "invest"do pacote "investr". Porém não estou conseguindo fazer dar certo para o modelo em questão. Commented 17/11/2016 às 20:39
  • Boa noite Marcus! A avaliação da significância dos termos já foi feita com o comando 'summary (Model)' e o p-valor não foi significante para o termo quadrático (p=0.237). Por outro lado, devido às “6” repetições para cada nível de X, o teste de falta de ajuste indicou a falta de ajuste (p-valor de 0.0232) para o modelo polinômial de primeira ordem. Mas este não é o foco do problema. Commented 17/11/2016 às 23:32

2 Respostas 2

1

Desconheço alguma função pronta para fazer isso, no entanto este problema pode ser tratado como um problema de otimização.

Queremos encontrar x em um dado intervalo que vai minimizar uma função. A função que pretendo minimizar é:

objetivo <- function(x, k, model){
  df <- data.frame(x = x)
  (k - predict(model, df))^2
}

Dado um valor y = k desejo encontrar um valor x que mais aproxima a predição do modelo desse valor k. No fundo isso quer dizer: me encontre x que quando eu aplicar o modelo eu vou chegar mais perto do y que estou procurando.

Exemplo:

Suponha os seguintes dados:

x <- runif(100)
y <- 10 + x + x^2 + rnorm(n = 100, mean = 0, sd = 0.1)
plot(x, y)

model<-lm(formula = y ~ x + I(x^2))

Call:
lm(formula = y ~ x + I(x^2))

Residuals:
      Min        1Q    Median        3Q       Max 
-0.208921 -0.064506  0.001537  0.061107  0.276347 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.98818    0.03736 267.335  < 2e-16 ***
x            1.05912    0.15697   6.747 1.10e-09 ***
I(x^2)       0.95822    0.14229   6.734 1.17e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09685 on 97 degrees of freedom
Multiple R-squared:  0.9728,    Adjusted R-squared:  0.9723 
F-statistic:  1735 on 2 and 97 DF,  p-value: < 2.2e-16

Agora, usando a função optim podemos obter o que desejamos.

v <- optim(0, objetivo, k = 11, model = model, method = "Brent", lower = min(x), upper = max(x))
v$par
[1] 0.6240072

Assim, obtivemos o valor de x que é o mais provável dado y = k = 11.

Note: a função optim por padrão usa o método Nelder-Mead, no entanto ela solta um warning quando a usei neste exemplo. Por isso troquei para o método "Brent" que exige valores mínimos e máximos para o valor estimado de x. Isso pode ser bom, uma vez que estimando polinômios, é possível que existam outros valores que minimizem esta função.

2
  • Daniel. Obrigado pela sua contribuição. Mas não era exatamente o que pretendia. Depois faça uma busca no google pelo pacote 'investr' e veja o que consegue fazer com estes dados... abraços. Commented 18/11/2016 às 11:11
  • 1
    @WeidsonC.deSouza o que ficou faltando? Com o investr vc vai conseguir os mesmos valores do que com a solução acima. Enfim, sugiro que você especifique melhor esses detalhes na sua pergunta para que alguém possa responder adequadamente. Commented 18/11/2016 às 11:21
1

Adaptando as respostas anteriores com dados mais parecidos com o da pergunta e entendendo que o resultado de x deve estar nas 7 categorias ("níveis") de valores descritos:

v <- c(25,50,75,100,125)
valores <- sort(rep.int(v,6))

set.seed(13)
x <- sample(valores, 36, replace = T)
y <- -26500 + 97000*x + -14*x^2 + rnorm(n = 36, mean = 0, sd = 100000)
plot(x, y, pch = 21)

O Modelo pode se dar da forma abaixo:

model <-lm(formula = y ~ x + I(x^2))

Bem como a função objetivo para encontrar o valor predito:

objetivo <- function(x, k, model){
  df <- data.frame(x = x)
  (k - predict(model, df))^2
}

Dessa forma fazendo um loop com 'optim':

    aux_1 <- numeric(length(y))

    for(i in seq_along(y)){ 
    aux_1[i] <- optim(0, objetivo, k = y[i], model = model, method = "Brent",
 lower = min(x), upper = max(x))$par
    }

    # Tomando os resultados mais próximos:

    result <- round(aux_1,0)
    result
[1] 100  51  50  27 125  25  76 100 124  25  99 124 125  76  75  50  50  74
[19] 125 100  25  74 100  77  26 102  25  74  49  98  75  76 125 125  75  25

Se quisermos apenas os valores presentes em (25, 50, 75, 100, 125) pode se fazer a função abaixo:

extract <- function(x,y){
y <- sort(y, decreasing = T)

for(i in seq_along(y)){
x[x/y[i] < 1.25 & x/y[i] >= 0.98] = y[i]}
return(x)
}

result <- extract(result,v)
result
[1] 100  50  50  25 125  25  75 100 125  25 100 125 125  75  75  50  50  75
[19] 125 100  25  75 100  75  25 100  25  75  50 100  75  75 125 125  75  25

Para valores de 'x' mais próximos uns dos outros pode ser preciso refinar mais a função que busca os valores mais próximos.

EDIT_1 (27/11/17, 23:35): estava gerando um vetor aux_2 que foi desnecessário na versão final publicada.

EDIT_2 (15/06/18, 20:20): Refiz o vetor de valores para ficar mais similar ao do data.frame apresentado como exemplo e acrescentei as boas sugestões do comentário do Rui Barradas abaixo. Ainda, como essa resposta carecia de aproximação para ficar apenas com valores do vetor (25, 50, 75, 100, 125), criei a função 'extract'. No entanto, é uma função específica para os valores desse exemplo, pode se ajustar os para casos mais gerais.

2
  • Porque não aux_1 <- numeric(length(y))? É mais eficaz. E ainda for(i in seq_along(y)), para evitar possíveis erros quando y é um vetor de comprimento zero. Eu sei que não é este o caso mas é um bom hábito substituir 1:length(.) por seq_along(.). Commented 14/06/2018 às 20:40
  • Muito obrigado @RuiBarradas, são sugestões de boas práticas, fiz edição para incorporá-las.
    – VituMaia
    Commented 15/06/2018 às 22:21

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .