Форум врачей-аспирантов

Здравствуйте, гость ( Вход | Регистрация )

2 страниц V   1 2 >  
Добавить ответ в эту темуОткрыть тему
> Помогите, пожалуйста, написать алгоритм для R
Olga_
сообщение 29.08.2012 - 12:02
Сообщение #1





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



Добрый день.

Помогите, пожалуйста, написать алгоритм для R.
У меня две базы данных, одна из них включает номер пациента, дату каждого измерения (как правило, измерение проводилось в среднем каждуjу неделю в течении 1 года ) и номер недели (1,2,3...54), вторая база включает информацию о ежедневной температуре воздуха.

Задача: нужно посчитать среднюю температуру, предшествующую каждому измерению. Например, для первой недели нужно рассчитать среднюю темп за 7 дней до измерения (28apr1999-7дней), для второй недели за период 28apr1999-07may1999 и т.д.
FID studydate week
1 28apr1999 1
1 07may1999 2
1 14may1999 3
1 21may1999 4


Second dataset
temp date
5.37 1999-04-20
2.13 1999-04-21
1.6 1999-04-22
-0.17 1999-04-23
2.53 1999-04-24


Спасибо!
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 29.08.2012 - 16:37
Сообщение #2





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(Olga_ @ 29.08.2012 - 12:02) *
Добрый день.

Помогите, пожалуйста, написать алгоритм для R.
У меня две базы данных, одна из них включает номер пациента, дату каждого измерения (как правило, измерение проводилось в среднем каждуjу неделю в течении 1 года ) и номер недели (1,2,3...54), вторая база включает информацию о ежедневной температуре воздуха.

Задача: нужно посчитать среднюю температуру, предшествующую каждому измерению. Например, для первой недели нужно рассчитать среднюю темп за 7 дней до измерения (28apr1999-7дней), для второй недели за период 28apr1999-07may1999 и т.д.
FID studydate week
1 28apr1999 1
1 07may1999 2
1 14may1999 3
1 21may1999 4


Second dataset
temp date
5.37 1999-04-20
2.13 1999-04-21
1.6 1999-04-22
-0.17 1999-04-23
2.53 1999-04-24


Спасибо!


Пожалуйста! smile.gif

Я предлагаю посчитать скользящим окном среднее для всех дней в базе. При импорте формата даты первого файла на некоторых системах возможно понадобится переключить локаль (мне в utf8 понадобилось).

Примеры слишком короткие для демонстрации работоспособности.

Код
> temp.time<-read.table("time.txt", header=TRUE)
> temp.time$date<-as.Date(temp.time$date)
> str(temp.time)
'data.frame':   5 obs. of  2 variables:
$ temp: num  5.37 2.13 1.6 -0.17 2.53
$ date:Class 'Date'  num [1:5] 10701 10702 10703 10704 10705
>
> library(caTools)
> temp.time <- data.frame(temp.time, avgtemp=runmean(temp.time$temp, 3, alg= "exact", endrule="mean", align =  "right"))
> temp.time
   temp       date  avgtemp
1  5.37 1999-04-20 3.750000
2  2.13 1999-04-21 3.750000
3  1.60 1999-04-22 3.033333
4 -0.17 1999-04-23 1.186667
5  2.53 1999-04-24 1.320000

> lct <- Sys.getlocale("LC_TIME"); Sys.setlocale("LC_TIME", "C")
> data<-read.table("data.txt", header=TRUE)
> data$studydate<-as.Date(data$studydate, "%d%b%Y")
> data
  FID  studydate week
1   1 1999-04-28    1
2   1 1999-05-07    2
3   1 1999-05-14    3
4   1 1999-05-21    4
> str(data)
'data.frame':   4 obs. of  3 variables:
$ FID      : int  1 1 1 1
$ studydate:Class 'Date'  num [1:4] 10709 10718 10725 10732
$ week     : int  1 2 3 4

> Sys.setlocale("LC_TIME", lct)

> temp.time$avgtemp[temp.time$date==data[1,2]]


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
Olga_
сообщение 29.08.2012 - 18:43
Сообщение #3





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



Спасибо!

У меня вопрос, почему к(width of moving window)=3?
Не совсем поняла, как работает предложенный вами метод..


