Interpolation des données sur la grippe qui conserve la moyenne hebdomadaire

13

Éditer

J'ai trouvé un document décrivant exactement la procédure dont j'ai besoin. La seule différence est que le papier interpole les données moyennes mensuelles au quotidien, tout en préservant les moyennes mensuelles. J'ai du mal à mettre en œuvre l'approche en R. Tous les indices sont appréciés.

Original

Pour chaque semaine, j'ai les données de comptage suivantes (une valeur par semaine):

  • Nombre de consultations de médecins
  • Nombre de cas de grippe

Mon objectif est d'obtenir des données quotidiennes par interpolation (j'ai pensé aux splines linéaires ou tronquées). L'important est que je souhaite conserver la moyenne hebdomadaire, c'est-à-dire que la moyenne des données interpolées quotidiennes doit être égale à la valeur enregistrée de cette semaine. De plus, l'interpolation doit être fluide. Un problème qui pourrait survenir est qu'une certaine semaine compte moins de 7 jours (par exemple au début ou à la fin d'une année).

Je vous serais reconnaissant de me donner des conseils à ce sujet.

Merci beaucoup.

Voici un exemple d'ensemble de données pour l'année 1995 ( mis à jour ):

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
la source
4
Cette question demande une version unidimensionnelle de l' interpolation zone à point , qui est assez bien étudiée dans l'industrie minière. L'abrégé référencé note explicitement que les méthodes géostatistiques produisent des "prévisions cohérentes (préservant la masse ...)". Je crois que ces approches surmontent les objections formulées par @Nick Cox.
whuber
@whuber Merci pour la référence, je n'étais pas au courant que ce genre de problème est bien connu en géostatistique. Êtes-vous au courant de l'implémentation de telles méthodes dans Rou d'autres packages statistiques (je n'ai pas accès à ArcGIS)? Sans implémentation concrètement disponible, je suis toujours bloqué, je le crains.
COOLSerdash
2
Je crois que cela pourrait être fait en utilisant le code dans geoRglm, à condition que vous ayez une très bonne compréhension de la variographie et du changement de support (qui est nécessaire pour développer le modèle de corrélation spatiale). Le manuel est publié par Springer Verlag sous le nom de Model-Based Geostatistics, Diggle & Ribeiro Jr.
whuber
3
La subdivision des données groupées est une procédure courante en démographie. Un terme de recherche est "interpolation Sprague"; cela vous conduira à de nombreuses variations. En ajustant une spline du cinquième degré aux valeurs cumulatives de manière à assurer une courbe monotone, cette méthode et ses variantes redivisent efficacement les données groupées. (Elle existe depuis 1880.) Le terme générique est «interpolation osculatoire». Rob Hyndman, entre autres, a écrit sur ce sujet: voir Smith, Hyndman et Wood, Spline Interpolation for Demographic Variables: the Monotonicity Problem, J. Pop. Res. 21 n ° 1 (2004), 95-98.
whuber
2
Votre question peut également être considérée comme une cartographie dasymétrique en une dimension. Il s'agit d'une procédure pour produire des cartes détaillées des quantités qui ont été mesurées à un niveau agrégé, comme les unités de recensement standard. (Il peut être retracé au moins jusqu'en 1936: voir John K. Wright, Une méthode de cartographie des densités de population: avec Cape Cod comme exemple. Revue géographique 26: 1 (janvier 1936), pp 103-110.) approche récente (quelque peu ad hoc , mais avec une courte bibliographie utile) voir giscience.org/proceedings/abstracts/giscience2012_paper_179.pdf .
whuber

Réponses:

8

J'ai réussi à créer une Rfonction qui interpole des points régulièrement espacés de manière linéaire et avec des splines tout en préservant les moyens (par exemple hebdomadaire, mensuel, etc.). Il utilise les fonctions na.approxet à na.splinepartir du zoopackage et calcule de manière itérative les splines avec les propriétés souhaitées. L'algorithme est décrit dans cet article .

Voici le code:

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

Appliquons la fonction à l'exemple de jeu de données donné dans la question:

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

Interpolations

Les interpolations linéaires et splines semblent bonnes. Vérifions si les moyennes hebdomadaires sont préservées (sortie tronquée):

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
la source
1
Vous devriez trouver un package approprié pour cela et demander au responsable s'il souhaite l'inclure.
Spacedman
4

Toute ligne droite qui passe par la moyenne au milieu de la plage produira des valeurs quotidiennes qui ont la moyenne requise. Le dernier commentaire de Nick Cox sur «diviser les comptes hebdomadaires par nombre de jours» en est un cas particulier avec gradient = 0.

Nous pouvons donc ajuster cela et choisir le gradient pour rendre les choses peut-être un peu plus lisses. Voici trois fonctions R pour faire quelque chose comme ça:

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

Ajoutez une mesure journalière à vos données, puis tracez, puis tracez l'interpolateur:

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

interpolateur linéaire à conservation moyenne

Une autre possibilité est de contraindre la continuité le week-end, mais cela vous donne un système avec un seul degré de liberté - c'est-à-dire qu'il est complètement défini par la pente de la première section (car alors toutes les autres sections doivent se rejoindre). Je n'ai pas codé cela - vous avez un essai!

[Apole pour le code R légèrement minable, il devrait vraiment retourner les points plutôt que de les tracer]

Spacedman
la source
+1, merci. Le problème est que les valeurs interpolées ne sont pas lisses et qu'il y a des étapes assez abruptes entre les semaines. J'ai édité ma question, y compris un document qui explique essentiellement l'approche dont j'ai besoin.
COOLSerdash
Quel est le but ici? Pourquoi présumer que les cas de grippe varient en douceur? Plus vous mettez de structure dans ces données par interpolation, plus la structure introduite devra simplement être démêlée à un certain stade de la modélisation. Je ne pense pas que vous ayez répondu à mon commentaire du 19 mai. "Puffer des données hebdomadaires aux données quotidiennes crée simplement des problèmes avec la dépendance introduite et des degrés de liberté extrêmement optimistes qui entraveront l'ajustement et l'évaluation du modèle."
Nick Cox
Contraindre à la moyenne est cependant une erreur. La moyenne que vous voyez ici est une moyenne d'échantillon et est sujette à des variations statistiques d'une certaine manière. Conjurez un modèle, puis utilisez un interpolateur dont la moyenne est attendue, puis effectuez plusieurs imputations de données quotidiennes et effectuez votre analyse cent fois ou plus pour comprendre comment cette incertitude affecte vos conclusions.
Spacedman
1
@Spacedman Les méthodes de l'API géostatistique auxquelles j'ai fait référence (dans un commentaire à la question) traiteront votre objection (tout à fait valide) avec aplomb, au moyen d'une composante non nulle dans le paramètre de pépite de variogramme. Les simulations conditionnelles géostatistiques sont une méthode contrôlée pour effectuer les imputations multiples auxquelles vous faites référence.
whuber
2
Absolument. Vous semblez avoir une situation unidimensionnelle qui est presque exactement comme un exemple courant dans le manuel Diggle & Ribeiro pour geoRglm (cas de paludisme en Gambie, avec proximité de marais, etc., comme covariables). La principale complication est de gérer le changement de support, mais cela n'affecterait pas vraiment la prédiction: cela affecterait principalement l'estimation du variogramme. Voir ncbi.nlm.nih.gov/pmc/articles/PMC2995922 pour une théorie et des exemples similaires ("krigeage binomial" des cas de maladie).
whuber
3

n

(Si les données avaient été des mesures plutôt que des nombres, je pencherais pour la modélisation des proportions via un modèle de Dirichlet, mais c'est un peu plus compliqué.)

Le fait que le nombre de jours ne soit pas toujours le même ne devrait pas être un problème particulier, tant que vous savez de quoi il s'agit - tant que vous utilisez un décalage pour mettre les choses au même «niveau».

Glen_b -Reinstate Monica
la source
1
Corrigez-moi si je me trompe, mais je pense que la question est inversée. Ce n'est pas comment lisser les comptes quotidiens; c'est comment deviner les comptes quotidiens à partir des données hebdomadaires. (On peut supposer que l'affiche contient des données quotidiennes sur autre chose, par exemple les températures.) À part cela, comment est-ce multinomial ou Dirichlet? Ça ressemble plus à un poisson pour moi.
Nick Cox
@NickCox Vous avez tout à fait raison, merci pour la clarification: j'ai des données hebdomadaires et je veux des données quotidiennes parce que j'ai d'autres données quotidiennes (c.-à-d. Variables météorologiques, mortalité, pollution de l'air, etc.).
COOLSerdash
3
Ma propre interprétation de la question est de vous demander pourquoi vous voulez faire cela. Je suppose, comme ci-dessus, que vous avez des données quotidiennes et que vous voulez tout sur la même base. Dans l'affirmative, envisagez une certaine réduction des données quotidiennes à min, moyenne, médiane, max sur plusieurs semaines ou tout ce qui a un sens scientifique. Le fait de gonfler les données hebdomadaires aux données quotidiennes ne fait que créer des problèmes avec la dépendance introduite et les degrés de liberté extrêmement optimistes qui entraveront l'ajustement et l'évaluation du modèle.
Nick Cox
@ Nick Cox, c'est absolument "deviner", mais d'après les informations fournies, il semble que ce soit ce que le PO recherchait.
Glen_b -Reinstate Monica
2
Une autre approche prudente consiste simplement à diviser les comptes hebdomadaires par le nombre de jours. Je sais qu'il y a un présupposé que le véritable processus sera plus fluide que cela, mais il préservera la moyenne.
Nick Cox
3

Je vais regrouper quelques commentaires supplémentaires comme autre réponse.

Il a fallu un certain temps pour que la structure de ce projet devienne plus claire. Étant donné que la grippe est maintenant révélée comme une covariable parmi plusieurs, tout ce que vous faites ne semble pas si crucial, ou du moins ne mérite pas le scepticisme exprimé dans certains de mes commentaires précédents. Comme tout le reste est sur une base quotidienne, réduire tout le reste à des semaines jetterait trop de détails.

L'orientation initiale de la question demeure, sur l'interpolation qui préserve la moyenne hebdomadaire à laquelle une réponse (extrême) est que la moyenne hebdomadaire préserve la moyenne hebdomadaire. Comme cela semble sans surprise peu attrayant ou irréaliste, d'autres méthodes d'interpolation semblent plus attrayantes et / ou méthodes d'imputation comme proposé par @Spacedman. (Je ne sais pas exactement si ce serait une imputation avec une saveur temporelle ou une interpolation avec une saveur stochastique ajoutée.)

Deux autres réflexions spécifiques:

  • Le fait de prendre les valeurs hebdomadaires (divisées par le nombre de jours) puis de les lisser avec des moyennes pondérées serait en pratique susceptible de préserver la moyenne à une bonne approximation.

  • Étant donné que les cas de grippe sont des dénombrements, le lissage des dénombrements des racines ou des journaux, puis la rétrotransformation pourraient mieux fonctionner que le lissage des dénombrements.

Nick Cox
la source