R/apc_plot_data.R

In apc: Age-Period-Cohort Analysis

#######################################################
#	apc package
#	Bent Nielsen, 15 Aug 2018, version 1.3.2
#	Bent Nielsen, 2 May 2017, version 1.3.1
#	Bent Nielsen, 17 Nov 2016, version 1.3
#	Bent Nielsen, 26 Apr 2015, version 1.1.1
#	functions to plot data
#######################################################
#	Copyright 2014-2017 Bent Nielsen
#	Nuffield College, OX1 1NF, UK
#	bent.nielsen@nuffield.ox.ac.uk
#
#	This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#######################################################

apc.data.sums	<- 	function(apc.data.list,data.type="r",average=FALSE,keep.incomplete=TRUE,apc.index=NULL)
#	BN 15 Aug 2018: added possibility for sums/averages over incomplete sequences.
#	BN 17 Nov 2016
#	BN 15 Dec 2013
#	input:	apc.data.list		
#			data.type			optional. Character.  "r" response only, "d" dose only, "m" mortality rates.
#			average				optional. Logical. sums if TRUE, averages if FALSE. Default=FALSE
#			keep.incomplete		optional. Logical. If true perform calculation for incomplete sequences by removing NAs.
#								If false incomplete sequences are NA.  Default=TRUE
#			apc.index			generated by apc.index	
#	output:				sums.age	sums of data by age
#						sums.per	sums of data by period
#						sums.coh	sums of data by cohort
#	note:	if apc.index supplied then it suffices to input
#			apc.data.list = list(response=response)	  ... if data.type="r"
#			apc.index does not need to be a full apc.index list. Sufficient entries are
#						age.max
#						per.max
#						coh.max
#						index.trap
#						index.data
#						per.zero
{ 	#	apc.data.sums
	######################
	#	get index
	if(is.null(apc.index)==TRUE)
		apc.index	<- apc.get.index(apc.data.list)
	######################
	#	set plot array dimensions
	if(data.type %in% c("d","m") & is.null(apc.data.list$dose)==TRUE) 	return(cat("apc.error: Doses are not available \n"))
	if(data.type == "r")	data.matrix	<- apc.data.list$response
	if(data.type == "d") 	data.matrix	<- apc.data.list$dose
	if(data.type == "m") 	data.matrix	<- apc.data.list$response / apc.data.list$dose
	#################
	#	construct trapezoid matrix in age/cohort coordinates
	trap		<- matrix(data=NA,nrow=apc.index$age.max,ncol=apc.index$coh.max)
	trap[apc.index$index.trap]	<- data.matrix[apc.index$index.data]
#	sums.age	<- rowSums(trap,na.rm=TRUE)
#	sums.coh	<- colSums(trap,na.rm=TRUE)
	#	construct trapezoid matrix in age/period coordinates
	trap.ap		<- matrix(data=NA,nrow=apc.index$age.max,ncol=apc.index$per.max)
	for(row in 1:apc.index$age.max)
	{
		col.lower	<- max(1,apc.index$per.zero+2-row)
		col.upper	<- min(apc.index$coh.max,apc.index$per.zero+1-row+apc.index$per.max)
		per.lower	<- max(1,row-apc.index$per.zero)
		per.upper	<- col.upper-col.lower+per.lower	
		trap.ap[row,per.lower:per.upper]	<- trap[row,col.lower:col.upper]
	}
	#	take sums
	if(average==TRUE)	
	{	sums.age	<- rowMeans(trap.ap,na.rm=keep.incomplete)
		sums.coh	<- colMeans(trap,	na.rm=keep.incomplete)
		sums.per	<- colMeans(trap.ap,na.rm=keep.incomplete)
	}
	if(average==FALSE)
	{	sums.age	<- rowSums(trap.ap,na.rm=keep.incomplete)	#ifelse(apply(is.na(trap   ),1,all),NA,rowSums(trap   ,na.rm=TRUE))
		sums.coh	<- colSums(trap   ,na.rm=keep.incomplete)	#ifelse(apply(is.na(trap   ),1,all),NA,colSums(trap   ,na.rm=TRUE))
		sums.per	<- colSums(trap.ap,na.rm=keep.incomplete)	#ifelse(apply(is.na(trap.ap),1,all),NA,colSums(trap.ap,na.rm=TRUE))
	}	
	#	return				
	return(list(sums.age=sums.age,
				sums.per=sums.per,
				sums.coh=sums.coh))
}	#	apc.data.sums

