3

Estoy intentando ajustar un polinomio cuadrático a mis datos usando el siguiente código:

#Polynimial regression for y ~ x
model <- lm(y ~ x + I(x^2))
summary(model)

#Box and whisker plot + polynomial
boxplot(y ~ x,
        col=c("white","lightgray"), ylab= "y", xlab= "x", dat)
means <- tapply(y,x,mean)
points(means,col="red",pch=18)
predicted.intervals <- predict(model,data.frame(x=x),interval='confidence',
                               level=0.99)
lines(x,predicted.intervals[,1],col='green',lwd=3)
lines(x,predicted.intervals[,2],col='black',lwd=1)
lines(x,predicted.intervals[,3],col='black',lwd=1)

Enlace a datos

Enlace datos B

El caso es que cuando corro el programa me aparece el gráfico de cajas, los puntos rojos representando las medias y la línea verde del polinomio ajustado a los datos. No obstante, también aparece una extraña línea recta que une las medias de los niveles 1 y 11 que no tengo ni idea de dónde proviene. Aquí va el gráfico:

introducir la descripción de la imagen aquí

He ajustado más veces polinomios a mis datos en regresiones no lineales, pero esto nunca me había pasado.

¿Alguna solución?

Edición 1:

Gráfico obtenido con datos B.

introducir la descripción de la imagen aquí

Finalmente, he conseguido ajustar el polinomio a los datos. El código utilizado es el siguiente:

#Factorizo la variable x
x <- as.factor(x)
#Vuelvo a transformar la variable a numérica
x <- as.numeric(x)

#Regresión cuadrática
model <- lm(y ~ x + I(x^2))
summary(model)

#Ajuste del polinomio
boxplot(y ~ x,
        col=c("white","lightgray"), ylab= "y", xlab= "x", dat)
means <- tapply(y,x,mean)
points(means,col="red",pch=18)
predicted.intervals <- predict(model,data.frame(x=x),interval='confidence',
                               level=0.99)
lines(x,predicted.intervals[,1],col='green',lwd=3)

El resultado para los datos B es este:

introducir la descripción de la imagen aquí

No obstante, me quedan algunas dudas:

