Interpolación de datos de influenza que conserva la media semanal

13

Editar

He encontrado un documento que describe exactamente el procedimiento que necesito. La única diferencia es que el documento interpola datos medios mensuales a diarios, al tiempo que conserva los medios mensuales. Tengo problemas para implementar el enfoque R. Cualquier pista es apreciada.

Original

Para cada semana, tengo los siguientes datos de conteo (un valor por semana):

  • Número de consultas médicas.
  • Número de casos de influenza.

Mi objetivo es obtener datos diarios por interpolación (pensé en splines lineales o truncadas). Lo importante es que quiero conservar la media semanal, es decir, la media de los datos interpolados diarios debería ser igual al valor registrado de esta semana. Además, la interpolación debe ser suave. Un problema que podría surgir es que cierta semana tiene menos de 7 días (por ejemplo, al comienzo o al final de un año).

Le agradecería consejos sobre este asunto.

Muchas gracias.

Aquí hay un conjunto de datos de muestra para el año 1995 ( actualizado ):

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
fuente
44
Esta pregunta solicita una versión unidimensional de la interpolación de área a punto , que está bastante bien estudiada en la industria minera. El resumen al que se hace referencia señala explícitamente que los métodos geoestadísticos arrojan "predicciones coherentes (preservando la masa ...)". Creo que estos enfoques superan las objeciones hechas por @Nick Cox.
whuber
@whuber Gracias por la referencia, no sabía que este tipo de problema es bien conocido en geoestadística. ¿Conoce alguna implementación de tales métodos en Ru otros paquetes estadísticos (no tengo acceso a ArcGIS)? Sin una implementación concretamente disponible, todavía estoy atascado, me temo.
COOLSerdash
2
Creo que esto podría hacerse utilizando el código geoRglm, siempre que comprenda muy bien la variografía y el cambio de soporte (que es necesario para desarrollar el modelo de correlación espacial). Springer Verlag publica el manual como Geoestadística basada en modelos, Diggle & Ribeiro Jr.
whuber
3
La subdivisión de datos agrupados es un procedimiento común en demografía. Un término de búsqueda es "interpolación Sprague"; te llevará a muchas variaciones. Al ajustar una spline de quinto grado a los valores acumulativos de una manera que asegure una curva monotónica, este método y sus variantes redividen efectivamente los datos agrupados. (Ha existido desde 1880). El término genérico es "interpolación osculatoria". Rob Hyndman, entre otros, ha escrito sobre este tema: véase Smith, Hyndman y Wood, Spline Interpolation for Demographic Variables: the Monotonicity Problem, J. Pop. Res. 21 No. 1 (2004), 95-98.
whuber
2
Su pregunta también se puede ver como mapeo dasimétrico en una dimensión. Este es un procedimiento para producir mapas detallados de cantidades que se han medido en algún nivel agregado, como las unidades censales estándar. (Se remonta al menos a 1936: ver John K. Wright, Un método para mapear las densidades de población: con Cape Cod como ejemplo. Revisión geográfica 26: 1 (enero de 1936), pp 103-110.) Para un enfoque reciente (algo ad hoc , pero con una breve bibliografía útil) ver giscience.org/proceedings/abstracts/giscience2012_paper_179.pdf .
whuber

Respuestas:

8

Me las arreglé para crear una Rfunción que interpola puntos espaciados de manera lineal y con splines mientras conserva los medios (por ejemplo, semanalmente, mensualmente, etc.). Utiliza las funciones na.approxy na.splinedel zoopaquete y calcula iterativamente las splines con las propiedades deseadas. El algoritmo se describe en este documento .

Aquí está el código:

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)
}

Apliquemos la función al conjunto de datos de ejemplo dado en la pregunta:

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

Interpolaciones

Tanto las interpolaciones lineales como las de spline parecen estar bien. Verifiquemos si se conservan las medias semanales (salida truncada):

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
fuente
1
Debería encontrar un paquete apropiado para eso y preguntarle al responsable si desea incluirlo.
Spacedman
4

Cualquier línea recta que atraviese la media en el punto medio del rango producirá valores diarios que tienen la media requerida. El último comentario de Nick Cox sobre 'dividir recuentos semanales por número de días' es un caso especial de eso con gradiente = 0.

Entonces podemos ajustar esto y elegir el gradiente para que las cosas sean un poco más suaves. Aquí hay tres funciones R para hacer algo así:

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)
  }
}

Agregue una medida de día a sus datos, luego trace, y luego trace el interpolador:

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

interpolador lineal de preservación media

Otra posibilidad es restringir la continuidad los fines de semana, pero esto le brinda un sistema con un solo grado de libertad, es decir, está completamente definido por la pendiente de la primera sección (porque entonces todas las otras secciones tienen que unirse). No he codificado esto, ¡tienes una oportunidad!

[Apols por el código R ligeramente en mal estado, realmente debería devolver los puntos en lugar de trazarlos]