apc.plot.data.sums	<- function(apc.data.list,data.type="a",average=FALSE,keep.incomplete=TRUE,apc.index=NULL,type="o",log="",main.outer=NULL,main.sub=NULL)
#	BN 15 Aug 2018: added possibility for sums/averages over incomplete sequences.
#	BN 15 Aug 2018:	updated so that one can get averages 
#	BN 15 Dec 2013
#	input:	apc.data.list	list. Data
#			data.type
#			average			logical. Sums if FALSE, Averages if TRUE. Default=FALSE
#			keep.incomplete		optional. Logical. If true perform calculation for incomplete sequences by removing NAs.
#								If false incomplete sequences are NA.  Default=TRUE
#			apc.index		optional list. Generated by apc.index	
#			type			optional plot parameter. Character. "o" if overlaid points and lines. "l" if lines. "p" if points.  Default is "o".
#			log				optional plot parameter. Character. "y" if y-scale is logarithmic, otherwise ""
#			data.type		optional. Character.  "r" response only, "d" dose only, "m" mortality rates, "a" all (default).
#	output:				plots of age sums, cohort sums, period sums
{	#	apc.plot.data.sums
	######################
	#	get index
	if(is.null(apc.index)==TRUE)
		apc.index	<- apc.get.index(apc.data.list)
	######################
	#	get title
	if(is.null(main.outer)==TRUE)
		main.outer	<- "Data sums by age/period/cohort index"			
	######################
	#	set plot array dimensions
	if(data.type %in% c("r","d","m") | is.null(apc.data.list$dose)==TRUE) 	par(mfrow=c(1,3))
	if(data.type == "a" & is.null(apc.data.list$dose)==FALSE) 				par(mfrow=c(3,3))
	if(data.type %in% c("d","m") & is.null(apc.data.list$dose)==TRUE) 		return(cat("apc.plot.data.sums error: Doses are not available \n"))
	par(mar=c(5,5,2,0),oma=c(0,0,1,1))
	######################
	#	String for plot
	#	BN 15 Aug 2018
	if(average==TRUE)	s.ylab="averages of data"	
	if(average==FALSE)	s.ylab="sums of data"	
	######################
	#	plot response
	if(data.type %in% c("r","a"))
	{	sums.response	<- apc.data.sums(apc.data.list,data.type="r",average,keep.incomplete,apc.index)
		if(is.null(main.sub)==TRUE)	main	<- "response"
		plot(seq(from=apc.index$age1,by=apc.index$unit,length=apc.index$age.max),sums.response$sums.age,main=main,xlab="age"   ,ylab=s.ylab,type=type,log=log)
		plot(seq(from=apc.index$per1,by=apc.index$unit,length=apc.index$per.max),sums.response$sums.per,main=main,xlab="period",ylab=s.ylab,type=type,log=log)
		plot(seq(from=apc.index$coh1,by=apc.index$unit,length=apc.index$coh.max),sums.response$sums.coh,main=main,xlab="cohort",ylab=s.ylab,type=type,log=log)
	}
	######################
	#	plot dose
	if(data.type %in% c("d","a") & is.null(apc.data.list$dose)==FALSE)
	{
		sums.dose		<- apc.data.sums(apc.data.list,data.type="d",average,keep.incomplete,apc.index)
		if(is.null(main.sub)==TRUE)	main	<- "dose"
		plot(seq(from=apc.index$age1,by=apc.index$unit,length=apc.index$age.max),sums.dose$sums.age,main=main,xlab="age"   ,ylab=s.ylab,type=type,log=log)
		plot(seq(from=apc.index$per1,by=apc.index$unit,length=apc.index$per.max),sums.dose$sums.per,main=main,xlab="period",ylab=s.ylab,type=type,log=log)
		plot(seq(from=apc.index$coh1,by=apc.index$unit,length=apc.index$coh.max),sums.dose$sums.coh,main=main,xlab="cohort",ylab=s.ylab,type=type,log=log)
	}
	######################
	#	plot mortality rates
	if(data.type %in% c("m","a") & is.null(apc.data.list$dose)==FALSE)
	{
		sums.dose		<- apc.data.sums(apc.data.list,data.type="m",average,keep.incomplete,apc.index)
		if(is.null(main.sub)==TRUE)	main	<- "rates"
		plot(seq(from=apc.index$age1,by=apc.index$unit,length=apc.index$age.max),sums.dose$sums.age,main=main,xlab="age"   ,ylab=s.ylab,type=type,log=log)
		plot(seq(from=apc.index$per1,by=apc.index$unit,length=apc.index$per.max),sums.dose$sums.per,main=main,xlab="period",ylab=s.ylab,type=type,log=log)
		plot(seq(from=apc.index$coh1,by=apc.index$unit,length=apc.index$coh.max),sums.dose$sums.coh,main=main,xlab="cohort",ylab=s.ylab,type=type,log=log)
	}
	title(main.outer,outer=TRUE)
}	#	apc.plot.data.sums

