Nothing
inst/script/PLOS_influenza2fig.R
In fda.usc: Functional Data Analysis and Utilities for Statistical Computing
# rm(list=ls())
completa.fdata<-function(y){
n<-nrow(y)
p<-ncol(y)
x<-y$argvals
ynew<-y
aa<-na.omit(y)
ii<-which(is.na.fdata(y))
for (i in ii){
app=approx(x,y$data[i,],n=p)
x=app$x
ynew$data[i,]<-app$y
}
ynew
}
library(fda.usc)
################################################################################
load("influenza.rda")
nt<-53*52*2
nrep<-156
a1<-209+ 52
a2<-a1+nrep-10
nmodel<-12
flu2<-flu
# selection response variable (n+1, to n+5)
ii1<-c(12,21,20,19)
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<-1
a2<-a1+128
iispace<-3500
#iispace<-0
indspace<-flu$df$ispace>iispace
indspace2<-sum(unique(flu$df$ispace)>iispace)
lev<-names(table(flu$df$ispace))
indi3<-flu$df$ispace<16000000
######epidemic<-1:156 # in the paper, high time-consuming
## Moderate Time-consuming (only 17 rolling iterations (last 17 weeks)
p1<-2
q1<-0
p2<-1
q2<-0
ii1<-c(12,21,20,19)
## High Time-consuming
levcomar<-names(table(flu$df$ispace[indspace]))
#levcomar <- setdiff(levcomar,c(levcomar[23],"3207","3209", "1506"))
ncomar<-length(levcomar)
epidemic<-c(145:164+52,145:164+2*52)
flu$irra$names<-list("main"= "RS" ,"xlab"= "Day","ylab"= "W/m2")
flu$irra<-completa.fdata(flu$irra)
iresp <- 25 # donde se guarda la ytest
n.ahead<-4
# repetir para cor.ARMA(1,0) o ARMA(2,0)
p1 <- 1
p2 <- 1
par.AR<-list("cor.AR"=list("order.max"=4))
par.ARMA1<-list("cor.ARMA"=list("p"=1))
par.ARMA2<-list("cor.ARMA"=list("p"=2))
################################################################################
# Figuras
################################################################################
nt<-53*52*2
#nrep<-156
a1<-209+ 52
a2<-a1+nrep-10
nmodel<-12
## High Time-consuming
levcomar<-names(table(flu$df$ispace[indspace]))
#levcomar <- setdiff(levcomar,c(levcomar[23],"3207","3209", "1506"))
ncomar<-length(levcomar)
epidemic<-c(145:164+52,145:164+2*52)
icomar<-2
nam.var<-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)")
nam.mod<-c("FLM","GLS","iGLS")
ij<-1 #horizonete
j<-1
flu2<-flu
flu$df$y<-flu$df[,ii1[j]]
ind<-which(flu$df$itime>=(epidemic[1]+a1) & flu$df$itime<=616
& 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")
#a plicar logartibmo en la tasa !!!!!
################# INICIO #####################
add_legend <- function(...) {
opar <- par(fig=c(0, 1, 0, 1), oma=c(1, 2,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)
# se guarda manualmente como pdf
windows(1000,700)
par(mfrow=c(1,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,cex=.9,
main=mn2,xlab="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,cex=.9,
main=mn2,xlab=" Time t: 14 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 Grafico #####################
################################################################################
################################################################################
# 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 6 PLOS paper
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.
# rm(list=ls())
completa.fdata<-function(y){
n<-nrow(y)
p<-ncol(y)
x<-y$argvals
ynew<-y
aa<-na.omit(y)
ii<-which(is.na.fdata(y))
for (i in ii){
app=approx(x,y$data[i,],n=p)
x=app$x
ynew$data[i,]<-app$y
}
ynew
}
library(fda.usc)
################################################################################
load("influenza.rda")
nt<-53*52*2
nrep<-156
a1<-209+ 52
a2<-a1+nrep-10
nmodel<-12
flu2<-flu
# selection response variable (n+1, to n+5)
ii1<-c(12,21,20,19)
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<-1
a2<-a1+128
iispace<-3500
#iispace<-0
indspace<-flu$df$ispace>iispace
indspace2<-sum(unique(flu$df$ispace)>iispace)
lev<-names(table(flu$df$ispace))
indi3<-flu$df$ispace<16000000
######epidemic<-1:156 # in the paper, high time-consuming
## Moderate Time-consuming (only 17 rolling iterations (last 17 weeks)
p1<-2
q1<-0
p2<-1
q2<-0
ii1<-c(12,21,20,19)
## High Time-consuming
levcomar<-names(table(flu$df$ispace[indspace]))
#levcomar <- setdiff(levcomar,c(levcomar[23],"3207","3209", "1506"))
ncomar<-length(levcomar)
epidemic<-c(145:164+52,145:164+2*52)
flu$irra$names<-list("main"= "RS" ,"xlab"= "Day","ylab"= "W/m2")
flu$irra<-completa.fdata(flu$irra)
iresp <- 25 # donde se guarda la ytest
n.ahead<-4
# repetir para cor.ARMA(1,0) o ARMA(2,0)
p1 <- 1
p2 <- 1
par.AR<-list("cor.AR"=list("order.max"=4))
par.ARMA1<-list("cor.ARMA"=list("p"=1))
par.ARMA2<-list("cor.ARMA"=list("p"=2))
################################################################################
# Figuras
################################################################################
nt<-53*52*2
#nrep<-156
a1<-209+ 52
a2<-a1+nrep-10
nmodel<-12
## High Time-consuming
levcomar<-names(table(flu$df$ispace[indspace]))
#levcomar <- setdiff(levcomar,c(levcomar[23],"3207","3209", "1506"))
ncomar<-length(levcomar)
epidemic<-c(145:164+52,145:164+2*52)
icomar<-2
nam.var<-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)")
nam.mod<-c("FLM","GLS","iGLS")
ij<-1 #horizonete
j<-1
flu2<-flu
flu$df$y<-flu$df[,ii1[j]]
ind<-which(flu$df$itime>=(epidemic[1]+a1) & flu$df$itime<=616
& 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")
#a plicar logartibmo en la tasa !!!!!
################# INICIO #####################
add_legend <- function(...) {
opar <- par(fig=c(0, 1, 0, 1), oma=c(1, 2,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)
# se guarda manualmente como pdf
windows(1000,700)
par(mfrow=c(1,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,cex=.9,
main=mn2,xlab="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,cex=.9,
main=mn2,xlab=" Time t: 14 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 Grafico #####################
################################################################################
################################################################################
# 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 6 PLOS paper
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.