What statistic should I use when I have two devices measuring 20 different samples over three days? I have two measuring devices that I use to measure 20 samples over the course of three days. Each device produce 3 measurements per day. As a result, I have nine values per device per sample.
Please let me know what statistical analysis I should use if I want to "prove" that there is no statistically significant difference between these two measuring devices.
EDIT
To be clear -- each "sample" (specimen) is measured three times each day (in the morning, at noon, and in the evening).
EDIT2
All measurements are continuous. The following is an example of how the data look like from two samples: (and I am interested in finding out if both devices are equivalent)

 A: Design of experiment: repeat the same time point.
A problem with comparing these time series is that you have variations between time points that may not be due to the random variation within the measurement device but due to random changes in the samples.
This makes it difficult to assess whether a variation in the measurements is due to the device or due to something else(like the sample).
For this reason it is better to have the measurements repeated for single time points, in which case you know that the variation in the sample is minimal. In that case you could ignore the variations in the sample and compare the measurements as if they come from a single distribution. For instance apply a Kolmogornov-Smirnov test to compare any difference a t-test to compare difference in means or an F-test to compare difference in variance (there are many more variants and they depend on how you want to compare the difference and what assumptions you can make about the distribution).

Paired measurements
In your case you have paired measurements and you can use that to still make some comparisons (although I believe your life would have been easier with repetitions of the same time point).
This is especially easy when you are only looking to compare whether the devices have a difference in their means. In order to test this you can use the difference in the measurements and test whether it's mean is zero. This is a paired difference test. The exact computation depends on the assumption about the distribution. It could be a z-test, t-test or a non-parametric test like a sign test or rank test.
To make other comparisons, like for instance whether the devices have a different variance, you could apply an F-test, but potentially you might make used of the paired data as well and consider a correlation between the measurements between the device 1 and device 2. The variations of the measurements of a sample around it's mean will be correlated for the two devices when there is a variation in the sample (it is correlated for the two devices because they measure the same sample, the same source of variation). I imagine that this should increase the power of the statistical test (less likely to make a type II error if the two devices are different). As a statistical test one could model the data with equal variance and model the data with different variance, then use a likelihood ratio test.(I am not sure whether there is a standard description, method or software for such comparison)
Example with code
The example below shows how you can model this as a mixed effects model in R. Here the package nlme is a bit easier to model the variable variance as a function of the device, on the other hand lme4 is a bit easier to model crossed mixed effects (it might be that you whish to not model crossed mixed effects).
Currently the code just performs the model and I wished to still extend it (I did it on my phone in an online R terminal but I need to do it on my computer for the next steps). What might be instructive is the step
y = z_sample[sample] + z_time[time] + z_id[id] + error

which formulates the model and describes the variable $y$ as a sum of factors that depend on the 'sample' the 'time' (and this makes that observations in the same sample or same time correlate) as well as the 'id' which is an interaction term, for each combination of sample and time a draw is made from a random normal distribution. Finally there is a noise term and this is allowed to be different for the two devices, here we used sd = device which means that device numbered 1 will have noise with a standard deviation of 1 and device numbered 2 will have noise with a standard deviation of 2.
The first three terms z_sample[sample], z_time[time] and z_id[id] are the same for both devices. This means that the paired measurements will correlate, and that is something that regular hypothesis test (e.g. F-test for testing equivalence of variance) do not take into account.
What I still need to add is graphs that help to see the relationships. Also I wanted to show how to test for the equivalence of variance and difference in means by using a likelihood ratio test with this mixed effects model. Then show with a simulation how the incorporation of the mixed effects will be more powerful and not overly conservative (ignoring the variance of the sample, time and id, will make the computed, nominal, p-value higher than the true p-value really is).
### settings
set.seed(1)
library(nlme)
library (lme4)
ns = 20
nt = 9

### generate data

### regressors 
device = c(rep(1,nt*ns),rep(2,nt*ns))
time = rep(1:nt,ns*2)
sample = rep(rep(1:ns, each = nt),2) 
id = rep(1:(ns*nt),2)

### dummy variables to trick lme in modeling different variance for device 1 and device 2 
### this is done with the term
### (0+d1+d2|id)
### it models a slope effect for d1 and d2
### because thes are dummy variables, only value 0 or 1, these are two intercept terms depending on the device and wil be able to have different variance

d1 = ifelse(device == 1,1,0)
d2 = ifelse(device == 2,1,0)

### vectors with the random noise and random effects
z_sample = rnorm(ns, mean =0, sd = 0.5)
z_time = rnorm(nt, mean = 0, sd = 1)
z_id = rnorm(nt*ns, mean = 0, sd = 1.5)
error = rnorm(nt*ns*2, mean = 0, sd = device)

### an observation as a sum of the noise terms
y = z_sample[sample] + z_time[time] + z_id[id] + error

### turn the numerical indices into a factor (such that the model will treat them as a factor instead of a numerical variable)
sample = as.factor(sample)
time = as.factor(time)
device = as.factor(device)
id = as.factor(id)
dummy = as.factor(rep(1,nt*ns*2))

### model with lme4
lmer(y ~ 1 + (1 | time) + (1 | sample) + (0+d1+d2|id), control = lmerControl(check.nobs.vs.nRE= "ignore"))

### model with nlmer
lme(y ~ 1, random = list(dummy = pdBlocked(list(pdIdent(~0+sample),
               pdIdent(~0+time), 
               pdIdent(~0+id)))), weights = varIdent(form = ~ 1|device))

