Quiero leer...
Introducción de los datos en R commander
Los datos con los que trabajaremos se encuentran a continuación y están organizados en dos columnas con el nombre de X e Y.
X Y 97.6470588235 5.98290598291 92.9411764706 14.5299145299 90 22.2222222222 85.2941176471 29.9145299145 80 37.6068376068 74.1176470588 46.1538461538 68.8235294118 52.9914529915 64.1176470588 58.9743589744 58.2352941176 64.9572649573 52.3529411765 70.9401709402 47.0588235294 76.0683760684 41.7647058824 81.1965811966 36.4705882353 84.6153846154 31.1764705882 88.8888888889 25.8823529412 92.3076923077 19.4117647059 95.7264957265 77.6470588235 41.0256410256 83.5294117647 33.3333333333 96.4705882353 9.40170940171 88.8235294118 24.7863247863 100 0 96.186440678 7.27272727273 91.9491525424 15.7575757576 88.9830508475 21.8181818182 85.1694915254 28.4848484848 80.9322033898 36.3636363636 78.813559322 39.3939393939 74.5762711864 46.0606060606 70.3389830508 52.1212121212 64.406779661 61.2121212121 56.7796610169 70.303030303 49.5762711864 78.1818181818 42.7966101695 84.8484848485 36.8644067797 89.696969697 31.7796610169 93.3333333333 26.2711864407 97.5757575758 21.186440678 100 61.4406779661 64.8484848485 68.2203389831 55.7575757576 53.813559322 73.3333333333
Seleccionamos toda la tabla anterior, la copiamos (Ctrl+C), y la cargamos en R commander siguiendo la siguiente ruta:
Y seleccionamos todo como está en la imagen anterior, dando por nombre a la tabla como Datos, y seleccionando que viene del Portapapeles:
Si no ha habido problemas, tendremos una tabla con 40 filas y 2 columnas.
Generación del modelo polinomial
En R commander seguimos la siguiente ruta:
A continuación ponemos el nombre del modelo (en nuestro caso EJEMPLO), y añadimos la ecuación polinomial, en la que nuestra variable dependiente (en nuestro caso Y) depende de la variable independiente (en nuestro caso X). Tendréis que poner en el modelo el nombre de las variables dependiente e independiente, ya que nosotros hemos usado X e Y para simplificar. La fórmula, en nuestro caso cuadrática o de segundo grado se expresa de la siguiente manera:
Y~X+I(X^2)
Es muy importante poner i mayúscula (I) y paréntesis para los grados de magnitud superiores a 1.
[sourcecode lang=»r»] EJEMPLO <- lm(Y ~ X + I(X^2), data=Datos)> summary(EJEMPLO)</span>
Call:
lm(formula = Y ~ X + I(X^2), data = Datos)
Residuals:
Min 1Q Median 3Q Max
-3.2630 -1.3456 0.1563 1.3596 3.2895
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 101.318990 1.971593 51.389 <2e-16 ***
X -0.007965 0.070725 -0.113 0.911
I(X^2) -0.009886 0.000573 -17.255 <2e-16 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Residual standard error: 1.825 on 37 degrees of freedom
Multiple R-squared: 0.9964, Adjusted R-squared: 0.9962
F-statistic: 5127 on 2 and 37 DF, p-value: < 2.2e-16
[/sourcecode]
Observamos que la regresión cuadrática es altamente significativa (p < 0.05) y que el valor ajustado de R^2 es 0.9962.
Intervalos de confianza del modelo
En R Commander seguimos la siguiente ruta, asegurándonos tener seleccionado el modelo EJEMPLO:
Y aquí seleccionamos el nivel de confianza, que por defecto viene en 0.95. Quien desee hacerlo al 99%, simplemente escriban 0.99.
> Confint(EJEMPLO, level=0.95) Estimate 2.5 % 97.5 % (Intercept) 101.318989988 97.32416422 105.313815751 X -0.007965000 -0.15126740 0.135337399 I(X^2) -0.009886546 -0.01104748 -0.008725612
Automáticamente te aparecen los coeficientes y los intervalos para la función cuadrática (columna de la izquierda, con nombre Estimate), el intervalo de confianza inferior al 95% (en la columna del medio, con nombre 2.5%) y el intervalo de confianza superior al 95% (columna de la derecha, con nombre 97.5%)
Representación gráfica de la función cuadrática y los intervalos de confianza al 95%
Paquete necesario: ggplot2
¿Cómo instalar y cargar nuevos paquetes en R Commander de un modo muy sencillo?
A continuación, para obtener el gráfico deseado, introducimos y ejecutamos el siguiente comando:
qplot(X, Y, data=Datos, geom=c("point", "smooth"), method="lm", formula= y ~ poly(x, 2))
Lo resaltado en negrita son los nombres de las variables (en nuestro caso X e Y, y la tabla de datos, en nuestro caso Datos). El «lm» es linear model y la fórmula es polinomial de segundo grado (por eso aparece un 2). Atención que la x (minúscula) de la fórmula no es la X de nuestros datos.
El gráfico resultante es el siguiente:
La línea azul es la función cuadrática mientras que el sombreado gris es el intervalo de confianza al 95%.
Ahora bien, ¿cómo hacemos para cambiar el intervalo de confianza? Pongamos que por ejemplo lo queremos al 99%. Habría que añadir lo siguiente al código previo:
level=0.99
qplot(X, Y, data=Datos, geom=c("point", "smooth"), method="lm", formula= y ~ poly(x, 2),level=0.99)
Predicción de valores de Y y su intervalo de confianza
Ahora bien, conociendo ya las funciones cuadráticas de la estimación y de los intervalos de confianza, si X = 85 ¿cuánto vale Y y su 95% de intervalo de confianza?
Para ello, ejecutamos el siguiente comando para crear un data.frame con el valor que nos interesa. En nuestro caso lo hemos llamado valoracalcular, pero podéis ponerle el que queráis. Es importante que el nombre de la variable esté escrito (en nuestro caso, el valor de la variable X se llama X):
valoracalcular=data.frame(X=85)
Y por fin, ejecutar el siguiente comando para ver los valores de Y y sus intervalos de confianza cuando X=85:
predict(EJEMPLO,valoracalcular,interval="confidence")
Obteniendo el siguiente resultado:
> predict(EJEMPLO,valoracalcular,interval="confidence") fit lwr upr 1 29.21167 28.44465 29.97869
En este caso, el valor de Y estimado (fit) es 29.21, mientras que el límite inferior del 95% del intervalo de confianza (lwr) es 28.44 y el superior (upr) es de 29.98.
¡Viva el Software Libre!
Referencias
- Para gráficos de gran calidad basados en regresiones y usando la librería ggplot2:
https://seriousstats.wordpress.com/2012/09/05/highqualitygraphs/ - Para el cálculo de predicción de valores de Y en la ecuación de regresión con sus 95% intervalos de confianza:
http://www.r-tutor.com/elementary-statistics/simple-linear-regression/prediction-interval-linear-regression
¡Viva el Software Libre!
I’m Yameng Zhang. This kind of post is great. It will be great If you can write it in English. I also use R for analysis and plot now…
Dear Yameng Zhang. I am planning to translate into English by request. As you can see in recent posts, I’m trying to compress all the information in one single picture, instead of printing screen several times and placing them along the post. The more simplier graphically the more understandable. So, I will rewrite that post following this phylosophy, and doing it in English. However, I will do it in the near future, so if you have any particular question, I would invite you to ask!
Great! I’m learning to use package «geomorph» now. I would like to ask some questions about the «pls» analysis. It’ll be great if you can post some essays on it.