Revision fd92910f232786c2531c1c74e02ece6b49a2f76f authored by J. O. Ramsay on 03 November 2009, 00:00:00 UTC, committed by Gabor Csardi on 03 November 2009, 00:00:00 UTC
1 parent 8ce7a35
fdarm-ch10.m
%
% Ramsay, Hooker & Graves (2009)
% Functional Data Analysis with R and Matlab (Springer)
%
% Remarks and disclaimers
% These R commands are either those in this book, or designed to
% otherwise illustrate how R can be used in the analysis of functional
% data.
% We do not claim to reproduce the results in the book exactly by these
% commands for various reasons, including:
% -- the analyses used to produce the book may not have been
% entirely correct, possibly due to coding and accuracy issues
% in the functions themselves
% -- we may have changed our minds about how these analyses should be
% done since, and we want to suggest better ways
% -- the R language changes with each release of the base system, and
% certainly the functional data analysis functions change as well
% -- we might choose to offer new analyses from time to time by
% augmenting those in the book
% -- many illustrations in the book were produced using Matlab, which
% inevitably can imply slightly different results and graphical
% displays
% -- we may have changed our minds about variable names. For example,
% we now prefer "yearRng" to "dayrange" for the weather data.
% -- three of us wrote the book, and the person preparing these scripts
% might not be the person who wrote the text
% Moreover, we expect to augment and modify these command scripts from time
% to time as we get new data illustrating new things, add functionality
% to the package, or just for fun.
%
% ch. 10. Linear Models for Functional Responses
%
% Set up some strings for constructing paths to folders.
% These strings should be modified so as to provided access
% to the specified folders on your computer.
% Path to the folder containing the Matlab functional data analysis
% software
fdaMPath = 'c:/Program Files/MATLAB/R2009a/fdaM';
addpath(fdaMPath)
% Path to the folder containing the examples
examplesPath = [fdaMPath,'/examples'];
addpath(examplesPath)
%
% Section 10.1 Functional Responses and an Analysis of Variance Model
%
% ----------------- Climate zone effects for temperature -----------------
% Section 10.1.1 Climate Region Effects on Temperature
% path to the daily weather data example folder
weatherPath = [examplesPath,'/weather'];
addpath(weatherPath)
load daily
% set up the data for the analysis
tempav = daily.tempav;
precav = daily.precav;
yearRng = daily.rng;
time = daily.time;
place = daily.place;
dayperiod = daily.period;
region_names = daily.region_names;
regions = daily.regions;
% tempfd from chapter 9
% define the harmonic acceleration operator
harmaccelLfd = vec2Lfd([0,(2*pi/dayperiod)^2,0], dayperiod);
% organize data to have winter = the center of the plot
dayOfYearShifted = [182:365, 1:181];
tempavShifted = tempav(dayOfYearShifted,:);
tempbasis65 = create_fourier_basis(yearRng,65);
lambda = 1e6;
tempfdPar65 = fdPar(tempbasis65, harmaccelLfd, lambda);
tempfd = smooth_basis(time, tempavShifted, tempfdPar65);
station{1} = 'Weather Station';
station{2} = place;
fdnames{1} = 'Day (July 1 to June 30)';
fdnames{2} = station;
fdnames{3} = 'Temperature (deg C)';
tempfd = putnames(tempfd, fdnames);
% augment tempfd by adding a 36th observation with temp(t) = 0
coef = getcoef(tempfd);
coef36 = [coef,zeros(65,1)];
temp36fd = fd(coef36,tempbasis65,fdnames);
% set up the covariates
p = 5;
regionCell = cell(p,1);
regionCell{1} = ones(36,1);
for j=2:p
xj = zeros(36,1);
xj(regions == j-1) = 1;
xj(36) = 1;
regionCell{j} = xj;
end
% set up the regression coefficient list
betabasis = create_fourier_basis(yearRng, 11);
betafdPar = fdPar(betabasis);
betaCell = cell(p,1);
for j = 1:p
betaCell{j} = betafdPar;
end
% carry out the functional analysis of variance
fRegressStr = fRegress(temp36fd, regionCell, betaCell);
% extract the estimated regression coefficients and y-values
betaestCell = fRegressStr.betahat;
regionFit = fRegressStr.yhat;
% Figure 10.1
for j = 1:p
subplot(2,3,j)
plot(getfd(betaestCell{j}))
xlabel('')
ylabel(['\fontsize{13} ',region_names(j,:)])
end
subplot(2,3,6)
plot(regionFit)
xlabel('')
ylabel('\fontsize{13} Prediction')
% ------ 10.1.2 Trends in Sea Bird Populations on Kodiak Island ---------
% 10.1.2 Trends = Sea Bird Populations on Kodiak Island
seabirdPath = [examplesPath,'/seabirds'];
addpath(seabirdPath)
load seabirds
% select only the data for sites Uyak and Uganik, which have data
% from 1986 to 2005, except for 1998
sites = seabirds.Bay;
birdlabels = seabirds.birdlabels;
birdcountmat = seabirds.birdcountmat;
% year indices for which there are data.
selYear = [1:12, 14:20];
yearObs = selYear + 1985;
yearRng = [1986,2005];
% keep only rows for which there are data
birdcountmat = birdcountmat(selYear,:);
% Compute mean counts taken over both sites and transects
meanCounts = (birdcountmat(:, 1:13) + ...
birdcountmat(:,14:26))/2;
% Compute log (base 10) mean counts
logCounts = log10(meanCounts);
% Figure 10.2
shellfishindex = [1,2,5,6,12,13];
fishindex = [3,4,7,8,9,10,11];
logcountlim = [min(min(logCounts)), max(max(logCounts))];
meanShellfish = mean(meanCounts(:, shellfishindex), 2);
subplot(2,1,1)
plot(yearObs, logCounts(:,shellfishindex), 'bo-', ...
yearObs, log10(meanShellfish), 'b--', ...
yearRng, [0,0], 'b:')
xlabel('\fontsize{13} year')
ylabel('\fontsize{13} Log_{10} count')
title('\fontsize{13} Shellfish Diet')
axis([yearRng, logcountlim])
meanfish = mean(meanCounts(:,fishindex), 2);
subplot(2,1,2)
plot(yearObs, logCounts(:,fishindex), 'bo-', ...
yearObs, log10(meanfish), 'b--', ...
yearRng, [0,0], 'b:')
xlabel('\fontsize{13} year')
ylabel('\fontsize{13} Log_{10} count')
title('\fontsize{13} fish Diet')
axis([yearRng, logcountlim])
% Compute mean counts taken over transects only within sites
% so we have 2 observations for each bird species each year.
% Two of these counts are zero, and are replaced by 1/(2*n)
meanCounts2 = birdcountmat;
for j = 1:13
for k = 1:2
jcol = j + (k-1)*13;
n = sum(meanCounts2(:,jcol));
zeroind = find(meanCounts2(:,jcol) == 0);
if ~isempty(zeroind)
meanCounts2(zeroind,jcol) = 1/(2*n);
end
end
end
logCounts2 = log10(meanCounts2);
% Represent log mean counts exactly with a polygonal basis
birdbasis = create_polygonal_basis([1,20],selYear);
[birdfd2, df, gcv, SSE, penmat, y2cMap] = ...
smooth_basis(selYear, logCounts2, birdbasis);
% -----------------------------------------------------------------
% After some preliminary analyses we determined that there was no
% contribution from either site or food*site interaction.
% Now we use a reduced model with only a feed effect,
% but we add bird effects, which were seen = the plot to be
% strong. Birds are nested within feed groups, and either their
% effects must sum to zero within each group, or we must designate
% a bird = each group as a baseline, and provide dummy variables
% for the remainder. We opt for the latter strategy.
% -----------------------------------------------------------------
% The design matrix contains an intercept dummy variable, a
% feed dummy variable, and dummy variables for birds, excluding
% the second bird = each group, which turns out to be the each
% group's most abundant species, and which is designated as the
% baseline bird for that group.
Zmat0 = zeros(26,15);
% Intercept or baseline effect
Intercept = ones(26,1);
% Crustacean/Mollusc feeding effect: a contrast between the two groups
foodindex = [1,2,5,6,12,13];
fooddummy = zeros(26,1);
fooddummy([foodindex,foodindex+13]) = 1;
% Bird effect, one for each species
birddummy = diag(ones(13,1));
birdvarbl = [birddummy;birddummy];
% fill the columns of the design matrix
Zmat0(:,1) = Intercept;
Zmat0(:,2) = fooddummy;
Zmat0(:,3:15) = birdvarbl;
% Two extra dummy observations are added to the functional data
% object for log counts, and two additional rows are added to
% the design matrix to force the bird effects within each diet
% group to equal 0.
birdcoef3 = [getcoef(birdfd2), zeros(19,2)];
birdfd3 = fd(birdcoef3, birdbasis);
Zmat = [Zmat0; zeros(2,15)];
Zmat(28,fishindex+2) = 1;
Zmat(27,shellfishindex+2) = 1;
p = 15;
xfdcell = cell(p,1);
betacell = xfdcell;
for j = 1:p
xfdcell{j} = Zmat(:,j);
end
% set up the functional parameter object for the regression fns.
% use cubic b-spline basis for intercept and food coefficients
betabasis1 = create_bspline_basis([1,20],21,4,selYear);
lambda = 10;
betafdPar1 = fdPar(betabasis1,2,lambda);
betacell{1} = betafdPar1;
betacell{2} = betafdPar1;
betabasis2 = create_constant_basis([1,20]);
betafdPar2 = fdPar(betabasis2);
for j = 3:15
betacell{j} = betafdPar2;
end
birdRegressStr = fRegress(birdfd3, xfdcell, betacell);
betaestcell = birdRegressStr{4};
% Figure 10.3 is produced = Section 10.2.2 below
% after estimating the smoothing parameter = Section 10.1.3
%
% Here we plot the regression parameters
% without the confidence intervals.
subplot(2,1,1)
plot(getfd(betaestcell{1}))
subplot(2,1,2)
plot(getfd(betaestcell{2}))
%
% Section 10.1.3 Choosing Smoothing Parameters
%
% Choose the level of smoothing by minimizing cross-validated
% error sums of squares.
loglam = -2:0.25:4;
SSE_CV = zeros(length(loglam),1);
betafdPari = betafdPar1;
CV_obs = 1:26;
wt = ones(28,1);
for i = 1:length(loglam)
disp(loglam(i))
betafdPari = putlambda(betafdPari, 10^loglam(i));
betacelli = betacell;
for j = 1:2
betacelli{j} = betafdPari;
end
CVi = fRegress_CV(birdfd3, xfdcell, betacelli, wt, CV_obs);
SSE_CV(i) = CVi;
end
% Figure 10.4
phdl = plot(loglam, SSE_CV, 'bo-');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} log smoothing parameter')
ylabel('\fontsize{13} cross validated sum of squares')
% Cross-validation is minimized at something like lambda = sqrt(10),
% although the discontinous nature of the CV function is disquieting.
betafdPar1 = putlambda(betafdPar1, 10^0.5);
for j = 1:2
betacell{j} = betafdPar1;
end
% carry out the functional regression analysis
fRegressStr = fRegress(birdfd3, xfdcell, betacell);
% plot regression functions
betanames{1} = 'Intercept';
betanames{2} = 'Food Effect';
birdBetaestcell = fRegressStr{4};
for j = 1:2
subplot(2,1,j)
betaestParfdj = birdBetaestcell{j};
betaestfdj = getfd(betaestParfdj);
betaestvecj = eval_fd(selYear, betaestfdj);
plot(selYear, betaestvecj, 'b-')
xlabel('\fontsize{13} Year')
ylabel(['\fontsize{13} ',betanames{j}])
end
% plot predicted functions
yhatfdobj = fRegressStr{5};
subplot(1,1,1)
plotfit_fd(logCounts2, selYear, yhatfdobj(1:26))
%
% Section 10.2 Functional Responses with Functional Predictors:
% The Concurrent Model
%
% Section 10.2.2 Confidence Intervals for Regression Functions
birdYhatmat = eval_fd(selYear, yhatfdobj(1:26));
rmatb = logCounts2 - birdYhatmat;
SigmaEb = cov(rmatb');
birdbetastderrCell = fRegress_stderr(fRegressStr, y2cMap, SigmaEb);
plotbeta(birdBetaestcell(1:2), birdbetastderrCell(1:2))
% Section 10.2.3 Knee Angle Predicted from Hip Angle
% path to the gait data example folder
gaitPath = [examplesPath,'/gait'];
addpath(gaitPath)
% load the data from gait.mat
load gait
hip = gait.hip;
knee = gait.knee;
gaitarray = gait.gaitarray;
gaittime = 0.5:1:19.5;
gaitrange = [0,20];
gaitfine = linspace(0,20,101)';
harmaccelLfd20 = vec2Lfd([0, (2*pi/20).^2, 0], gaitrange);
nbasis = 21;
gaitbasis = create_fourier_basis(gaitrange, nbasis);
gaitLoglam = -4:0.25:0;
nglam = length(gaitLoglam);
% First select smoothing for the raw data
gaitSmoothStats = zeros(nglam, 3);
gaitSmoothStats(:,1) = gaitLoglam;
% loop through smoothing parameters
for ilam = 1:nglam
lambdai = 10^gaitLoglam(ilam);
gaitfdPari = fdPar(gaitbasis, harmaccelLfd20, lambdai);
[gaitfdi, dfi, gcvi] = smooth_basis(gaittime, gaitarray, gaitfdPari);
gaitSmoothStats(ilam, 2) = dfi;
gaitSmoothStats(ilam, 3) = sum(sum(gcvi));
end
% display and plot GCV criterion and degrees of freedom
disp(gaitSmoothStats)
subplot(1,1,1)
plot(gaitSmoothStats(:,1), gaitSmoothStats(:,3), 'bo-')
xlabel('\fontsize{13} log_{10} lambda')
ylabel('\fontsize{13} GCV')
% two panel display of df and GCV
subplot(2,1,1)
plot(gaitSmoothStats(:,1), gaitSmoothStats(:,3), 'bo-')
xlabel('\fontsize{13} log_{10} lambda')
ylabel('\fontsize{13} GCV')
subplot(2,1,2)
plot(gaitSmoothStats(:,1), gaitSmoothStats(:,2), 'bo-')
xlabel('\fontsize{13} log_{10} lambda')
ylabel('\fontsize{13} degrees of freedom')
% GCV is minimized with lambda = 10^(-1.5).
gaitfdPar = fdPar(gaitbasis, harmaccelLfd20, 10^(-1.5));
[gaitfd, df, gcv, SSE, penmat, y2cMap] = ...
smooth_basis(gaittime, gaitarray, gaitfdPar);
gaitfdnames{1} = 'Normalized time';
gaitfdnames{2} = 'Child';
gaitfdnames{3} = cell(2,1);
gaitfdnames{3,1} = 'Angle';
gaitfdnames{3,2} = ['Hip ', 'Knee'];
gaitfd = putnames(gaitfd, gaitfdnames);
hipfd = gaitfd(:,1);
kneefd = gaitfd(:,2);
% Figure 10.5
% compute mean knee function and its first two derivatives
kneefdMean = mean(kneefd);
D1kneevecMean = eval_fd(gaitfine, kneefdMean, 1);
D2kneevecMean = eval_fd(gaitfine, kneefdMean, 2);
% plot the mean knee function and its first two derivatives in a
% 3-panel plot
% plot mean knee function
subplot(3,1,1)
kneeMeanvec = eval_fd(gaitfine, kneefdMean);
plot(gaitfine, kneeMeanvec, 'b-', ...
[7.5, 7.5], [0, 80], 'b--', [14.7, 14.7], [0, 80], 'b--')
xlabel('')
title('\fontsize{13} Mean Knee Angle')
% plot first derivative of mean knee function
subplot(3,1,2)
D1kneeRng = [min(D1kneevecMean), max(D1kneevecMean)];
plot(gaitfine, D1kneevecMean, 'b-', [0,20], [0,0], 'b:', ...
[7.5, 7.5], D1kneeRng, 'b--', [14.7, 14.7], D1kneeRng, 'b--')
xlabel('')
title('\fontsize{13} Mean Knee Angle Velocity')
% plot second derivative of mean knee function
subplot(3,1,3)
D2kneeRng = [min(D2kneevecMean), max(D2kneevecMean)];
plot(gaitfine, D2kneevecMean, 'b-', [0,20], [0,0], 'b:', ...
[7.5, 7.5], D2kneeRng, 'b--', [14.7, 14.7], D2kneeRng, 'b--')
xlabel('')
title('\fontsize{13} Mean Knee Angle Acceleration')
% Figure 10.6
% plot a phase/plane diagram for the mean knee function
% values of the first and second derivative at observation times
D1vec = eval_fd(gaittime, kneefdMean, 1);
D2vec = eval_fd(gaittime, kneefdMean, 2);
% labels for these points
plotlab = cell(20,1);
for i=1:20, plotlab{i} = [' ',num2str(i)]; end
% plot the phase/plane diagram using function phaseplanePlot
subplot(1,1,1)
[phdl, thdl] = phaseplanePlot(kneefdMean, 1, 2, plotlab, gaittime);
set(phdl, 'LineWidth', 2)
set(thdl, 'fontsize', 12)
xlabel('\fontsize{13} Knee Velocity')
ylabel('\fontsize{13} Knee Acceleration')
axis([-25,25,-15,15])
% Set up a functional linear regression to predict knee angle
% from hip angle
% It is essential to remove previous cell arrays with the same names
% define new covariate cell array with the first covariate being
% the constant function and the second being hip angle
hipxfdCell{1} = ones(39,1);
hipxfdCell{2} = hipfd;
% define regression function cell array kneebetaCell
kneebetafdPar = fdPar(gaitbasis, harmaccelLfd20);
kneebetaCell = cell(2,1);
kneebetaCell{1} = kneebetafdPar;
kneebetaCell{2} = kneebetafdPar;
% the regression analysis
kneefRegressStr = fRegress(kneefd, hipxfdCell, kneebetaCell);
% Figure 10.7
% cell array for estimated regression coefficients
kneebetaestCell = kneefRegressStr.betahat;
% evaluate regression functions
kneeIntercept = eval_fd(gaitfine, getfd(kneebetaestCell{1}));
kneeHipCoef = eval_fd(gaitfine, getfd(kneebetaestCell{2}));
% Plot intercept along with the mean knee function
subplot(2,1,1)
phdl=plot(gaitfine, kneeIntercept, 'b-', gaitfine, kneeMeanvec, 'b--', ...
[0,20], [0,0], 'b:', ...
[7.5,7.5], [0,80], 'b--', [14.7,14.7], [0,80], 'b--');
set(phdl, 'LineWidth', 2)
xlabel('')
ylabel('')
title('\fontsize{16} Intercept and Mean Knee Angle')
% Plot hip coefficient
ylim = [0, max(kneeHipCoef)];
subplot(2,1,2)
phdl=plot(gaitfine, kneeHipCoef, 'b-', [0,20], [0,0], 'b:', ...
[7.5,7.5], ylim, 'b--', [14.7,14.7], ylim, 'b--');
set(phdl, 'LineWidth', 2)
xlabel('')
ylabel('')
title('\fontsize{16} Hip Coefficient')
% Squared multiple correlation
% fitted knee functions
kneehatfd = kneefRegressStr.yhat;
% matrix of residual values computed at observation times
kneehatmat = eval_fd(gaitfine, kneehatfd);
kneemat = eval_fd(gaitfine, kneefd);
resmat1 = kneemat - kneehatmat;
% matrix of residuals for fitting the data with the mean knee angle
kneemeanvec = eval_fd(gaitfine, mean(kneefd));
resmat0 = kneemat - kneemeanvec * ones(1,39);
SSE0 = sum(resmat0.^2, 2);
SSE1 = sum(resmat1.^2, 2);
RSQRD = (SSE0-SSE1)./SSE0;
% Plot Hip Coefficient & Squared Multiple Correlation
subplot(1,1,1)
phdl=plot(gaitfine, kneeHipCoef, 'b-', gaitfine, RSQRD, 'b--', ...
[7.5,7.5], ylim, 'b--', [14.7,14.7], ylim, 'b--');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} ')
ylabel('\fontsize{13} ')
title('\fontsize{16} Hip Coefficient and Squared Multiple Correlation')
legend('\fontsize{12} Hip Coef.', '\fontsize{12} R^2')
% Figure 10.8
% covariance matrix for residuals
resmat = eval_fd(gaittime, kneefd) - eval_fd(gaittime, kneehatfd);
SigmaE = cov(resmat');
% compute cell array of standard error functions
kneebetastderrCell = fRegress_stderr(kneefRegressStr, y2cMap, SigmaE);
% plot regression functions with 95% confidence intervals
plotbeta(kneebetaestCell, kneebetastderrCell)
% Figure 10.9
% predict knee acceleration from hip acceleration
clear hipxfdcell % important to clear earlier version of cell array
hipxfdCell2 = cell(2,1);
hipxfdCell2{1} = ones(39,1);
hipxfdCell2{2} = deriv(hipfd, 2);
knee_accelfd = deriv(kneefd, 2);
kneefRegressStr2 = fRegress(knee_accelfd, hipxfdCell2, kneebetaCell);
% cell array for estimated regression coefficients
kneebetaestCell2 = kneefRegressStr2.betahat;
% evaluate regression functions
kneeIntercept2 = eval_fd(gaitfine, getfd(kneebetaestCell2{1}));
kneeHipCoef2 = eval_fd(gaitfine, getfd(kneebetaestCell2{2}));
% Squared multiple correlation
% fitted knee functions
knee_accelhatfd = kneefRegressStr2.yhat;
% matrix of residual values computed at observation times
knee_accelmat = eval_fd(gaitfine, knee_accelfd);
resmat1 = knee_accelmat - eval_fd(gaitfine, knee_accelhatfd);
% matrix of residuals for fitting the data with the mean knee angle
knee_accelmeanvec = eval_fd(gaitfine, mean(knee_accelfd), 2);
resmat0 = knee_accelmat - knee_accelmeanvec * ones(1,39);
SSE0 = sum(resmat0.^2, 2);
SSE1 = sum(resmat1.^2, 2);
RSQRD = (SSE0-SSE1)./SSE0;
ylim=[0, max(kneeHipCoef2)];
phdl = plot(gaitfine, kneeHipCoef2, 'b-', ...
[7.5,7.5], ylim, 'b--', [14.7,14.7], ylim, 'b--');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} ')
ylabel('\fontsize{13} ')
title('\fontsize{13} Hip acceleration')
% Squared multiple correlation
phdl = plot(gaitfine, RSQRD, 'b-', ...
[7.5,7.5], [0,1], 'b--', [14.7,14.7], [0,1], 'b--');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} ')
ylabel('\fontsize{13} ')
title('\fontsize{13} R^2')
%
% Section 10.3 Beyond the Concurrent Model
%
% (no computations = this section)
%
% Section 10.4 A Functional Linear Model for Swedish Mortality
%
% SwedeLogHazardPath = [examplesPath,'/SwedeLogHazard'];
%
% addpath(SwedeLogHazardPath)
%
% load SwedeLogHazard
% load SwedeLogHazard1914
% The Swedish
% Figure 10.10
% indices of years 1751, 1810 and 1860
threeyears = [1751, 1810, 1860] - 1750;
tempLogHazard = [SwedeLogHazard(:,threeyears), SwedeLogHazard1914];
SwedeTime = (0:80)';
phdl = plot(SwedeTime, tempLogHazard, '-');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} age')
ylabel('\fontsize{13} log hazard')
legend('\fontsize{12} 1751', '\fontsize{12} 1810', ...
'\fontsize{12} 1860', '\fontsize{12} 1914', ...
'location', 'north')
% plot the data for 1751 to 1884
plot(SwedeLogHazard)
% set up a saturated order six basis for smoothing the data
SwedeRng = [0,80];
nbasis = 85;
norder = 6;
SwedeBasis = create_bspline_basis(SwedeRng,nbasis,norder);
% set the smoothing parameters
lambda = 1e-7;
SwedefdPar = fdPar(SwedeBasis, 4, lambda);
% smooth the data
Swedefd = smooth_basis(SwedeTime,SwedeLogHazard,SwedefdPar);
plotfit_fd(SwedeLogHazard, SwedeTime, Swedefd)
% set up coefficient functions
nbasis = 23;
SwedeBetaBasis = create_bspline_basis(SwedeRng, nbasis);
SwedeBeta0fdPar = fdPar(SwedeBetaBasis, 2, 1e-5);
SwedeBeta1fd = bifd(zeros(23), SwedeBetaBasis, SwedeBetaBasis);
SwedeBeta1fdPar = bifdPar(SwedeBeta1fd, 2, 2, 1e3, 1e3);
p = 2;
SwedeBetaCell = cell(p, 1);
SwedeBetaCell{1} = SwedeBeta0fdPar;
SwedeBetaCell{2} = SwedeBeta1fdPar;
% Define the dependent and independent variables
NextYear = Swedefd(2:144);
LastYear = Swedefd(1:143);
Swede_linmodStr = linmod(NextYear, LastYear, SwedeBetaCell);
Swede_beta1fd = eval_bifd(SwedeTime, SwedeTime, ...
Swede_linmodStr.beta);
% Figure 10.11
surfc(SwedeTime, SwedeTime, Swede_beta1fd)
xlabel('\fontsize{13} age')
ylabel('\fontsize{13} age')
zlabel('\fontsize{13} \beta(s,t)')
%
% Section 10.5 Permutation Tests of Functional Hypotheses
%
% Section 10.5.1 Functional t-Tests
% --------------- Comparing male and female growth data ----------------
% path to the Berkeley growth data example folder
growthPath = [examplesPath,'/growth'];
addpath(growthPath)
% load the growth data struct object from the .mat file growth.mat
load growth
% Figure 10.12
phdl = plot(age, hgtmmat(:, 1:10), 'b--', ...
age, hgtfmat(:, 1:10), 'b-');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} age')
ylabel('\fontsize{13} height (cm)')
axis([1,18, 60,200])
nbasis = 35;
norder = 6;
growthbasis = create_bspline_basis([1,18], nbasis, norder, age);
growfdPar = fdPar(growthbasis, 3, 10^(-0.5));
hgtffd = smooth_basis(age,hgtfmat,growfdPar);
hgtmfd = smooth_basis(age,hgtmmat,growfdPar);
tresStr = tperm_fd(hgtffd, hgtmfd);
% Figure 10.13
% Section 10.5.2 Functional F-Tests
% --------------- Testing for no effect of climate zone ----------------
% temp36fd, regionCell, betaCell from Section 10.1.1 above
F_res = Fperm_fd(temp36fd, regionCell, betaCell);
% Figure 10.14
qval = 0.95;
phdl = plot(argvals, F_res.Fvals, 'b-', ...
argvals, F_res.qvals_pts, 'b--', ...
[min(argvals), max(argvals)], [qval,qval], 'b:');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} day')
ylabel('\fontsize{13} F-statistic')
legend('\fontsize{13} Observed Statistic', ...
['\fontsize{13} pointwise ', 1-q, ' critical value'], ...
['\fontsize{13} maximum ', 1-q, ' critical value'])
%
% 10.6 Details for R Functions fRegress, fRegress_CV and fRegress_stderr
%
help fRegress
help fRegress_CV
help fRegress_stderr
%
% 10.7 Details for Function plotbeta
%
help plotbeta
%
% 10.8 Details for Function linmod
%
help linmod
%
% 10.9 Details for Functions Fperm_fd and tperm_fd
%
help(Fperm_fd)
help(tperm_fd)
%
% Section 10.10 Some Things to Try
%
% (exercises for the reader)
%
% Section 10.11 More to Read
%
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...