Is the population of blue-eyed Martians decreasing? Suppose that we want to test the hypothesis that the proportion of blue-eyed martians has been decreasing throughout the 20th century. Unfortunately, the Martian population fluctuates highly, so every decade there is large difference in the total population [update: consider the Martian population to be constant at one billion Martians. The data below are random samples at each year]. The data set (made up as I'm writing this) could look something like this:
Year | Total martian population | Blue-eyed martians | Proportion
1910 | 400                      | 250                | 0.625
1920 | 2000                     | 1000               | 0.500
1930 | 70                       | 40                 | 0.571
1940 | 30                       | 14                 | 0.467
1950 | 10                       | 4                  | 0.400
1960 | 140                      | 52                 | 0.371
1970 | 50 000                   | 15 400             | 0.308
1980 | 70 000                   | 22 000             | 0.314
1990 | 1500                     | 80                 | 0.053
2000 | 5000                     | 800                | 0.160

Analyzing years when the Martian population is under 100 is clearly not as meaningful statistically as when the population is over 10 000, since in the latter case we have a larger data set. Still, we would like to use all available data to verify our hypothesis with a conventional 95% significance level.
How do we proceed? Do we weight the importance of each year according to the size of the sample at the time?
More edit to fit concerns: the concern here is how do we adequately weight each dataset keeping in mind that they are of such different sizes. There is no sample bias as data is randomly selected.
 A: This answer describes three ways to handle the varying sample sizes appropriately: a Generalized Linear Model and two weighted Ordinary Least Squares regressions.  In this case all three work well.  In general, when some proportions are near $0$ or $1$, the GLM is better.

Because the sample sizes are so small compared to the populations (less than ten percent of them), to an excellent approximation the distribution of blue-eyed and non-blue-eyed results in a sample of size $n$ is Binomial (because the samples are random).  The other Binomial parameter, $p$, is the true (but unknown) proportion of blue-eyed subjects in the population.  Thus, the chance of observing $k$ blue-eyed people is
$$\binom{n}{k}p^k(1-p)^{n-k}.\tag{1}$$
Each decade we know $n$ and $k$--those are given by the data--but we don't know $p$.  We may estimate it by assuming that the log odds corresponding to $p$ varies by year linearly (at least to a good approximation).  This means we assume there are numbers $\beta_0$ and $\beta_1$ such that
$$\log(p) - \log(1-p) = \beta_0 + \beta_1 \times \text{Year}.$$
Equivalently,
$$p = \frac{1}{1 + e^{-\beta_0-\beta_1\text{Year}}};\ 1-p = \frac{ e^{-\beta_0-\beta_1\text{Year}}}{1 +  e^{-\beta_0-\beta_1\text{Year}}}.$$
Plugging this into (1) gives the chance of observing $k$ out of $n$ during a given year $t$ as
$$\binom{n}{k} \frac{e^{-(\beta_0+\beta_1t)(n-k)}}{\left(1 + e^{-(\beta_0+\beta_1t)}\right)^n}.\tag{2}$$
Assuming the samples are independently obtained at years $t_1, t_2,$ etc and writing the corresponding sample sizes and counts of blue-eyed subjects as $n_i$ and $k_i$, the probability of the data is the product of the probabilities of the individual results.  This product is (by definition) the likelihood of $(\beta_0, \beta_1)$.  We may estimate these parameters as the values $(\hat\beta_0, \hat\beta_1)$ that maximize the likelihood; equivalently, they maximize the log likelihood
$$\Lambda(\beta_0,\beta_1) = \sum_t \log\left(\binom{n}{k} \frac{e^{-(\beta_0+\beta_1t)(n-k)}}{\left(1 + e^{-(\beta_0+\beta_1t)}\right)^n}\right)\tag{3}$$
obtained from $(2)$. 
(This simplifies considerably, using rules of logarithms, which is one reason to express the time-proportion relationship in terms of log odds.  When all proportions are between $0.2$ and $0.8$, approximately, there is little qualitative difference between using probabilities $p$ or their log odds: the fitted curve will be linear or close to linear, respectively.)
$(3)$ is a Binomial Generalized Linear Model.  It must be fitted by numerically minimizing $\Lambda$.  The glm procedure in R (shown at the end of this post) gives the solution
$$(\hat\beta_0, \hat\beta_1)_\text{GLM}  = (31.498711, -0.0163568).$$

The data in this figure are plotted with disks whose areas are proportional to the sample sizes.  The GLM fit is curvilinear.  Shown for comparison, in gray, is the line we would get just by dumping the $(\text{Year},\text{Proportion})$ data shown in the question into an Ordinary Least Squares solver.  Both fits are influenced by the greater proportions in earlier years, despite the small sample sizes then.  However, the GLM fit does a better job of approximating the proportions in the largest samples obtained in 1970 and 1980.  The dotted blue line is described below.
By adding a quadratic term we can test the goodness of fit.  It significantly improves the GLM fit (although visually the difference is not great), providing evidence that this model does not describe the variation in results well.  Looking at the plot indicates the result in 1990 was much lower than the model predicts.

An alternative, but comparable, approach is to estimate $p$ individually for each year $t_i$, perhaps as $k_i / n_i$ (although other estimators are possible).  A linear regression of the log odds of these estimates against the year, weighted by the sample sizes $n_i$, or Weighted Least Squares regression, yields
$$(\hat\beta_0, \hat\beta_1)_\text{WLS} = (36.12744, -0.018706).$$
The standard errors of these estimates are $15.55$ and $0.00787$, respectively, indicating that the WLS estimates are not significantly different from the Binomial GLM.  (The GLM's standard errors are considerably smaller, though: it "knows" these sample sizes are pretty large whereas the linear regression "knows" nothing about the sample sizes at all: it only has a sequence of ten separate observations.)  Note that this alternative might not be available if $k_i=n_i$ or $k_i=0$, unless a different estimator of the probabilities is used (which doesn't produce values of $0$ or $1$).
Finally, we might simply perform a weighted least squares regression of the raw probability estimates $k/n$ against the year, inversely weighted by an estimate of sample variance.  The variance of a Binomial$(n,p)$ variable $X$, re-expressed as a proportion $X/n$ is $p(1-p)/n$.  That may be estimated from a sample as
$$p(1-p)n \approx \frac{k}{n}\frac{n-k}{n}/n = \frac{k(n-k)}{n^3}.$$
Its result appears in the figure as a dotted blue line.  In this case it appears to compromise between the GLM and OLS fits.

The following R code performed the analyses and produced the figure.
year <- seq(1910, 2000, by=10)
total <- c(40, 200, 7, 3, 1, 14, 5000, 7000, 150, 500) * 10
blue <- c(250, 1000, 40, 14, 4, 52, 15400, 22000, 80, 800)
X <- data.frame(Year=year, Success=blue, Failure=total-blue,
                Proportion=blue/total, Total=total)
#
# GLM
#
fit <- glm(cbind(Success, Failure) ~ Year, X, family="binomial")
summary(fit)
#
# WLS of the log odds (an alternative)
#
fit.WLS <- lm(log(Success/Failure) ~ Year, X, weights=Total)
summary(fit.WLS)
#
# Plot the results.
#
X.more <- data.frame(Year=1901:2010)
X.more$Prediction <- predict(fit, X.more, type="response")
plot(X$Year, X$Proportion, ylim=0:1,
     type="p", pch=21, bg="Red", cex=sqrt(X$Total/2000),
     xlab="Year", ylab="Proportion",
     main="GLM and OLS Fits", sub="GLM: solid line; OLS: dotted line")
lines(X.more, lwd=2)
abline(lm(Proportion ~ Year, X), 
       lty=3, lwd=3, col="Gray") #The OLS fit
abline(lm(Proportion ~ Year, X, weights=Total^3/(Success*Failure)), 
       lty=3, lwd=3, col="Blue") #The weighted OLS fit to the proportions

