Отображение доверительных интервалов на графике

[Igg]

Здравствуйте!
В продолжении задачи: нужно построить в R линейный график с нанесенным диапазоном ошибки среднего (error bar)

[Igg]

+ файл с примером данных

[Максим Дубинин]

вот код, разумеется, так как в данных смешены наборы типа 2000 и 65000 получается на картинке полная фигня, но забавно Код снабдил кое-какими комментариями. Он не идеален, по мелочи можно еще поправлять и универсализировать, но работает хорошо.

Код: Выделить всё

#загружаем данные
data = read.table("D:\\Programming\\R\\error-bars\\ndvi2007source.csv",header=T)
#определим номер последнего поля ID
IDfieldpos = dim(data)[2]
#как обсуждалось в другом посте, сделаем массив средних и стандартных отклонений
means = aggregate(data[,-IDfieldpos],list(ID = data$ID),mean)
sds = aggregate(data[,-IDfieldpos],list(ID = data$ID),sd)
#определим сколько у нас уникальных ID, по ним будем итерировать
numID = nlevels(as.factor(data$ID))
#цикл
for (i in 1:numID) {
        #тут немного глючит у меня, конвертирую выборку из массива средних и СТД в простые единичные вектора
	amean = as.vector(means[i,-1],mode="integer")
	asds = as.vector(sds[i,-1],mode="integer")
        #это вектор 1-кол-во полей данных для error bars (это не numID)
	xs = seq(length(amean))
        #для каждого вектора значений строим график своего цвета
	plot(amean,type="l",xlab="Julian Day",xaxt='n',ylab="NDVI",col=i,ylim=c(0,max(data)))
	axis(1,at=seq(length(amean)),labels=substr(names(means)[-1],2,4),las=3)
        #стрелки - error bars вверх-вниз на стандартное отклонение
	arrows(xs,amean-asds,xs,amean+asds,code=3,angle=90,length=0.1,col=i)
	par(new=T)
}

вот что у меня получилось с картинкой

[Igg]

Sim, спасибо!
все работает)

правда у меня почему то выражение
means = aggregate(data[,-IDfieldpos],list(ID = data$ID),mean)
работает только в таком виде
means = aggregate(data[,-IDfieldpos],list(ID = data$ID),FUN = "mean")
(иначе выдет ошибку...)

вот что получилось(см. вложение)

meadows_steppe_pine_birch.pdf

(19.91 КБ) 775 скачиваний

[Максим Дубинин]

правда у меня почему то выражение
means = aggregate(data[,-IDfieldpos],list(ID = data$ID),mean)
работает только в таком виде
means = aggregate(data[,-IDfieldpos],list(ID = data$ID),FUN = "mean")
(иначе выдет ошибку...)

Странно, у меня 2.8.1 и все работает и так и так без ошибок.