Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
605 lines (449 sloc) 20.3 KB
title author date output
Timeseries plots and trend analysis
Russ Poldrack
November 2, 2014
html_document

Timeseries plots for MyConnectome dataset

Code available at: https://github.com/poldrack/myconnectome/blob/master/myconnectome/timeseries/Make_timeseries_plots.Rmd

This code generates timeseries plots for each variable, and performs a test for linear and polynomial trends.

library(forecast)
library(knitr)
basedir=Sys.getenv('MYCONNECTOME_DIR')


datamat=load_behav_data() 
varname='behav'
sample_spacing=1

xvars=c('panas.positive','panas.negative','panas.fatigue','afterscan.Anxietyduringscan','afterscan.diastolic','afterscan.pulse','afterscan.systolic','morning.Sleepquality','morning.Soreness','prevevening.Alcohol','prevevening.Guthealth','prevevening.Psoriasisseverity','prevevening.Stress', 'prevevening.Timespentoutdoors','temp.mean',"email.LIWCcdi","email.LIWCnegemo","email.LIWCposemo",'zeo.zq')

datamat=subset(datamat,select=c(xvars,'date'))

save_latex=TRUE


Behavioral variables

FDR-corrected p-values for trends

pvals=c()
# first get fdr p-vals on trends
for (varnum in 1:(dim(datamat)[2]-1)) {
  alldays = seq(min(datamat$date), max(datamat$date), by='1 day')
  alldays_weekly=seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=na.approx(zoo(x[,varnum]),maxgap=7)
	
	x_interp=x_ts
	x_interp=na.approx(zoo(x[,varnum]))
	
	s_poly=(as.numeric(alldays) - mean(as.numeric(alldays)))**2
	s_poly=(s_poly - mean(s_poly))
	s_poly=s_poly/max(s_poly)
	
	s_lin=as.numeric(alldays) - mean(as.numeric(alldays))
	s_lin=s_lin/max(s_lin)

	# take weekly mean for auto.arima analysis, due to NAs
	a=auto.arima(x[,varnum],xreg=cbind(s_poly,s_lin),allowdrift=FALSE)
	p=cwp(a)
	poly.loc=which(names(data.frame(p)) == "s_poly")
	lin.loc=which(names(data.frame(p)) == "s_lin")
	pvals=rbind(pvals,c(mean(x_ts,na.rm=TRUE), var(a$residuals,na.rm=TRUE),p[4,lin.loc],p[4,poly.loc]))
}
p_fdr=p.adjust(pvals[,3:4],method='BH')
p_fdr=matrix(p_fdr,length(p_fdr)/2,2)

pvals=as.data.frame(pvals)
pvals[,3:4]=p_fdr
names(pvals)=c('mean','var(resid)','Linear trend (p-val)','Poly trend (p-val)')
row.names(pvals)=names(datamat)[1:19]

write.table(pvals,file=sprintf('%s/timeseries/behav_timeseries_stats.txt',basedir),row.names=FALSE,col.names=FALSE)

kable(pvals)
if (save_latex) {
  k=kable(pvals,format='latex',digits=3)
  write(k,file=sprintf('%s/timeseries/tables/behav_timeseries_stats.tex',basedir))
}
for (varnum in 1:(dim(datamat)[2]-1)) {
  alldays = seq(min(datamat$date), max(datamat$date), by='1 day')
	alldays_weekly=seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=na.approx(zoo(x[,varnum]),maxgap=7,na.rm=FALSE)
	
	x_interp=x_ts
	x_interp=na.approx(zoo(x[,varnum]),na.rm=FALSE)
	x_interp_nonan=na.approx(zoo(x[,varnum]),na.rm=TRUE)

	lo <- loess(x_interp ~ seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),span=0.5)
  
	plot(x_ts,xaxt='n',xlab='',ylab=varname,col='black')
	points(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),x[,varnum],pch=18)
	lines(seq(from=start(x_interp_nonan)[1],to=end(x_interp_nonan)[1],by=sample_spacing),lo$fitted,col='blue',lwd=3)
	axis(2)
	axis(1,labels=FALSE,at=seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing),xlab=FALSE)
	offset=0.1*(max(x_ts,na.rm=TRUE)-min(x_ts,na.rm=TRUE))
	text(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing), par("usr")[3] - offset, labels = labels, srt = 45,adj=0.9, xpd = TRUE)
	title(sprintf('%s - lin p = %0.3f, poly p=%0.3f',names(datamat)[varnum],p_fdr[varnum,1],p_fdr[varnum,2]))	
}

Network measures

FDR-corrected p-values for trends