Я думала, что нужно написать loop для решения этой задачи. Но мои попытки не дали желаемого результата..
airpollution <- read.table("C:.../Pay.csv",blank.lines.skip = TRUE,fill=T,header=T,sep=";")
summary(airpollution)
dat <- paste(airpollution[,3],"-",airpollution[,2],"-",airpollution[,1],sep="")
dat <- as.Date(dat, format="%Y-%m-%d")
pay <- cbind(airpollution,dat)



x <- read.dta("C:/.../symptoms.dta")

x$stday <- as.Date(x$stday, origin="1960-01-01",format="%Y-%m-%d")
temp<- rep(NA, length(x$FID))
index <- 1:length(pay$dat)
pm10<- rep(NA, length(x$FID))


for (i in 1:length(x$FID)) {
start <- x$stday[i]==pay$dat
s <- index[start]
end <- x$stday[i]+1==pay$dat
e <- index[end]
pm10[i] <- mean(pay$pm10[s:e],na.rm = T)
paytemp[i] <- mean(pay$temp[s:e],na.rm = T)

}

Не уверена насчет этой части
end <- x$stday[i]+1==pay$dat

+1 это значит, что следующая строка, или я ошибаюсь?
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 29.08.2012 - 20:57
Сообщение #4





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(Olga_ @ 29.08.2012 - 18:43) *
Спасибо!

У меня вопрос, почему к(width of moving window)=3?
Не совсем поняла, как работает предложенный вами метод..


Я думала, что нужно написать loop для решения этой задачи. Но мои попытки не дали желаемого результата..
airpollution <- read.table("C:.../Pay.csv",blank.lines.skip = TRUE,fill=T,header=T,sep=";")
summary(airpollution)
dat <- paste(airpollution[,3],"-",airpollution[,2],"-",airpollution[,1],sep="")
dat <- as.Date(dat, format="%Y-%m-%d")
pay <- cbind(airpollution,dat)



x <- read.dta("C:/.../symptoms.dta")

x$stday <- as.Date(x$stday, origin="1960-01-01",format="%Y-%m-%d")
temp<- rep(NA, length(x$FID))
index <- 1:length(pay$dat)
pm10<- rep(NA, length(x$FID))


for (i in 1:length(x$FID)) {
start <- x$stday[i]==pay$dat
s <- index[start]
end <- x$stday[i]+1==pay$dat
e <- index[end]
pm10[i] <- mean(pay$pm10[s:e],na.rm = T)
paytemp[i] <- mean(pay$temp[s:e],na.rm = T)

}

Не уверена насчет этой части
end <- x$stday[i]+1==pay$dat

+1 это значит, что следующая строка, или я ошибаюсь?


1) на таких коротких данных я выбрал усреднение по 3 дням (или меньше за счет выбранного алгоритма) чтобы хоть что то работало для отладки.

Программа работает просто --- для всех дней справочника температуры рассчитывается среднее за предыдущие 7 дней (у меня для отладки на Ваших лапидарных примерах таблиц -- за 3 дня).

После этого Вы просто запрашиваете последней строкой из справочника среднюю температуру для даты инцидента.

Скажите куда Вам это среднее надо подставить? И приведите более длинные куски данных. Иначе мне непонятно чем Вас не устраивает предложенное мной решение.

2) Ну что могу сказать, как говорится --- " Хороший программист на фортране, на любом языке напишет программу на фортране" smile.gif. По поводу указанной строчки кода -- это какая то невозможная конструкция

Я не использую циклов в R, и Вам не советую. В R все делается через *apply() операторы. Для использование for() нужна крайне исключительная причина.

Прочитайте пожалуйста руководство "Введение в R" (или на английском, или на русском со страницы перевода http://m7876.wiki.zoho.com/Introduction-to-R.html , а то у меня какое то чувство что я его зря переводил)


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
Olga_
сообщение 29.08.2012 - 21:17
Сообщение #5





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



Цитата(p2004r @ 29.08.2012 - 19:57) *
1) на таких коротких данных я выбрал усреднение по 3 дням (или меньше за счет выбранного алгоритма) чтобы хоть что то работало для отладки.

Программа работает просто --- для всех дней справочника температуры рассчитывается среднее за предыдущие 7 дней (у меня для отладки на Ваших лапидарных примерах таблиц -- за 3 дня).

После этого Вы просто запрашиваете последней строкой из справочника среднюю температуру для даты инцидента.

Скажите куда Вам это среднее надо подставить? И приведите более длинные куски данных. Иначе мне непонятно чем Вас не устраивает предложенное мной решение.