#apc.plot.data.sparsity	<- 	function(apc.data.list,data.type="a",apc.index=NULL,sparsity.limits=c(1,2),legend=TRUE,cex=NULL,pch=c(15,0),col.gray=c(0,0),main.outer=NULL)
##	BN 25 Apr 2015 (25 nov 2013)
##	input:	apc.data.list		list. Data
##			apc.index			optional. List. Generated by apc.index	
##			sparsity.limits		optional. Vector with two entries giving sparsity thresholds, default c(1,2)
##			data.type			optional. Character.  "r" response only, "d" dose only, "a" all (default).
##	output:	plots indicating where data are sparse
#{	#	apc.plot.data.sparsity	
#	######################
#	#	get index
#	if(is.null(apc.index)==TRUE)	apc.index	<- apc.get.index(apc.data.list)
#	######################
#	#	get title
#	if(is.null(main.outer)==TRUE)	main.outer	<- paste("Sparsity")			
##		main.outer	<- paste("Sparsity plots \n","(black <",as.character(sparsity.limits[1]),",grey <",as.character(sparsity.limits[2]),")")			
#	######################
#	#	legend plot dimension correction
#											legend.cor	<- 0
#	if(legend && max(sparsity.limits<10))	legend.cor	<- 5		
#	if(legend && max(sparsity.limits>=10))	legend.cor	<- 5.5		
#	######################
#	#	get variables
#	data.format	<- apc.index$data.format
#	xlab		<- apc.index$data.xlab	
#	ylab		<- apc.index$data.ylab	
#	x1			<- apc.index$data.x1		
#	y1			<- apc.index$data.y1
#	xmax		<- apc.index$data.xmax
#	ymax		<- apc.index$data.ymax
#	unit		<- apc.index$unit
#	######################
#	#	set limits	
#	xlim	<- c(x1,x1 +(xmax-1)*unit)
#	if(data.format=="CL")	ylim	<- c(y1 +(ymax-1)*unit,y1)	
#	else					ylim	<- c(y1,y1 +(ymax-1)*unit)
#	######################
#	#	set plot array dimensions
#	if(data.type %in% c("r","d") | is.null(apc.data.list$dose)==TRUE) 	par(mfrow=c(1,1))
#	if(data.type == "a" & is.null(apc.data.list$dose)==FALSE) 			par(mfrow=c(1,2))
#	if(data.type == "d" & is.null(apc.data.list$dose)==TRUE) 			return(cat("apc.error: Doses are not available \n"))
#	par(mar=c(5,5,2,legend.cor),oma=c(0,0,1,1),xpd=TRUE)
#	######################
#	#	plot response
#	function.sparsity.plot	<- function(data.matrix,main)
#	{	if(is.null(cex)==TRUE)
#		{	nmax	<- max(xmax,ymax)
#						cex	<- 0.5
#			if(nmax<20) cex	<- 1
#			if(nmax<10) cex	<- 2
#			if(nmax<5 ) cex	<- 5
#		}
#		plot(1,1,pch=NA,xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,main=main)
#		for(row in 1:xmax)
#			for(col in 1:ymax)
#			{
#				x	<- x1+(row-1)*unit
#				y	<- y1+(col-1)*unit
#				z	<- data.matrix[row,col]
#				if(is.na(z)==TRUE)				points(x,y,cex=cex,pch="x")
#				else
#				{
#					if(z<sparsity.limits[2] && z>=sparsity.limits[1])
#						points(x,y,cex=cex,pch=pch[2],col=gray(col.gray[2]))	
#					if(z<sparsity.limits[1])
#						points(x,y,cex=cex,pch=pch[1],col=gray(col.gray[1]))
#				}	
#			}
#		if(legend)	
#			legend("topright",
##					inset=c(-0.2*legend.cor/5,0),
#					inset=c(-0.2,0),
#					pch=pch,col=gray(col.gray),
#					legend=paste(c("<","<"),as.character(sparsity.limits)) )
#	}
#	if(data.type %in% c("r","a"))
#		function.sparsity.plot(apc.data.list$response,"responses")
#	#	plot dose
#	if(data.type %in% c("d","a") & is.null(apc.data.list$dose)==FALSE)
#		function.sparsity.plot(apc.data.list$dose,"doses")	
#	title(main.outer,outer=TRUE)
#}	#	apc.plot.data.sparsity

