Skip to content

Commit 25572df

Browse files
JamesRamsay5JamesRamsay5
committed
fRegress made more readable
1 parent 14d76be commit 25572df

File tree

7 files changed

+383
-280
lines changed

7 files changed

+383
-280
lines changed

‎NAMESPACE

+1
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,7 @@ export(AmpPhaseDecomp,
6767
fRegress.fd,
6868
fRegress.CV,
6969
fRegress.stderr,
70+
fRegress.double,
7071
fRegress.formula,
7172
Fstat.fd,
7273
geigen,

‎R/Fperm.fd.R

+114-73
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,5 @@
1+
#Fperm.fd(annualprec, xfdlist, betalist, nperm=10)
2+
13
Fperm.fd <- function(yfdPar, xfdlist, betalist, wt=NULL,
24
nperm=200,argvals=NULL,q=0.05,plotres=TRUE,...)
35
{
@@ -7,96 +9,135 @@ Fperm.fd <- function(yfdPar, xfdlist, betalist, wt=NULL,
79
# q = quantile to compare
810
# plotres: Do we plot the results?
911

12+
# set up containers for null results
13+
1014
Fnull = rep(0,nperm)
1115
Fnullvals = c()
12-
13-
q = 1-q
14-
16+
17+
# switch q from residual to main body
18+
19+
q <- 1-q
20+
21+
# keep track of processing time
22+
1523
begin <- proc.time()
1624
fRegressList <- fRegress(yfdPar, xfdlist, betalist)
1725
elapsed.time <- max(proc.time()-begin,na.rm=TRUE)
18-
19-
if( elapsed.time > 30/nperm ){
26+
if( elapsed.time > 30/nperm ) {
2027
print(paste('Estimated Computing time =',
2128
round(nperm*elapsed.time),'seconds.'))
2229
}
23-
24-
yhat <- fRegressList$yhatfdobj
25-
if(is.null(yhat)){ yhat = fRegressList$yhat }
30+
31+
#. obtain fit vector yhat from fRegressList
32+
33+
yhat <- fRegressList$yhatfdobj. #. the fitted values
34+
if(is.null(yhat)) yhat = fRegressList$yhat
2635
if(is.list(yhat) && ('fd' %in% names(yhat))) yhat <- yhat$fd
27-
36+
37+
#. initial invokation of Fstat.fd
38+
2839
tFstat <- Fstat.fd(yfdPar,yhat,argvals)
29-
40+
41+
# two sets of results
42+
3043
Fvals <- tFstat$F
31-
Fobs = max(Fvals)
32-
33-
argvals = tFstat$argvals
34-
44+
Fobs <- max(Fvals)
45+
46+
#. obtain argument values
47+
48+
argvals <- tFstat$argvals
49+
50+
#. determine length of data vector
51+
3552
if(is.vector(yfdPar)){
36-
n = length(yfdPar)
53+
n <- length(yfdPar)
3754
} else {
38-
n = ncol(yfdPar$coefs)
55+
n <- ncol(yfdPar$coefs)
3956
}
40-
57+
58+
#. --------------------------------------------------------------------------
59+
# Loop through permutations
60+
#. --------------------------------------------------------------------------
61+
4162
for(i in 1:nperm){
42-
43-
tyfdPar = yfdPar[sample(n)]
44-
63+
64+
tyfdPar <- yfdPar[sample(n)]
65+
66+
# analyze of a vector of observations
67+
4568
fRegressList <- fRegress(tyfdPar, xfdlist, betalist)
46-
47-
if(is.fd(yhat)){
48-
yhat <- fRegressList$yhatfdobj
49-
if(is.list(yhat) && ('fd' %in% names(yhat))) yhat <- yhat$fd
50-
} else{ yhat = fRegressList$yhat }
51-
52-
tFstat = Fstat.fd(tyfdPar,yhat,argvals)
53-
69+
70+
# extract fit vector from fRegressList
71+
72+
if(is.fd(yhat)) {
73+
yhat <- fRegressList$yhatfdobj
74+
if(is.list(yhat) && ('fd' %in% names(yhat))) yhat <- yhat$fd
75+
} else {
76+
yhat <- fRegressList$yhat
77+
}
78+
79+
# ith analysis by Fstat.fd
80+
81+
tFstat <- Fstat.fd(tyfdPar,yhat,argvals)
82+
83+
# add Fnullvals to container Fnullvals
84+
5485
Fnullvals <- cbind(Fnullvals,tFstat$F)
55-
56-
Fnull[i] = max(Fnullvals[,i])
86+
87+
#.
88+
Fnull[i] <- max(Fnullvals[,i])
5789
}
58-
59-
pval = mean( Fobs < Fnull )
60-
qval = quantile(Fnull,q)
61-
62-
pvals.pts = apply(Fvals<Fnullvals,1,mean)
63-
qvals.pts = apply(Fnullvals,1,quantile,q)
64-
65-
if(plotres){
66-
if(is.fd(yfdPar)){
67-
ylims = c(min(c(Fvals,qval,qvals.pts)),max(c(Fobs,qval)))
68-
69-
if( is.null(names(yhat$fdnames)) ){ xlab = 'argvals' }
70-
else{ xlab = names(yhat$fdnames)[1] }
71-
72-
plot(argvals,Fvals,type="l",ylim=ylims,col=2,lwd=2,
73-
xlab=xlab,ylab='F-statistic',main='Permutation F-Test',...)
74-
lines(argvals,qvals.pts,lty=3,col=4,lwd=2)
75-
abline(h=qval,lty=2,col=4,lwd=2)
76-
77-
legendstr = c('Observed Statistic',
78-
paste('pointwise',1-q,'critical value'),
79-
paste('maximum',1-q,'critical value'))
80-
81-
legend(argvals[1],ylims[2],legend=legendstr,col=c(2,4,4),
82-
lty=c(1,3,2),lwd=c(2,2,2))
83-
}
84-
else{
85-
xlims = c(min(c(Fnull,Fobs)),max(c(Fnull,Fobs)))
86-
hstat = hist(Fnull,xlim=xlims,lwd=2,xlab='F-value',
87-
main = 'Permutation F-Test')
88-
abline(v = Fobs,col=2,lwd=2)
89-
abline(v = qval,col=4,lty=2,lwd=2)
90-
91-
legendstr = c('Observed Statistic',
92-
paste('Permutation',1-q,'critical value'))
93-
94-
legend(xlims[1],max(hstat$counts),legend=legendstr,col=c(2,4),
95-
lty=c(1,2),lwd=c(2,2))
96-
}
90+
91+
#. --------------------------------------------------------------------------
92+
# Display results
93+
#. --------------------------------------------------------------------------
94+
95+
pval <- mean( Fobs < Fnull )
96+
qval <- quantile(Fnull,q)
97+
98+
pvals.pts <- apply(Fvals<Fnullvals,1,mean)
99+
qvals.pts <- apply(Fnullvals,1,quantile,q)
100+
101+
if(plotres) {
102+
#. dispslay results using scattter plots
103+
if(is.fd(yfdPar)){
104+
ylims <- c(min(c(Fvals,qval,qvals.pts)),max(c(Fobs,qval)))
105+
106+
if( is.null(names(yhat$fdnames)) ){ xlab <- 'argvals' }
107+
else{ xlab <- names(yhat$fdnames)[1] }
108+
109+
plot(argvals,Fvals,type="l",ylim=ylims,col=2,lwd=2,
110+
xlab=xlab,ylab='F-statistic',main='Permutation F-Test',...)
111+
lines(argvals,qvals.pts,lty=3,col=4,lwd=2)
112+
abline(h=qval,lty=2,col=4,lwd=2)
113+
114+
legendstr <- c('Observed Statistic',
115+
paste('pointwise',1-q,'critical value'),
116+
paste('maximum',1-q,'critical value'))
117+
118+
legend(argvals[1],ylims[2],legend=legendstr,col=c(2,4,4),
119+
lty=c(1,3,2),lwd=c(2,2,2))
120+
} else {
121+
#. display results with histogram
122+
123+
xlims <- c(min(c(Fnull,Fobs)),max(c(Fnull,Fobs)))
124+
hstat <- hist(Fnull,xlim=xlims,lwd=2,xlab='F-value',
125+
main = 'Permutation F-Test')
126+
abline(v = Fobs,col=2,lwd=2)
127+
abline(v = qval,col=4,lty=2,lwd=2)
128+
129+
legendstr <- c('Observed Statistic',
130+
paste('Permutation',1-q,'critical value'))
131+
132+
legend(xlims[1],max(hstat$counts),legend=legendstr,col=c(2,4),
133+
lty=c(1,2),lwd=c(2,2))
97134
}
98-
99-
return(list(pval=pval,qval=qval,Fobs=Fobs,Fnull=Fnull,
100-
Fvals=Fvals,Fnullvals=Fnullvals,pvals.pts=pvals.pts,qvals.pts=qvals.pts,
101-
fRegressList=fRegressList,argvals=argvals))
135+
}
136+
137+
#. return results
138+
139+
return(list(pval=pval, qval=qval, Fobs=Fobs, Fnull=Fnull,
140+
Fvals=Fvals, Fnullvals=Fnullvals, pvals.pts=pvals.pts,
141+
qvals.pts=qvals.pts,
142+
fRegressList=fRegressList, argvals=argvals))
102143
}

‎R/Fstat.fd.R

+33-34
Original file line numberDiff line numberDiff line change
@@ -1,38 +1,37 @@
11
Fstat.fd <- function(y,yhat,argvals=NULL) {
2-
3-
# observed, predicted and where to evaluate
4-
5-
if( is.numeric(yhat) ){ yhat = as.vector(yhat) }
6-
7-
if( (is.vector(y) & !is.vector(yhat)) | (is.fd(y) &!is.fd(yhat)) ){
8-
stop("y and yhat must both be either scalars or functional data objects.")
2+
3+
# observed, predicted and where to evaluate
4+
5+
if( is.numeric(yhat) ) yhat <- as.vector(yhat)
6+
7+
if( (is.vector(y) & !is.vector(yhat)) | (is.fd(y) &!is.fd(yhat)) ) {
8+
stop("y and yhat must both be either scalars or functional data objects.")
9+
}
10+
11+
12+
if( is.fd(y) ) {
13+
rangeobs <- y$basis$range
14+
rangehat <- yhat$basis$range
15+
16+
if( !prod(rangeobs == rangehat) ){
17+
stop("y and yhat do not have the same range")
918
}
10-
11-
12-
if( is.fd(y) ){
13-
rangeobs = y$basis$range
14-
rangehat = yhat$basis$range
15-
16-
if( !prod(rangeobs == rangehat) ){
17-
stop("y and yhat do not have the same range")
18-
}
19-
20-
21-
if(is.null(argvals)){
22-
argvals = seq(rangeobs[1],rangeobs[2],length.out=101)
23-
}
24-
25-
yvec = eval.fd(argvals,y)
26-
yhatvec = eval.fd(argvals,yhat)
27-
28-
F = apply(yhatvec,1,var)/apply( (yvec-yhatvec)^2,1,mean)
29-
19+
20+
21+
if(is.null(argvals)){
22+
argvals <- seq(rangeobs[1],rangeobs[2],length.out=101)
3023
}
31-
else{
32-
yvec = y
33-
yhatvec = yhat
34-
F = var(yhatvec)/mean( (yvec-yhatvec)^2 )
35-
}
36-
37-
return( list(F=F,argvals=argvals) )
24+
25+
yvec <- eval.fd(argvals,y)
26+
yhatvec <- eval.fd(argvals,yhat)
27+
28+
F <- apply(yhatvec,1,var)/apply( (yvec-yhatvec)^2,1,mean)
29+
30+
} else {
31+
yvec <- y
32+
yhatvec <- yhat
33+
F <- var(yhatvec)/mean( (yvec-yhatvec)^2 )
34+
}
35+
36+
return( list(F=F,argvals=argvals) )
3837
}

0 commit comments

Comments
 (0)