Mi variable x comprende valores de 0 a 1 (11 niveles en pasos de 0.1)

  1. ¿Por qué he tenido que factorizar mi variable original x, después volver a transformar la variable a numérica (adoptando esta valores discretos entre 0 y 11? Solo así consigo ajustar el polinomio a los datos, pero la regresión está ejecutada sobre valores numéricos de 1 a 11.

  2. ¿Por qué en el caso del "conjunto de datos B" debo utilizar

    lines(x,predicted.intervals[,1],col='green',lwd=3)

... mientras que en el caso del "conjunto de datos A" debo utilizar

lines(predicted.intervals[,1],col='green',lwd=3)

?

10
  • pyring, podrías añadir los datos que generan la gráfica? si no son demasiado grandes puedes hacer dput(x) y dput(y) y agregar esa salida a tu pregunta. Saludos. el 21 feb. 2018 a las 18:44
  • hola @Patricio Moracho. El problema es que el conjunto de datos es enorme (son datos de simulaciones): 100000 valores cada variable.
    – pyring
    el 22 feb. 2018 a las 17:36
  • Te entiendo, pero al menos yo, sin ver los datos me cuesta entender que es lo que puede estar pasando. Si se me ocurre algo te comento. el 22 feb. 2018 a las 17:43
  • No tendría inconveniente en compartir un archivo csv. Creo que una posibilidad sería por medio del chat de stackoverflow. Un saludo.
    – pyring
    el 22 feb. 2018 a las 23:24
  • No tengo mucho experiencia en el chat, pero si quieres lo intentamos, aunque sin no tienes algún problema de confidencialidad con los datos, tal vez lo mejor es algún dropbox o similar y un enlace en la pregunta así incluso tienes mas chances de respuesta. Saludos el 23 feb. 2018 a las 1:31

1 respuesta 1

4

Vamos al primer problema, los datos del conjunto A. Voy a trabajar solo con la regresión, veamos:

datA <- read.table(file = "datA.csv", header = T, sep=",", stringsAsFactors = F, dec = ".")

model <- lm(datA$y ~ datA$x + I(datA$x^2))
datA$fitted <- predict(model,data.frame(x=datA$x))

Ahora vamos a graficar solo los puntos y la curva para la regresión:

plot(datA$x, datA$fitted)
lines(datA$x, datA$fitted)

El resultado seguramente te es conocido:

introducir la descripción de la imagen aquí

Esto es explicable por que datA no está ordenada por x, al dibujar las lineas a partir de los puntos x y fitted eventualmente podríamos tener un punto (1, ?) y luego un (0, ?), por lo que volveríamos al origen haciendo que la gráfica sea "circular". Para resolver esto, simplemente ordenamos por x:

datA.ordered <- datA[order(datA$x),]

plot(datA.ordered$x, datA.ordered$fitted)
lines(datA.ordered$x, datA.ordered$fitted)

Ahora el resultado es más acorde a lo buscado:

introducir la descripción de la imagen aquí

¿Es esta la solución y explicación del problema? Sí y No. Veamos, si incorporamos está gráfica al boxplot

datA <- read.table(file = "C:/Users/pmoracho/Downloads/datA.csv", header = T, sep=",", stringsAsFactors = F, dec = ".")

model <- lm(datA$y ~ datA$x + I(datA$x^2))
boxplot(datA$y ~ datA$x,
        col=c("white","lightgray"), ylab= "y", xlab= "x", datA)

means <- tapply(datA$y,datA$x,mean)
points(means,col="red",pch=18)
datA$fitted <- predict(model,data.frame(x=datA$x))
datA.ordered <- datA[order(datA$x),]
lines(datA.ordered$x, datA.ordered$fitted,col='green',lwd=3)

podemos ver este resultado:

introducir la descripción de la imagen aquí

¿Qué podemos notar? la curva original se "comprimió" entre los valores de 0 y 1, la explicación es que no son compatibles los sistemas de coordenadas, esto por que el boxplot considera los valores de x como variables discretas, sin embargo nuestros valores de datA$x no lo son, y aquí es dónde entra el uso de factor(), así:

lines(factor(datA.ordered$x), datA.ordered$fitted,col='green',lwd=3)

Ahora sí los datos de x son consistentes con las x del boxplot:

introducir la descripción de la imagen aquí

Esto explica y me ha funcionado para los dos juegos de datos, no publico los resultados para no hacer más larga la respuesta. De todas formas te preguntarás por que originalmente han sido distintos los comportamientos de los dos juegos de datos, la explicación es sencilla, el conjunto A está desordenado y el B está ordenado (siempre hablando del valor de x).

Te recomiendo además que uses spline() para dibujar estas curvas, te evitaría tener que ordenar previamente y la otra ventaja más importante es que no le pasa a line() el juego completo de datos sino los puntos mínimos para interporlar la línea en la gráfica:

lines(spline(factor(datA$x), datA$fitted),col='green',lwd=3)
4
  • Gracias @Patricio Moracho. ¿Llegas a ver por qué con unos datos diferentes y utilizando el mismo código, no se consigue una línea ajustada a los datos sino una extraña línea recta? He añadido enlace a datos B.
    – pyring
    el 23 feb. 2018 a las 19:59
  • He añadido también el gráfico obtenido con el mismo código para datos B.
    – pyring
    el 23 feb. 2018 a las 23:07
  • @pyring, rehíce completamente la respuesta, ahora creo que está la explicación del problema, antes solo teníamos soluciones parciales. el 28 feb. 2018 a las 14:25
  • Fenómeno, de gran ayuda.
    – pyring
    el 28 feb. 2018 a las 16:36

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.