gistfile1.txt

m1_gam_comp1 <- gamm4(comp1_provided~
                        s(curve_day)+
                        s(animals_actual)+
                        s(feed_provided_cum)+
                        s(comp2_provided)+
                        ogn_id+
                        newField,
                      data=training,
                      random=~(1+curve_day|lcn_id))
m1_gam_comp2 <- gamm4(comp2_provided~
                        s(curve_day)+
                        s(animals_actual)+
                        s(feed_provided_cum)+
                        s(comp1_provided)+
                        ogn_id+
                        newField,
                      data=training,
                      random=~(1+curve_day|lcn_id))
m1_gam_comp3 <- gamm4(comp3_provided~
                        s(curve_day)+
                        s(animals_actual)+
                        s(feed_provided_cum)+
                        s(comp1_provided)+
                        s(comp2_provided),
                      data=training,
                      random=~(1+curve_day|lcn_id))
testing$pred1<-predict(m1_gam_comp1$gam, 
                      newdata=testing)
testing$comp2_provided<-ifelse(is.na(testing$comp2_provided),0, 
                               testing$comp2_provided)
testing$pred2<-predict(m1_gam_comp2$gam, 
                      newdata=testing)
testing$comp3_provided<-ifelse(is.na(testing$comp3_provided),0, 
                               testing$comp3_provided)
testing$pred3<-predict(m1_gam_comp3$gam, 
                       newdata=testing)
testing%>%
  group_by(ogn_id, curve_day)%>%
  summarize(prov    = sum(feed_provided)*100,
            pred1   = sum(pred1),
            pred2   = sum(pred2),
            pred3   = sum(pred3),
            pred12  = pred1 + pred2,
            pred123 = pred1 + pred2 + pred3
            )%>%
  ggplot(., aes(x=curve_day))+
  geom_area(aes(y=prov, 
                fill=as.factor(ogn_id)), 
            alpha=0.5)+
  geom_point(aes(y=pred12), alpha=0.5)+
  geom_line(aes(y=pred12))+
  theme_bw()+
  labs(x="Day in Curve", 
       y="comp2", 
       fill="Farm", 
       title="Predicted vs Observed on Component 1+2 - TEST data")+
  facet_wrap(~ogn_id)