Nothing
inst/script/influenza.R
In fda.usc: Functional Data Analysis and Utilities for Statistical Computing
library(fda.usc)
################################################################################
load("influenza.rda")
nt<-53*52*2
nrep<-156
a1<-209+ 52
a2<-a1+nrep
nmodel<-12
flu2<-flu
# selection response variable (n+1, to n+5)
ii0<-c(12,21,20,19,18)
tres <- flu$temp$data
tres[tres>10] <- 0
tres[flu$temp$data<=10] <- flu$temp$data[flu$temp$data<=10]-10
ff1 <- y~ratet13
ff2 <- y~temp
ff3 <- y~tres
ff4 <- y~irra
ff5 <- y~ratet13 + temp
ff6 <- y~ratet13 + tres
ff7 <- y~ratet13 + irra
ff8 <- y~temp + irra
nmodel<-24
ahead<-nahead<-4
a1<-209+ 52
a2<-a1+53
a1<-209+ 26
a2<-a1+53+26
nrep<-500
indi<-1:nrep
iispace<-3500
indspace<-flu$df$ispace>iispace
indspace2<-sum(unique(flu$df$ispace)>iispace)
pred.res<-array(NA,dim=c(nahead,nrep,nmodel+1,indspace2)) #nmodel +1 (la y)
pred0<-pred<-array(NA,dim=c(nahead,nrep,nmodel))
lev<-names(table(flu$df$ispace))
indi3<-flu$df$ispace<16000000
# 10 o 53 comarcas
nrep<-500
######epidemic<-1:156 # in the paper, high time-consuming
## Moderate Time-consuming (only 17 rolling iterations (last 17 weeks)
p1<-1
q1<-0
p2<-1
q2<-0
ii1<-c(12,21,20,19)
## High Time-consuming
levcomar<-names(table(flu$df$ispace[indspace]))
ncomar<-length(levcomar)
epidemic<-c(145:164+52,145:164+2*52)
for (i in epidemic) {
for (icomar in 1:ncomar) {
j<-1
flu2<-flu
flu$df$y<-flu$df[,ii1[j]]
ind<-which(flu$df$itime>=(i+a1) & flu$df$itime<=(i+a2)
& flu$df$ispace==levcomar[icomar])
flu2$df<-flu$df[ind,]
flu2$temp<-flu$temp[ind]
flu2$tmp.d1<-flu$tmp.d1[ind]
flu2$temp.d2<-flu$temp.d2[ind]
flu2$temt10<-flu$tempt10[ind]
flu2$ratet13<-flu$ratet13[ind]
flu2$temp.max<-flu$temp.max[ind]
flu2$temp.min<-flu$temp.min[ind]
flu2$vtemp<-flu$rtemp[ind]
flu2$hum<-flu$hum[ind,8:14]
flu2$hum$rangeval<-c(8,14)
flu2$tempt10<-flu$tempt10[ind,]
flu2$tres<-flu$tres[ind,]
flu2$trs.d1<-flu$trs.d1[ind,]
flu2$irra<-flu$irra[ind]
flu2$humt10<-flu$humt10[ind,]
res.lm1<-fregre.lm(ff1,data=flu2)
res.lm2<-fregre.lm(ff2,data=flu2)
res.lm3<-fregre.lm(ff3,data=flu2)
res.lm4<-fregre.lm(ff4,data=flu2)
res.lm5<-fregre.lm(ff5,data=flu2)
res.lm6<-fregre.lm(ff6,data=flu2)
res.lm7<-fregre.lm(ff7,data=flu2)
res.lm8<-fregre.lm(ff8,data=flu2)
# par.cor1<-list("cor.ARMA"=list("index"=c("itime"),"group"="ispace",method="lm","p"=1)r)
# res.fou5<-fregre.igls(ff2,data=flu2,correlation=par.cor1)
res1<-fregre.gls(ff1,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res2<-fregre.gls(ff2,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res3<-fregre.gls(ff3,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res4<-fregre.gls(ff4,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res5<-fregre.gls(ff5,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res6<-fregre.gls(ff6,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res7<-fregre.gls(ff7,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res8<-fregre.gls(ff8,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res.gls1<-fregre.gls(ff1,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls2<-fregre.gls(ff2,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls3<-fregre.gls(ff3,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls4<-fregre.gls(ff4,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls5<-fregre.gls(ff5,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls6<-fregre.gls(ff6,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls7<-fregre.gls(ff7,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls8<-fregre.gls(ff8,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
ind2<-which(flu$df$itime==(i+a2+1)& flu$df$ispace==levcomar[icomar])
flu3<-flu
flu3$df$y<-flu$df[,ii0[j]]
flu3$df<-flu$df[ind2,,drop=FALSE]
flu3$temp<-flu$temp[ind2]
flu3$tmp.d1<-flu$tmp.d1[ind2]
flu3$tempt10<-flu$tempt10[ind2]
flu3$tres<-flu$tres[ind2]
flu3$trs.d1<-flu$trs.d1[ind2]
flu3$temp.min<-flu$temp.min[ind2]
flu3$temp.max<-flu$temp.max[ind2]
flu3$ratet13<-flu$ratet13[ind2]
flu3$hum<-flu$hum[ind2,8:14]
flu3$hum$rangeval<-c(8,14)
flu3$irra<-flu$irra[ind2]
# prediction
pred.lm1<-predict.fregre.lm(res.lm1,flu3)
pred.lm2<-predict.fregre.lm(res.lm2,flu3)
pred.lm3<-predict.fregre.lm(res.lm3,flu3)
pred.lm4<-predict.fregre.lm(res.lm4,flu3)
pred.lm5<-predict.fregre.lm(res.lm5,flu3)
pred.lm6<-predict.fregre.lm(res.lm6,flu3)
pred.lm7<-predict.fregre.lm(res.lm7,flu3)
pred.lm8<-predict.fregre.lm(res.lm8,flu3)
pred1<-predict.fregre.gls(res1,flu3,n.ahead=ahead)
pred2<-predict.fregre.gls(res2,flu3,n.ahead=ahead)
pred3<-predict.fregre.gls(res3,flu3,n.ahead=ahead)
pred4<-predict.fregre.gls(res4,flu3,n.ahead=ahead)
pred5<-predict.fregre.gls(res5,flu3,n.ahead=ahead)
pred6<-predict.fregre.gls(res6,flu3,n.ahead=ahead)
pred7<-predict.fregre.gls(res7,flu3,n.ahead=ahead)
pred8<-predict.fregre.gls(res8,flu3,n.ahead=ahead)
yresp<-t(flu3$df[,c("ratet1","ratet2","ratet3","ratet4")])
pred.res[,i,1,icomar]<-rbind(pred.lm1,pred.lm1,pred.lm1,pred.lm1)
pred.res[,i,2,icomar]<-rbind(pred.lm2,pred.lm2,pred.lm2,pred.lm2)
pred.res[,i,3,icomar]<-rbind(pred.lm3,pred.lm3,pred.lm3,pred.lm3)
pred.res[,i,4,icomar]<-rbind(pred.lm4,pred.lm4,pred.lm4,pred.lm4)
pred.res[,i,5,icomar]<-rbind(pred.lm5,pred.lm5,pred.lm5,pred.lm5)
pred.res[,i,6,icomar]<-rbind(pred.lm6,pred.lm6,pred.lm6,pred.lm6)
pred.res[,i,7,icomar]<-rbind(pred.lm7,pred.lm7,pred.lm7,pred.lm7)
pred.res[,i,8,icomar]<-rbind(pred.lm8,pred.lm8,pred.lm8,pred.lm8)
pred.res[,i,9,icomar]<-t(pred1)
pred.res[,i,10,icomar]<-t(pred2 )
pred.res[,i,11,icomar]<-t(pred3)
pred.res[,i,12,icomar]<-t(pred4)
pred.res[,i,13,icomar]<-t(pred5)
pred.res[,i,14,icomar]<-t(pred6 )
pred.res[,i,15,icomar]<-t(pred7 )
pred.res[,i,16,icomar]<-t(pred8 )
pred.res[,i,17,icomar]<-yresp
}
iresp <- 17
pred[,i,1]<-apply((pred.res[,i,iresp,]-pred.res[,i,1,])^2,1,mean)
pred[,i,2]<-apply((pred.res[,i,iresp,]-pred.res[,i,2,])^2,1,mean)
pred[,i,3]<-apply((pred.res[,i,iresp,]-pred.res[,i,3,])^2,1,mean)
pred[,i,4]<-apply((pred.res[,i,iresp,]-pred.res[,i,4,])^2,1,mean)
pred[,i,5]<-apply((pred.res[,i,iresp,]-pred.res[,i,5,])^2,1,mean)
pred[,i,6]<-apply((pred.res[,i,iresp,]-pred.res[,i,6,])^2,1,mean)
pred[,i,7]<-apply((pred.res[,i,iresp,]-pred.res[,i,7,])^2,1,mean)
pred[,i,8]<-apply((pred.res[,i,iresp,]-pred.res[,i,8,])^2,1,mean)
pred[,i,9]<-apply((pred.res[,i,iresp,]-pred.res[,i,9,])^2,1,mean)
pred[,i,10]<-apply((pred.res[,i,iresp,]-pred.res[,i,10,])^2,1,mean)
pred[,i,11]<-apply((pred.res[,i,iresp,]-pred.res[,i,11,])^2,1,mean)
pred[,i,12]<-apply((pred.res[,i,iresp,]-pred.res[,i,12,])^2,1,mean)
pred[,i,13]<-apply((pred.res[,i,iresp,]-pred.res[,i,13,])^2,1,mean)
pred[,i,14]<-apply((pred.res[,i,iresp,]-pred.res[,i,14,])^2,1,mean)
pred[,i,15]<-apply((pred.res[,i,iresp,]-pred.res[,i,15,])^2,1,mean)
pred[,i,16]<-apply((pred.res[,i,iresp,]-pred.res[,i,16,])^2,1,mean)
print( i)
pred[,i,]
}
save.image(file="fluExample.RData")
###############################################################################
load(file="fluExample.RData")
table6<-t(round(apply(pred[,epidemic,1:16],c(1,3),mean,na.rm=T),2))
table6<-cbind(table6[c(1:8),1],table6[c(9:16),1],table6[c(1:8),2],table6[c(9:16),2])
colnames(table6)<-c("FLM (t+1)","FGLS (t+1)","FLM (t+2)","FGLS (t+2)")
rownames(table6)<-c("Rate(w)","Temp(t)","Temp.thres(t)","SR(t)",
"Rate(w),Temp(t)","Rate(w),Temp(t)",
"Rate(w)Temp.thres(t)","Temp(t),SR(t)")
table6
###############################################################################
j<-1
flu2<-flu
flu$df$y<-flu$df[,ii1[j]]
ind<-which(flu$df$itime>=(epidemic[1]+a1) & flu$df$itime<=(epidemic[40]+a2)
& flu$df$ispace %in% levcomar)
length(ind)
flu2$df<-flu$df[ind,]
flu2$temp<-flu$temp[ind]
flu2$tmp.d1<-flu$tmp.d1[ind]
flu2$temp.d2<-flu$temp.d2[ind]
flu2$temt10<-flu$tempt10[ind]
flu2$ratet13<-flu$ratet13[ind]
flu2$temp.max<-flu$temp.max[ind]
flu2$temp.min<-flu$temp.min[ind]
flu2$vtemp<-flu$rtemp[ind]
flu2$hum<-flu$hum[ind,8:14]
flu2$hum$rangeval<-c(8,14)
flu2$tempt10<-flu$tempt10[ind,]
flu2$tres<-flu$tres[ind,]
flu2$trs.d1<-flu$trs.d1[ind,]
flu2$irra<-flu$irra[ind]
flu2$humt10<-flu$humt10[ind,]
res1<-fregre.gls(ff1,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res2<-fregre.gls(ff2,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res3<-fregre.gls(ff3,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res4<-fregre.gls(ff4,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res5<-fregre.gls(ff5,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res6<-fregre.gls(ff6,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res7<-fregre.gls(ff7,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res8<-fregre.gls(ff8,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
################# INICIO #####################
add_legend <- function(...) {
opar <- par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0),
mar=c(0, 0, 0, 0), new=TRUE)
on.exit(par(opar))
plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
legend(...)
}
pdf(file = "S2_fig.pdf",width = 10.24, height =7.68, pointsize = 10)
par(mfrow=c(2,2))
cols=c("blue","sky blue","red","dark red")
# Rate
covar<- "ratet13"
res <-res1
x<-flu2[[covar]]
a<-fdata(res$beta.l[[covar]],flu2[[covar]]$argval,flu2[[covar]]$rangeval)
b1<-b<- +inprod.fdata(x,a)
qq<-quantile(b,prob=c(.25,.5,.75,1))
ind1<-ifelse(b[,1]<qq[1],TRUE,FALSE)
ind2<-ifelse(b[,1]>=qq[1],TRUE,FALSE)
ind3<-ifelse(b[,1]>=qq[2],TRUE,FALSE)
ind4<-ifelse(b[,1]>=qq[3],TRUE,FALSE)
mn0<-expression(paste( beta[rate]))
ylab <- "Rate(w)"
ylab <- expression(Rate["s,n"](w))
mn2<-expression(paste(symbol("<"), list(Rate["s,n"](w), beta[Rate]),symbol(">")))
plot(c((func.mean(x[ind1,])),(func.mean(x[ind2,]))
,(func.mean(x[ind3,])),(func.mean(x[ind4,])))
,lwd=2,lty=1
,col=cols,
main=mn2,xlab="Range time w: 13 weeks",ylab=ylab)
# Temp
covar<- "temp"
res <-res2
x<-flu2[[covar]]
a<-fdata(res$beta.l[[covar]],flu2[[covar]]$argval,flu2[[covar]]$rangeval)
b1<-b<- +inprod.fdata(x,a)
qq<-quantile(b,prob=c(.25,.5,.75,1))
ind1<-ifelse(b[,1]<qq[1],TRUE,FALSE)
ind2<-ifelse(b[,1]>=qq[1],TRUE,FALSE)
ind3<-ifelse(b[,1]>=qq[2],TRUE,FALSE)
ind4<-ifelse(b[,1]>=qq[3],TRUE,FALSE)
mn0<-expression(paste( beta[rate]))
ylab <- expression(Temp["s,n"](t))
mn2<-expression(paste(symbol("<"), list(Temp["s,n"](t), beta[Temp]),symbol(">")))
plot(c((func.mean(x[ind1,])),(func.mean(x[ind2,]))
,(func.mean(x[ind3,])),(func.mean(x[ind4,])))
,lwd=2,lty=1
,col=cols,
main=mn2,xlab="Range time t: 14 days",ylab=ylab)
# Temp
covar<- "tres"
res <-res3
x<-flu2[[covar]]
a<-fdata(res$beta.l[[covar]],flu2[[covar]]$argval,flu2[[covar]]$rangeval)
b1<-b<- +inprod.fdata(x,a)
qq<-quantile(b,prob=c(.25,.5,.75,1))
ind1<-ifelse(b[,1]<qq[1],TRUE,FALSE)
ind2<-ifelse(b[,1]>=qq[1],TRUE,FALSE)
ind3<-ifelse(b[,1]>=qq[2],TRUE,FALSE)
ind4<-ifelse(b[,1]>=qq[3],TRUE,FALSE)
ylab <- expression(Temp.thres["s,n"](t))
mn2<-expression(paste(symbol("<"), list(Temp.thres["s,n"](t), beta[Temp.thres]),symbol(">")))
plot(c((func.mean(x[ind1,])),(func.mean(x[ind2,]))
,(func.mean(x[ind3,])),(func.mean(x[ind4,])))
,lwd=2,lty=1
,col=cols,
# main=mn2,xlab="t: previous 14 days",ylab=ylab)
main=mn2,xlab="t in days",ylab=ylab)
# Solar Radiation
covar<- "irra"
res <-res4
x<-flu2[[covar]]
a<-fdata(res$beta.l[[covar]],flu2[[covar]]$argval,flu2[[covar]]$rangeval)
b1<-b<- +inprod.fdata(x,a)
qq<-quantile(b,prob=c(.25,.5,.75,1))
ind1<-ifelse(b[,1]<qq[1],TRUE,FALSE)
ind2<-ifelse(b[,1]>=qq[1],TRUE,FALSE)
ind3<-ifelse(b[,1]>=qq[2],TRUE,FALSE)
ind4<-ifelse(b[,1]>=qq[3],TRUE,FALSE)
ylab <- expression(SR["s,n"](t))
mn2<-expression(paste(symbol("<"), list(SR["s,n"](t), beta[SR]),symbol(">")))
plot(c((func.mean(x[ind1,])),(func.mean(x[ind2,]))
,(func.mean(x[ind3,])),(func.mean(x[ind4,])))
,lwd=2,lty=1
,col=cols,
main=mn2,xlab="t (in days)",ylab=ylab)
leg <- c(expression(v[q[1]]), expression(v[q[2]]), expression(v[q[3]]) ,expression(v[q[4]]))
#legend(8,2500,inset=c(-0.2,0),col=cols,legend=leg,cex=.8,bty="n",lty=1,lwd=1)
add_legend("bottom", legend=leg,lty=1, #pch=20,
col=cols, lwd=2,
horiz=TRUE, bty='n',
cex=0.95)
dev.off()
################# fin Gráfico #####################
################################################################################
################################################################################
# Figure 1 Series Temporales
flu4<-aggregate(flu2$df,list(flu2$df[,"itime"]),mean)
#cambiar: utilizar un ggplot o algo así más chulo (ver links/linkedIN)
png(file = "S1_fig.png")
par(mfrow=c(4,1))
FluRate<-ts(exp(flu4[,"ratet0"]),2001,2011,53)
Temperature<-ts(flu4[,"tmed"] ,2001,2011,53)
SolarRadiation<-ts(flu4[,"imed"] ,2001,2011,53)
RelativeHumidity<-ts(flu4[,"hmed"] ,2001,2011,53)
TemperatureThreshold<-ts(flu4[,"thold"] ,2001,2011,53)
aa<-cbind(FluRate,Temperature,SolarRadiation,RelativeHumidity)
plot(aa,main="",xlab="Years")
dev.off()
################################################################################
################################################################################
################################################################################
################################################################################
# Pearson and distance correlation
################################################################################
################################################################################
#Upper diagonal matrix of Table 1: Pearson correlation
# no se muestran
A<-cor(flu2$df[,c(19,12,22,13,23,26,29)])
################################################################################
#Lower diagonal matrix of Table 1: Distance correlation
################################################################################
# Distance computation for Multivariate-Multivariate data
names(flu)
A1<-metric.lp(flu2$ratet13)
A2<-metric.lp(flu2$temp)
A3<-metric.lp(flu2$tres)
A4<-metric.lp(flu2$irra)
A5<-metric.lp(flu2$hum)
A6<-metric.dist(data.matrix(flu2$df$ratet0))
A7<-metric.dist(data.matrix(flu2$df$ratet1))
A8<-metric.dist(data.matrix(flu2$df$ratet2))
A9<-metric.dist(data.matrix(flu2$df$tmed))
A10<-metric.dist(data.matrix(flu2$df$thold))
A11<-metric.dist(data.matrix(flu2$df$imed))
A12<-metric.dist(data.matrix(flu2$df$hmed))
# Distance correlation
A<-cor(flu2$df[,c(21,12,22,13,23,26,29)])
names(flu2$df[,c(21,12,22,13,23,26,29)])
A[2,1]<-dcor.dist(A8,A7)
A[3,1]<-dcor.dist(A8,A6)
A[4,1]<-dcor.dist(A8,A9)
A[5,1]<-dcor.dist(A8,A10)
A[6,1]<-dcor.dist(A8,A11)
A[7,1]<-dcor.dist(A8,A12)
A[3,2]<-dcor.dist(A7,A6)
A[4,2]<-dcor.dist(A7,A9)
A[5,2]<-dcor.dist(A7,A10)
A[6,2]<-dcor.dist(A7,A11)
A[7,2]<-dcor.dist(A7,A12)
A[4,3]<-dcor.dist(A6,A9)
A[5,3]<-dcor.dist(A6,A10)
A[6,3]<-dcor.dist(A6,A11)
A[7,3]<-dcor.dist(A6,A12)
A[5,4]<-dcor.dist(A9,A10)
A[6,4]<-dcor.dist(A9,A11)
A[7,4]<-dcor.dist(A9,A12)
A[6,5]<-dcor.dist(A10,A11)
A[7,5]<-dcor.dist(A10,A12)
A[7,6]<-dcor.dist(A11,A12)
round(A,2)
################################################################################
################################################################################
B<-matrix(NA,5,12)
# Distance computation for Functional-Multivariate data
B[1,1]<-dcor.dist(A8,A1)
B[2,1]<-dcor.dist(A8,A2)
B[3,1]<-dcor.dist(A8,A3)
B[4,1]<-dcor.dist(A8,A4)
B[5,1]<-dcor.dist(A8,A5)
B[1,2]<-dcor.dist(A7,A1)
B[2,2]<-dcor.dist(A7,A2)
B[3,2]<-dcor.dist(A7,A3)
B[4,2]<-dcor.dist(A7,A4)
B[5,2]<-dcor.dist(A7,A5)
B[1,3]<-dcor.dist(A6,A1)
B[2,3]<-dcor.dist(A6,A2)
B[3,3]<-dcor.dist(A6,A3)
B[4,3]<-dcor.dist(A6,A4)
B[5,3]<-dcor.dist(A6,A5)
B[1,4]<-dcor.dist(A9,A1)
B[2,4]<-dcor.dist(A9,A2)
B[3,4]<-dcor.dist(A9,A3)
B[4,4]<-dcor.dist(A9,A4)
B[5,4]<-dcor.dist(A9,A5)
B[1,5]<-dcor.dist(A10,A1)
B[2,5]<-dcor.dist(A10,A2)
B[3,5]<-dcor.dist(A10,A3)
B[4,5]<-dcor.dist(A10,A4)
B[5,5]<-dcor.dist(A10,A5)
B[1,6]<-dcor.dist(A11,A1)
B[2,6]<-dcor.dist(A11,A2)
B[3,6]<-dcor.dist(A11,A3)
B[4,6]<-dcor.dist(A11,A4)
B[5,6]<-dcor.dist(A11,A5)
B[1,7]<-dcor.dist(A12,A1)
B[2,7]<-dcor.dist(A12,A2)
B[3,7]<-dcor.dist(A12,A3)
B[4,7]<-dcor.dist(A12,A4)
B[5,7]<-dcor.dist(A12,A5)
# Distance computation for Functional-Functional data
B[1,8]<- B[2,9]<- B[3,10]<- B[4,11]<- B[5,12]<-1
B[2,8]<-dcor.dist(A1,A2)
B[3,8]<-dcor.dist(A1,A3)
B[4,8]<-dcor.dist(A1,A4)
B[5,8]<-dcor.dist(A1,A5)
B[3,9]<-dcor.dist(A2,A3)
B[4,9]<-dcor.dist(A2,A4)
B[5,9]<-dcor.dist(A2,A5)
B[4,10]<-dcor.dist(A3,A4)
B[5,10]<-dcor.dist(A3,A5)
B[5,11]<-dcor.dist(A4,A5)
B[1,8:12]<-B[,8]
B[2,9:12]<-B[2:5,9]
B[3,10:12]<-B[3:5,10]
B[4,11:12]<-B[4:5,11]
rnames<-c("raten+2","raten+1","rate","Temp","Thold","Irra","Hum")
fnames<-c("rate(w)","Temp(t)","Thold(t)","Irra(t)","Hum(t)")
colnames(B)<-c(rnames,fnames)
rownames(B)<-fnames
# TABLE 5
round(B[,c(8:12,2,1)],2)
save.image(file="fluExample.RData")
Try the fda.usc package in your browser
Any scripts or data that you put into this service are public.
fda.usc documentation built on Oct. 17, 2022, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(fda.usc)
################################################################################
load("influenza.rda")
nt<-53*52*2
nrep<-156
a1<-209+ 52
a2<-a1+nrep
nmodel<-12
flu2<-flu
# selection response variable (n+1, to n+5)
ii0<-c(12,21,20,19,18)
tres <- flu$temp$data
tres[tres>10] <- 0
tres[flu$temp$data<=10] <- flu$temp$data[flu$temp$data<=10]-10
ff1 <- y~ratet13
ff2 <- y~temp
ff3 <- y~tres
ff4 <- y~irra
ff5 <- y~ratet13 + temp
ff6 <- y~ratet13 + tres
ff7 <- y~ratet13 + irra
ff8 <- y~temp + irra
nmodel<-24
ahead<-nahead<-4
a1<-209+ 52
a2<-a1+53
a1<-209+ 26
a2<-a1+53+26
nrep<-500
indi<-1:nrep
iispace<-3500
indspace<-flu$df$ispace>iispace
indspace2<-sum(unique(flu$df$ispace)>iispace)
pred.res<-array(NA,dim=c(nahead,nrep,nmodel+1,indspace2)) #nmodel +1 (la y)
pred0<-pred<-array(NA,dim=c(nahead,nrep,nmodel))
lev<-names(table(flu$df$ispace))
indi3<-flu$df$ispace<16000000
# 10 o 53 comarcas
nrep<-500
######epidemic<-1:156 # in the paper, high time-consuming
## Moderate Time-consuming (only 17 rolling iterations (last 17 weeks)
p1<-1
q1<-0
p2<-1
q2<-0
ii1<-c(12,21,20,19)
## High Time-consuming
levcomar<-names(table(flu$df$ispace[indspace]))
ncomar<-length(levcomar)
epidemic<-c(145:164+52,145:164+2*52)
for (i in epidemic) {
for (icomar in 1:ncomar) {
j<-1
flu2<-flu
flu$df$y<-flu$df[,ii1[j]]
ind<-which(flu$df$itime>=(i+a1) & flu$df$itime<=(i+a2)
& flu$df$ispace==levcomar[icomar])
flu2$df<-flu$df[ind,]
flu2$temp<-flu$temp[ind]
flu2$tmp.d1<-flu$tmp.d1[ind]
flu2$temp.d2<-flu$temp.d2[ind]
flu2$temt10<-flu$tempt10[ind]
flu2$ratet13<-flu$ratet13[ind]
flu2$temp.max<-flu$temp.max[ind]
flu2$temp.min<-flu$temp.min[ind]
flu2$vtemp<-flu$rtemp[ind]
flu2$hum<-flu$hum[ind,8:14]
flu2$hum$rangeval<-c(8,14)
flu2$tempt10<-flu$tempt10[ind,]
flu2$tres<-flu$tres[ind,]
flu2$trs.d1<-flu$trs.d1[ind,]
flu2$irra<-flu$irra[ind]
flu2$humt10<-flu$humt10[ind,]
res.lm1<-fregre.lm(ff1,data=flu2)
res.lm2<-fregre.lm(ff2,data=flu2)
res.lm3<-fregre.lm(ff3,data=flu2)
res.lm4<-fregre.lm(ff4,data=flu2)
res.lm5<-fregre.lm(ff5,data=flu2)
res.lm6<-fregre.lm(ff6,data=flu2)
res.lm7<-fregre.lm(ff7,data=flu2)
res.lm8<-fregre.lm(ff8,data=flu2)
# par.cor1<-list("cor.ARMA"=list("index"=c("itime"),"group"="ispace",method="lm","p"=1)r)
# res.fou5<-fregre.igls(ff2,data=flu2,correlation=par.cor1)
res1<-fregre.gls(ff1,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res2<-fregre.gls(ff2,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res3<-fregre.gls(ff3,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res4<-fregre.gls(ff4,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res5<-fregre.gls(ff5,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res6<-fregre.gls(ff6,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res7<-fregre.gls(ff7,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res8<-fregre.gls(ff8,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime),method="ML")
res.gls1<-fregre.gls(ff1,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls2<-fregre.gls(ff2,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls3<-fregre.gls(ff3,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls4<-fregre.gls(ff4,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls5<-fregre.gls(ff5,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls6<-fregre.gls(ff6,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls7<-fregre.gls(ff7,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
res.gls8<-fregre.gls(ff8,data=flu2,correlation=corARMA(p=p2,q=q2,form=~itime),method="ML")
ind2<-which(flu$df$itime==(i+a2+1)& flu$df$ispace==levcomar[icomar])
flu3<-flu
flu3$df$y<-flu$df[,ii0[j]]
flu3$df<-flu$df[ind2,,drop=FALSE]
flu3$temp<-flu$temp[ind2]
flu3$tmp.d1<-flu$tmp.d1[ind2]
flu3$tempt10<-flu$tempt10[ind2]
flu3$tres<-flu$tres[ind2]
flu3$trs.d1<-flu$trs.d1[ind2]
flu3$temp.min<-flu$temp.min[ind2]
flu3$temp.max<-flu$temp.max[ind2]
flu3$ratet13<-flu$ratet13[ind2]
flu3$hum<-flu$hum[ind2,8:14]
flu3$hum$rangeval<-c(8,14)
flu3$irra<-flu$irra[ind2]
# prediction
pred.lm1<-predict.fregre.lm(res.lm1,flu3)
pred.lm2<-predict.fregre.lm(res.lm2,flu3)
pred.lm3<-predict.fregre.lm(res.lm3,flu3)
pred.lm4<-predict.fregre.lm(res.lm4,flu3)
pred.lm5<-predict.fregre.lm(res.lm5,flu3)
pred.lm6<-predict.fregre.lm(res.lm6,flu3)
pred.lm7<-predict.fregre.lm(res.lm7,flu3)
pred.lm8<-predict.fregre.lm(res.lm8,flu3)
pred1<-predict.fregre.gls(res1,flu3,n.ahead=ahead)
pred2<-predict.fregre.gls(res2,flu3,n.ahead=ahead)
pred3<-predict.fregre.gls(res3,flu3,n.ahead=ahead)
pred4<-predict.fregre.gls(res4,flu3,n.ahead=ahead)
pred5<-predict.fregre.gls(res5,flu3,n.ahead=ahead)
pred6<-predict.fregre.gls(res6,flu3,n.ahead=ahead)
pred7<-predict.fregre.gls(res7,flu3,n.ahead=ahead)
pred8<-predict.fregre.gls(res8,flu3,n.ahead=ahead)
yresp<-t(flu3$df[,c("ratet1","ratet2","ratet3","ratet4")])
pred.res[,i,1,icomar]<-rbind(pred.lm1,pred.lm1,pred.lm1,pred.lm1)
pred.res[,i,2,icomar]<-rbind(pred.lm2,pred.lm2,pred.lm2,pred.lm2)
pred.res[,i,3,icomar]<-rbind(pred.lm3,pred.lm3,pred.lm3,pred.lm3)
pred.res[,i,4,icomar]<-rbind(pred.lm4,pred.lm4,pred.lm4,pred.lm4)
pred.res[,i,5,icomar]<-rbind(pred.lm5,pred.lm5,pred.lm5,pred.lm5)
pred.res[,i,6,icomar]<-rbind(pred.lm6,pred.lm6,pred.lm6,pred.lm6)
pred.res[,i,7,icomar]<-rbind(pred.lm7,pred.lm7,pred.lm7,pred.lm7)
pred.res[,i,8,icomar]<-rbind(pred.lm8,pred.lm8,pred.lm8,pred.lm8)
pred.res[,i,9,icomar]<-t(pred1)
pred.res[,i,10,icomar]<-t(pred2 )
pred.res[,i,11,icomar]<-t(pred3)
pred.res[,i,12,icomar]<-t(pred4)
pred.res[,i,13,icomar]<-t(pred5)
pred.res[,i,14,icomar]<-t(pred6 )
pred.res[,i,15,icomar]<-t(pred7 )
pred.res[,i,16,icomar]<-t(pred8 )
pred.res[,i,17,icomar]<-yresp
}
iresp <- 17
pred[,i,1]<-apply((pred.res[,i,iresp,]-pred.res[,i,1,])^2,1,mean)
pred[,i,2]<-apply((pred.res[,i,iresp,]-pred.res[,i,2,])^2,1,mean)
pred[,i,3]<-apply((pred.res[,i,iresp,]-pred.res[,i,3,])^2,1,mean)
pred[,i,4]<-apply((pred.res[,i,iresp,]-pred.res[,i,4,])^2,1,mean)
pred[,i,5]<-apply((pred.res[,i,iresp,]-pred.res[,i,5,])^2,1,mean)
pred[,i,6]<-apply((pred.res[,i,iresp,]-pred.res[,i,6,])^2,1,mean)
pred[,i,7]<-apply((pred.res[,i,iresp,]-pred.res[,i,7,])^2,1,mean)
pred[,i,8]<-apply((pred.res[,i,iresp,]-pred.res[,i,8,])^2,1,mean)
pred[,i,9]<-apply((pred.res[,i,iresp,]-pred.res[,i,9,])^2,1,mean)
pred[,i,10]<-apply((pred.res[,i,iresp,]-pred.res[,i,10,])^2,1,mean)
pred[,i,11]<-apply((pred.res[,i,iresp,]-pred.res[,i,11,])^2,1,mean)
pred[,i,12]<-apply((pred.res[,i,iresp,]-pred.res[,i,12,])^2,1,mean)
pred[,i,13]<-apply((pred.res[,i,iresp,]-pred.res[,i,13,])^2,1,mean)
pred[,i,14]<-apply((pred.res[,i,iresp,]-pred.res[,i,14,])^2,1,mean)
pred[,i,15]<-apply((pred.res[,i,iresp,]-pred.res[,i,15,])^2,1,mean)
pred[,i,16]<-apply((pred.res[,i,iresp,]-pred.res[,i,16,])^2,1,mean)
print( i)
pred[,i,]
}
save.image(file="fluExample.RData")
###############################################################################
load(file="fluExample.RData")
table6<-t(round(apply(pred[,epidemic,1:16],c(1,3),mean,na.rm=T),2))
table6<-cbind(table6[c(1:8),1],table6[c(9:16),1],table6[c(1:8),2],table6[c(9:16),2])
colnames(table6)<-c("FLM (t+1)","FGLS (t+1)","FLM (t+2)","FGLS (t+2)")
rownames(table6)<-c("Rate(w)","Temp(t)","Temp.thres(t)","SR(t)",
"Rate(w),Temp(t)","Rate(w),Temp(t)",
"Rate(w)Temp.thres(t)","Temp(t),SR(t)")
table6
###############################################################################
j<-1
flu2<-flu
flu$df$y<-flu$df[,ii1[j]]
ind<-which(flu$df$itime>=(epidemic[1]+a1) & flu$df$itime<=(epidemic[40]+a2)
& flu$df$ispace %in% levcomar)
length(ind)
flu2$df<-flu$df[ind,]
flu2$temp<-flu$temp[ind]
flu2$tmp.d1<-flu$tmp.d1[ind]
flu2$temp.d2<-flu$temp.d2[ind]
flu2$temt10<-flu$tempt10[ind]
flu2$ratet13<-flu$ratet13[ind]
flu2$temp.max<-flu$temp.max[ind]
flu2$temp.min<-flu$temp.min[ind]
flu2$vtemp<-flu$rtemp[ind]
flu2$hum<-flu$hum[ind,8:14]
flu2$hum$rangeval<-c(8,14)
flu2$tempt10<-flu$tempt10[ind,]
flu2$tres<-flu$tres[ind,]
flu2$trs.d1<-flu$trs.d1[ind,]
flu2$irra<-flu$irra[ind]
flu2$humt10<-flu$humt10[ind,]
res1<-fregre.gls(ff1,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res2<-fregre.gls(ff2,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res3<-fregre.gls(ff3,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res4<-fregre.gls(ff4,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res5<-fregre.gls(ff5,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res6<-fregre.gls(ff6,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res7<-fregre.gls(ff7,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
res8<-fregre.gls(ff8,data=flu2,correlation=corARMA(p=p1,q=q1,form=~itime|ispace),method="ML")
################# INICIO #####################
add_legend <- function(...) {
opar <- par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0),
mar=c(0, 0, 0, 0), new=TRUE)
on.exit(par(opar))
plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
legend(...)
}
pdf(file = "S2_fig.pdf",width = 10.24, height =7.68, pointsize = 10)
par(mfrow=c(2,2))
cols=c("blue","sky blue","red","dark red")
# Rate
covar<- "ratet13"
res <-res1
x<-flu2[[covar]]
a<-fdata(res$beta.l[[covar]],flu2[[covar]]$argval,flu2[[covar]]$rangeval)
b1<-b<- +inprod.fdata(x,a)
qq<-quantile(b,prob=c(.25,.5,.75,1))
ind1<-ifelse(b[,1]<qq[1],TRUE,FALSE)
ind2<-ifelse(b[,1]>=qq[1],TRUE,FALSE)
ind3<-ifelse(b[,1]>=qq[2],TRUE,FALSE)
ind4<-ifelse(b[,1]>=qq[3],TRUE,FALSE)
mn0<-expression(paste( beta[rate]))
ylab <- "Rate(w)"
ylab <- expression(Rate["s,n"](w))
mn2<-expression(paste(symbol("<"), list(Rate["s,n"](w), beta[Rate]),symbol(">")))
plot(c((func.mean(x[ind1,])),(func.mean(x[ind2,]))
,(func.mean(x[ind3,])),(func.mean(x[ind4,])))
,lwd=2,lty=1
,col=cols,
main=mn2,xlab="Range time w: 13 weeks",ylab=ylab)
# Temp
covar<- "temp"
res <-res2
x<-flu2[[covar]]
a<-fdata(res$beta.l[[covar]],flu2[[covar]]$argval,flu2[[covar]]$rangeval)
b1<-b<- +inprod.fdata(x,a)
qq<-quantile(b,prob=c(.25,.5,.75,1))
ind1<-ifelse(b[,1]<qq[1],TRUE,FALSE)
ind2<-ifelse(b[,1]>=qq[1],TRUE,FALSE)
ind3<-ifelse(b[,1]>=qq[2],TRUE,FALSE)
ind4<-ifelse(b[,1]>=qq[3],TRUE,FALSE)
mn0<-expression(paste( beta[rate]))
ylab <- expression(Temp["s,n"](t))
mn2<-expression(paste(symbol("<"), list(Temp["s,n"](t), beta[Temp]),symbol(">")))
plot(c((func.mean(x[ind1,])),(func.mean(x[ind2,]))
,(func.mean(x[ind3,])),(func.mean(x[ind4,])))
,lwd=2,lty=1
,col=cols,
main=mn2,xlab="Range time t: 14 days",ylab=ylab)
# Temp
covar<- "tres"
res <-res3
x<-flu2[[covar]]
a<-fdata(res$beta.l[[covar]],flu2[[covar]]$argval,flu2[[covar]]$rangeval)
b1<-b<- +inprod.fdata(x,a)
qq<-quantile(b,prob=c(.25,.5,.75,1))
ind1<-ifelse(b[,1]<qq[1],TRUE,FALSE)
ind2<-ifelse(b[,1]>=qq[1],TRUE,FALSE)
ind3<-ifelse(b[,1]>=qq[2],TRUE,FALSE)
ind4<-ifelse(b[,1]>=qq[3],TRUE,FALSE)
ylab <- expression(Temp.thres["s,n"](t))
mn2<-expression(paste(symbol("<"), list(Temp.thres["s,n"](t), beta[Temp.thres]),symbol(">")))
plot(c((func.mean(x[ind1,])),(func.mean(x[ind2,]))
,(func.mean(x[ind3,])),(func.mean(x[ind4,])))
,lwd=2,lty=1
,col=cols,
# main=mn2,xlab="t: previous 14 days",ylab=ylab)
main=mn2,xlab="t in days",ylab=ylab)
# Solar Radiation
covar<- "irra"
res <-res4
x<-flu2[[covar]]
a<-fdata(res$beta.l[[covar]],flu2[[covar]]$argval,flu2[[covar]]$rangeval)
b1<-b<- +inprod.fdata(x,a)
qq<-quantile(b,prob=c(.25,.5,.75,1))
ind1<-ifelse(b[,1]<qq[1],TRUE,FALSE)
ind2<-ifelse(b[,1]>=qq[1],TRUE,FALSE)
ind3<-ifelse(b[,1]>=qq[2],TRUE,FALSE)
ind4<-ifelse(b[,1]>=qq[3],TRUE,FALSE)
ylab <- expression(SR["s,n"](t))
mn2<-expression(paste(symbol("<"), list(SR["s,n"](t), beta[SR]),symbol(">")))
plot(c((func.mean(x[ind1,])),(func.mean(x[ind2,]))
,(func.mean(x[ind3,])),(func.mean(x[ind4,])))
,lwd=2,lty=1
,col=cols,
main=mn2,xlab="t (in days)",ylab=ylab)
leg <- c(expression(v[q[1]]), expression(v[q[2]]), expression(v[q[3]]) ,expression(v[q[4]]))
#legend(8,2500,inset=c(-0.2,0),col=cols,legend=leg,cex=.8,bty="n",lty=1,lwd=1)
add_legend("bottom", legend=leg,lty=1, #pch=20,
col=cols, lwd=2,
horiz=TRUE, bty='n',
cex=0.95)
dev.off()
################# fin Gráfico #####################
################################################################################
################################################################################
# Figure 1 Series Temporales
flu4<-aggregate(flu2$df,list(flu2$df[,"itime"]),mean)
#cambiar: utilizar un ggplot o algo así más chulo (ver links/linkedIN)
png(file = "S1_fig.png")
par(mfrow=c(4,1))
FluRate<-ts(exp(flu4[,"ratet0"]),2001,2011,53)
Temperature<-ts(flu4[,"tmed"] ,2001,2011,53)
SolarRadiation<-ts(flu4[,"imed"] ,2001,2011,53)
RelativeHumidity<-ts(flu4[,"hmed"] ,2001,2011,53)
TemperatureThreshold<-ts(flu4[,"thold"] ,2001,2011,53)
aa<-cbind(FluRate,Temperature,SolarRadiation,RelativeHumidity)
plot(aa,main="",xlab="Years")
dev.off()
################################################################################
################################################################################
################################################################################
################################################################################
# Pearson and distance correlation
################################################################################
################################################################################
#Upper diagonal matrix of Table 1: Pearson correlation
# no se muestran
A<-cor(flu2$df[,c(19,12,22,13,23,26,29)])
################################################################################
#Lower diagonal matrix of Table 1: Distance correlation
################################################################################
# Distance computation for Multivariate-Multivariate data
names(flu)
A1<-metric.lp(flu2$ratet13)
A2<-metric.lp(flu2$temp)
A3<-metric.lp(flu2$tres)
A4<-metric.lp(flu2$irra)
A5<-metric.lp(flu2$hum)
A6<-metric.dist(data.matrix(flu2$df$ratet0))
A7<-metric.dist(data.matrix(flu2$df$ratet1))
A8<-metric.dist(data.matrix(flu2$df$ratet2))
A9<-metric.dist(data.matrix(flu2$df$tmed))
A10<-metric.dist(data.matrix(flu2$df$thold))
A11<-metric.dist(data.matrix(flu2$df$imed))
A12<-metric.dist(data.matrix(flu2$df$hmed))
# Distance correlation
A<-cor(flu2$df[,c(21,12,22,13,23,26,29)])
names(flu2$df[,c(21,12,22,13,23,26,29)])
A[2,1]<-dcor.dist(A8,A7)
A[3,1]<-dcor.dist(A8,A6)
A[4,1]<-dcor.dist(A8,A9)
A[5,1]<-dcor.dist(A8,A10)
A[6,1]<-dcor.dist(A8,A11)
A[7,1]<-dcor.dist(A8,A12)
A[3,2]<-dcor.dist(A7,A6)
A[4,2]<-dcor.dist(A7,A9)
A[5,2]<-dcor.dist(A7,A10)
A[6,2]<-dcor.dist(A7,A11)
A[7,2]<-dcor.dist(A7,A12)
A[4,3]<-dcor.dist(A6,A9)
A[5,3]<-dcor.dist(A6,A10)
A[6,3]<-dcor.dist(A6,A11)
A[7,3]<-dcor.dist(A6,A12)
A[5,4]<-dcor.dist(A9,A10)
A[6,4]<-dcor.dist(A9,A11)
A[7,4]<-dcor.dist(A9,A12)
A[6,5]<-dcor.dist(A10,A11)
A[7,5]<-dcor.dist(A10,A12)
A[7,6]<-dcor.dist(A11,A12)
round(A,2)
################################################################################
################################################################################
B<-matrix(NA,5,12)
# Distance computation for Functional-Multivariate data
B[1,1]<-dcor.dist(A8,A1)
B[2,1]<-dcor.dist(A8,A2)
B[3,1]<-dcor.dist(A8,A3)
B[4,1]<-dcor.dist(A8,A4)
B[5,1]<-dcor.dist(A8,A5)
B[1,2]<-dcor.dist(A7,A1)
B[2,2]<-dcor.dist(A7,A2)
B[3,2]<-dcor.dist(A7,A3)
B[4,2]<-dcor.dist(A7,A4)
B[5,2]<-dcor.dist(A7,A5)
B[1,3]<-dcor.dist(A6,A1)
B[2,3]<-dcor.dist(A6,A2)
B[3,3]<-dcor.dist(A6,A3)
B[4,3]<-dcor.dist(A6,A4)
B[5,3]<-dcor.dist(A6,A5)
B[1,4]<-dcor.dist(A9,A1)
B[2,4]<-dcor.dist(A9,A2)
B[3,4]<-dcor.dist(A9,A3)
B[4,4]<-dcor.dist(A9,A4)
B[5,4]<-dcor.dist(A9,A5)
B[1,5]<-dcor.dist(A10,A1)
B[2,5]<-dcor.dist(A10,A2)
B[3,5]<-dcor.dist(A10,A3)
B[4,5]<-dcor.dist(A10,A4)
B[5,5]<-dcor.dist(A10,A5)
B[1,6]<-dcor.dist(A11,A1)
B[2,6]<-dcor.dist(A11,A2)
B[3,6]<-dcor.dist(A11,A3)
B[4,6]<-dcor.dist(A11,A4)
B[5,6]<-dcor.dist(A11,A5)
B[1,7]<-dcor.dist(A12,A1)
B[2,7]<-dcor.dist(A12,A2)
B[3,7]<-dcor.dist(A12,A3)
B[4,7]<-dcor.dist(A12,A4)
B[5,7]<-dcor.dist(A12,A5)
# Distance computation for Functional-Functional data
B[1,8]<- B[2,9]<- B[3,10]<- B[4,11]<- B[5,12]<-1
B[2,8]<-dcor.dist(A1,A2)
B[3,8]<-dcor.dist(A1,A3)
B[4,8]<-dcor.dist(A1,A4)
B[5,8]<-dcor.dist(A1,A5)
B[3,9]<-dcor.dist(A2,A3)
B[4,9]<-dcor.dist(A2,A4)
B[5,9]<-dcor.dist(A2,A5)
B[4,10]<-dcor.dist(A3,A4)
B[5,10]<-dcor.dist(A3,A5)
B[5,11]<-dcor.dist(A4,A5)
B[1,8:12]<-B[,8]
B[2,9:12]<-B[2:5,9]
B[3,10:12]<-B[3:5,10]
B[4,11:12]<-B[4:5,11]
rnames<-c("raten+2","raten+1","rate","Temp","Thold","Irra","Hum")
fnames<-c("rate(w)","Temp(t)","Thold(t)","Irra(t)","Hum(t)")
colnames(B)<-c(rnames,fnames)
rownames(B)<-fnames
# TABLE 5
round(B[,c(8:12,2,1)],2)
save.image(file="fluExample.RData")
Try the fda.usc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.