Интерполяция данных по гриппу, сохраняющих среднее значение за неделю

13

редактировать

Я нашел статью с описанием именно той процедуры, которая мне нужна. Единственное отличие состоит в том, что документ интерполирует среднемесячные данные на ежедневные, сохраняя при этом среднемесячные значения. У меня есть проблемы, чтобы реализовать подход в R. Любые намеки приветствуются.

оригинал

Для каждой недели у меня есть следующие данные подсчета (одно значение в неделю):

  • Количество консультаций врачей
  • Количество случаев гриппа

Моя цель - получать ежедневные данные путем интерполяции (я думал о линейных или усеченных сплайнах). Важно то, что я хочу сохранить среднее недельное значение, то есть среднее значение ежедневных интерполированных данных должно равняться зарегистрированному значению этой недели. Кроме того, интерполяция должна быть плавной. Одна проблема, которая может возникнуть, состоит в том, что определенная неделя имеет менее 7 дней (например, в начале или в конце года).

Буду признателен за совет по этому вопросу.

Большое спасибо.

Вот примерный набор данных за 1995 год ( обновлен ):

structure(list(daily.ts = structure(c(9131, 9132, 9133, 9134, 
9135, 9136, 9137, 9138, 9139, 9140, 9141, 9142, 9143, 9144, 9145, 
9146, 9147, 9148, 9149, 9150, 9151, 9152, 9153, 9154, 9155, 9156, 
9157, 9158, 9159, 9160, 9161, 9162, 9163, 9164, 9165, 9166, 9167, 
9168, 9169, 9170, 9171, 9172, 9173, 9174, 9175, 9176, 9177, 9178, 
9179, 9180, 9181, 9182, 9183, 9184, 9185, 9186, 9187, 9188, 9189, 
9190, 9191, 9192, 9193, 9194, 9195, 9196, 9197, 9198, 9199, 9200, 
9201, 9202, 9203, 9204, 9205, 9206, 9207, 9208, 9209, 9210, 9211, 
9212, 9213, 9214, 9215, 9216, 9217, 9218, 9219, 9220, 9221, 9222, 
9223, 9224, 9225, 9226, 9227, 9228, 9229, 9230, 9231, 9232, 9233, 
9234, 9235, 9236, 9237, 9238, 9239, 9240, 9241, 9242, 9243, 9244, 
9245, 9246, 9247, 9248, 9249, 9250, 9251, 9252, 9253, 9254, 9255, 
9256, 9257, 9258, 9259, 9260, 9261, 9262, 9263, 9264, 9265, 9266, 
9267, 9268, 9269, 9270, 9271, 9272, 9273, 9274, 9275, 9276, 9277, 
9278, 9279, 9280, 9281, 9282, 9283, 9284, 9285, 9286, 9287, 9288, 
9289, 9290, 9291, 9292, 9293, 9294, 9295, 9296, 9297, 9298, 9299, 
9300, 9301, 9302, 9303, 9304, 9305, 9306, 9307, 9308, 9309, 9310, 
9311, 9312, 9313, 9314, 9315, 9316, 9317, 9318, 9319, 9320, 9321, 
9322, 9323, 9324, 9325, 9326, 9327, 9328, 9329, 9330, 9331, 9332, 
9333, 9334, 9335, 9336, 9337, 9338, 9339, 9340, 9341, 9342, 9343, 
9344, 9345, 9346, 9347, 9348, 9349, 9350, 9351, 9352, 9353, 9354, 
9355, 9356, 9357, 9358, 9359, 9360, 9361, 9362, 9363, 9364, 9365, 
9366, 9367, 9368, 9369, 9370, 9371, 9372, 9373, 9374, 9375, 9376, 
9377, 9378, 9379, 9380, 9381, 9382, 9383, 9384, 9385, 9386, 9387, 
9388, 9389, 9390, 9391, 9392, 9393, 9394, 9395, 9396, 9397, 9398, 
9399, 9400, 9401, 9402, 9403, 9404, 9405, 9406, 9407, 9408, 9409, 
9410, 9411, 9412, 9413, 9414, 9415, 9416, 9417, 9418, 9419, 9420, 
9421, 9422, 9423, 9424, 9425, 9426, 9427, 9428, 9429, 9430, 9431, 
9432, 9433, 9434, 9435, 9436, 9437, 9438, 9439, 9440, 9441, 9442, 
9443, 9444, 9445, 9446, 9447, 9448, 9449, 9450, 9451, 9452, 9453, 
9454, 9455, 9456, 9457, 9458, 9459, 9460, 9461, 9462, 9463, 9464, 
9465, 9466, 9467, 9468, 9469, 9470, 9471, 9472, 9473, 9474, 9475, 
9476, 9477, 9478, 9479, 9480, 9481, 9482, 9483, 9484, 9485, 9486, 
9487, 9488, 9489, 9490, 9491, 9492, 9493, 9494, 9495), class = "Date"), 
    wdayno = c(0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 
    5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 
    6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 
    0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 
    1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 
    2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 
    3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 
    4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 
    5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 
    6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 
    0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 
    1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 
    2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 
    3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 
    4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 
    5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 
    6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 
    0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 
    1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 
    2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 
    3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 
    4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 
    5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 
    6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 
    0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L), month = c(1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 
    3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
    3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
    4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
    5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 
    6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 
    6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 
    8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
    8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 
    9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 
    9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 
    10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 
    10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 
    11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 
    11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
    12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
    12, 12, 12, 12), year = c(1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 
    1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995), yearday = 0:364, 
    no.influ.cases = c(NA, NA, NA, 168L, NA, NA, NA, NA, NA, 
    NA, 199L, NA, NA, NA, NA, NA, NA, 214L, NA, NA, NA, NA, NA, 
    NA, 230L, NA, NA, NA, NA, NA, NA, 267L, NA, NA, NA, NA, NA, 
    NA, 373L, NA, NA, NA, NA, NA, NA, 387L, NA, NA, NA, NA, NA, 
    NA, 443L, NA, NA, NA, NA, NA, NA, 579L, NA, NA, NA, NA, NA, 
    NA, 821L, NA, NA, NA, NA, NA, NA, 1229L, NA, NA, NA, NA, 
    NA, NA, 1014L, NA, NA, NA, NA, NA, NA, 831L, NA, NA, NA, 
    NA, NA, NA, 648L, NA, NA, NA, NA, NA, NA, 257L, NA, NA, NA, 
    NA, NA, NA, 203L, NA, NA, NA, NA, NA, NA, 137L, NA, NA, NA, 
    NA, NA, NA, 78L, NA, NA, NA, NA, NA, NA, 82L, NA, NA, NA, 
    NA, NA, NA, 69L, NA, NA, NA, NA, NA, NA, 45L, NA, NA, NA, 
    NA, NA, NA, 51L, NA, NA, NA, NA, NA, NA, 45L, NA, NA, NA, 
    NA, NA, NA, 63L, NA, NA, NA, NA, NA, NA, 55L, NA, NA, NA, 
    NA, NA, NA, 54L, NA, NA, NA, NA, NA, NA, 52L, NA, NA, NA, 
    NA, NA, NA, 27L, NA, NA, NA, NA, NA, NA, 24L, NA, NA, NA, 
    NA, NA, NA, 12L, NA, NA, NA, NA, NA, NA, 10L, NA, NA, NA, 
    NA, NA, NA, 22L, NA, NA, NA, NA, NA, NA, 42L, NA, NA, NA, 
    NA, NA, NA, 32L, NA, NA, NA, NA, NA, NA, 52L, NA, NA, NA, 
    NA, NA, NA, 82L, NA, NA, NA, NA, NA, NA, 95L, NA, NA, NA, 
    NA, NA, NA, 91L, NA, NA, NA, NA, NA, NA, 104L, NA, NA, NA, 
    NA, NA, NA, 143L, NA, NA, NA, NA, NA, NA, 114L, NA, NA, NA, 
    NA, NA, NA, 100L, NA, NA, NA, NA, NA, NA, 83L, NA, NA, NA, 
    NA, NA, NA, 113L, NA, NA, NA, NA, NA, NA, 145L, NA, NA, NA, 
    NA, NA, NA, 175L, NA, NA, NA, NA, NA, NA, 222L, NA, NA, NA, 
    NA, NA, NA, 258L, NA, NA, NA, NA, NA, NA, 384L, NA, NA, NA, 
    NA, NA, NA, 755L, NA, NA, NA, NA, NA, NA, 976L, NA, NA, NA, 
    NA, NA, NA, 879L, NA, NA, NA, NA), no.consultations = c(NA, 
    NA, NA, 15093L, NA, NA, NA, NA, NA, NA, 20336L, NA, NA, NA, 
    NA, NA, NA, 20777L, NA, NA, NA, NA, NA, NA, 21108L, NA, NA, 
    NA, NA, NA, NA, 20967L, NA, NA, NA, NA, NA, NA, 20753L, NA, 
    NA, NA, NA, NA, NA, 18782L, NA, NA, NA, NA, NA, NA, 19778L, 
    NA, NA, NA, NA, NA, NA, 19223L, NA, NA, NA, NA, NA, NA, 21188L, 
    NA, NA, NA, NA, NA, NA, 22172L, NA, NA, NA, NA, NA, NA, 21965L, 
    NA, NA, NA, NA, NA, NA, 21768L, NA, NA, NA, NA, NA, NA, 21277L, 
    NA, NA, NA, NA, NA, NA, 16383L, NA, NA, NA, NA, NA, NA, 15337L, 
    NA, NA, NA, NA, NA, NA, 19179L, NA, NA, NA, NA, NA, NA, 18705L, 
    NA, NA, NA, NA, NA, NA, 19623L, NA, NA, NA, NA, NA, NA, 19363L, 
    NA, NA, NA, NA, NA, NA, 16257L, NA, NA, NA, NA, NA, NA, 19219L, 
    NA, NA, NA, NA, NA, NA, 17048L, NA, NA, NA, NA, NA, NA, 19231L, 
    NA, NA, NA, NA, NA, NA, 20023L, NA, NA, NA, NA, NA, NA, 19331L, 
    NA, NA, NA, NA, NA, NA, 18995L, NA, NA, NA, NA, NA, NA, 16571L, 
    NA, NA, NA, NA, NA, NA, 15010L, NA, NA, NA, NA, NA, NA, 13714L, 
    NA, NA, NA, NA, NA, NA, 10451L, NA, NA, NA, NA, NA, NA, 14216L, 
    NA, NA, NA, NA, NA, NA, 16800L, NA, NA, NA, NA, NA, NA, 18305L, 
    NA, NA, NA, NA, NA, NA, 18911L, NA, NA, NA, NA, NA, NA, 17812L, 
    NA, NA, NA, NA, NA, NA, 18665L, NA, NA, NA, NA, NA, NA, 18977L, 
    NA, NA, NA, NA, NA, NA, 19512L, NA, NA, NA, NA, NA, NA, 17424L, 
    NA, NA, NA, NA, NA, NA, 14464L, NA, NA, NA, NA, NA, NA, 16383L, 
    NA, NA, NA, NA, NA, NA, 19916L, NA, NA, NA, NA, NA, NA, 18255L, 
    NA, NA, NA, NA, NA, NA, 20113L, NA, NA, NA, NA, NA, NA, 20084L, 
    NA, NA, NA, NA, NA, NA, 20196L, NA, NA, NA, NA, NA, NA, 20184L, 
    NA, NA, NA, NA, NA, NA, 20261L, NA, NA, NA, NA, NA, NA, 22246L, 
    NA, NA, NA, NA, NA, NA, 23030L, NA, NA, NA, NA, NA, NA, 10487L, 
    NA, NA, NA, NA)), .Names = c("daily.ts", "wdayno", "month", 
"year", "yearday", "no.influ.cases", "no.consultations"), row.names = c(NA, 
-365L), class = "data.frame")
r time-series interpolation COOLSerdash
источник
4
Этот вопрос требует одномерной версии интерполяции между точками , которая достаточно хорошо изучена в горнодобывающей промышленности. В указанном реферате прямо отмечается, что геостатистические методы дают «согласованные (сохраняющие массу ...) прогнозы». Я полагаю, что эти подходы преодолевают возражения, высказанные @Nick Cox.
whuber
@whuber Спасибо за ссылку, я не знал, что такого рода проблемы хорошо известны в геостатистике. Вам известно о какой-либо реализации таких методов в Rдругих статистических пакетах (у меня нет доступа к ArcGIS)? Боюсь, что без конкретной реализации я все еще застрял.
COOLSerdash
2
Я считаю, что это можно сделать с помощью кода geoRglm, если у вас есть очень хорошее понимание вариографии и изменения поддержки (что необходимо для разработки модели пространственной корреляции). Руководство издано Springer Verlag как Геостатистика на основе моделей, Diggle & Ribeiro Jr.
whuber
3
Разделение сгруппированных данных является обычной процедурой в демографии. Термин поиска "Sprague интерполяция"; это приведет вас ко многим вариантам. Подгоняя сплайн пятой степени к совокупным значениям таким образом, который обеспечивает монотонную кривую, этот метод и его варианты эффективно разделяют сгруппированные данные. (Он существует примерно с 1880 года.) Общий термин «осцилляционная интерполяция». Роб Хиндман, в частности, писал на эту тему: см. Смит, Хиндман и Вуд, Сплайн-интерполяция для демографических переменных: проблема монотонности, J. Pop. Местожительство 21 № 1 (2004), 95-98.
whuber
2
Ваш вопрос также можно рассматривать как дазиметрическое отображение в одном измерении. Это процедура для составления подробных карт величин, которые были измерены на некотором агрегированном уровне, таких как стандартные единицы переписи. (Это можно проследить, по крайней мере, до 1936 года: см. Джон К. Райт, «Метод картирования плотности населения: на примере Кейп-Код». Географический обзор 26: 1 (январь 1936), с. 103–110.) недавний подход (несколько специальный , но с краткой полезной библиографией) см. giscience.org/proceedings/abstracts/giscience2012_paper_179.pdf .
whuber

