Simple linear regression:

X={1.0, 1.0, 3.3, 3.3, 4.0, 4.0, 4.0, 5.6, 5.6, 5.6, 6.0, 6.0, 6.5, 6.5};
Y = {1.6, 1.8, 1.8, 2.7, 2.6, 2.6, 2.2, 3.5, 2.8, 2.1, 3.4, 3.2, 3.4, 3.9};
Data=Transpose[{X,Y}];
model = LinearModelFit[Data,d,d]

SSE = Total[model["FitResiduals"]^2]
Ssquare = model["EstimatedVariance"]
Rsquare = model["RSquared"]
paraconf = model["ParameterConfidenceIntervals", ConfidenceLevel->0.95]
model["ParameterTable"]
conf=model["MeanPredictionBands",ConfidenceLevel->0.95]
pred=model["SinglePredictionBands",ConfidenceLevel->0.95]

s = 3.5;
{1.2396331223356167` +0.3244412866554706` s-2.1788128296672284` √(0.09200313867789942` -0.03561853737126769` s+0.00399566925639221` (s^2)),1.2396331223356167` +0.3244412866554706` s+2.1788128296672284` √(0.09200313867789942` -0.03561853737126769` s+0.00399566925639221` (s^2))}

Hand Calculations:

n = Length[X];
sumx = Mean[X]*n
sumy = Mean[Y]*n
sumx2 = Mean[X^2]*n
sumy2 = Mean[Y^2]*n
sumxy = Mean[X*Y]*n

sxx = 1/n*(n*sumx2-sumx^2)
syy = 1/n*(n*sumy2-sumy^2)
sxy = 1/n*(n*sumxy-sumx*sumy)
R2 = sxy^2/(sxx*syy)

b1 = sxy/sxx
b0 = Mean[Y]-b1*Mean[X]

sse = syy  - b1*sxy
ssepe = sumy2 - ((1.6+1.8)^2/2+(1.8+2.7)^2/2+(2.6+2.6+2.2)^2/3+(3.5+2.8+2.1)^2/3+(3.4+3.2)^2/2+(3.4+3.9)^2/2)
sself = sse - ssepe

Lack of Fit

k = 6;
n = 14;
a2 = 0.05;
Fkn = (sself/(k-2))/(ssepe/(n-k))
InverseCDF[FRatioDistribution[k-2, n-k],1-a2]