Hombre espacial
fuente
+1, gracias. El problema es que los valores interpolados no son uniformes y hay pasos bastante abruptos entre las semanas. He editado mi pregunta, incluido un documento que básicamente explica exactamente el enfoque que necesito.
COOLSerdash
¿Cuál es el propósito aquí? ¿Por qué se presume que los casos de influenza varían sin problemas? Cuanta más estructura coloque en estos datos por interpolación, más tendrá que desenredarse la estructura introducida en alguna etapa de modelado. No creo que haya abordado mi comentario del 19 de mayo: "Recopilar datos semanales en datos diarios solo crea problemas con la dependencia introducida y grados de libertad excesivamente optimistas que perjudicarán el ajuste y la evaluación del modelo".
Nick Cox
Sin embargo, limitarse a la media está mal. La media que ve aquí es una media de muestra, y está sujeta a variación estadística de alguna manera. Elabore un modelo, luego use un interpolador que tenga la media como su expectativa, luego haga múltiples imputaciones de datos diarios y haga su análisis cien o más veces para descubrir cómo esta incertidumbre afecta sus conclusiones.
Spacedman 01 de
1
@Spacedman Los métodos de API geoestadística a los que me referí (en un comentario a la pregunta) manejarán su objeción (bastante válida) con aplomo, por medio de un componente distinto de cero en el parámetro nugget de variograma. Las simulaciones condicionales geoestadísticas son un método controlado para realizar las imputaciones múltiples a las que se refiere.
whuber
2
Absolutamente. Parece tener una situación unidimensional que es casi exactamente como un ejemplo corriente en el manual de Diggle & Ribeiro para geoRglm (casos de malaria en Gambia, con proximidad a pantanos, etc., como covariables). La principal complicación es manejar el cambio de soporte, pero eso realmente no afectaría la predicción: afectaría principalmente la estimación del variograma. Consulte ncbi.nlm.nih.gov/pmc/articles/PMC2995922 para obtener algunos ejemplos de teoría y similares ("kriging binomial" de casos de enfermedades).
whuber
3

norte

(Si los datos hubieran sido mediciones en lugar de conteos, me inclinaría por modelar las proporciones a través de un modelo de Dirichlet, pero eso es un poco más complicado).

El hecho de que el número de días no siempre sea el mismo no debería ser un problema particular, siempre y cuando sepas lo que es, siempre y cuando uses un desplazamiento para poner las cosas en el mismo 'nivel'.

Glen_b -Reinstate a Monica
fuente
1
Corrígeme si me equivoco, pero creo que esto tiene la pregunta al revés. No es cómo suavizar los recuentos diarios; es cómo adivinar recuentos diarios a partir de datos semanales. (Presumiblemente, el póster tiene datos diarios para otra cosa, por ejemplo, temperaturas). Aparte de eso, ¿cómo es este multinomial o Dirichlet? Me parece más un Poisson.
Nick Cox
@NickCox Tienes toda la razón, gracias por aclarar: tengo datos semanales y quiero datos diarios porque tengo otros datos diarios (es decir, variables meteorológicas, mortalidad, contaminación del aire, etc.).
COOLSerdash
3
Mi propia opinión sobre la pregunta es preguntar por qué quieres hacer esto. Supongo, como anteriormente, que tiene algunos datos diarios y quiere todo sobre la misma base. Si es así, considere alguna reducción de los datos diarios a mínimo, promedio, mediana, máximo durante semanas o lo que tenga sentido científico. Recopilar datos semanales en datos diarios solo crea problemas con la dependencia introducida y grados de libertad excesivamente optimistas que perjudicarán el ajuste y la evaluación del modelo.
Nick Cox
@Nick Cox es absolutamente "adivinar", pero en la información dada parece ser lo que buscaba el OP.
Glen_b -Reinstale a Monica
2
Otro enfoque conservador es dividir los recuentos semanales por la cantidad de días. Sé que hay una presuposición de que el proceso real será más fluido que eso, pero preservará la media.
Nick Cox
3

Agruparé algunos comentarios adicionales como otra respuesta.

La estructura de este proyecto ha tardado un poco en aclararse. Dado que la influenza ahora se revela como una covariable entre varias, lo que haces no parece tan crucial, o al menos no merece el escepticismo expresado en algunos de mis comentarios anteriores. Como todo lo demás es diario, reducir todo lo demás a semanas arrojaría demasiados detalles.

El enfoque original de la pregunta permanece, en la interpolación que conserva la media semanal a la que una respuesta (extrema) es que la media semanal conserva la media semanal. Como eso, como era de esperar, parece poco atractivo o poco realista, otros métodos de interpolación parecen más atractivos y / o métodos de imputación propuestos por @Spacedman. (No estoy claro si eso sería imputación con un sabor temporal o interpolación con sabor estocástico agregado).

Dos pensamientos específicos adicionales:

  • Tomar los valores semanales (divididos por el número de días) y luego suavizarlos con promedios ponderados probablemente en la práctica preservaría la media a una buena aproximación.

  • Como los casos de influenza son recuentos, suavizar los recuentos de raíz o registro y luego la transformación inversa podría funcionar mejor que simplemente suavizar los recuentos.

Nick Cox
fuente