Ответы:

8

Мне удалось создать Rфункцию, которая интерполирует четные точки линейно и со сплайнами, сохраняя при этом средние значения (например, еженедельно, ежемесячно и т. Д.). Он использует функции na.approxи na.splineиз zooпакета и итеративно вычисляет сплайны с требуемыми свойствами. Алгоритм описан в этой статье .

Вот код:

interpol.consmean <- function(y, period=7, max.iter=100, tol=1e-4, plot=FALSE) {

  require(zoo)

  if( plot == TRUE ) {
    require(ggplot2)
  }

  y.temp.linear <- matrix(NA, ncol=length(y), nrow=max.iter+1)
  y.temp.linear[1, ] <- y

  y.temp.spline <- y.temp.linear

  y.temp.pred.spline <- matrix(NA, ncol=length(y), nrow=max.iter)
  y.temp.pred.linear <- matrix(NA, ncol=length(y), nrow=max.iter)

  ind.actual <- which(!is.na(y))

  if ( !all(diff(ind.actual)[1]== diff(ind.actual)) ) {
    stop("\"y\" must contain an evenly spaced time series")
  }

  partial <- ifelse((length(y) - ind.actual[length(ind.actual)]) < period/2,
                    TRUE, FALSE)

  for(k in 1:max.iter) {

    y.temp.pred.linear[k,] <- na.approx(y.temp.linear[k, ], na.rm=FALSE, rule=2)
    y.temp.pred.spline[k,] <- na.spline(y.temp.spline[k, ], method="fmm")

    interpol.means.linear <- rollapply(y.temp.pred.linear[k,], width=period, mean,
                                       by=period, align="left", partial=partial) 
    interpol.means.splines <- rollapply(y.temp.pred.spline[k,], width=period, mean,
                                        by=period, align="left", partial=partial) 

    resid.linear <- y.temp.linear[k, ][ ind.actual ] - interpol.means.linear
    resid.spline <- y.temp.spline[k, ][ ind.actual ] - interpol.means.splines

    if ( max(resid.linear, na.rm=TRUE) < tol & max(resid.spline, na.rm=TRUE) < tol ){
      cat("Converged after", k, "iterations with tolerance of", tol, sep=" ")
      break
    }

    y.temp.linear[k+1, ][!is.na(y.temp.linear[k, ])] <-  resid.linear
    y.temp.spline[k+1, ][!is.na(y.temp.spline[k, ])] <-  resid.spline

  }  

  interpol.linear.final <- colSums(y.temp.pred.linear, na.rm=TRUE)
  interpol.spline.final <- colSums(y.temp.pred.spline, na.rm=TRUE)

  if ( plot == TRUE ) {

    plot.frame <- data.frame(
      y=rep(y,2)/7,
      x=rep(1:length(y),2),
      inter.values=c(interpol.linear.final, interpol.spline.final)/7,
      method=c(rep("Linear", length(y)), rep("Spline", length(y)))
    )

    p <- ggplot(data=plot.frame, aes(x=x)) +
      geom_point(aes(y=y, x=x), size=4) +
      geom_line(aes(y=inter.values, color=method), size=1) +
      ylab("y") +
      xlab("x") +
      theme(axis.title.y =element_text(vjust=0.4, size=20, angle=90)) +
      theme(axis.title.x =element_text(vjust=0, size=20, angle=0)) +
      theme(axis.text.x =element_text(size=15, colour = "black")) +
      theme(axis.text.y =element_text(size=17, colour = "black")) +
      theme(panel.background =  element_rect(fill = "grey85", colour = NA),
            panel.grid.major =  element_line(colour = "white"),
            panel.grid.minor =  element_line(colour = "grey90", size = 0.25))+
      scale_color_manual(values=c("#377EB8", "#E41A1C"), 
                         name="Interpolation method",
                         breaks=c("Linear", "Spline"),
                         labels=c("Linear", "Spline")) +
      theme(legend.position="none") +
      theme(strip.text.x = element_text(size=16)) +
      facet_wrap(~ method)

    suppressWarnings(print(p))

  }
  list(linear=interpol.linear.final, spline=interpol.spline.final)
}

