Оценка точек пересечения и наклонов для трех обработок

Предположим, что мои данные выглядят так (это просто данные, иллюстрирующие мой вопрос):

set.seed(1234)
x=runif(12, 22,28);x
y=runif(12,88,120);y

y1=y[1:4];y1
y2=y[5:8];y2
y3=y[9:12];y3

x1=x[1:4];x1
x2=x[5:8];x2
x3=x[9:12];x3

dat= cbind(y1,x1,y2,x2,y3,x3);dat

            y1       x1        y2       x2        y3       x3
    [1,] 118.72718 27.79948 100.56455 23.59740  95.37221 25.51741
    [2,]  90.46189 26.12509  94.38618 27.96536 103.42347 24.44862
    [3,]  92.48808 25.67335 108.15713 24.08918 108.85686 26.52390
    [4,] 103.06760 25.84996 108.88237 26.37924 103.26130 22.05003

Предположим, что (y1,x1) представляет лечение 1, (y2,x2) представляет лечение 2, а (y3,x3) представляет лечение 3. Моя цель состоит в том, чтобы оценить точки пересечения и наклоны (используя пошаговый план). Затем определите матрицу гипотезы (нулевая гипотеза о том, что все наклоны одинаковы).

Я думал, что могу использовать индикаторы. Например, для x1= 1 при лечении и 0 в противном случае. Тогда x2= 1, если лечение 2, и 0 в противном случае. Вот иначе относится к третьему обращению.

Вот что я придумал:

    X< as.matrix(cbind(rep(1,12),c(rep(1,4),rep(0,8)),c(rep(0,4),rep(1,4),rep(0,4))))

Во всяком случае, я не знаю, куда идти отсюда. Мне было интересно, если кто-нибудь даст мне что-нибудь, чтобы начать. Пожалуйста, все обозначения должны быть в матричных формах. Благодарю вас!

Редактировать: я пытаюсь сделать это в матричной нотации, вот что я получил до сих пор, хотя у меня есть некоторые коэффициенты как NA:

set.seed(1234)
x=runif(12, 22,28);x
y=runif(12,88,120);y

y1=y[1:4];y1
y2=y[5:8];y2
y3=y[9:12];y3

x1=x[1:4];x1
x2=x[5:8];x2
x3=x[9:12];x3

dat= cbind(y1,x1,y2,x2,y3,x3);dat



 X1<- cbind( dat[,2],rep(0,4),rep(0,4), dat[,2],rep(1,4),rep(1,4));X1
 X2 <- cbind(dat[,4],rep(1,4),dat[,4], rep(0,4),rep(1,4),rep(1,4));X2
 X3 <- cbind(dat[,6],rep(1,4),dat[,6], rep(0,4),rep(0,4),rep(0,4));X3

 X <- rbind(X1,X2,X3)

 model <- lm(formula = y ~ X);summary(model)

person Smith Pay    schedule 19.10.2018    source источник


Ответы (2)


Если вы измените вторую часть, вы можете добиться того, что хотите сделать, связав свои значения y и x с третьей «категориальной» переменной, которая сообщает испытание, к которому оно принадлежит.

Затем сложите их в один фрейм данных из трех столбцов, пометьте столбцы, чтобы ваши переменные были прямыми, и свяжите их вместе.

Убедитесь, что ваши номера не были классифицированы в привязке (как мои).

Затем вы можете сделать линейную регрессию

one<- cbind(y1,x1,'one')
two<- cbind(y2, x2, 'two')
three<- cbind(y3, x3,'three')

dat<- data.frame(rbind(one, two, three))
colnames(dat)<- c('y', 'x', 'treatment')
dat$y<-as.numeric(dat$y) #binding tweaked numerics to categories
dat$x<-as.numeric(dat$x) #I converted them back
model <- lm(formula = y ~ 0 + x + treatment, data = dat)

Тогда ваш вывод выглядит следующим образом:

