Adaptando a resposta acima 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:

    valores <- c(25,50,75,100,125)
    
    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 <- NULL
        
        for(i in 1:length(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

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 ao invés do round.

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