apc.plot.data.sparsity	<- 	function(apc.data.list,data.type="a",swap.axes=FALSE,apc.index=NULL,sparsity.limits=c(1,2),cex=NULL,pch=15,main.outer=NULL)
#	BN 7 feb 2015
#	BN 27 apr 2017: swap.axes added
#	input:	apc.data.list		list. Data
#			data.type			"r": response only
#								"d": dose only
#								"a": all, so both response and dose
#			swap.axes				logical. If true swap data matrix so plot axes and matrix axes match.  Default is FALSE unless data.format=CL
#			apc.index			optional. List. Generated by apc.index	
#			sparsity.limits		optional. Vector with two entries giving sparsity thresholds, default c(1,2)
#			data.type				optional. Character.  "r" response only, "d" dose only, "a" all (default).
#	output:	plots indicating where data are sparse
{	#	apc.plot.data.sparsity	
	######################
	#	get index
	if(is.null(apc.index)==TRUE)
		apc.index	<- apc.get.index(apc.data.list)
	######################
	#	get title
	if(is.null(main.outer)==TRUE)
		main.outer	<- paste("Sparsity plots \n","(black <",as.character(sparsity.limits[1]),",grey <",as.character(sparsity.limits[2]),")")			
	######################
	#	get variables
	data.format	<- apc.index$data.format
	xlab		<- apc.index$data.xlab	
	ylab		<- apc.index$data.ylab	
	x1			<- apc.index$data.x1		
	y1			<- apc.index$data.y1
	xmax		<- apc.index$data.xmax
	ymax		<- apc.index$data.ymax
	unit		<- apc.index$unit
	######################
	#	swap.axes for CL
	if(data.format=="CL")	swap.axes=1-swap.axes
	######################
	#	set limits	
	xlim	<- c(x1,x1 +(xmax-1)*unit)
	ylim	<- c(y1,y1 +(ymax-1)*unit)
	if(data.format=="CL")
		if(swap.axes)	xlim	<- c(x1 +(xmax-1)*unit,x1)	
		else			ylim	<- c(y1 +(ymax-1)*unit,y1)	
	######################
	#	set plot array dimensions
	if(data.type %in% c("r","d") | is.null(apc.data.list$dose)==TRUE) 	par(mfrow=c(1,1))
	if(data.type == "a" & is.null(apc.data.list$dose)==FALSE) 			par(mfrow=c(1,2))
	if(data.type == "d" & is.null(apc.data.list$dose)==TRUE) 			return(cat("apc.error: Doses are not available \n"))
	par(mar=c(5,5,2,0),oma=c(0,0,2,1))
	######################
	#	plot response
	function.sparsity.plot	<- function(data.matrix,main,swap.axes)
	{	if(is.null(cex)==TRUE)
		{	nmax	<- max(xmax,ymax)
						cex	<- 0.5
			if(nmax<20) cex	<- 1
			if(nmax<10) cex	<- 2
			if(nmax<5 ) cex	<- 5
		}
		if(swap.axes){	l.x1 <- y1;	l.y1 <- x1; l.xmax <- ymax; l.ymax <- xmax }
		else		 {	l.x1 <- x1;	l.y1 <- y1;	l.xmax <- xmax; l.ymax <- ymax }
		if(swap.axes){
				data.matrix <- t(data.matrix)
				plot(1,1,pch=NA,xlim=ylim,ylim=xlim,xlab=ylab,ylab=xlab,main=main)
		}
		else 	plot(1,1,pch=NA,xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,main=main)
		for(row in 1:l.xmax)
			for(col in 1:l.ymax)
			{
				x	<- l.x1+(row-1)*unit
				y	<- l.y1+(col-1)*unit
				z	<- data.matrix[row,col]
				if(is.na(z)==TRUE)				points(x,y,cex=cex,pch="x")
				else
				{
					if(z<sparsity.limits[2])	points(x,y,cex=cex,pch=pch,col=gray(0.66))	
					if(z<sparsity.limits[1]) 	points(x,y,cex=cex,pch=pch,col=gray(0.0))
				}	
			}
	}
	if(data.type %in% c("r","a"))
		function.sparsity.plot(apc.data.list$response,"responses",swap.axes)
	#	plot dose
	if(data.type %in% c("d","a") & is.null(apc.data.list$dose)==FALSE)
		function.sparsity.plot(apc.data.list$dose,"doses",swap.axes)	
	title(main.outer,outer=TRUE)
}	#	apc.plot.data.sparsity