2) Ну что могу сказать, как говорится --- " Хороший программист на фортране, на любом языке напишет программу на фортране" smile.gif. По поводу указанной строчки кода -- это какая то невозможная конструкция

Я не использую циклов в R, и Вам не советую. В R все делается через *apply() операторы. Для использование for() нужна крайне исключительная причина.

Прочитайте пожалуйста руководство "Введение в R" (или на английском, или на русском со страницы перевода http://m7876.wiki.zoho.com/Introduction-to-R.html , а то у меня какое то чувство что я его зря переводил)


Проблема в том, что это в среднем нужно рассчитать температуру за предыдущие семь дней, но это не всегда так. Есть периоды более 7 дней и менее 7 дней. Только для первого измерения всегда нужна средняя температура за предыдущие 7 дней. Меня то как раз рассчитать среднее за 7 дней очень устраивает:)

Я не программист и в R только начала работать. Мои коллеги рекомендовали использовать for для того, чтобы рассчитать среднее...
Почему вы пишите, что для for нужна крайне исключительная причина?

Данные сейчас загружу
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
Olga_
сообщение 29.08.2012 - 21:17
Сообщение #6





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



---

Сообщение отредактировал Olga_ - 8.11.2012 - 22:35
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 30.08.2012 - 10:44
Сообщение #7





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(Olga_ @ 29.08.2012 - 21:17) *
Проблема в том, что это в среднем нужно рассчитать температуру за предыдущие семь дней, но это не всегда так. Есть периоды более 7 дней и менее 7 дней. Только для первого измерения всегда нужна средняя температура за предыдущие 7 дней. Меня то как раз рассчитать среднее за 7 дней очень устраивает:)

Я не программист и в R только начала работать. Мои коллеги рекомендовали использовать for для того, чтобы рассчитать среднее...
Почему вы пишите, что для for нужна крайне исключительная причина?

Данные сейчас загружу


Данные это хорошо. Но:

1) Что значит "менее-более семи дней"?! Вы в первом сообщении написали а) что у Вас база данных температуры именно за _каждый_ день, б) что для каждого инциндента из таблицы случаев надо получить среднюю температуру за предыдущие 7 дней.

2) Объясните структуру своего эксперимента. Типа: "Для каждой строки из таблицы такой то получить то то и то то...". Ну или просто раскажите что надо получитью Иначе мне не понятно как надо группировать данные, и я не смогу Вам помочь.

3) Привет коллегам, они хорошие "программисты на фортране" smile.gif В 95% случаев анализа данных в R for() не нужен (а следовательно , по заветам медицинской статистики не существует smile.gif, именно для этого R и разрабатывали собственно --- что бы избежать низкоуровневого программирования.


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
Olga_
сообщение 30.08.2012 - 12:02
Сообщение #8





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



Цитата(p2004r @ 30.08.2012 - 09:44) *
Данные это хорошо. Но:

1) Что значит "менее-более семи дней"?! Вы в первом сообщении написали а) что у Вас база данных температуры именно за _каждый_ день, б) что для каждого инциндента из таблицы случаев надо получить среднюю температуру за предыдущие 7 дней.

2) Объясните структуру своего эксперимента. Типа: "Для каждой строки из таблицы такой то получить то то и то то...". Ну или просто раскажите что надо получитью Иначе мне не понятно как надо группировать данные, и я не смогу Вам помочь.