Call:
lm(formula = y ~ 0 + x + treatment, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4579 -2.2156 -0.4115  1.5442  6.6039 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)  
x               -0.1461     0.4283  -0.341   0.7418  
treatmentone     7.9185     3.9775   1.991   0.0817 .
treatmentthree   8.0112     2.5156   3.185   0.0129 *
treatmenttwo     6.4185     3.9775   1.614   0.1453  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.04 on 8 degrees of freedom
Multiple R-squared:  0.7991,    Adjusted R-squared:  0.6986 
F-statistic: 7.954 on 4 and 8 DF,  p-value: 0.006839

Расписывая формулу с нулем, вы получаете влияние x плюс эффект для каждого из коэффициентов.

Каждая бета представляет собой поправку пересечения для обработки, а коэффициент х представляет собой общий наклон.

РЕДАКТИРОВАТЬ в соответствии с предложением Роланда, это даст вам индивидуальные склоны.

Call:
lm(formula = y ~ x + treatment, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4579 -2.2156 -0.4115  1.5442  6.6039 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)  
(Intercept)      7.9185     3.9775   1.991   0.0817 .
x               -0.1461     0.4283  -0.341   0.7418  
treatmentthree   0.0927     3.4463   0.027   0.9792  
treatmenttwo    -1.5000     2.8570  -0.525   0.6138  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.04 on 8 degrees of freedom
Multiple R-squared:  0.08671,   Adjusted R-squared:  -0.2558 
F-statistic: 0.2532 on 3 and 8 DF,  p-value: 0.857
person sconfluentus    schedule 19.10.2018
comment
Не подавляйте перехват и включайте взаимодействие, чтобы разрешить разные наклоны: y ~ x * treatment - person Roland; 19.10.2018
comment
Это означает, что, включив точку пересечения, вы получите точку пересечения, как указано выше, а затем наклон для X и коэффициент для каждого УРОВНЯ категории. Поскольку они являются взаимоисключающими категориями, значение бета для двух других становится равным нулю при умножении на двоичную переменную. два других. В итоге вы получите коэффициент X и постоянную корректировку категории для точки пересечения y. - person sconfluentus; 20.10.2018
comment
Вы можете сложить переменные X, как я сделал выше, и все будет работать нормально. Вы не можете использовать традиционную матрицу с категориями, вам нужно преобразовать их в дамми, два столбца с 0, 1, чтобы сделать это. Лучший способ - использовать `lm() и позволить R управлять матричными математическими вычислениями. - person sconfluentus; 20.10.2018

Создайте каркас модели d (из данных в примечании в конце):

d <- with(DF, data.frame(y = c(y1, y2, y3), x = c(x1, x2, x3), g = gl(3, 4)))

а затем запустите регрессии для отдельных точек пересечения и наклонов по сравнению с отдельными точками пересечения с одинаковым наклоном. Наконец, используйте anova для выполнения теста.

fm1 <- lm(y ~ g + x:g + 0, d) # separate intercepts & slopes
coef(fm1)
##           g1           g2           g3         g1:x         g2:x         g3:x 
## -201.0983992  141.2889141   94.7617449   11.4666919   -1.5011629    0.3233902 

fm2 <- lm(y ~ g + x + 0, d) # separate intercepts, same slope
coef(fm2)
##         g1         g2         g3          x 
## 83.5468017 85.9297193 86.2446353  0.6691224 

anova(fm1, fm2)

давая:

Analysis of Variance Table

Model 1: y ~ g + x/g + 0
Model 2: y ~ g + x + 0
  Res.Df    RSS Df Sum of Sq      F Pr(>F)  
1      6 330.52                             
2      8 723.85 -2   -393.33 3.5701 0.0952 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

поэтому p = 0,0952, и две модели существенно не отличаются, т. е. наклоны существенно не отличаются друг от друга на уровне 5%.

Примечание

Предполагаемый ввод:

Lines <- "
       y1       x1        y2       x2        y3       x3
118.72718 27.79948 100.56455 23.59740  95.37221 25.51741
 90.46189 26.12509  94.38618 27.96536 103.42347 24.44862
 92.48808 25.67335 108.15713 24.08918 108.85686 26.52390
103.06760 25.84996 108.88237 26.37924 103.26130 22.05003"
DF <- read.table(text = Lines, header = TRUE)
person G. Grothendieck    schedule 19.10.2018