Петля R: выполняет функцию для каждых 3 строк - PullRequest
Новая
       21

Петля R: выполняет функцию для каждых 3 строк

3 голосов
/ 06 июля 2019

У меня 2000 растений пшеницы, растущих в течение 40 дней.Я хотел бы выполнить функцию коэффициентов на каждом заводе, чтобы найти коэффициенты квадратного уравнения, которые составляют 3 момента времени.(a, b и c)

(1) Функция coef(lm(y~poly(x,2,raw=TRUE)) работает именно так, как я хочу.

(2) Однако способ представления моих данных требует от меня ручного задания x и y.

(3) Таким образом, я растопил свои данные и заказал их.

(4) Я хотел бы сделать цикл, который возьмет первые три в столбце «День» и установитэто как х.Затем я хотел бы взять первые три в столбце «Высота» и установить его как y.

Затем я хотел бы выполнить функцию коэффициентов.

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

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

1) Эта функция работает, давая мне коэффициенты: a, b, c 

x<-c(1,2,3)
y<-c(1,10,4)
coef(lm(y~poly(x,2,raw=TRUE)))

2) Это то, как мои данные изначально выглядели 

A = matrix(c(5, 4, 2, 10, 10, 4, 5, 15, 6),nrow=3, ncol=3)
colnames(A)<-c("10", "25", "40")
rownames(A)<-c("Wheat 1", "Wheat 2", "Wheat 3")
A

3) Это мой растопленный формат 

A.melted<-as.data.frame(melt(A, id.vars="ID"))
A.melted<-A.melted[with(A.melted,order(Var1)),]
colnames(A.melted) <- c("WheatID", "Day", "Height")
A.melted$Day<-as.numeric(as.character(A.melted$Day))
A.melted
#

4) Это то, что я пытаюсь сделать с моим циклом ....

  • для каждых 3 строк,
  • x<-A.melted[,2]
  • y<-A.melted[,3]
  • coef(lm(y~poly(x,2,raw=TRUE)))
  • что-то для компиляции коэффициентов: a, b, c

Я просто не знаком с синтаксисом циклов, и я хотел бы любые советы и предложения.Изучая Google говорит мне, что не следует делать циклы, если это не является абсолютно необходимым, так как я могу столкнуться с большим количеством проблем - таким образом, я также открыт для техник без циклов.

Ответы [ 4 ]

0 голосов
/ 06 июля 2019

Или вы можете выполнить apply() непосредственно на исходной матрице: 

x <- as.numeric(colnames(A))
apply(A, 1, function(y) coef(lm(y~poly(x,2,raw=TRUE))))

                            Wheat 1      Wheat 2       Wheat 3
(Intercept)             -3.88888889 -0.555555556  6.666667e-01
poly(x, 2, raw = TRUE)1  1.11111111  0.477777778  1.333333e-01
poly(x, 2, raw = TRUE)2 -0.02222222 -0.002222222 -2.417315e-18

Или вы можете транспонировать данные и напрямую использовать вызов coef(...): 

x <- as.numeric(colnames(A))
coef(lm(t(A) ~ poly(x, 2, raw = TRUE)))
0 голосов
/ 06 июля 2019

Если вы хотите сделать это в цикле, попробуйте это. Важнейшей частью является использование seq вместе с аргументом by =, чтобы индекс мог выполнять необходимые шаги. 

library(tibble)

df <- tibble(
  WheatID = rep(NA_character_, nrow(A)),
  Intercept = rep(NA_real_, nrow(A)),
  poly1 = rep(NA_real_, nrow(A)),
  poly2 = rep(NA_real_, nrow(A))
)
cnt <- 1
for (i in seq(1, nrow(A.melted), by = 3)) {
  x <- A.melted$Day[i + 0:2]
  y <- A.melted$Height[i + 0:2]
  df$WheatID[cnt] <- as.character(A.melted$WheatID[i])
  df[cnt, 2:4] <- coef(lm(y~poly(x,2,raw=TRUE)))
  cnt <- cnt + 1
}
df

Примечание: я не data.table парень. Поэтому я представляю вам tibble.

0 голосов
/ 06 июля 2019

Это то, что вы можете сделать с пакетом data.table. 

data.list <- split(A.melted, f = (1:nrow(A.melted) - 1) %/% 3)
coefs <- lapply(data.list, function(x) {
  coefs <- coef(lm(Day ~ poly(Height, raw=TRUE), data = x))
  data.table(
    intercept = coefs[1],
    poly.height = coefs[2]
  )
})

coefs <- rbindlist(coefs)
0 голосов
/ 06 июля 2019

Мы можем сделать это с помощью data.table, см. ?data.table: 

library(data.table)

A.models = A.melted[, model := list(.(lm(Height ~ poly(Day, 2), 
                                         data = list(.(.SD[WheatID == .BY[[1]]]))))), 
                    by = WheatID] 

A.models[, coefs := list(.(coefficients(model[[1]]))), 
         by = WheatID]

Вы можете получить доступ к каждой модели следующим образом: 

A.models[WheatID == "Wheat 1", model[[1]]]

и даже 

A.models[WheatID == "Wheat 1", summary(model[[1]])]

Волшебство здесь происходит, потому что data.table принимает J выражения , а не только функции.

...