datamat=load_network_data() 
meancorsim=read.table(sprintf("%s/rsfmri/corrsim_mean.txt",basedir))$V1
datamat=cbind(datamat[,1:2],meancorsim,datamat$date)
names(datamat)=c('modularity_weighted','eff_weighted','meancorrsim','date')

varname='netdat'
sample_spacing=1

pvals=c()
# first get fdr p-vals on trends
for (varnum in 1:(dim(datamat)[2]-1)) {
  alldays = seq(min(datamat$date), max(datamat$date), by='1 day')
	alldays_weekly=seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=na.approx(zoo(x[,varnum]),maxgap=7)
	
	x_interp=x_ts
	x_interp=na.approx(zoo(x[,varnum]))
	
	s_poly=(seq(length(x[,varnum])) - mean(seq(length(x[,varnum]))))**2
	s_poly=(s_poly - mean(s_poly))
	s_poly=s_poly/max(s_poly)
	s_lin=seq(length(x[,varnum])) - mean(seq(length(x[,varnum])))
	s_lin=s_lin/max(s_lin)
	
	# take weekly mean for auto.arima analysis, due to NAs
	a=auto.arima(x[,varnum],xreg=cbind(s_poly,s_lin),allowdrift=FALSE)
	p=cwp(a)
	poly.loc=which(names(data.frame(p)) == "s_poly")
	lin.loc=which(names(data.frame(p)) == "s_lin")
	pvals=rbind(pvals,c(mean(x_ts,na.rm=TRUE), var(a$residuals,na.rm=TRUE),p[4,lin.loc],p[4,poly.loc]))
}
p_fdr=p.adjust(pvals[,3:4],method='BH')
p_fdr=matrix(p_fdr,length(p_fdr)/2,2)

pvals=as.data.frame(pvals)
pvals[,3:4]=p_fdr
names(pvals)=c('mean','var(resid)','Linear trend (p-val)','Poly trend (p-val)')
row.names(pvals)=names(datamat)[1:3]

kable(pvals)
if (save_latex) {
  k=kable(pvals,format='latex',digits=4)
  write(k,file=sprintf('%s/timeseries/tables/netdat_timeseries_stats.tex',basedir))
}


for (varnum in 1:(dim(datamat)[2]-1)) {
	alldays = seq(min(datamat$date), max(datamat$date), by='1 day')
	alldays_weekly=seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=na.approx(zoo(x[,varnum]),maxgap=7)
	
	x_interp=x_ts
	x_interp=na.approx(zoo(x[,varnum]))
	lo <- loess(x_interp ~ seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),span=0.5)
	plot(x_ts,xaxt='n',xlab='',ylab=varname,col='black')
	points(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),x[,varnum],pch=18)
	lines(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),lo$fitted,col='blue',lwd=3)
	axis(2)
	axis(1,labels=FALSE,at=seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing),xlab=FALSE)
	offset=0.1*(max(x_ts,na.rm=TRUE)-min(x_ts,na.rm=TRUE))
	text(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing), par("usr")[3] - offset, labels = labels, srt = 45,adj=0.9, xpd = TRUE)
	title(sprintf('%s - lin p = %0.3f, poly p=%0.3f',names(datamat)[varnum],p_fdr[varnum,1],p_fdr[varnum,2]))
}

Within-network connectivity

FDR-corrected p-values for trends

datamat=load_fmri_data() 
varname='wincorr'
sample_spacing=1

pvals=c()
# first get fdr p-vals on trends
for (varnum in 1:(dim(datamat)[2]-1)) {
  alldays = seq(min(datamat$date), max(datamat$date), by='1 day')
	alldays_weekly=seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=na.approx(zoo(x[,varnum]),maxgap=7)
	
	x_interp=x_ts
	x_interp=na.approx(zoo(x[,varnum]))
	
	s_poly=(seq(length(x[,varnum])) - mean(seq(length(x[,varnum]))))**2
	s_poly=(s_poly - mean(s_poly))
	s_poly=s_poly/max(s_poly)
	s_lin=seq(length(x[,varnum])) - mean(seq(length(x[,varnum])))
	s_lin=s_lin/max(s_lin)
	
	# take weekly mean for auto.arima analysis, due to NAs
	a=auto.arima(x[,varnum],xreg=cbind(s_poly,s_lin),allowdrift=FALSE)
	p=cwp(a)
	poly.loc=which(names(data.frame(p)) == "s_poly")
	lin.loc=which(names(data.frame(p)) == "s_lin")
	pvals=rbind(pvals,c(mean(x_ts,na.rm=TRUE), var(a$residuals,na.rm=TRUE),p[4,lin.loc],p[4,poly.loc]))
}
p_fdr=p.adjust(pvals[,3:4],method='BH')
p_fdr=matrix(p_fdr,length(p_fdr)/2,2)