3) Привет коллегам, они хорошие "программисты на фортране" smile.gif В 95% случаев анализа данных в R for() не нужен (а следовательно , по заветам медицинской статистики не существует smile.gif, именно для этого R и разрабатывали собственно --- что бы избежать низкоуровневого программирования.

Используя базу данных температуры за каждый день. (файл temp.txt ) нужно рассчитать среднюю температуру для каждого участника исследования (data.txt) за период с момента предыдущего наблюдения до следующего. Обычно разница между наблюдениями составляет 7 дней, но может быть и больше, например, для 1 пациента период между первым наблюдением и вторым
stday
28apr1999
07may1999
Конечно, такие случаи встречаются редко и их можно просто игнорировать, используя предложенную вами методику рассчета (средняя температура за предыдущие 7 дней).
Для первого измерения всегда нужна средняя температура за предыдущие 7 дней.

Одна из целей исследования - оценить влияние факторов окружающей среды (темп., РМ10, О3 и т.д.) на частоту респираторных симптомов в течении 1 года.

Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 30.08.2012 - 15:37
Сообщение #9





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(Olga_ @ 30.08.2012 - 12:02) *
Используя базу данных температуры за каждый день. (файл temp.txt ) нужно рассчитать среднюю температуру для каждого участника исследования (data.txt) за период с момента предыдущего наблюдения до следующего. Обычно разница между наблюдениями составляет 7 дней, но может быть и больше, например, для 1 пациента период между первым наблюдением и вторым
stday
28apr1999
07may1999
Конечно, такие случаи встречаются редко и их можно просто игнорировать, используя предложенную вами методику рассчета (средняя температура за предыдущие 7 дней).
Для первого измерения всегда нужна средняя температура за предыдущие 7 дней.

Одна из целей исследования - оценить влияние факторов окружающей среды (темп., РМ10, О3 и т.д.) на частоту респираторных симптомов в течении 1 года.



Ольга, Вам не кажется что эти данные немного отличаются от первоначальных? Вы не ошиблись с файлом?


Код
day     month   year    temp    temp_min        temp_max                                                                
1       1       1998    1.15                                                                            
2       1       1998    6.64                                                                            
3       1       1998    6.01                                                                            
4       1       1998    6.16                                                                            
5       1       1998    7.55                                                                            
6       1       1998    6.36                                                                            
7       1       1998    8.13                                                                            
8       1       1998    7.66                                                                            
9       1       1998    2.42



Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
Olga_
сообщение 30.08.2012 - 16:10
Сообщение #10





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



Цитата(p2004r @ 30.08.2012 - 14:37) *
Ольга, Вам не кажется что эти данные немного отличаются от первоначальных? Вы не ошиблись с файлом?


Код
day     month   year    temp    temp_min        temp_max                                                                
1       1       1998    1.15                                                                            
2       1       1998    6.64                                                                            
3       1       1998    6.01                                                                            
4       1       1998    6.16                                                                            
5       1       1998    7.55                                                                            
6       1       1998    6.36                                                                            
7       1       1998    8.13                                                                            
8       1       1998    7.66                                                                            
9       1       1998    2.42


Нет, сначала , как пример, я привела модифицированные данные, где я обьединила в одну коллонку время исследования.
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 30.08.2012 - 19:50
Сообщение #11





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(Olga_ @ 30.08.2012 - 16:10) *
Нет, сначала , как пример, я привела модифицированные данные, где я обьединила в одну коллонку время исследования.


ну так почему бы и не выкладывать именно рабочий файл, раз уж Вы время привели к вменяемому формату?

ну тогда уж продолжу прошлый "алгоритм":

1) получаем массив средних за 3-14 последних дней на текущую дату (колонку даты возьмете свою smile.gif

Код
tmp<-sapply(1:10, function(i) avgtemp=runmean(temp.time$temp, i+3, alg= "exact", endrule="mean", align =  "right") )


2) разницу между текущим обследованием и предыдущим получаем

Код
>  data<-read.table("data.txt", header=TRUE)
>  data$stday<-as.Date(data$stday, "%d%b%Y")
> head(data)
  FID      stday woche
1   1 1999-04-28     1
2   1 1999-05-07     2
3   1 1999-05-14     3
4   1 1999-05-21     4
5   1 1999-05-27     5
6   1 1999-06-04     6
> head(diff(data$stday))
Time differences in days
[1] 9 7 7 6 8 7


Эти данные имеют на один день меньше, и их надо дополнить NA вначале если присоединять в data.

3) Берем из diff(data$stday)[i+1] где i номер строки из data имеющий смысл и получаем номер столбца скользящей средней температуры: tmp[,diff(data$stday)[i+1]]

4) Теперь вставляем выборку даты data$stday[i]==temp$day, получаем итоговую конструкцию

tmp[ data$stday[i]==temp$day , diff(data$stday)[i+1] ]

которая по i тому номеру строки из первой таблицы возвращает среднюю температуру за предыдущий период

5) осмысленные строки -- по красивому надо написать агрегирующую функцию в plyr smile.gif но можно и проще
Код
sapply(2:nrow(data), function(i) if(data[i,1]==data[i-1,1]) tmp[ data$stday[i]==temp$day , diff(data$stday)[i+1] ] else NA )


Сообщение отредактировал p2004r - 30.08.2012 - 19:53


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
Olga_
сообщение 30.08.2012 - 21:16
Сообщение #12





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



Цитата(p2004r @ 30.08.2012 - 18:50) *
ну так почему бы и не выкладывать именно рабочий файл, раз уж Вы время привели к вменяемому формату?