Давайте применим функцию к набору данных, приведенному в вопросе:

interpolations <- interpol.consmean(y=dat.frame$no.influ.cases, period=7,
                                    max.iter = 100, tol=1e-6, plot=TRUE)

Вставки

И линейная, и сплайн-интерполяция кажутся хорошими. Давайте проверим, сохранены ли еженедельные средние (усеченный вывод):

cbind(dat.frame$no.influ.cases[!is.na(dat.frame$no.influ.cases)],
      rollapply(interpolations$linear, 7, mean, by=7, align="left", partial=F))

      [,1] [,2]
 [1,]  168  168
 [2,]  199  199
 [3,]  214  214
 [4,]  230  230
 [5,]  267  267
 [6,]  373  373
 [7,]  387  387
 [8,]  443  443
 [9,]  579  579
[10,]  821  821
[11,] 1229 1229
COOLSerdash
источник
1
Вы должны найти подходящий пакет для этого и спросить сопровождающего, хотят ли они его включить.
Космонавт
4

Любая прямая линия, которая проходит через среднее значение в средней точке диапазона, будет давать дневные значения, которые имеют требуемое среднее значение. Последний комментарий Ника Кокса о «делении еженедельных отсчетов по количеству дней» является частным случаем с градиентом = 0.

Таким образом, мы можем отрегулировать это и выбрать градиент, чтобы сделать вещи немного более плавными. Вот три функции R, чтобы сделать что-то подобное:

