¿Qué pasa si mis datos de regresión lineal contienen varias relaciones lineales combinadas?

Question

Digamos que yo soy el estudio de cómo los narcisos responder a diferentes condiciones de suelo. Me han recopilado datos sobre el pH del suelo frente a la altura madura de el narciso. Estoy esperando a una relación lineal, por lo que voy sobre la ejecución de una regresión lineal.

Sin embargo, no me di cuenta cuando empecé mi estudio que la población en realidad contiene dos variedades de narcisos, cada uno de los cuales responde de manera muy diferente a que el pH del suelo. De manera que la gráfica contiene dos diferentes relaciones lineales:

Puedo globo ocular y separar manualmente, por supuesto. Pero me pregunto si hay un enfoque más riguroso.

Preguntas:

1. Hay una prueba estadística para determinar si un conjunto de datos sería mejor ajuste por una sola línea o N líneas?

2. ¿Cómo puedo ejecutar una regresión lineal para adaptarse a las N líneas? En otras palabras, ¿cómo puedo separar la mezcla de datos?

No puedo pensar de algunos enfoques combinatorios, pero parecen costosas computacionalmente.

---

Aclaraciones:

1. La existencia de dos variedades era desconocido en el momento de la recolección de datos. La variedad de cada narciso no se observó, no se ha notado, y no se registran.

2. Es imposible recuperar esta información. Los narcisos han muerto desde el momento de la recolección de datos.

Tengo la impresión de que este problema es algo similar a la aplicación de algoritmos de clustering, en que casi se necesita saber el número de racimos antes de empezar. Yo creo que con CUALQUIER conjunto de datos, incrementar el número de líneas de disminuir el total de r.m.s. error. En el extremo, usted puede dividir el conjunto de datos en arbitraria pares y simplemente dibujar una línea a través de cada par. (E. g., si tuviera 1000 puntos de datos, se podría dividir en 500 arbitraria pares y dibujar una línea a través de cada par.) El ajuste sería exacto y la r.m.s. el error sería exactamente cero. Pero eso no es lo que queremos. Queremos que el "derecho" número de líneas.

Answer 1

Creo que Demetri la respuesta es un gran si asumimos que usted tiene las etiquetas de las distintas variedades. Cuando leí tu pregunta que no parece ser el caso para mí. Podemos utilizar un enfoque basado en el algoritmo EM, básicamente, de ajustarse al modelo que Demetri sugiere, pero sin saber las etiquetas para la variedad. Por suerte, el mixtools paquete en R proporciona esta funcionalidad para nosotros. Porque la información es bastante separados y parece que tienes un poco debe ser bastante exitosa.

library(mixtools)

# Generate some fake data that looks kind of like yours
n1 <- 150
ph1 = runif(n1, 5.1, 7.8)
y1 <- 41.55 + 5.185*ph1 + rnorm(n1, 0, .25)

n2 <- 150
ph2 <- runif(n2, 5.3, 8)
y2 <- 65.14 + 1.48148*ph2 + rnorm(n2, 0, 0.25)

# There are definitely better ways to do all of this but oh well
dat <- data.frame(ph = c(ph1, ph2), 
                  y = c(y1, y2), 
                  group = rep(c(1,2), times = c(n1, n2)))

# Looks about right
plot(dat

Podemos examinar los resultados

$ph, dat$

Por lo que caben dos regresiones y se estima que el 49,7% de las observaciones cayó en la regresión para el componente 1 y el 50,2% cayó en la regresión para el componente 2. La forma en que me simulada de los datos se estableció una división de 50-50 por lo que esta es buena.

La 'verdadera' de los valores que he usado para la simulación debe dar a las líneas:

y = 41.55 + 5.185*ph y y = 65.14 + 1.48148*ph

(que yo estima "a mano" de su parcela, de manera que los datos que crear una apariencia similar a la suya) y las líneas que el algoritmo EM dio en este caso fueron:

y = 41.514 + 5.19*ph y y = 64.655 + 1.55*ph

Muy cerca de los valores reales.

Podemos trazar el conjunto de líneas, junto con los datos

y)