apc.plot.data.within	<- 	function(apc.data.list,data.type="r",plot.type="awc",average=FALSE,thin=NULL,apc.index=NULL,ylab=NULL,type="o",log="",legend=TRUE,lty=1:5,col=1:6,bty="n",main=NULL,x="topleft",return=FALSE)
#	BN 17 Nov 2016: NA dealt with better.
#	BN 25 apr 2015
#	input:	apc.data.list	list. Data
#			plot.type		Character. "awp", "pwa" "awc", "cwa, "cwp", "pwc": for example: "awp" gives time series in age within each period level: for an AP data-array these are the column sums.  
#			apc.index		optional. List. Generated by apc.index	
#			data.type		optional. Character.  "r" or "response": response, "d" of "dose": dose, "m" or "mortality" or "rates": mortality.  Default is "r".
#			average			optional. If TRUE/FALSE reports averages/sums. Default is FALSE.
#			thin			Optional. Numerical.  age/periods/cohorts are grouped in groups of size thin.  Default=NULL
#			ylab			optional plot parameter. Character: common label on y axes, for instance "log mortality rates". Default ""
#			type			optional plot parameter. Character. "o" if overlaid points and lines. "l" if lines. "p" if points.  Default is "o".
#			log				optional plot parameter. Character: common scale for y axes, Default: "" for normal scale. Use "y" if log scale
#			legend			optional plot parameter. Logical: should legends be drawn. Default: "TRUE"
#			lty				optional plot parameter. Vector of line types.
#							The first element is for the first column, the second element for the second column, etc.,
#							even if lines are not plotted for all columns. Line types will be used cyclically
#							until all plots are drawn.  Default: 1:5
#			col				optional plot parameter. Vector of colors.
#							The first element is for the first column, the second element for the second column, etc.,
#							even if lines are not plotted for all columns. Colors will be used cyclically
#							until all plots are drawn.  Default is 1:6.
#			bty				optional plot parameter. Character: the type of box to be drawn around the legend.
#							The allowed values are "n" and "o".	Default is "n".
#			x				optional plot parameter: x legend. Location of legend, Default is "topleft"
#			return			Logical. If TRUE returns
#							x	vector.	time coordinates
#							m	matrix. columns are the outcomes for different within groups
#	output:	6 plots of data: age vs period, cohort vs age, cohort vs period
{ 	#	apc.plot.data.within
	######################
	#	check input
	if((plot.type %in% c("awp","awc","cwp","pwa","cwa","pwc"))==FALSE)
		return(cat("apc.plot.data.within error: plot.type not recognised \n"))
	l.data.type <- data.type
	if(data.type == "response") 	l.data.type <- "r"
	if(data.type == "dose")			l.data.type <- "d"
	if(data.type ==	"rates")		l.data.type <- "m"
	if(data.type == "mortality")	l.data.type <- "m"
	if(l.data.type %in% c("d","m") & is.null(apc.data.list$dose)==TRUE)
		return(cat("apc.plot.data.within error: Doses are not available \n"))
	######################
	#	get index
	if(is.null(apc.index)==TRUE)
		apc.index	<- apc.get.index(apc.data.list)
	######################
	#	get data type
	if(l.data.type == "r")	
		{ title.add.to <- "response"; 	m <- apc.data.list$response;	}
	if(l.data.type == "d")
		{ title.add.to <- "dose"; 		m <- apc.data.list$dose;	}
	if(l.data.type == "m")
		{ title.add.to <- "rates"; 		m <- apc.data.list$response/apc.data.list$dose;	}
	######################
	#	get trapezoid in age/cohort coordinate system
	trap.m	<- matrix(data=NA,nrow=apc.index$age.max,ncol=apc.index$coh.max)
	trap.m[apc.index$index.trap]	<- m[apc.index$index.data]
	######################
	#	get ylab
	if(is.null(ylab))
	{
		if(log=="y")	ylab <- "log "
		if(l.data.type=="r")	ylab <- paste(ylab,"response",sep="")
		if(l.data.type=="d")	ylab <- paste(ylab,"dose",sep="")
		if(l.data.type=="m")	ylab <- paste(ylab,"rate",sep="")
	}	
	######################
	function.trapezoid.to.ap	<- function(trapezoid,transpose=FALSE)
	#	BN, 17 Sep 2013
	#	transforms a trapezoid in AC format to a trapezoid in AP format
	{	#	function.trapezoid.to.ap
		######################
		#	get values
		age.max		<- apc.index$age.max
		per.max		<- apc.index$per.max
		coh.max		<- apc.index$coh.max
		per.zero	<- apc.index$per.zero
		nrow		<- age.max
		ncol		<- coh.max
		if(transpose==TRUE)
		{
			trapezoid	<- t(trapezoid)
			nrow		<- coh.max
			ncol		<- age.max
		}	
		#	create ap matrix
		m			<- matrix(data=NA,nrow=nrow,ncol=per.max)
		for(row in 1:nrow)
		{
			col.lower	<- max(1,per.zero+2-row)
			col.upper	<- min(ncol,per.zero+1-row+per.max)
			per.lower	<- max(1,row-per.zero)
			per.upper	<- col.upper-col.lower+per.lower	
			m[row,per.lower:per.upper]	<- trapezoid[row,col.lower:col.upper]
		}
		return(m)
	}	#	function.trapezoid.to.ap
	######################
	#	choose two time scales
	#	1. "awp": age within period
	if(plot.type=="awp")
	{	m.to.plot	<- function.trapezoid.to.ap(trap.m);
		x1		<- apc.index$age1;	x.max	<- apc.index$age.max;	xlab	<- "age";
		w1		<- apc.index$per1;	w.max	<- apc.index$per.max;	title.sub	<- "within period"	}
	#	2. "awc": age within cohort	
	if(plot.type=="awc")
	{	m.to.plot	<- trap.m;
		x1		<- apc.index$age1;	x.max	<- apc.index$age.max;	xlab	<- "age";
		w1		<- apc.index$coh1;	w.max	<- apc.index$coh.max;	title.sub	<- "within cohort"	}
	#	3. "cwp": cohort within period
	if(plot.type=="cwp")
	{	m.to.plot	<- function.trapezoid.to.ap(trap.m,transpose=TRUE);
		x1		<- apc.index$coh1;	x.max	<- apc.index$coh.max;	xlab	<- "cohort";
		w1		<- apc.index$per1;	w.max	<- apc.index$per.max;	title.sub	<- "within period"	}
	#	4. "pwa": period within age
	if(plot.type=="pwa")
	{	m.to.plot	<- t(function.trapezoid.to.ap(trap.m));
		x1		<- apc.index$per1;	x.max	<- apc.index$per.max;	xlab	<- "period";
		w1		<- apc.index$age1;	w.max	<- apc.index$age.max;	title.sub	<- "within age"	}
	#	5. "cwa": cohort within age
	if(plot.type=="cwa")
	{	m.to.plot	<- t(trap.m);
		x1		<- apc.index$coh1;	x.max	<- apc.index$coh.max;	xlab	<- "cohort";
		w1		<- apc.index$age1;	w.max	<- apc.index$age.max;	title.sub	<- "within age"	}
	#	6. "pwc": period within cohort
	if(plot.type=="pwc")
	{	m.to.plot	<- t(function.trapezoid.to.ap(trap.m,transpose=TRUE));
		x1		<- apc.index$per1;	x.max	<- apc.index$per.max;	xlab	<- "period";
		w1		<- apc.index$coh1;	w.max	<- apc.index$coh.max;	title.sub	<- "within cohort"	}
	######################
	#	get thinning value function
	function.thin.value		<- function(m,thin=NULL)
	#	BN 27 Jun 2014
	#	function for getting default thin value for grouping
	{	#	function.thin.value
		l.thin	<- thin
		if(is.null(l.thin)==TRUE)
			{	ncol	<- ncol(m)
							l.thin <- 1
				if(ncol>10) l.thin <- 2
				if(ncol>20) l.thin <- 4
				if(ncol>40) l.thin <- 8
				if(ncol>80) l.thin <- 16
			}
		return(l.thin)
	}	#	function.thin.value
	######################
	#	thinning function
	function.thin.matrix	<- function(m,l.thin)
	#	BN 17 Nov 2016: averages added.
	#	BN 27 Jun 2013
	#	function for grouping columns: it takes row sums within each group
	{ 	#	function.thin.matrix
		if(l.thin==1)	mm <- m
		else
		{
			nrow	<- nrow(m)
			ncol	<- ncol(m)
			ngroup	<- ceiling(ncol/l.thin)
			mm		<- matrix(data=NA,nrow=nrow,ncol=ngroup)
			if(average==TRUE)
			{	for(group in 1:(ngroup-1))                                                                   
				{	mmm	<- as.matrix(m[,((group-1)*l.thin+1):(group*l.thin)])                                           
					mm[,group]	<- as.matrix(	rowMeans(mmm,na.rm=TRUE))
				}                                                                                            
				mmm	<- as.matrix(m[,((ngroup-1)*l.thin+1):ncol])                                                        
				mm[,ngroup]	<- as.matrix(	rowMeans(mmm,na.rm=TRUE))
			}
			if(average==FALSE)
			{	for(group in 1:(ngroup-1))                                                                    
				{	mmm	<- as.matrix(m[,((group-1)*l.thin+1):(group*l.thin)])                                           
					mm[,group]	<- as.matrix(	rowSums(mmm,na.rm=TRUE))
				#	mm[,group]	<- as.matrix(	ifelse(apply(is.na(mmm),1,all),NA,rowSums(mmm,na.rm=TRUE)))   
				}                                                                                             
				mmm	<- as.matrix(m[,((ngroup-1)*l.thin+1):ncol])                                                         
				mm[,ngroup]	<- as.matrix(	rowSums(mmm,na.rm=TRUE))       
				#	mm[,ngroup]	<- as.matrix(	ifelse(apply(is.na(mmm),1,all),NA,rowSums(mmm,na.rm=TRUE)))       
			}
			if(ngroup != ncol/l.thin)
				print("apc.plot.data.within warning: maximal index not divisible by thin, so last group smaller than other groups")
		}
		return(mm)	
	}	#	function.thin.matrix
	######################
	#	thin matrix
	l.thin	<- function.thin.value(m.to.plot,thin)
	m.to.plot	<- function.thin.matrix(m.to.plot,l.thin=l.thin)
	######################
	#	set plot array dimensions
	old.par <- par()
	par(mar=c(5,5,3,1)+0.1)
	######################
	#	plot
	v.x		<- seq(from=x1,length=x.max,by=apc.index$unit)
	if(is.null(main))	main <- title.sub
	matplot(v.x,m.to.plot,type=type,pch=20,log=log,lty=lty,col=col,xlab=xlab,ylab=ylab,main=main)
	within	<- seq(from=w1,length=ceiling(w.max/l.thin),by=apc.index$unit*l.thin)
	if(legend==TRUE)	legend(x=x,legend=as.character(within),lty=lty,col=col,bty=bty)
	######################
	#	return to old par
	par	<- old.par
	######################
	#	return data sums
	if(return==TRUE)
	{	dimnames(m.to.plot)	<- list(as.character(v.x),as.character(within))
		return(m.to.plot)	}
}	#	apc.plot.data.within