interpwk <- function(x,y,delta){
  offset=-3:3
  yout=y+delta*offset
  xout=x+offset
  cbind(xout,yout)
}

get_delta <- function(x,y,pos){
  (y[pos+1]-y[pos-1])/(x[pos+1]-x[pos-1])
}

#' get slope from neighbours
interpall <- function(x,y,delta1,f=1){
  for(i in 2:(length(x)-1)){
    delta=get_delta(x,y,i)
    xyout=interpwk(x[i],y[i],delta/f)
    points(xyout)
  }
}

Добавьте меру дня к вашим данным, затем построите график, а затем построите интерполятор:

> data$day=data$week*7
> plot(data$day,data$no.influ.cases,type="l")
> interpall(data$day,data$no.influ.cases,f=1)

линейный средне-сохраняющий интерполятор

Другой возможностью является ограничение непрерывности в выходные дни, но это дает вам систему с одной степенью свободы - то есть она полностью определяется наклоном первого раздела (потому что тогда все остальные разделы должны объединиться). Я не закодировал это - ты можешь!

[Извините за слегка потертый код R, он должен действительно возвращать точки, а не строить их]

Spacedman
источник
+1, спасибо. Проблема в том, что интерполированные значения не являются гладкими, и между неделями происходят довольно резкие шаги. Я отредактировал свой вопрос, включая статью, которая в основном объясняет именно тот подход, который мне нужен.
COOLSerdash
Какова цель здесь? Почему предполагаемые случаи заболевания гриппом меняются плавно? Чем больше структуры вы вкладываете в эти данные путем интерполяции, тем больше придется вводить структуру на каком-то этапе моделирования. Я не думаю, что вы обратили внимание на мой комментарий от 19 мая «Увеличение еженедельных данных до ежедневных данных просто создает проблемы с введенной зависимостью и чрезмерно чрезмерно оптимистичными степенями свободы, которые будут мешать подгонке и оценке модели».
Ник Кокс
Ограничение среднего значения - это неправильно. Среднее значение, которое вы видите здесь, является средним для выборки и в некоторой степени подвержено статистическим изменениям. Придумайте модель, затем используйте интерполятор, который имеет среднее значение в качестве ожидаемого значения, затем выполните многократные вычисления ежедневных данных и проведите анализ сто или более раз, чтобы выяснить, как эта неопределенность влияет на ваши выводы.
Spacedman
1
@Spacedman Геостатистические методы API, на которые я ссылался (в комментарии к вопросу), будут обрабатывать ваше (вполне допустимое) возражение с помощью aplomb с помощью ненулевого компонента в параметре слепка вариограммы. Геостатистическое условное моделирование - это контролируемый метод выполнения нескольких импутаций, на которые вы ссылаетесь.
whuber
2
Абсолютно. Похоже, вы столкнулись с одномерной ситуацией, которая почти точно напоминает работающий пример в руководстве Diggle & Ribeiro для geoRglm (случаи малярии в Гамбии, с близостью к болотам и т. Д. В качестве ковариат). Основным осложнением является обработка смены поддержки, но это не повлияет на прогноз: в первую очередь это повлияет на оценку вариограммы. См. Ncbi.nlm.nih.gov/pmc/articles/PMC2995922 для ознакомления с некоторыми теориями и аналогичными примерами («биноминальное кригинг» случаев заболевания).
whuber
3