pvals=as.data.frame(pvals)
pvals[,3:4]=p_fdr
names(pvals)=c('mean','var(resid)','Linear trend (p-val)','Poly trend (p-val)')
row.names(pvals)=names(datamat)[1:12]
write.table(pvals,file=sprintf('%s/timeseries/wincorr_timeseries_stats.txt',basedir),row.name=FALSE,col.names=FALSE)

kable(pvals)
if (save_latex) {
  k=kable(pvals,format='latex',digits=4)
  write(k,file=sprintf('%s/timeseries/tables/wincorr_timeseries_stats.tex',basedir))
}


for (varnum in 1:(dim(datamat)[2]-1)) {
	alldays = seq(min(datamat$date), max(datamat$date), by='1 day')
	alldays_weekly=seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=na.approx(zoo(x[,varnum]),maxgap=7)
	
	x_interp=x_ts
	x_interp=na.approx(zoo(x[,varnum]))

	lo <- loess(x_interp ~ seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),span=0.5)
	plot(x_ts,xaxt='n',xlab='',ylab=varname,col='black')
	points(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),x[,varnum],pch=18)
	lines(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),lo$fitted,col='blue',lwd=3)
	axis(2)
	axis(1,labels=FALSE,at=seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing),xlab=FALSE)
	offset=0.1*(max(x_ts,na.rm=TRUE)-min(x_ts,na.rm=TRUE))
	text(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing), par("usr")[3] - offset, labels = labels, srt = 45,adj=0.9, xpd = TRUE)
	title(sprintf('%s - lin p = %0.3f, poly p=%0.3f',names(datamat)[varnum],p_fdr[varnum,1],p_fdr[varnum,2]))
}

Between-module correlations

FDR-corrected p-values for trends

datamat=load_fmri_data(type='bwcorr') 
varname='bwcorr'
sample_spacing=1

pvals=c()
# first get fdr p-vals on trends
for (varnum in 1:(dim(datamat)[2]-1)) {
  alldays = seq(min(datamat$date), max(datamat$date), by='1 day')
  alldays_weekly=seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=na.approx(zoo(x[,varnum]),maxgap=7)
	
	x_interp=x_ts
	x_interp=na.approx(zoo(x[,varnum]))
	
	s_poly=(seq(length(x[,varnum])) - mean(seq(length(x[,varnum]))))**2
	s_poly=(s_poly - mean(s_poly))
	s_poly=s_poly/max(s_poly)
	s_lin=seq(length(x[,varnum])) - mean(seq(length(x[,varnum])))
	s_lin=s_lin/max(s_lin)
	
	# take weekly mean for auto.arima analysis, due to NAs
	a=auto.arima(x[,varnum],xreg=cbind(s_poly,s_lin),allowdrift=FALSE)
	p=cwp(a)
	poly.loc=which(names(data.frame(p)) == "s_poly")
	lin.loc=which(names(data.frame(p)) == "s_lin")
	pvals=rbind(pvals,c(mean(x_ts,na.rm=TRUE), var(a$residuals,na.rm=TRUE),p[4,lin.loc],p[4,poly.loc]))
}
p_fdr=p.adjust(pvals[,3:4],method='BH')
p_fdr=matrix(p_fdr,length(p_fdr)/2,2)

pvals=as.data.frame(pvals)
pvals[,3:4]=p_fdr
names(pvals)=c('mean','var(resid)','Linear trend (p-val)','Poly trend (p-val)')
row.names(pvals)=names(datamat)[1:66]
write.table(pvals,file=sprintf('%s/timeseries/bwcorr_timeseries_stats.txt',basedir),row.names=FALSE,col.names=FALSE)
kable(pvals)
if (save_latex) {
  k=kable(pvals,format='latex',digits=4)
  write(k,file=sprintf('%s/timeseries/tables/bwcorr_timeseries_stats.tex',basedir))
}


Skipping plots for between-module correlations, too may variables

RNA-seq

FDR-corrected p-values for trends

datamat=load_rnaseq_data()
mod_des=read.table(sprintf('%s/rna-seq/WGCNA/module_descriptions.txt',basedir),sep='\t')
data_names=c()
for (i in 1:(dim(datamat)[2]-1)) {
  	data_names=rbind(data_names,sprintf('ME%d:%s',i,mod_des$V2[i]))
		}