apc.plot.data.within.all.six	<- 	function(apc.data.list,data.type="r",average=FALSE,thin=NULL,apc.index=NULL,ylab=NULL,type="o",log="",legend=TRUE,lty=1:5,col=1:6,bty="n",main.outer=NULL,x="topleft")
#	BN 25 apr 2015
{	#	apc.plot.data.within.all.six
	######################
	#	get index
	if(is.null(apc.index)==TRUE)
		apc.index	<- apc.get.index(apc.data.list)
	######################
	if(is.null(main.outer))	main.outer <- "plots of data using two indices"
	old.par <- par()
	par(mfrow=c(2,3))
	par(oma=c(0,0,2,0))
	l.data.type = c("awp","awc","cwp","pwa","cwa","pwc")
	for(i in 1:6)
		apc.plot.data.within(apc.data.list,data.type=data.type,plot.type=l.data.type[i],average=average,thin=thin,apc.index=apc.index,ylab=ylab,type=type,log=log,legend=legend,lty=lty,col=col,bty=bty,x=x)
	title(main.outer,outer=TRUE)	
	par	<- old.par	
}	#	apc.plot.data.within.all.six

apc.plot.data.level	<- 	function(apc.data.list,data.type="r",rotate=FALSE,apc.index=NULL,main=NULL,lab=NULL,contour=FALSE,colorkey=TRUE)
#	BN 25 apr 2015
#	input:	apc.data.list	list. Data
#			data.type		optional. Character.
#							"r" or "response": response,
#							"d" of "dose": dose,
#							"m" or "mortality" or "rates": mortality.  
#							"residual": residual
#							"fitted.values": fitted values
#							"linear.predictors": linear predictors
#							Default is "r".
#			main			optional. Character. Main title
#			lab				optional plot parameter.	A numerical vector of the form c(x, y, len)
#							which modifies the default way that axes are annotated.
#							The values of x and y give the (approximate) number of tickmarks on
#							the x and y axes. len is not implemented.
#			contour			optional levelplot (lattice) parameter. Logical. Contour lines drawn if TRUE. Default FALSE.
#			colorkey		optional levelplot (lattice) parameter. Logical or list. Determines color key. Default TRUE.
{	#	apc.plot.data.level
	######################
	#	get index
	if(is.null(apc.index)==TRUE)
		apc.index	<- apc.get.index(apc.data.list)
	######################
	#	check input
	l.data.type <- data.type
	if(data.type == "response") 	l.data.type <- "r"
	if(data.type == "dose")			l.data.type <- "d"
	if(data.type ==	"rates")		l.data.type <- "m"
	if(data.type == "mortality")	l.data.type <- "m"
	if(l.data.type %in% c("d","m") & is.null(apc.data.list$dose)==TRUE)
		return(cat("apc.plot.data.within error: Doses are not available \n"))
	######################
	#	get data type
	if(l.data.type == "r")	
		{ l.main <- "response"; 	m <- apc.data.list$response;	}
	if(l.data.type == "d")
		{ l.main <- "dose"; 		m <- apc.data.list$dose;	}
	if(l.data.type == "m")
		{ l.main <- "rates"; 		m <- apc.data.list$response/apc.data.list$dose;	}
	if(l.data.type == "residual")
		{ l.main <- "residual"; 	m <- apc.data.list$m.residual;	}
	if(l.data.type == "fitted.values")
		{ l.main <- "fitted values"; 	m <- apc.data.list$m.fitted.values;	}
	if(l.data.type == "linear.predictors")
		{ l.main <- "linear predictors"; 	m <- apc.data.list$m.linear.predictors;	}
	if(is.null(main))	main <- l.main	
	######################
	#	Rotate
	l.rotate	<- rotate
	if(apc.data.list$data.format=="CL")	l.rotate <- 1-l.rotate
	######################
	#	Generate nice labels
	x1	<- apc.index$data.x1;	xmax <- apc.index$data.xmax;	xlab <- apc.index$data.xlab;
	y1	<- apc.index$data.y1;	ymax <- apc.index$data.ymax;	ylab <- apc.index$data.ylab;
	unit<- apc.index$unit	
	if(is.null(lab))	lab[1:2] <- c(min(5,xmax),min(5,ymax)) 	
	x.at <- seq(from=1,to=xmax,length=lab[1])
	x.lab<- as.character(x1+(x.at-1)*unit)
	if(l.rotate)	x.lab <- x.lab[length(x.lab):1]
	x	 <- list(at=x.at,labels=x.lab)
	y.at <- seq(from=1,to=ymax,length=lab[2])
	y.lab<- as.character(y1+(y.at-1)*unit)
	y	 <- list(at=y.at,labels=y.lab)
	######################
	# 	Level plot
	if(l.rotate)
		print(levelplot(t(m[nrow(m):1,]),xlab=ylab,ylab=xlab,scales=list(x=y,y=x),main=main,contour=contour,colorkey=colorkey))
	else	
		print(levelplot(m,xlab=xlab,ylab=ylab,scales=list(x=x,y=y),main=main,contour=contour,colorkey=colorkey))
}	#	apc.plot.data.level

