Regresión lineal en R (datos normales y logarítmicos)

Question

Deseo realizar una regresión lineal en R para datos en una gráfica logarítmica normal y en una doble.

Para datos normales el conjunto de datos podría ser la siguientes aparatos:

lin <- data.frame(x = c(0:6), y = c(0.3, 0.1, 0.9, 3.1, 5, 4.9, 6.2)) 
plot (lin$x, lin$y) 

No quiero calcular dibujar una línea para la regresión lineal de los puntos de datos solamente 2, 3 y 4.

Para datos logarítmica doble el conjunto de datos podrían ser los siguientes:

data = data.frame(
    x=c(1:15), 
    y=c(
     1.000, 0.742, 0.623, 0.550, 0.500, 0.462, 0.433, 
     0.051, 0.043, 0.037, 0.032, 0.028, 0.025, 0.022, 0.020 
    ) 
    ) 
plot (data$x, data$y, log="xy") 

Aquí quiero dibujar la línea de regresión para los conjuntos de datos 1: 7 y para 8:15.

Ho puedo calcular la pendiente y los Y-offset ALS así como parámetros para el ajuste (R^2, p-valor)?

¿Cómo se hace para datos normales y logarítmicos?

Gracias por su ayuda,

Sven

Answer 1

En R, lineales modelos de mínimos cuadrados están equipados con la función de lm(). Usando la interfaz de fórmula podemos utilizar el argumento subset para seleccionar los puntos de datos utilizados para ajustar el modelo real, por ejemplo:

lin <- data.frame(x = c(0:6), y = c(0.3, 0.1, 0.9, 3.1, 5, 4.9, 6.2)) 
linm <- lm(y ~ x, data = lin, subset = 2:4) 

entrega:

R> linm 

Call: 
lm(formula = y ~ x, data = lin, subset = 2:4) 

Coefficients: 
(Intercept)   x 
    -1.633  1.500 

R> fitted(linm) 
     2   3   4 
-0.1333333 1.3666667 2.8666667 

En cuanto a la doble registro, tiene dos elecciones, supongo; i) estimar dos modelos separados como hicimos anteriormente, o ii) estimar a través de ANCOVA. La transformación del registro se realiza en la fórmula usando log().

a través de dos modelos distintos:

logm1 <- lm(log(y) ~ log(x), data = dat, subset = 1:7) 
logm2 <- lm(log(y) ~ log(x), data = dat, subset = 8:15) 

O a través de ANCOVA, en el que necesitamos un indicador variable de

dat <- transform(dat, ind = factor(1:15 <= 7)) 
logm3 <- lm(log(y) ~ log(x) * ind, data = dat) 

Usted podría preguntarse si estos dos enfoques son equivalentes? Bueno, lo son y podemos mostrar esto a través de los coeficientes del modelo.

R> coef(logm1) 
    (Intercept)  log(x) 
-0.0001487042 -0.4305802355 
R> coef(logm2) 
(Intercept)  log(x) 
    0.1428293 -1.4966954 

Entonces las dos pendientes son -0.4306 y -1.4967 para los modelos separados. Los coeficientes para el modelo ANCOVA son:

R> coef(logm3) 
    (Intercept)   log(x)  indTRUE log(x):indTRUE 
    0.1428293  -1.4966954  -0.1429780  1.0661152 

¿Cómo conciliar los dos? Bueno, la forma en que configuré ind, logm3 está parametrizada para proporcionar valores más directamente estimados a partir del logm2; las intersecciones de logm2 y logm3 son las mismas, como lo son los coeficientes para log(x). Para obtener los valores equivalentes a los coeficientes de logm1, tenemos que hacer una manipulación, primero para la intercepción:

R> coefs[1] + coefs[3] 
    (Intercept) 
-0.0001487042 

donde el coeficiente de indTRUE es la diferencia en la media del grupo 1 sobre la media del grupo 2. y por la pendiente:

R> coefs[2] + coefs[4] 
    log(x) 
-0.4305802 

que es el mismo que nos dieron por logm1 y se basa en la pendiente para el grupo 2 (coefs[2]) modificada por la diferencia en la pendiente para el grupo 1 (coefs[4]).

En cuanto al trazado, una forma fácil es a través de abline() para modelos simples. P.ej. para el ejemplo de datos normal:

plot(y ~ x, data = lin) 
abline(linm) 

Para los datos de registro que podría tener que ser un poco más creativo, y la solución general es de predecir en el rango de los datos y la trama de las predicciones:

pdat <- with(dat, data.frame(x = seq(from = head(x, 1), to = tail(x,1), 
            by = 0.1)) 
pdat <- transform(pdat, yhat = c(predict(logm1, pdat[1:70,, drop = FALSE]), 
           predict(logm2, pdat[71:141,, drop = FALSE]))) 

Cuál puede trazar en la escala original, por exponentiating yhat

plot(y ~ x, data = dat) 
lines(exp(yhat) ~ x, dat = pdat, subset = 1:70, col = "red") 
lines(exp(yhat) ~ x, dat = pdat, subset = 71:141, col = "blue") 

o en la escala logarítmica:

plot(log(y) ~ log(x), data = dat) 
lines(yhat ~ log(x), dat = pdat, subset = 1:70, col = "red") 
lines(yhat ~ log(x), dat = pdat, subset = 71:141, col = "blue") 

Por ejemplo ...

Esta solución general funciona bien para el modelo ANCOVA más complejo también. Aquí se crea una nueva PDAT como antes y añadir en un indicador

pdat <- with(dat, data.frame(x = seq(from = head(x, 1), to = tail(x,1), 
            by = 0.1)[1:140], 
          ind = factor(rep(c(TRUE, FALSE), each = 70)))) 
pdat <- transform(pdat, yhat = predict(logm3, pdat)) 

Aviso cómo obtenemos todas las predicciones que queremos de la única llamada a predict() debido al uso de ANCOVA para adaptarse logm3. Ahora podemos trazar como antes:

plot(y ~ x, data = dat) 
lines(exp(yhat) ~ x, dat = pdat, subset = 1:70, col = "red") 
lines(exp(yhat) ~ x, dat = pdat, subset = 71:141, col = "blue") 

Answer 2

#Split the data into two groups 
data1 <- data[1:7, ] 
data2 <- data[8:15, ] 

#Perform the regression 
model1 <- lm(log(y) ~ log(x), data1) 
model2 <- lm(log(y) ~ log(x), data2) 
summary(model1) 
summary(model2) 

#Plot it 
with(data, plot(x, y, log="xy")) 
lines(1:7, exp(predict(model1, data.frame(x = 1:7)))) 
lines(8:15, exp(predict(model2, data.frame(x = 8:15)))) 

En general, la división de los datos en diferentes grupos y que ejecutan diferentes modelos en diferentes subconjuntos es inusual, y, probablemente, la forma mala. Es posible que desee considerar la adición de una variable de agrupación

data$group <- factor(rep(letters[1:2], times = 7:8)) 

y corriendo algún tipo de modelo en el conjunto de datos, por ejemplo,

model_all <- lm(log(y) ~ log(x) * group, data) 
summary(model_all)