data_names=rbind(data_names,'date')
names(datamat)=data_names

varname='Eigengene expression'
sample_spacing=7

# first get pvals
pvals=c()

for (varnum in 1:(dim(datamat)[2]-1)) {
	alldays = seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=zoo(x[,varnum])
	
	x_interp=x_ts
	x_interp=na.approx(x_interp)
	
	
	s_poly=(as.numeric(alldays) - mean(as.numeric(alldays)))**2
	s_poly=(s_poly - mean(s_poly))
	s_poly=s_poly/max(s_poly)
	
	s_lin=as.numeric(alldays) - mean(as.numeric(alldays))
	s_lin=s_lin/max(s_lin)
	
	a=auto.arima(x[,varnum],xreg=cbind(s_poly,s_lin),allowdrift=FALSE)
	p=cwp(a)
	poly.loc=which(names(data.frame(p)) == "s_poly")
	lin.loc=which(names(data.frame(p)) == "s_lin")
	pvals=rbind(pvals,c(mean(x_ts,na.rm=TRUE), var(a$residuals,na.rm=TRUE),p[4,lin.loc],p[4,poly.loc]))

}

p_fdr=p.adjust(pvals[,3:4],method='BH')
p_fdr=matrix(p_fdr,length(p_fdr)/2,2)

pvals=as.data.frame(pvals)
pvals[,3:4]=p_fdr
names(pvals)=c('mean','var(resid)','Linear trend (p-val)','Poly trend (p-val)')


row.names(pvals)=names(datamat)[1:dim(pvals)[1]]


kable(pvals)
if (save_latex) {
  k=kable(pvals,format='latex',digits=3)
  write(k,file=sprintf('%s/timeseries/tables/wgcna_timeseries_stats.tex',basedir))
}



for (varnum in 1:(dim(datamat)[2]-1)) {
	alldays = seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=zoo(x[,varnum])
	
	x_interp=x_ts
	x_interp=na.approx(x_interp)
	

	lo <- loess(x_interp ~ seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),span=0.5)
	pred=predict(lo,se=TRUE)
	plot(x_ts,xaxt='n',xlab='',ylab=varname,col='black')
	lines(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),lo$fitted,col='blue')
	points(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),x[,varnum],pch=18)
	axis(2)
	axis(1,labels=FALSE,at=seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing),xlab=FALSE)
	offset=0.1*(max(x_ts,na.rm=TRUE)-min(x_ts,na.rm=TRUE))
	text(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing), par("usr")[3] - offset, labels = labels, srt = 45,adj=0.9, xpd = TRUE)
	title(sprintf('%s - lin p = %0.3f, poly p=%0.3f',names(datamat)[varnum],p_fdr[varnum,1],p_fdr[varnum,2]))

  }

Metabolite clusters

FDR-corrected p-values for trends

datamat=load_metab_data() 
varname='Metabolite eigenconcentration'
sample_spacing=7

data_names=c()
for (i in 1:(dim(datamat)[2]-1)) {
    data_names=rbind(data_names,sprintf('C%d:%s',i,names(datamat)[i]))
		}
data_names=rbind(data_names,'date')
names(datamat)=data_names



pvals=c()

# first get pvals

for (varnum in 1:(dim(datamat)[2]-1)) {
	if (sample_spacing==7) {
		alldays = seq(min(datamat$date), max(datamat$date), by='1 week')
	} else {
		alldays = seq(min(datamat$date), max(datamat$date), by='1 day')
	}

	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=zoo(x[,varnum])
	
	x_interp=x_ts
	x_interp=na.approx(x_interp,,maxgap=7)
	
	
	s_poly=(as.numeric(alldays) - mean(as.numeric(alldays)))**2
	s_poly=(s_poly - mean(s_poly))
	s_poly=s_poly/max(s_poly)
	
	s_lin=as.numeric(alldays) - mean(as.numeric(alldays))
	s_lin=s_lin/max(s_lin)
	
	a=auto.arima(x[,varnum],allowdrift=FALSE, xreg=cbind(s_poly,s_lin))
	p=cwp(a)
	poly.loc=which(names(data.frame(p)) == "s_poly")
	lin.loc=which(names(data.frame(p)) == "s_lin")
	pvals=rbind(pvals,c(mean(x_ts,na.rm=TRUE), var(a$residuals,na.rm=TRUE),p[4,lin.loc],p[4,poly.loc]))
}

p_fdr=p.adjust(pvals[,3:4],method='BH')
p_fdr=matrix(p_fdr,length(p_fdr)/2,2)