N

(Если бы данные были скорее измерениями, чем подсчетами, я бы склонялся к моделированию пропорций с помощью модели Дирихле, но это немного сложнее.)

Тот факт, что количество дней не всегда будет одинаковым, не должно быть особой проблемой, если вы знаете, что это такое - если вы используете смещение, чтобы поместить вещи на один и тот же «уровень».

Glen_b - Восстановить Монику
источник
1
Поправь меня, если я ошибаюсь, но я думаю, что это вопрос задом наперед. Это не то, как сгладить ежедневные счета; это как угадать ежедневные подсчеты по недельным данным. (Предположительно, на плакате есть ежедневные данные для чего-то другого, например, температуры.) Кроме этого, как это многочлен или Дирихле? Для меня это больше похоже на Пуассона.
Ник Кокс
@NickCox Вы абсолютно правы, спасибо за разъяснения: у меня есть еженедельные данные, и я хочу получать данные за день, потому что у меня есть другие данные за день (например, метеорологические переменные, смертность, загрязнение воздуха и т. Д.).
COOLSerdash
3
Моя собственная задача - спросить, почему вы хотите это сделать. Я полагаю, как и выше, что у вас есть некоторые ежедневные данные и вы хотите, чтобы все было одинаково. Если это так, рассмотрите некоторое сокращение ежедневных данных до минимума, среднего значения, медианы, максимума за недели или любого другого, что имеет научный смысл Наращивание еженедельных данных на ежедневные данные только создает проблемы с введенной зависимостью и чрезмерно чрезмерно оптимистичными степенями свободы, которые будут мешать подгонке и оценке модели.
Ник Кокс
@ Ник Кокс, это абсолютно «догадка», но на основании предоставленной информации, похоже, то, что было после ОП.
Glen_b
2
Другой консервативный подход состоит в том, чтобы просто делить еженедельные подсчеты на количество дней. Я знаю, что есть предположение, что реальный процесс будет более плавным, но он сохранит среднее значение.
Ник Кокс
3

Я соберу несколько дополнительных комментариев в качестве другого ответа.

Потребовалось время, чтобы структура этого проекта стала понятнее. Учитывая, что грипп теперь выявляется как один ковариат среди нескольких, вполне то, что вы делаете с этим, не кажется столь важным, или, по крайней мере, не заслуживает скептицизма, выраженного в некоторых из моих предыдущих комментариев. Как и все остальное на ежедневной основе, сокращение всего остального до недель приведет к выбрасыванию слишком большого количества деталей.

Первоначальный фокус вопроса остается на интерполяции, которая сохраняет недельное среднее значение, на которое один (крайний) ответ заключается в том, что недельное среднее значение сохраняет недельное среднее значение. Поскольку это неудивительно, что кажется непривлекательным или нереалистичным, другие методы интерполяции кажутся более привлекательными и / или методами вменения, предложенными @Spacedman. (Вполне понятно, будет ли это вменением с временным ароматом или интерполяцией с добавленным стохастическим ароматом.)

Еще две конкретные мысли:

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

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

Ник Кокс
источник