apc.plot.data.all	<- 	function(apc.data.list,log="",rotate=FALSE)
#	BN 25 Apr 2015 (25 nov 2013)
#	input:	apc.data.list	list. Data
#			log				optional plot parameter. Character: common scale for y axes, Default: "y" for log scale. Use "" if normal scale
#	output:	all descriptive plots.
{	#	apc.plot.data.all
	apc.index	<- apc.get.index(apc.data.list)
					apc.plot.data.sums(apc.data.list,log=log,apc.index=apc.index)
		dev.new(); 	apc.plot.data.sparsity(apc.data.list,apc.index=apc.index)
		dev.new(); 	apc.plot.data.within.all.six(apc.data.list,"r",log=log,apc.index=apc.index)
	if(is.null(apc.data.list$dose)==FALSE)
	{	dev.new();	apc.plot.data.within.all.six(apc.data.list,"d",log=log,apc.index=apc.index) 	
		dev.new();	apc.plot.data.within.all.six(apc.data.list,"m",log=log,apc.index=apc.index) 			
	}
		dev.new();	apc.plot.data.level(apc.data.list,"r",rotate=rotate,apc.index=apc.index)
	if(is.null(apc.data.list$dose)==FALSE)
	{	dev.new();	apc.plot.data.level(apc.data.list,"d",rotate=rotate,apc.index=apc.index)
		dev.new();	apc.plot.data.level(apc.data.list,"m",rotate=rotate,apc.index=apc.index)
	}
}	#	apc.plot.data.all