# Fit the regression. One line for each component. This defaults
# to assuming there are two underlying groups/components in the data
out <- regmixEM(y = dat

Answer 2

EDIT: yo pensaba originalmente OP sabía que las observaciones de vino a partir de la cual las especies. OP editar deja en claro que mi enfoque original no es factible. Lo voy a dejar para la posteridad, pero la otra respuesta es mucho mejor. Como consuelo, he codificado de un modelo de mezcla de Stan. No estoy diciendo que un enfoque Bayesiano es particularmente bueno en este caso, pero es solo algo limpio que puedo aportar.

Stan Código

data{

  //Number of data points
  int N; 

  real y[N];
  real x[N];
}
parameters{
  //mixing parameter
  real<lower=0, upper =1>  theta;

  //Regression intercepts
  real beta_0[2];

  //Regression slopes.
  ordered[2] beta_1;

  //Regression noise
  real<lower=0> sigma[2];
}
model{

  //priors
  theta ~ beta(5,5);
  beta_0 ~ normal(0,1);
  beta_1 ~ normal(0,1);
  sigma ~ cauchy(0,2.5);

  //mixture likelihood
  for (n in 1:N){
    target+=log_mix(theta,
                     normal_lpdf(y[n] | beta_0[1] + beta_1[1]*x[n], sigma[1]),
                     normal_lpdf(y[n] | beta_0[2] + beta_1[2]*x[n], sigma[2]));
  }
}
generated quantities {
  //posterior predictive distribution
  //will allow us to see what points belong are assigned
  //to which mixture 
  matrix[N,2] p;
  matrix[N,2] ps;
  for (n in 1:N){
    p[n,1] = log_mix(theta,
                     normal_lpdf(y[n] | beta_0[1] + beta_1[1]*x[n], sigma[1]),
                     normal_lpdf(y[n] | beta_0[2] + beta_1[2]*x[n], sigma[2]));

    p[n,2]= log_mix(1-theta,
                     normal_lpdf(y[n] | beta_0[1] + beta_1[1]*x[n], sigma[1]),
                     normal_lpdf(y[n] | beta_0[2] + beta_1[2]*x[n], sigma[2]));

    ps[n,]= p[n,]/sum(p[n,]);
  }
}

Ejecutar El Modelo Stan De R

library(tidyverse)
library(rstan)


#Simulate the data
N = 100
x = rnorm(N, 0, 3)
group = factor(sample(c('a','b'),size = N, replace = T))

y = model.matrix(~x*group)%*% c(0,1,0,2) 
y = as.numeric(y) + rnorm(N)

d = data_frame(x = x, y = y)

d %>% 
  ggplot(aes(x,y))+
  geom_point()

#Fit the model
N = length(x)
x = as.numeric(x)
y = y

fit = stan('mixmodel.stan', 
           data = list(N= N, x = x, y = y),
           chains = 8,
           iter = 4000)

Resultados

Las líneas discontinuas son terreno de la verdad, las líneas continuas son estimadas.

---

Respuesta Original

Si usted sabe que la muestra proviene de la variedad de narciso, se puede estimar una interacción entre variedad y el PH del suelo.

Su modelo se verá como

$$ y = \beta_0 + \beta_1 \text{variety} + \beta_2\text{PH} + \beta_3\text{variety}\cdot\text{PH} $$

Here is an example in R. I've generated some data that looks like this:

Clearly two different lines, and the lines correspond to two species. Here is how to estimate the lines using linear regression.

library(tidyverse)

#Simulate the data
N = 1000
ph = runif(N,5,8)
species = rbinom(N,1,0.5)

y = model.matrix(~ph*species)%*% c(20,1,20,-3) + rnorm(N, 0, 0.5)
y = as.numeric(y)

d = data_frame(ph = ph, species = species, y = y)

#Estimate the model
model = lm(y~species*ph, data = d)
summary(model)

And the result is

> summary(model)

Call:
lm(formula = y ~ species * ph, data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.61884 -0.31976 -0.00226  0.33521  1.46428 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 19.85850    0.17484  113.58   <2e-16 ***
species     20.31363    0.24626   82.49   <2e-16 ***
ph           1.01599    0.02671   38.04   <2e-16 ***
species:ph  -3.03174    0.03756  -80.72   <2e-16 ***
---
Signif. codes:  
0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1

Residual standard error: 0.4997 on 996 degrees of freedom
Multiple R-squared:  0.8844,    Adjusted R-squared:  0.8841 
F-statistic:  2541 on 3 and 996 DF,  p-value: < 2.2e-16

For species labeled 0, the line is approximately

$$ y = 19 + 1\cdot \text{PH}$$

For species labeled 1, the line is approximately

$$ y = 40 - 2 \cdot \text{PH} $$

Answer 3

El enfoque estadístico es muy similar a las dos de la respuesta anterior, pero se trata de un poco más con cómo escoger el número de clases latentes si la falta de conocimientos previos. Puede utilizar los criterios de información o de parsimonia como guía para elegir el número de clases latentes.

Aquí está una Stata ejemplo mediante una secuencia finita modelos de mezcla (FMMs) con 2-4 latente de clases y/o componentes. La primera tabla de los coeficientes de las latente de pertenencia a una clase. Estas son un poco difíciles de interpretar, pero pueden ser convertidos a probabilidades más tarde con estat lcprob. Para cada clase, usted también consigue una intercepción y un ph pendiente parámetro, seguido por latentes de la clase marginal probabilidades, y dos en la muestra de ICs. Estos coeficientes estimados se interpretan como los coeficientes de un modelo de regresión lineal. Aquí los más pequeños en la muestra de BIC dice que usted escoja el modelo de dos componentes como la mejor. AIC extrañamente selecciona el 3 de modelo de componente. También se puede utilizar fuera de la muestra de ICs para recoger o utilizar la validación cruzada.

Otra manera de medir que están empujando a los datos demasiado lejos es si la última participación de la clase es muy pequeño, ya que un adicional de los componentes, puede simplemente reflejar la presencia de valores atípicos en los datos. En ese caso, la parsimonia de los favores de simplificar el modelo y la extracción de los componentes. Sin embargo, si usted piensa que las clases son pequeñas posible en su configuración, esto puede no ser el canario en la mina de carbón. Aquí la parsimonia de los favores de la 2 modelo de componente desde la tercera clase sólo contiene $.0143313 \cdot 300 \approx 4$ observaciones.

La FMM enfoque no siempre funcionan bien en la práctica si las clases son menos duras. Usted puede ejecutar en las dificultades de cálculo, con muchas clases latentes, especialmente si usted no tiene la cantidad suficiente de datos, o de la probabilidad de que la función de múltiples locales de maxima.

. clear

. /* Fake Data */
. set seed 10011979

. set obs 300
number of observations (_N) was 0, now 300

. gen     ph = runiform(5.1, 7.8) in 1/150
(150 missing values generated)

. replace ph = runiform(5.3, 8)   in 151/300
(150 real changes made)

. gen y      = 41.55 + 5.185*ph   + rnormal(0, .25)  in 1/150
(150 missing values generated)

. replace y  = 65.14 + 1.48148*ph + rnormal(0, 0.25) in 151/300
(150 real changes made)

. 
. /* 2 Component FMM */
. fmm 2, nolog: regress y ph

Finite mixture model                            Number of obs     =        300
Log likelihood =  -194.5215

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.Class      |  (base outcome)
-------------+----------------------------------------------------------------
2.Class      |
       _cons |   .0034359   .1220066     0.03   0.978    -.2356927    .2425645
------------------------------------------------------------------------------

Class          : 1
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   5.173137   .0251922   205.35   0.000     5.123761    5.222513
       _cons |     41.654   .1622011   256.80   0.000      41.3361    41.97191
-------------+----------------------------------------------------------------
     var(e.y)|   .0619599   .0076322                      .0486698     .078879
------------------------------------------------------------------------------

Class          : 2
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   1.486062    .026488    56.10   0.000     1.434147    1.537978
       _cons |   65.10664   .1789922   363.74   0.000     64.75582    65.45746
-------------+----------------------------------------------------------------
     var(e.y)|   .0630583   .0075271                      .0499042    .0796797
------------------------------------------------------------------------------

. estat lcprob

Latent class marginal probabilities             Number of obs     =        300

--------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
       Class |
          1  |    .499141   .0305016      .4396545    .5586519
          2  |    .500859   .0305016      .4413481    .5603455
--------------------------------------------------------------

. estat ic

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        300         .  -194.5215       7     403.043   428.9695
-----------------------------------------------------------------------------
               Note: N=Obs used in calculating BIC; see [R] BIC note.

. 
. /* 3 Component FMM */
. fmm 3, nolog: regress y ph

Finite mixture model                            Number of obs     =        300
Log likelihood =  -187.4824

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.Class      |  (base outcome)
-------------+----------------------------------------------------------------
2.Class      |
       _cons |  -.0312504    .123099    -0.25   0.800    -.2725199    .2100192
-------------+----------------------------------------------------------------
3.Class      |
       _cons |  -3.553227   .5246159    -6.77   0.000    -4.581456   -2.524999
------------------------------------------------------------------------------

Class          : 1
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   5.173077   .0252246   205.08   0.000     5.123637    5.222516
       _cons |   41.65412     .16241   256.48   0.000      41.3358    41.97243
-------------+----------------------------------------------------------------
     var(e.y)|   .0621157   .0076595                      .0487797    .0790975
------------------------------------------------------------------------------

Class          : 2
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   1.476049   .0257958    57.22   0.000      1.42549    1.526608
       _cons |   65.18698   .1745018   373.56   0.000     64.84496    65.52899
-------------+----------------------------------------------------------------
     var(e.y)|   .0578413   .0070774                      .0455078    .0735173
------------------------------------------------------------------------------

Class          : 3
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   1.776746   .0020074   885.09   0.000     1.772811     1.78068
       _cons |   62.76633   .0134072  4681.54   0.000     62.74005    62.79261
-------------+----------------------------------------------------------------
     var(e.y)|   9.36e-06   6.85e-06                      2.23e-06    .0000392
------------------------------------------------------------------------------

. estat lcprob

Latent class marginal probabilities             Number of obs     =        300

--------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
       Class |
          1  |   .5005343   .0304855      .4410591    .5599944
          2  |   .4851343   .0306119      .4256343    .5450587
          3  |   .0143313   .0073775      .0051968     .038894
--------------------------------------------------------------

. estat ic

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        300         .  -187.4824      11    396.9648   437.7064
-----------------------------------------------------------------------------
               Note: N=Obs used in calculating BIC; see [R] BIC note.

. 
. /* 4 Component FMM */
. fmm 4, nolog: regress y ph

Finite mixture model                            Number of obs     =        300
Log likelihood = -188.06042

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.Class      |  (base outcome)
-------------+----------------------------------------------------------------
2.Class      |
       _cons |  -.6450345   .5853396    -1.10   0.270    -1.792279      .50221
-------------+----------------------------------------------------------------
3.Class      |
       _cons |  -.8026907   .6794755    -1.18   0.237    -2.134438    .5290568
-------------+----------------------------------------------------------------
4.Class      |
       _cons |  -3.484714   .5548643    -6.28   0.000    -4.572229     -2.3972
------------------------------------------------------------------------------

Class          : 1
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   5.173031   .0251474   205.71   0.000     5.123743    5.222319
       _cons |   41.65574    .161938   257.23   0.000     41.33835    41.97313
-------------+----------------------------------------------------------------
     var(e.y)|   .0617238   .0076596                      .0483975    .0787195
------------------------------------------------------------------------------

Class          : 2
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   1.503764   .0371216    40.51   0.000     1.431007    1.576521
       _cons |   65.13498   .2666049   244.31   0.000     64.61244    65.65751
-------------+----------------------------------------------------------------
     var(e.y)|   .0387473   .0188853                      .0149062    .1007195
------------------------------------------------------------------------------

Class          : 3
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   1.441334   .0443892    32.47   0.000     1.354333    1.528335
       _cons |   65.26791   .2765801   235.98   0.000     64.72582       65.81
-------------+----------------------------------------------------------------
     var(e.y)|   .0307352    .010982                      .0152578    .0619127
------------------------------------------------------------------------------

Class          : 4
Response       : y
Model          : regress

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
y            |
          ph |   1.665207   .0079194   210.27   0.000     1.649685    1.680728
       _cons |   63.42577   .0510052  1243.52   0.000      63.3258    63.52573
-------------+----------------------------------------------------------------
     var(e.y)|    .000096   .0000769                        .00002    .0004611
------------------------------------------------------------------------------

. estat lcprob

Latent class marginal probabilities             Number of obs     =        300

--------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
       Class |
          1  |   .4991443   .0304808      .4396979     .558615
          2  |   .2618733   .1506066      .0715338    .6203076
          3  |   .2236773    .150279      .0501835    .6110804
          4  |    .015305    .008329       .005234    .0438994
--------------------------------------------------------------

. estat ic

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        300         .  -188.0604      15    406.1208   461.6776
-----------------------------------------------------------------------------
               Note: N=Obs used in calculating BIC; see [R] BIC note.

Answer 4

Me centraré en la cuestión de la significación estadística desde Dason cubierto ya la parte de modelado.

Estoy familiarizado con alguno de pruebas formales para esto (que estoy seguro que existe), así que voy a lanzar algunas ideas que hay (y probablemente voy a agregar código R y detalles técnicos más tarde).

En primer lugar, es conveniente para inferir las clases. Supongamos que usted tiene dos líneas de ajuste a los datos, se puede reconstruir aproximadamente las dos clases asignando a cada punto a la clase de la línea más cercano a él. Para puntos cerca de la intersección, que se ejecutará en problemas, pero por ahora simplemente ignorar a aquellos (puede ser una manera de conseguir alrededor de esto, pero por ahora solo espero que esto no cambia mucho).

La manera de hacer esto es elegir a $x_{l}$ e $x_{r}$ (suelos con valores de pH) con $x_{l} \leq x_{r}$ tal que las partes izquierda de $x_{l}$ están suficientemente separadas y las partes derecha de $x_{r}$ están suficientemente separadas (el punto más cercano donde las distribuciones no se superponen).

Entonces hay dos maneras naturales veo a hacerlo.

La menos divertida forma es simplemente ejecutar el conjunto de datos original, combinado con el inferirse etiquetas de clase a través de una regresión lineal como en Demetri la respuesta.

Una más que interesante manera de hacerlo sería a través de una versión modificada de ANOVA. El punto es crear una artificial conjunto de datos que representa las dos líneas (con similar se extendió entre ellos) y, a continuación, aplicar el ANOVA. Técnicamente, usted necesita hacer esto una vez para el lado izquierdo, y una vez por la derecha (es decir, usted tendrá dos conjuntos de datos artificiales).

Empezamos con la izquierda, y aplicar un simple promedio de enfoque para obtener dos grupos. Básicamente, cada punto de decir que la primera clase es de la forma $$ y^{(i)}_{1} = \beta_{1,1} x_{1}^{(i)} + \beta_{1,0} + e_{1}^{(i)}$$ así que vamos a sustituir la expresión lineal $\beta_{1,1} x_{1}^{(i)} + \beta_{1,0}$ por una constante, es decir, el valor promedio del término lineal o $$ \beta_{1,1} x^{\mathrm{avg}} + \beta_{1, 0}$$ donde $x^{\mathrm{avg}}_{l}$ es, literalmente, el promedio de $x$ valor para el lado izquierdo (es importante destacar que este es de más de dos clases, ya que hace las cosas más consistente). Que es, se reemplazan $y_{1}^{(i)}$con $$ \tilde{y}_{1}^{(i)} = \beta_{1,1} x^{\mathrm{avg}} + \beta_{1, 0} + e_{1}^{(i)},$$ y hacemos lo mismo para la segunda clase. Es decir, el nuevo conjunto de datos consiste en la recopilación de $\tilde{y}_{1}^{(i)}$ y de manera similar a $\tilde{y}_{2}^{(i)}$.

Tenga en cuenta que ambos enfoques natural generalizar a $N$ de las clases.

Answer 5

Es posible que la inclusión de ambos en el mismo gráfico es un error? Dado que las variedades se comportan completamente diferente hay ningún valor en la superposición de los datos? A mí me parece que usted está buscando para los impactos a una especie de narciso, no de los impactos de los entornos similares en diferentes narcisos. Si usted ha perdido los datos que ayuda a determinar la especie "a" de la especie "B" puede simplemente el comportamiento del grupo "a" y el comportamiento "B" e incluyen el descubrimiento de dos especies en su narrativa. O, si usted realmente desea un gráfico, basta con utilizar dos conjuntos de datos sobre el mismo eje. No tengo a ninguna parte cerca de la experiencia que veo en los otros las respuestas dadas, así que tengo que encontrar menos "expertos" de los métodos. Me gustaría ejecutar un análisis de los datos en una hoja de cálculo entorno donde las ecuaciones son más fáciles de desarrollar. Luego, una vez que los grupos se vuelven obvias, crear las dos tablas de datos, seguido por la conversión de ellos en tablas/gráficos. Yo trabajo con una gran cantidad de datos y a menudo me encuentro que mis suposiciones de las diferentes correlaciones resultan mal; que es lo que los datos se supone que es para que nos ayude a descubrir. Una vez que me entero de que mis suposiciones son incorrectas, mostrar los datos, basado en los comportamientos descubierto y discutir los comportamientos y las resultantes de los análisis estadísticos como parte de la narrativa.