pvals=as.data.frame(pvals)
pvals[,3:4]=p_fdr
names(pvals)=c('mean','var(resid)','Linear trend (p-val)','Poly trend (p-val)')


row.names(pvals)=names(datamat)[1:dim(pvals)[1]]

kable(pvals)
if (save_latex) {
  k=kable(pvals,format='latex',digits=3)
  write(k,file=sprintf('%s/timeseries/tables/metab_timeseries_stats.tex',basedir))
}


for (varnum in 1:(dim(datamat)[2]-1)) {
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=zoo(x[,varnum])
	
	x_interp=x_ts
	x_interp=na.approx(x_interp,maxgap=7)
	

	lo <- loess(x_interp ~ seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),span=0.5)
	pred=predict(lo,se=TRUE)
	plot(x_ts,xaxt='n',xlab='',ylab=varname,col='black')
	lines(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),lo$fitted,col='blue')
	points(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),x[,varnum],pch=18)
	axis(2)
	axis(1,labels=FALSE,at=seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing),xlab=FALSE)
	offset=0.1*(max(x_ts,na.rm=TRUE)-min(x_ts,na.rm=TRUE))
	text(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing), par("usr")[3] - offset, labels = labels, srt = 45,adj=0.9, xpd = TRUE)
	title(sprintf('%s - lin p = %0.3f, poly p=%0.3f',names(datamat)[varnum],p_fdr[varnum,1],p_fdr[varnum,2]))
}


ImmPort immune system pathways

FDR-corrected p-values for trends

datamat=load_ImmPort_data()

varname='Eigengene expression'
sample_spacing=7
varnum=1
pvals=c()
# first get pvals

for (varnum in 1:(dim(datamat)[2]-1)) {
	alldays = seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=zoo(x[,varnum])
	
	x_interp=x_ts
	x_interp=na.approx(x_interp)
	
	
	s_poly=(as.numeric(alldays) - mean(as.numeric(alldays)))**2
	s_poly=(s_poly - mean(s_poly))
	s_poly=s_poly/max(s_poly)
	
	s_lin=as.numeric(alldays) - mean(as.numeric(alldays))
	s_lin=s_lin/max(s_lin)
	
	a=auto.arima(x[,varnum],xreg=cbind(s_poly,s_lin),allowdrift=FALSE)
	p=cwp(a)
	poly.loc=which(names(data.frame(p)) == "s_poly")
	lin.loc=which(names(data.frame(p)) == "s_lin")
	pvals=rbind(pvals,c(mean(x_ts,na.rm=TRUE), var(a$residuals,na.rm=TRUE),p[4,lin.loc],p[4,poly.loc]))

}

p_fdr=p.adjust(pvals[,3:4],method='BH')
p_fdr=matrix(p_fdr,length(p_fdr)/2,2)

pvals=as.data.frame(pvals)
pvals[,3:4]=p_fdr
names(pvals)=c('mean','var(resid)','Linear trend (p-val)','Poly trend (p-val)')

row.names(pvals)=names(datamat)[1:15]

kable(pvals)
if (save_latex) {
  k=kable(pvals,format='latex',digits=3)
  write(k,file=sprintf('%s/timeseries/tables/immport_timeseries_stats.tex',basedir))
}



for (varnum in 1:(dim(datamat)[2]-1)) {
	alldays = seq(min(datamat$date), max(datamat$date), by='1 week')
	
	x=get_x_ts(datamat,alldays,scale=FALSE)
	spacing=60
	labels=seq(from=as.Date(datamat$date[1]),to=as.Date(datamat$date[length(datamat$date)]),by=spacing)
	x_ts=zoo(x[,varnum])
	
	x_interp=x_ts
	x_interp=na.approx(x_interp)
	

	lo <- loess(x_interp ~ seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),span=0.5)
	pred=predict(lo,se=TRUE)
	plot(x_ts,xaxt='n',xlab='',ylab=varname,col='black')
	lines(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),lo$fitted,col='blue')
	points(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=sample_spacing),x[,varnum],pch=18)
	axis(2)
	axis(1,labels=FALSE,at=seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing),xlab=FALSE)
	offset=0.1*(max(x_ts,na.rm=TRUE)-min(x_ts,na.rm=TRUE))
	text(seq(from=start(x_interp)[1],to=end(x_interp)[1],by=spacing), par("usr")[3] - offset, labels = labels, srt = 45,adj=0.9, xpd = TRUE)
	title(sprintf('%s - lin p = %0.3f, poly p=%0.3f',names(datamat)[varnum],p_fdr[varnum,1],p_fdr[varnum,2]))
}