ну тогда уж продолжу прошлый "алгоритм":

1) получаем массив средних за 3-14 последних дней на текущую дату (колонку даты возьмете свою smile.gif

Код
tmp<-sapply(1:10, function(i) avgtemp=runmean(temp.time$temp, i+3, alg= "exact", endrule="mean", align =  "right") )


2) разницу между текущим обследованием и предыдущим получаем

Код
>  data<-read.table("data.txt", header=TRUE)
>  data$stday<-as.Date(data$stday, "%d%b%Y")
> head(data)
  FID      stday woche
1   1 1999-04-28     1
2   1 1999-05-07     2
3   1 1999-05-14     3
4   1 1999-05-21     4
5   1 1999-05-27     5
6   1 1999-06-04     6
> head(diff(data$stday))
Time differences in days
[1] 9 7 7 6 8 7


Эти данные имеют на один день меньше, и их надо дополнить NA вначале если присоединять в data.

3) Берем из diff(data$stday)[i+1] где i номер строки из data имеющий смысл и получаем номер столбца скользящей средней температуры: tmp[,diff(data$stday)[i+1]]

4) Теперь вставляем выборку даты data$stday[i]==temp$day, получаем итоговую конструкцию

tmp[ data$stday[i]==temp$day , diff(data$stday)[i+1] ]

которая по i тому номеру строки из первой таблицы возвращает среднюю температуру за предыдущий период

5) осмысленные строки -- по красивому надо написать агрегирующую функцию в plyr smile.gif но можно и проще
Код
sapply(2:nrow(data), function(i) if(data[i,1]==data[i-1,1]) tmp[ data$stday[i]==temp$day , diff(data$stday)[i+1] ] else NA )



Огромное спасибо!


То что я хотела, но никак не могла выразить из-за недостатка знаний в R. Вы, конечно, правы, нужно сначала читать, а потом еще перечитывать.
Мне кажется, что вы меня уже не в первый раз выручаете. Как то регестрировалась я на одном форуме биологов н- ное время назад...

Вчера файлы никак загружались, выложила, те файлы, которые удалось загрузить.






Сообщение отредактировал Olga_ - 30.08.2012 - 21:17
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
Olga_
сообщение 11.09.2012 - 11:50
Сообщение #13





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



Что то у меня не получается..

Мой код

airpollution <- read.table("...Payerne.csv",blank.lines.skip = TRUE,fill=T,header=T,sep=";")
summary(airpollution)
day <- paste(airpollution[,3],"-",airpollution[,2],"-",airpollution[,1],sep="")
day <- as.Date(day, format="%Y-%m-%d")
payerne <- cbind(airpollution,day)

tmp<-sapply(1:10, function(i) avgtemp=runmean(payerne$temp, i+3, alg= "exact", endrule="mean", align = "right") )

data<- read.dta("...symptoms.dta")

data$stday<-as.Date(data$stday, "%d%b%Y")
head(diff(data$stday))
sapply(2:nrow(data), function(i) if(data[i,1]==data[i-1,1]) tmp[ data$stday[i]==payerne$day , diff(data$stday)[i+1] ] else NA )
Error in FUN(2:20426[[177L]], ...) : subscript out of bounds

Сообщение отредактировал Olga_ - 11.09.2012 - 11:51
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 12.09.2012 - 12:30
Сообщение #14





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



ну а файлы данных запаковать в rar и выложить?


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
Olga_
сообщение 6.11.2012 - 18:08
Сообщение #15





Группа: Пользователи
Сообщений: 43
Регистрация: 4.01.2012
Пользователь №: 23400



Решили оставить вариант с уровнем полютантов за предыдущие 7 дней.


Хотела бы использовать также для анализа время, в моем исследовании это 1 неделя. Но как было отмечено выше, время Исследований имеет неодинаковый промежуток.

1. Сначала, хотела бы создать одинаковые промежутки времени- недели начиная, к примеру с 01.01.1999. Какой это пакет в R? Где то кто то уже рекомендовал этот пакет, но не могу вспомнить

2. Затем интерполировать реальные даты в искусственные промежутки времени.

С помощью какого пакета/пакетов R можно это выполнить?

Спасибо заранее.

Сообщение отредактировал Olga_ - 6.11.2012 - 18:09
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 

2 страниц V   1 2 >
Добавить ответ в эту темуОткрыть тему