title: Correction due to rounding error

text: The height for 1000 students is approximately normal with a mean 174.5cm and a standard deviation of 6.9cm. If 200 random samples of size 25 are chosen from this population and the values of the mean are recorded to the nearest integer, determine the probability that the mean height for the students is more than 176cm. Since the samples were rounded to the nearest integer, I should find $P(X>176.5)$ instead of $P(X>176)$. Is this how we account for the effect of rounding the observations? EDIT: In light of whuber's answer: The answer given by my module (no workings were provided): $\hspace{1cm} n=25; Normal$ $\hspace{1cm} \mu_{\overline{x}}=174.5cm$ $\hspace{1cm} \sigma_{\overline{x}}=6.9/5=1.38$ The answer is 0.1379. Which I'm pretty sure was found using $1-\phi(\dfrac{176-174.5}{1.38})$ So,   1. Is this an acceptable answer?   2. Since $n$ was less than 30, would it be ok to find the probability using a t-distribution? Thank you.


title: Explanation of finite correction factor

text: I understand that when sampling from a finite population and our sample size is more than 5% of the population, we need to a correction on the sample's mean and standard error using this formula: $\hspace{10mm} FPC=\sqrt{\frac{N-n}{N-1}}$ Where N is the population size and n is the sample size. I have 3 questions about this formula:   1. Why is the threshold set at 5%?   2. How was the formula derived?   3. Are there other online resources that comprehensively explain this formula besides this paper? Thank you,


title: Finding independent "clusters" in a matrix

text: I've called my question "clustering" but I am not sure if that's the right term. Imagine my matrix looks like this:               [ 0.  ,  0.92,  0.  ,  0.85,  0.  ]     [ 0.92,  0.  ,  0.  ,  0.89,  0.  ]     [ 0.85,  0.  ,  0.  ,  0.89,  0.  ]     [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ]     [ 0.  ,  0.89,  0.  ,  0.89,  0.  ]      What I am after is not a single clustering/dendrogram of indices, I am after extracting largest possible rectangles from the matrix (by re-arranging X/Y indices). For example, from the above I am expecting to get:               [ 0.92,  0.89,  0.  ,  0.  ,  0.  ]     [ 0.85,  0.89,  0.  ,  0.  ,  0.  ]     [ 0.  ,  0.85,  0.92,  0.  ,  0.  ]     [ 0.  ,  0.89,  0.89,  0.  ,  0.  ]     [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ]      However, I expect each cluster to be identified, even if it is not possible to show them all in a single matrix. What would be an algorithm I should be after? I've tried various clustering algorithms and they give me best single clustering, not focusing on the best individuals (which are mutually incompatible).


title: Choosing the right threshold for a biometric trait authentication system

text: I have a biometric authentication system that is using a person's gait to authenticate them. I extract features from gait, run it through a comparison versus a template and produce a similarity score (where if this similarity score is below a certain threshold, then the user is authenticated). So, I have 72 trials total (36 trials containing a positive case and 36 that contain a negative case). What I want to do is graph the ability of this system to authenticate people by illustrating it with a ROC graph. Unfortunately, I don't quite understand how to choose a threshold. Is there some mathematical procedure involved for choosing a threshold for the similarity scores? Do I just choose a bunch of different thresholds, and graph the corresponding ROC curves for all these different threshold values? The resulting similarity scores vary from [0.6,1.2] where the positive cases tend to lie around 0.6. All my coding is being done in Matlab.


title: Missing values in survival analysis

text: I am conducting a study on a cohort of people with a follow-up period of 7 years. I wish to use Cox Proportional Hazard model to estimate HR between an exposure and the length of time of an event. One missing information is the date of birth for the all subjects, but month and year are available.This prevents the calculation of exact age at the time of the study. Any suggestions will be much appreciated? Any sensitivity analysis should be conducted? Thanks


title: What statistical blogs would you recommend?

text: What statistical research blogs would you recommend, and why?


title: What are some good frameworks for method selection?

text: I have been looking into theoretical frameworks for method selection (note: not model selection) and have found very little systematic, mathematically- motivated work. By 'method selection', I mean a framework for distinguishing the appropriate (or better, optimal) method with respect to a problem, or problem type. What I have found is substantial, if piecemeal, work on particular methods and their tuning (i.e. prior selection in Bayesian methods), and method selection via bias selection (e.g. Inductive Policy: The Pragmatics of Bias Selection). I may be unrealistic at this early stage of machine learning's development, but I was hoping to find something like what measurement theory does in prescribing admissible transformations and tests by scale type, only writ large in the arena of learning problems. Any suggestions?


title: Are the conclusions we draw from comparing two distributions correct?

text: We want to compare two distributions of ages (birth years) of individuals. Given a set of individuals ( _all_ ) and a subset of that set ( _subset_ ), we want to find out:   1. Is it valid to compare the age distribution of _all_ with that of _subset_ (when | _subset_ | is much smaller than | _all_ |) or can we only compare _all_ minus _subset_?    2. We did a Shapiro-Wilk test on both distributions and got W = 0.7456 (p-value = 7.499e-08) for _all_ and W = 0.7467 (p-value = 7.865e-08) for _subset_. This confirms that both distributions are not normal distributions, rather exponential or power-law. Therefore, we can not use Pearson correlation to compare them but should rather use Spearman?   3. Comparing both distributions with Spearman, R says              t = 93.9151, df = 47, p-value < 2.2e-16     alternative hypothesis: true correlation is not equal to 0      95 percent confidence interval:      0.9952748 0.9985103      sample estimates:           cor      0.9973462      which shows a high correlation. How meaningful is this result, given the large amount of young people?   4. Can we infer from the correlation test that both distributions are (very) similar?   5. If so, does this allow us to conclude that _subset_ is not more biased towards younger people than _all_ is?   6. Finally, if our assumption was that _subset_ contains more younger people than one would expect, can we now infer that this assumption is wrong? Since I can not upload images, yet: here is a plot of both distributions and here the raw data:               # year  age(all)  age(subset)     1938  368  1     1939  360  1     1941  394  1     1942  809  3     1943  964  1     1944  686  1     1945  701  1     1946  1301  1     1947  1228  5     1948  1565  2     1949  2019  1     1950  2146  3     1951  2202  6     1952  2343  5     1953  2061  7     1954  2313  3     1955  2963  11     1956  3157  8     1957  3676  16     1958  4051  12     1959  5024  18     1960  5282  19     1961  6849  29     1962  7376  21     1963  8951  38     1964  10314  29     1965  13052  60     1966  13601  65     1967  14606  68     1968  18571  97     1969  19796  101     1970  21248  101     1971  20997  101     1972  24852  135     1973  26648  145     1974  30310  170     1975  34124  190     1976  39344  232     1977  45367  240     1978  54883  303     1979  63037  302     1980  73844  390     1981  77437  377     1982  83985  428     1983  90484  482     1984  100153  500     1985  100139  526     1986  106011  529     1987  101197  472


title: Logistic regression and error terms

text: In logistic regression, if we considered residuals, could they only take on the values $0$ or $1$? The data points themselves take on only $1$ or $0$. The logistic curve can take on any value between $0$ and $1$. What would the distribution of the residuals look like?


title: Deriving PDF in AutoRegressive Model

text: I hope that this is a right place and way to ask this question. I am trying to understand how to derive the probability density function of x(t) in an AR model of order K given (t-k) past observations.   I am primarily referring to this paper (Section 2.2) but it cites a Japanese book for this derivation. I could not find similar derivation in the other time series books(Shumway and stoffer, Pourahmadi). Can someone please provide me appropriate references or explain it. Thanks for your time. iinception


title: randomForest vs. cforest; Can I get partial dependence plots and percent variance explained in package party?

text: I have a data set with 24 predictor variables, all continuous, but with different scales and potential collinearity. I’m trying to decide whether to use `randomForest` or `cforest` in party with conditional importance permutation. I recognize that I should probably use `cforest` if I want to overcome variable selection bias, but I find the ability to get partial dependence plots and percent variance explained from the `randomForest` package to be quite appealing. I was wondering if anyone knew if it were possible to get partial dependence plots and percent variance explained from `cforest`? Also, it appears that `ctree` uses a significance test to select variables; is this the same for `cforest`? And how might I get these significance values for each variable in cforest?


title: What does a chi-square test mean when no p-value is returned?

text: I did multinomial logistic regression using SPSS chi-square is .000 , df is 0 and significance =. So what does it mean significance = .?


title: Simulating Poisson Processes

text: My book outlines a procedure but a preliminary part of it is unclear to me. Let X be the number of occurences of an event over a unit of time and assume that it has a Poisson distribution with mean $m=\lambda $. Let $T_1, T_2 , T_3, \ldots $ be the interarrival times of the occurences and they are iid with an exponential $\lambda $ distribution. Note that $ X=k $ iff $$ \sum_{j=1}^k T_j \leq 1 \quad \text{and} \quad \sum_{j=1}^{k+1} T_j >1\. $$ This is precisely what I do not understand. Why does the total waiting tme until $k$ occurences have to be less than or equal to 1? Any help is greatly appreciated. Thank you.


title: What factor extraction method should be used in the FA with optimally scaled variables?

text: I have obtained optimally scaled variables from a highly mixed nature of data containing binary, nominal, ordinal and scale type variables. The optimal scaling was obtained in SPSS through a CATPCA procedure. Now I want to use these variables in a Factor Analysis and want to use a rotation to obtain meaningful loads. What I think is that I should avoid any normality assumption for these optimally scaled variables. So, maximum likelihood method of factor extraction is probably not applicable here (although I am not so good at FA, so not exactly sure). Should I use principal components method of factor extraction instead? What else method can be useful that avoids distributional assumption? If I want the factor scores in a further regression as IVs, will it be at all a good idea to use oblique rotation? I think I should try to keep the factors as uncorrelated as possible if I want them as IVs in a further regression. Is this concept right? Thanks for any kind of suggestion. :)


title: Thresholds for detection of significant change from baseline in time-series data

text: I have a data set containing a daily sensor data measurements recorded from 20 participants for 60 days (baseline data). I am trying to develop methods for predicting/estimating decline in long-term monitoring studies, i.e. can measurement of a parameter on a daily basis be used to detect/predict decline (i.e. significant change) by identification of negative trends or abberant measurements. I hope to be able to generate statistical thresholds to allow identification of decline by examining the trends (using some combination of sensor derived parameters and referencing these to baseline clinical data. What is the most appropriate method to define a threshold for a significant change or decline for each participant and how do I best predict/detect negative trending or aberrant behaviour in unseen data? I wondered if I could get some opinions on a good approach? (note I have been using ICC, ANOVA as well as examining std. dev. of the baseline but have not found any of these approaches particularly useful)


title: MLE of a function of a parameter

text: I am working on a problem where we are interested in finding the MLE for a function of two parameters. I am having problems with going about finding this. Intuitively, the idea makes sense. I am just wondering about the definition of the MLE of a function of two parameters (Google isn't turning up much). The question is as follows: **Question:** Suppose that $X_1,\ldots,\,X_n$ are iid $N(\mu, \sigma^2)$ with unknown $\mu,\sigma^2$. Find the MLE for $\frac{\mu}{\sigma}$. Note that this is not just a homework problem, but part of a take home final. I really am not looking for much of an answer, but more or less the idea for such problems. **Edit** Apparently MLE's are invariant under function. TY


title: If C wins B 80% of the time, and B wins A 80% of the time, how often would C beat A?

text: I am currently trying to better understand probabilistic skill ranking systems for games, but I find that I have trouble properly understanding the basic concept of how skill as a pairwise comparison can be generalized. For instance, if all you know is that player C wins player B 80% of the time, while that same player B wins player A 80% of the time, would this be enough data to determine how often C would win against A? How would those calculations work? Of course it might even be possible for a game to have different styles of play where A might win specifically against C, which would completely confuse the issue, but I am talking about general ranking systems such as ELO or Trueskill that only take winning into account.


title: Comparison of means for individual vs consensus

text: I have a question regarding how to evaluate agreement between an individual's rating, and that of a group of people (of which the individual was a part). The group score was achieved through consensus (i.e. they agreed on one score as a group). I was originally planning to use kappa to look at agreement between the two scores, but I am now questioning this approach. Any ideas on how I can evaluate either the difference or the agreement between the two scores? My main worry is independence. The data looks something like this:   ID IndividualScore GroupScore   1 5 3   2 4 2   etc


title: Computing a "concentration of event occurrence" index

text: I'm looking at the concentration of an event occurrences in a given interval of time. For example, we suppose that an event occurred 4 times in an interval of length 10. I can represent this as a string where `X` means _the event occurred_ whereas `o` means _the event didn't occur_. The distribution of events could be rather "uniform" in time :               X o o X o o X o o X      Or it can be "concentrated" like this :               X X X o o o o o o X      Or like this :               X o o o X X o o o X      I'm looking for an index or a distance that would allow me to measure and compare the concentration of event occurrence in my interval. Ideally it would allow to compare interval of different lengths, for example it would distinguish between :               X o X o X o X o X o X      And :               X o o o X X o X      Do you know of any index or measure that would allow such a thing ? Many thanks in advance.


title: Ranking undergrad students by their future income - Mixture distribution

text: I would be very grateful for some advice on how to model mixture distributions with R. Given a problem to create a ranking of graduate students by their yearly income after completing their education, what are some suited models for this task? Specifically, my data has a distribution with a point mass at 0 (the majority of graduates doesn't find, or start, a full-time job right away). The rest of the data is sort of nicely distributed. The data $x$ was transformed $\log(x+1)$. ![histogram](http://i.stack.imgur.com/4lUnO.png)   * My first approach was a simple regression model    * My second approach were two models (one for classifying whether they get a job - is very weak, and second to predict the income) Simply chaining these two models works much worse than the simple model. My next step would be a Bayesian mixture model to predict income. I was thinking about fitting a mixture with 2 Gaussians, where I would set the mean for one of them to be known as equal to 0. Would that make sense? Has anyone a good experience with some package? Another problem might be that I am always predicting the income using a regression and building a ranking from that, rather than running an ordinal regression. What is the best way handle this situation - if the target variable (income) that the ranking is based on is itself available for training data? Full disclosure: This is a fictitious scenario, as I cannot discuss the exact details of the real case.


title: How to statistically compare the performance of machine learning classifiers?

text: Based on estimated classification accuracy, I want to test whether one classifier is statistically better on a base set than another classifier . For each classifier, I select a training and testing sample randomly from the base set, train the model, and test the model. I do this ten times for each classifier. I therefore have ten estimate classification accuracy measurements for each classifier. How do I statistically test whether the $classifier 1$ is a better classifier than the $classifier 2$ on the base dataset. What t-test is appropriate to use?


title: Confidence interval

text: A random sample of **388** married couples found that **292** had two or more personality preferences in common. In another random sample of **552** married couples, it was found that only **24** had no preferences in common. Let p1 be the population proportion of all married couples who have two or more personality preferences in common. Let p2 be the population proportion of all married couples who have no personality preferences in common. Find a 99% confidence interval for p1 – p2. Lower Limit: Upper Limit:


title: Best way to display data with outliers?

text: I am fairly new to data analysis and visualization, and I'm trying to figure out the best model to show some data regarding page load time (in seconds). The current view is a line graph where the x-axis is the page, and the y-axis is the load time. This graph is visually misleading though, because in my mind, the psychological expectation of a line graph is to convey something over time. Because there are a couple of outliers in the data where the page load time is much longer than the average, someone looking at the graph might be turned off because their eyes jump to the outliers, so it looks like our pages take forever to load. In actuality, about 95% of the pages take fewer than two seconds to load, and 5% of the pages take longer. So what's the best method to display this data? I want people to understand that on average, our page load time is _great_ , and there are very few exceptions. I feel like depicting a bell curve might help to show that the majority of the pages fall within x-standard deviations of the mean, but A) I'm not sure how to do that in Excel, and B) I'm not sure that's the best option.


title: pdf of the product of two independent random variables, normal and chi-square

text: what is the pdf of the product of two independent random variables X and Y, if X and Y are independent? X is normal distributed and Y is chi-square distributed. Z = XY if $X$ has normal distribution $$X\sim N(\mu_x,\sigma_x^2)$$ $$f_X(x)={1\over\sigma_x\sqrt{2\pi}}e^{-{1\over2}({x-\mu_x\over\sigma_x})^2}$$ and $Y$ has Chi-square distribution with $k$ degree of freedom $$Y\sim \chi_k^2$$ $$f_Y(y)={y^{(k/2)-1}e^{-y/2}\over{2^{k/2}\Gamma({k\over2})}}u(y)$$ whre $u(y)$ is unit step function. Now, what is the pdf of $Z$ if $X$ and $Y$ are independent? One way to find the solution is to use Rohatgi's well known result (1976,p.141) if $f_{XY}(x,y)$ be the joint pdf of continuous RV's $X$ and $Y$, the pdf of $Z$ is $$f_Z(z) = \int_{-\infty}^{\infty}{{1\over|y|}f_{XY}({z\over y},y)dy} $$ since, $X$ and $Y$ are independent $f_{XY}(x,y)=f_X(x)f_Y(y)$ $$f_Z(z) = \int_{-\infty}^{\infty}{{1\over|y|}f_{X}({z\over y})f_{Y}(y)dy} $$ $$f_Z(z) = {1\over\sigma_x\sqrt{2\pi}}{1\over{2^{k/2}\Gamma({k\over2})}}\int_{0}^{\infty}{{1\over|y|}e^{-{1\over2}({{z\over y}-\mu_x\over\sigma_x})^2} {y^{(k/2)-1}e^{-y/2}}dy} $$ Where we face the problem of solving the integral $\int_{0}^{\infty}{{1\over|y|}e^{-{1\over2}({{z\over y}-\mu_x\over\sigma_x})^2} {y^{(k/2)-1}e^{-y/2}}dy}$. Can anyone help me with this problem. is there any alternative way to solve this?


title: Two-Parameter Exponential Family Conjugate Prior

text: A probability distribution is a member of the two-parameter exponential family if the distribution can be expresses in the following form: $$ h(\theta, \phi) \text{exp}\left[\sum t(x_{i})\psi(\theta, \phi) + \sum u(x_{i})\chi(\theta, \phi)\right] $$ for parameters: $\theta$ and $\phi$, data: $x_{1}, x_{2}, \ldots, x_{n}$ and functions: $h$, $t$, $u$, $\psi$ and $\chi$. Furthermore, the conjugate prior for the two-exponential family takes the following form: $$ p(\theta, \phi) \propto h(\theta, \phi)^{v}\text{exp}\left[\tau\psi(\theta, \phi) + \omega\chi(\theta, \phi)\right] $$ Whilst looking at some members of the two-parameter exponential family, I observed the following example: The gamma distribution ($\text{Ga}(\alpha, \beta)$) can be expressed in the above form, where $h = \frac{\phi^{\theta}}{\Gamma(\theta)}$, $\sum t = \sum \text{log}x_{i}$, $\sum u = \sum x_{i}, \psi = \theta, \chi = -\phi$. The example also states, following from the above, that the conjugate prior form for this gamma distribution is given by: $$ \left[\frac{\phi^{\theta}}{\Gamma(\theta)}\right]^{v} \text{exp}[\tau\theta - \omega\phi] $$ Could somebody explain how this gamma conjugate prior can be expressed in this form?


title: Polychoric PCA and component loadings in Stata

text: I’m using Stata 12.0, and I’ve downloaded the `polychoricpca` command written by Stas Kolenikov, which I wanted to use with data that includes a mix of categorical and continuous variables. Given the number of variables (around 25), my hunch is that I will need to generate more than 3 components. Ultimately, I would like to generate a handful of meaningful components (rather than dozens of variables) and use the components as independent variables in logistic regression. Using `polychoricpca`, I am able to generate a table showing the eigenvalues and the eigenvectors (loadings) for each variable for the first three (3) components only. `polychoricpca` appears to call these loadings “scoring coefficients” and produces these for every level of the variable, such that if a variable has three categories you’ll see three scoring coefficients (“loadings”) for that variable. Never having worked with polychoric PCA before, I’m used to only seeing one loading per variable/item. I want to examine these coefficients (“loadings”) to try to understand what the components are and how they might be labelled. My questions: (1) What if it looks as if I should generate 4 components? It seems as if I wouldn’t be able to examine and understand what that 4th component is because I can’t see how each of the items load on that 4th component, only the first 3. Is there a way to see how each item loads on more than the first three components? (2) Can I simply use the polychoric correlation matrix combined with Stata’s `pcamat` command to examine how each item loads on each component (the eigenvector table). I thought this might be a way of being able to examine loadings if I have more than 3 components. The idea came from this UCLA stats help post on using `factormat` with a polychoric correlation matrix. `pcamat` in Stata, however, produces only 1 loading (coefficient) per variable, not 1 loading for every level of the variable. Any thoughts on whether it would be appropriate just to report the single loading from `pcamat`?


title: Survival Analysis--All variables failed the proportionality test. What next?

text: I am building a cox survival analysis model with a retailer's transaction data. Almost all the variables have failed the proportionality test. Can I continue with the Cox model? Should I build a LIFEREG model instead of the Cox model?


Let $\mathbf{X}=(X_1,X_2,\ldots,X_{2k})$ be a sample from normal distribution $\mathcal{N}(0,1)$. Prove that distribution of sample median (for even sample) is symmetric around 0. $me(\mathbf{X}) = \frac{X_{k:2k}+X_{(k+1):2k}}{2}$, where $X_{i:2k}$ is $i$-th order statistic. I was trying to use formula $f_{X_{k:n}}(x) =\frac{n!}{(k-1)!(n-k)!}[F_X(x)]^{k-1}[1-F_X(x)]^{n-k} f_X(x)$, but without any success (maybe there is a way to do it differently than finding $me(\mathbf{X})$ distribution?


title: Centering data in multiple regression

text: In a multiple regression analysis (with 4 continuous predictors and 2 categorical factors), we mean centered the data (for each continuous variable) due to issues of multicollinearity when the interaction terms are included. My question is whether I can center the response variable too. More specifically, the response variable and the 4 continuous predictors are all averaged survey responses (using scales of 1 to 5). I originally thought that if I were centering the explanatory variables, I might as well center the response.. but I realize that most references to centering seem to only apply to the predictor variables. Any help is much appreciated. **Update** : as I mentioned in an earlier comment, I questioned the validity of centering my response variable due to having different ANOVA results using centered versus non-centered response. My interpretation is that the linear model for a non-centered response (using the lm() in R) uses the mean of the 'reference level' of my two factors for computing the intercept. When I centered the response variable, it's being subtracted from the grand mean of this variable, rather than the mean of the reference level. Now I've verified that by 'centering' using the reference level mean and it does indeed yield identical results for p-values as the non-centered response model. I hope I am interpreting this issue correctly. IF someone could further confirm/clarify this I would really appreciate that.


title: Kolmogorov-Smirnov test?

text: I am studying disturbance caused by ship traffic to a small seabird. I observed focal animals for a set amount of time and record whether or not they fly from the water during the observation. This particular bird does not fly at high probabilities when not disturbed (about 10% of the time). Post hoc, I have appended the distance to the nearest ship to every observation (ships of interest had GPS locators recording a point every 5 sec). I have plotted the cumulative distribution function for ALL observations and for observations where the bird flew from the water as a function of distance to nearest ship. As expected, the majority of the observations in which the bird flew were observed when the ship was close. ![ecdf plot, fly, did not fly, all obs](http://i.stack.imgur.com/0X3MY.jpg) Can I use the Kolmogorov-Smirnov test to test if there is a statistical difference in the distribution of the flight observations and the total observations? My thought is if these two distributions are different it would suggest that ship distance has an influence on flight. I worry as these distribution functions are not independent as the flight observations are a subset of the total observations. Thoughts? Having read a bit further on this site I think I can test the distribution of the observations in which flight occurred (F) against the distribution of the observations in which it did not (NF) as these are independent. If these distributions are the same F=NF, then we can assume that the distribution of (F) and (TOT = all observations) are the same as we know that the distribution of (F) is equal to itself and (F)+(T) = (TOT). Right? **UPDATE: 2/12/14** Following the suggestions from @Scortchi I investigated the relationship of the incidence of flight vs distance to the nearest ship in a logistic regression framework. There was slight relationship present (negative slope) but the p-value was non-significant suggesting that the true slope could be zero. Based on the desciptive statistics (including the ecdf plots) I suspected that the effect of close ships was being drowned out by the many observations when the ship was not affecting behavior. I then used the R package segmented (http://cran.r-project.org/web/packages/segmented/segmented.pdf) to try and find a break-point in the model. The program found that breaking the data at 2.6 km from the ship and fitting two separate coefficients was better than the single coefficient model. The coefficient for the slope of the close ship approaches was negative and suggests that ships affect the flight response until about 2.6 km (p-value < 0.001). The coefficient for the second slope was slightly positive but the p-value was non-significant at the 0.05 alpha level (p-value = 0.11). So in summary, the segmented regression line was able to detect a threshold difference at which flight probability increases. The estimate for the probability of flight when the ship is further than 2.6 km is 0.11. Fittingly, I observed 79 birds in when no ships were even in the study bay (> 100 km away) and 9 of these individuals flew from the water during the observation (0.113 probability). Thanks for all the suggestions. I hope this question along with the suggestions and answers helps others.


title: How do you check ergodicity of a stochastic processes from its sample path(s)?

text: How do you check ergodicity of a wide-sense stationary stochastic processes from its sample path(s)? Can we check ergodicity from a single sample path? Or do we need multiple sample paths? One motivation of checking ergodicity is in time series to ensure that you could safely use the average of a sample path over time as estimate for the population mean? Thanks!


title: confidence level query

text: A new treatment for strokes is put on trial. There are two equal size groups,one group is given a drug, and one a placebo. The 95% confidence interval for the difference between the two proportion of subjects having strokes was (0.07,0.12). What can we infer: 1) not enough info to infer 2)at most 12% of the people had strokes 3) at least 4% had strokes I'm not sure on this. I think it might be 1). Can anyone help out? Thanks


title: How to apply hurdle models to panel data (using Stata)?

text: Is it possible to apply hurdle models (like the Craggit, probit and truncated models) to panel data, preferably with fixed effects to control for unobserved heterogeneity? In Stata, the user-written command `craggit` only allows to use pooled panel data, but not to control for unobserved heterogeneity... In general, the hurdle model is quite simple to estimate: probit model on the binary outcome variable, truncated (or lognormal) on the (positive values of the) outcome variable. Is it possible to just estimate a RE probit (or even FE logit) in a first step, and a fixed effects model in the second – and to actually call that a hurdle model? Thanks.


title: Heteroskedasticity in a Linear Mixed Model SAS PROC MIXED

text: Asked a version of this question before but realized it needed some clarification. I have a dataset with identical twin pairs and fraternal twin pairs. I want to examine the relationship between an outcome variable ("DEPENDENT") and a number of other variables ("INDEPENDENT"). However, I can't run a normal OLS regression because each twin's dependent variable is correlated with their co- twin's (for my purposes, this is a nuisance). The way I have been dealing with this is to use SAS PROC MIXED and include a random intercept defined by twin cluster (FAMILYID). Here is my syntax:               proc mixed method=ml covtest noclprint;        class FAMILYID;        model DEPENDENT = INDEPENDENT/solution;        random intercept/sub=FAMILYID type=un gcorr;      run;       When I do this, I get a covariance parameter estimate for the "FAMILYID" as well as a residual term. My understanding is that this decomposes the variability into "within-cluster" variability (residual estimate) and "between-cluster" variability (FAMILYID estimate). I get the intraclass correlation coefficient for the family clusters by computing: (FAMILYID estimate) divided by (FAMILYID estimate + residual estimate). However! I realized that I probably have a heteroskedasticity issue, because fraternal twin pair members are likely to be markedly less similar to each other than identical twin pair members are. I revised my syntax to this:               proc mixed method=ml covtest noclprint;        class FAMILYID;        model DEPENDENT = INDEPENDENT/solution;        random intercept/sub=FAMILYID group=TWINTYPE type=un gcorr;      run;       Where TWINTYPE signifies whether the twin pair is identical or fraternal. When I do this, I get one FAMILYID estimate for fraternal twins, one FAMILYID estimate for identical twins, and one residual term. So -- this means that I have told my model to make different "between-cluster" estimates for identical and fraternal twins -- but isn't it still disregarding twin type when it calculates the residual term? And if the residual term reflects the "within- cluster" variability, isn't it a problem that it doesn't take into consideration that there is likely to be LESS variability within an identical twin pair on the outcome variable and MORE variability within a fraternal twin pair? So I guess my question is one of two things:   1. is there syntax that would leave me with a between-cluster estimate for identical twins and fraternal twins, and then a residual term for identical twins and fraternal twins? or    2. Am I fundamentally confused about how all of this works, and if so, what am I missing? Help/clarification much much appreciated.


title: Finding transformation function for a distribution that looks like exponential

text: Suppose that we have two data sets, `R` and `P`. `R` is larger than or equal to `P`. `R` can be negative or positive. `P` can be negative or positive. So in all cases `(R-P)` is positive or zero. the `R`and `P` are dollars and the scale of them is high likes `1,000,000,0` or `500,000,000`. Now I want create a criterion looks like :               E = 1 - ( 1 / 1+(R-P) )      So I have `E` that is between `0` and `1`. I want `E` scale between `0` and `1`. But the problem is this. `(R-P)` is has high values,so all outputs of this criterion is near `1`. I revise this function, Now i have this :               E = 1 - ( 1 / 1+ ( (R-P)/mean(R-P) ) )      Now because of `mean(R-P)` we have better scaling between `0` and `1`. This is distribution of `R-P` for more help : ![enter image description here](http://i.stack.imgur.com/2XOxM.jpg) After doing that I want classify data with this `E` . For example samples that have a `E` greater than `0.7` are in `class 1` and samples with `E` lower than `0.7` are in `class 2` (In a binary classification problem). For a better distribution in final version of this criterion i will remove outliers of `(R-P)` with this function to have more robust criterion :               abs(X-mean(X)) >=2*std(X)      `(X here is (R-P)` After removing outliers in `R-P` i will calculate `mean(R-P)`what is your recommendations for better transformation of this problem? I have a theoretical problem with using `mean` in this criterion so i need comments about this. **PS.** We should have this behavior in this criterion:   * if `P=constant` and `R=increase` >> `E=increase`   * if `R=constant` and `P=increase (towards negative)` >> `E=increase`   * if `P=constant` and `R=decrease` >> `E=decrease`   * if `R=constant` and `P=decrease(towards positive)` >> `E=decrease`   * a sample that have highest `(R-P)` should have highest `E` (`1` or near `1`)   * a sample that have lowest `(R-P)` should have lowest `E` (`1` or near `1`)


title: Multi-armed bandit in face of full reward information

text: I am new to this area of machine learning. I am just walking myself through `UCB1` algorithm which seems to assume that the payoff can be learnt only for action that is taken. What I am curious about is if there is already analysis of the multi-armed bandit problem to the cases where the payoff can also be learnt for actions not taken (may be the term for this is "full reward vector"?), and that information used in subsequent steps. I could modify `UCB1` algorithm to update payoff for actions not taken as well. But, I would like to know if there are any drawbacks/pitfalls for doing this. I haven't done closed-form analysis of the regret bounds when reward vector is fully available. I don't know what the terminology is for such cases. So, I have a hard time finding work on such cases. I understand that when full reward vector is available, it is a supervised reinforcement learning problem. What I am curious about is if multi-armed bandit work can be used well in supervised reinforcement learning problem or not. And whether it makes sense to do this or not. I will very much appreciate pointers to work/paper in this area.


title: What is the proper way to estimate the CDF for a distribution from samples taken from that distribution?

text: Given $n$ samples from a (continuous) distribution X, the obvious thing to do is sort them, and distribute them equally across $[0,1]$ by taking $(x_{(k)}, (k-1/2)/n)$ as estimates of particular points on the CDF, and doing some sort of interpolation between points, as necessary. Is this the "right" way to make this estimate? How do I get error bars for the estimated points? It doesn't seem like they'd necessarily be symmetric.


title: How to conduct a meta-analysis on studies that report results variously as odds ratios, hazards ratios, or rate ratios?

text: I'm doing a meta-analysis of some studies that report results variously as odds ratios, hazards ratios or rate ratios (all with confidence intervals). Is there any way to combine these together/convert between them so that I can do a meta-analysis of all the studies?


title: How do I test whether one predictor is significantly better than another? Is Hotelling's T the best option?

text: ### Overview: I want to test if "emotional numbing" is a significantly better predictor of "lower intimate relationship functioning" than "reexperiencing" or "hyperarousal" It was suggested to me to use stepwise regression to analyze the data. However, one of my dissertation committee members suggested Hotelling's T. ### My hypothesis:   * Hypothesis 1: Emotional numbing will be a better predictor of lower relationship functioning than re-experiencing or hyperarousal symptoms.    * Hypothesis 1a: Emotional numbing will be a better predictor of lower dyadic adjustment than re-experiencing or hyperarousal symptoms.    * Hypothesis 1b: Emotional numbing will be a better predictor of intimacy than re-experiencing or hyperarousal symptoms.    * Hypothesis 1c: Emotional numbing will be a better predictor of communication than re-experiencing or hyperarousal symptoms. ### Questions:   * Is Hotelling's T the most appropriate way to determine whether emotional numbing is a significantly better predictor than the other 2 variables?   * Or is there is a better statistical test for this task?


title: Applications of Data Mining in economics

text: I am looking for a thesis topic for my Applied Economics degree where I can use Data Mining/Machine Learning for some kind of prediction. Do you guys have any good ideas about what I could do? My first thought was credit scoring, but that seems a little boring considering the number of papers out there and the lack of available data. I don't want to write about stock market forecasting because that's what I did for my mathematician thesis. A good topic would be predicting some sort of an economical indicator but I don't really have any good ideas so far. Thank you!


title: What distribution does “arrival of people” at an event in time follow?

text: Suppose that there is a Cricket match scheduled on Sunday, this weekend. We know that people do not arrive at the stadium at constant rate. Few hours before scheduled start of the match, people start arriving at very slow rate. This rate then, in general, keeps on increasing. A large number of people arrive just before start of the match. I would like to know which distribution should I use to model this, so that I could tell **What is the probability of arrival of x number of people in a particular time interval?**


title: What is loc parameter in GPD distribution in POT package for R?

text: I fitted the Generalized Pareto distribution (GPD) using the `POT` package in `R`. The fitted object provides shape and scale parameters, but not the location parameter. I need it to evaluate the quantile function with `qgpd(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, lambda = 0)`. I am sure this is very basic question (maybe it is mean value, but from which values - all, or only that which are bigger than threshold?) but I am begginer in this area and help is highly appreciated.


title: Interaction Term Insignificant but Composite Variable Significant

text: I have 4 latent variables each with 3-7 indicators. I created a composite for each variable by summing the data and then mean-centering (variables have some collinearity issues, highly related). Thus I have Var1, Var2, Var3, Var4. I then wanted to create a master composite of all four variables. I thus summed all the individual indicators and mean-centered the total, for Var1234. I'm trying to test the theory of complementarity in which these four variables are greater together than individually (i.e. var1 + var2 + var3 + var4 = 4 ... but Var1234 >4) I ran two linear regressions.   MODEL 1: All the individual indicators + Individual interaction terms + 4 way interaction term = DV MODEL 2: Var1234 = DV In the 2nd model Var1234 is significant, however in the 1st model the 4 way interaction term is insignificant. I believe this is a problem, but I am not sure. Is it related to the way I created my composite variables, or the fact that there is a 2nd order variable that is a composite of other composites? Is there a way to explain why a 4 way interaction term is insignificant when the composite of those 4 terms IS significant?


title: How to calculate estimated proportions and their confidence intervals from a mixed model?

text: I have an experiment with two treatments. It is a split plot experiment, with the structure Block/Treatment1/Treatment2. Each treatment has 2 levels. The dependent variable is presence/absence data for functional speceis groups. I am analysing the data using mixed models via lme4 and afex. The model structure is as follows:               m <- lmer (DV ~ Treatment 1 * Treatment 2 + (1|Block/Treatment1),                family = binomial)       For graphical purposes however, I want to plot the treatment means for each of the 4 treatment combinations (Treatment 1 level A and Treatment 2 level A; Treatment 1 level B Treatment 2 level A; Treatment 1 level B and Treatment 2 level A; Treatment 1 level B and Treatment 2 level B). I couldn't find an intuitive way to extract these mean and standard error values in R for data with a binomial error distribution, thus I am calculating these values by hand. I appreciate that the standard errors I calculate will not take into account the random effect structure in my dataset, but I cannot find a way to calculate them to include this random effect structure (using a mcmc simulation in language R to create HPD confidence intervals only works for Gaussian data). As such, I want to calculate treatment means for the proportion of a particular functional group in a particular treatment combination, and then standard errors of these 4 treatment means. Consider that the proportions I have for a particular treatment are:               Plot 1: 24/65     Plot 2 26/64     Plot 3: 25/65     Plot 4: 22/62     Plot 5: 30/66     Plot 6: 29/65      I understand that to calculate the average proportion I need to sum the total number of hits (sum of all the numerators above) and divide by the total sum of hits and misses (sum of all the denominators above). This gives me 0.40. My question then is, how does one calculate the standard error of the resulting proportion? Calling the proportion of hits p and misses q, I have the following formula for the standard error of a proportion:               sqrt((p*q)/n)      Is this correct? In addition, what is n? Is n the number of proportions that make up the mean proportion (in this case 6) or the total number of hits and misses (i.e. sum of all the denominators above, which here is 387)? I'm sorry that this is a simple question, I just cannot find an answer anywhere! Also i'm sorry if this is not the correct forum to post this in... Many thanks.


title: Optimisation model paramaters for time series data using k-fold cross validation

text: I'm trying to optimise model parameters for times series data using k-fold cross-validation. I'm using the forward chaining strategy described here: Using k-fold cross-validation for time-series model selection i.e. fold 1 : optimise [1], test [2] fold 2 : optimise [1 2], test [3] fold 3 : optimise [1 2 3], test [4] fold 4 : optimise [1 2 3 4], test [5] fold 5 : optimise [1 2 3 4 5], test [6] Then averaging the resulting optimised parameters from all five folds to give the final model parameters. My question is: does this not put an overly strong weighting on fold 1 which would be included in every optimisation run? Is there a way to avoid this, or a better cross-validation strategy?


title: When an analytical Jacobian is available, is it better to approximate the Hessian by $J^TJ$, or by finite differences of the Jacobian?

text: Let's say I'm computing some model parameters my minimizing the sum squared residuals, and I'm assuming my errors are Gaussian. My model produces analytical derivatives, so the optimizer does not need to use finite differences. Once the fit is complete, I want to calculate standard errors of the fitted parameters. Generally, in this situation, the Hessian of the error function is taken to be related to the covariance matrix by: $$ \sigma^2 H^{-1} = C $$ where $\sigma^2$ is the variance of the residuals. When no analytical derivatives of the error are available, it is typically impractical to compute the Hessian, so $J^TJ$ is taken as a good approximation. However, in my case, I've got an analytical J, so it's relatively cheap for me to compute H by finite differencing J. So, my question is this: Would it be more accurate to approximate H using my exact J and applying the above approximation, or to approximate H by finite differencing J?


title: Find all correlations between two sets of variables

text: I have a list of records with fields that can be seen as two sets of variables: dependent and independent. There are say 100 fields in each record and we consider 10 of them as dependent variables and 90 as independent. What are the methods to find all the correlations that exist between the two sets? (At introductory statistics course that I've taken we have been dealing with simple correlations between two variables but the task like this was not covered.)


title: Light bulb color problem

text: Please first have a look at the following little problem: > There are two indistinguishable light bulbs A and B. A flashes red light > with prob .8 and blue with prob .2; B red with .2 and blue .8. Now with .5 > prob you are presented with either A or B. You're supposed to observe its > flash color to make a best guess (maximizing probability of correct > guessing) which bulb it is. Before you start to make observations, however, > you must decide how many times you want to observe it (say n times, then you > observe it flashing n times and make your guess). Suppose flashes are > independent. Intuitively, one would think the more observations one makes, the better one's chances will be. Curiously though, it is easy calculation to show that n=2 doesn't improve upon n=1, and n=4 doesn't improve upon n=3. I didn't go further but I speculate n=2k doesn't improve upon n=2k-1. I'm not able to prove it for the general case. But is it true? If so, how can one intuitively understand the result?


title: Individual patient data meta-analysis with SD's as outcome measure

text: I'm reading the results of an individual patient data (IPD) meta-analysis. The outcomes of the study is expressed in standard deviations (SD's). Comparisons between groups are stated as 0,23 (95% CI 0.15-0.33), SD 0.34, etc. I can't grasp how to interpretate this. Can somebody help?


title: How to use R gbm with distribution = "adaboost"?

text: Documentation states that R gbm with distribution = "adaboost" can be used for 0-1 classification problem. Consider the following code fragment:               gbm_algorithm <- gbm(y ~ ., data = train_dataset, distribution = "adaboost", n.trees = 5000)     gbm_predicted <- predict(gbm_algorithm, test_dataset, n.trees = 5000)      It can be found in the documentation that predict.gbm > Returns a vector of predictions. By default the predictions are on the scale > of f(x). However the particular scale is not clear for the case of distribution = "adaboost". **Could anyone help with the interpretation of predict.gbm return values and provide an idea of conversion to the 0-1 output?**


title: How to derive a distance function based on multiple variables for cluster analysis?

text: _I am not a statistician, so please excuse my lack of statistics knowledge/terminology._ I have bunch of network nodes that I want to run cluster analysis on and identify clusters. So as far as I understand, I can follow the following steps to run a hierarchical agglomerative analysis (HAC):   1. Identify variables   2. Define a distance function   3. Run the algorithm to join closer clusters and create one big cluster   4. Cut the dendrogram tree at the height that makes meaningful clusters based on context My question is related to the second step even though it is not yet clear for me how I am going to do the last step. The data I want to analyse is bunch of computer network nodes which are down (not responding). I have the following information for each network node:   * location (longitude, latitude)   * time node went down   * network provider These are the most relevant information I believe I have to take into consideration for my clustering. Basically I want to cluster the nodes that went down probably because of the same reason in a region. For example if bunch of nodes went down at about the same time and physically they are close to each other and they have the same provider, so probably they fall into the same cluster. Now the question is how do I derive my distance function and include all these variables in it such that it would make sense? In other words **what is the mechanism to derive a distance function based on multiple variables?** Also as you notice the variables are of different types. So, in the distance function should I take care of this by considering **Gower's coefficient of similarity**? How? Any examples or suggestion regarding whether I am in the right direction or not can be very helpful too.


title: How to calculate sum of squared errors (SSE or SSR) with Stata 12 software?

text: I have a techincal issue with Stata: I need to calculate the SSE of a regression model, but the automatic output just gives me RMSE; I need SSE of this constrained model because I have to test if this model is better than the completely-free model. SPSS gives me the SSE value in the regression output, but I don't know how to calculate a constrained linear regression with this model.


title: Clustering based on interactions between members of two groups

text: I have a large dataset from a survey that describes what web pages people use. So for each person I have a list of pages that they visit and how frequently they visit them. What methods can be used to cluster both the webpages (based on similar users who use them) and the users (who use similar webpages).               Example>     User A visited web pages a, b, a, c     User B visited a, b, d, b     User C visited page e, e      Conclusion from that would be that pages `a`,`b`,`c`,`d` are somewhat similar and belong to the same cluster. The other cluster is page `e`. Also, users a,b are similar, the other user cluster is the one containing `User C`, who visits different pages.


title: Updating linear regression efficiently when adding observations and/or predictors in R

text: I would be interested in finding ways in R for efficiently updating a linear model when an observation or a predictor is added. biglm has an updating capability when adding observations, but my data are small enough to reside in memory (although I do have a large number of instances to update). There are ways to do this with bare hands, e.g., to update the QR factorization (see "Updating the QR Factorization and the Least Squares Problem", by Hammarling and Lucas), but I am hoping for an existing implementation.


title: Hypothesis testing on the inverse covariance matrix

text: Suppose I observe _i.i.d._ $x_i \sim \mathcal{N}\left(\mu,\Sigma\right)$, and wish to test $H_0: A\ $vech$\left(\Sigma^{-1}\right) = a$ for a conformable matrix $A$ and vector $a$. Is there known work on this problem? The obvious (to me) attempt would be via a likelihood ratio test, but it seems like maximizing the likelihood subject to the constraints of $H_0$ would require a SDP solver and could be pretty hairy.


title: Queries regarding gaussian processes and kriging

text: I am confused how gaussian processes and kriging are related? Can anyone please give me some simple explanation. I tried to go through the wiki but I didn't get it


title: Role of n.minobsinnode parameter of GBM in R

text: I wanted to know what the _n.minobsinnode_ parameter means in the GBM package. I read the manual, but it is not clear what it does. Should that number be small or large to improve the results?


title: How accurate is F test in panel data

text: I heard that the F-test is to advice you whether to use fixed effects or pooled OLS. However, I didn't find any details about it in books. Only in a very few studies. What is the hypothesis of the test? After I run my fixed model I got               F test that all u_i=0:     F(100, 389) =     2.98                 Prob > F = 0.0000      As I understand `u_i=0` means I shall use fixed effect (am I right?), in general how accurate is this F-test compared to the Hausman test and LM test given the fact that we will end to knowing whether to use fixed effects or pooled OLS?


title: A list of common time series tests?

text: I am reviewing my time series knowledge and looking for a document that has the commonly-used time series tests, what they are used for, how to use them, etc. e.g. Augmented Dickey–Fuller test, PACF tests, etc. I found a wikipedia page of common statistical tests, but I am looking for a list of such that is specific for time series analysis. http://en.wikipedia.org/wiki/Statistical_hypothesis_testing#Common_test_statistics Thanks!


title: Confusion regarding kriging

text: I was reading this wikipedia article related to kriging . I didn't understand the part when it says that Kriging computes the best linear unbiased estimator, $\hat Z (x_0)$, of $Z(x_0)$ such that kriging variance of is minimized with the unbiasedness condition. I didn't get the derivation and also how the variance is minimized. Any suggestions? Specially, I didn't get the part where applies minimized subject to unbiasedness condition. I think it should have been E[Z'(x0)-Z(x0)] instead of E[Z'(x)-Z(x)] isn't it. ' is equivalent to hat in the wiki article. Also I didn't get how the kriging error is derived


title: Does anyone use the web based spreadsheet EditGrid.com?

text: I have been looking for an alternative to Google Docs as a web based spreadsheet application as my employer's corporate firewall now blocks this application. I kicked around a few alternatives such as Zoho.com. I found the latter practically unusable. All its applications took a long time to load and really did not work well if at all most of the time. When I came across EditGrid, I was delighted. It feels pretty much like Excel on a lot of counts. However, there is no apparent support besides a Help Center that is inactive and invariably does not address the problems I have. For instance, I have no idea how to include titles and legends within charts. I also don't see how one can create folders. I thought maybe someone would be familiar with EditGrid and either being able to help me out directly or point me in a direction of a support group somewhere. Thanks in advance for any info.


title: Methods for fitting a "simple" measurement error model

text: I am looking for methods which can be used to estimate the "OLS" measurement error model. $$y_{i}=Y_{i}+e_{y,i}$$ $$x_{i}=X_{i}+e_{x,i}$$ $$Y_{i}=\alpha + \beta X_{i}$$ Where the errors are independent normal with unknown variances $\sigma_{y}^{2}$ and $\sigma_{x}^{2}$. "Standard" OLS won't work in this case. Wikipedia has some unappealing solutions - the two given force you to assume that either the "variance ratio" $\delta=\frac{\sigma_{y}^{2}}{\sigma_{x}^{2}}$ or the "reliability ratio" $\lambda=\frac{\sigma_{X}^{2}}{\sigma_{x}^{2}+\sigma_{X}^{2}}$ is known, where $\sigma_{X}^2$ is the variance of the true regressor $X_i$. I am not satisfied by this, because how can someone who doesn't know the variances know their ratio? Anyways, are there any other solutions besides these two which don't require me to "know" anything about the parameters? Solutions for just the intercept and slope are fine.


title: Spam filtering using naive Bayesian classifiers with the e1071/klaR package on R

text: I'm looking at doing text classification/spam filtering using naive Bayesian classifiers with the e1071 or klaR package on R. Is there a good tutorial out there to describe this? I'm kind of stuck because I'm not sure what to use as the data to input into the NaiveBayes function. Some help very much appreciated, thanks!


title: Period effects in pooled time series data in R

text: This is closely related to a question I asked yesterday but I've now got a much more complete answer on which I was hoping to get feedback. The previous question was just looking for conceptual advice and was very helpful. You can find the relevant data and introduction here. I wanted to find period effects for each age group. I've run two regressions using dummies as part of an interaction term. I'm hoping to see if my method is flawed and if my interpretation of the results is correct or not. They first regression is as follows:               > ## Generate YearDummy and AgeGroupDummy using factor()     >      > YearDummy <- factor(YearVar)     > AgeGroupDummy <- factor(AgeGroup)     >      > ## Check to see that YearDummy and CohortDummy are indeed factor variables     >      > is.factor(YearDummy)     [1] TRUE     > is.factor(AgeGroupDummy)     [1] TRUE     > ## Regress on AgeGroup and include AgeGroup*YearDummy interaction terms     >      > PooledOLS1 <- lm(PPHPY ~ AgeGroup + AgeGroup*YearDummy + 0,      data=maildatapooled)     > summary(PooledOLS1)     Call:     lm(formula = PPHPY ~ AgeGroup + AgeGroup * YearDummy + 0, data =      maildatapooled)     Residuals:      Min 1Q Median 3Q Max      -38.852 -10.632 3.298 11.275 26.481      Coefficients:      Estimate Std. Error t value Pr(>|t|)      AgeGroup 26.2212 3.5070 7.477 3.84e-10 ***     YearDummy1 119.8836 15.6840 7.644 1.99e-10 ***     YearDummy2 123.7458 15.6840 7.890 7.55e-11 ***     YearDummy3 103.2660 15.6840 6.584 1.28e-08 ***     YearDummy4 97.7102 15.6840 6.230 5.06e-08 ***     YearDummy5 103.3295 15.6840 6.588 1.26e-08 ***     YearDummy6 103.2330 15.6840 6.582 1.29e-08 ***     YearDummy7 84.8291 15.6840 5.409 1.16e-06 ***     YearDummy8 70.7114 15.6840 4.509 3.09e-05 ***     YearDummy9 90.9566 15.6840 5.799 2.65e-07 ***     YearDummy10 50.0885 15.6840 3.194 0.00224 **      YearDummy11 37.7004 15.6840 2.404 0.01933 *      YearDummy12 33.1947 15.6840 2.116 0.03846 *      AgeGroup:YearDummy2 1.8066 4.9597 0.364 0.71695      AgeGroup:YearDummy3 -3.8022 4.9597 -0.767 0.44632      AgeGroup:YearDummy4 -1.7436 4.9597 -0.352 0.72640      AgeGroup:YearDummy5 -6.0494 4.9597 -1.220 0.22735      AgeGroup:YearDummy6 -6.7992 4.9597 -1.371 0.17552      AgeGroup:YearDummy7 -3.6752 4.9597 -0.741 0.46158      AgeGroup:YearDummy8 -0.4799 4.9597 -0.097 0.92323      AgeGroup:YearDummy9 -9.8190 4.9597 -1.980 0.05232 .      AgeGroup:YearDummy10 -2.2452 4.9597 -0.453 0.65241       My interpretation of the interaction term coefficients is that they represent the difference in slope of AgeGroup between the period of the corresponding YearDummy and the AgeGroup slope at the top of the results. This is kind of like the AgeGroup effect across different periods. My second regression is as follows:               > ## Regress YearVar and Include YearVar*AgeGroupDUmmy     >      > PooledOLS2 <- lm(PPHPY ~ YearVar + YearVar*AgeGroupDummy + 0,      data=maildatapooled)     > summary(PooledOLS2)     Call:     lm(formula = PPHPY ~ YearVar + YearVar * AgeGroupDummy + 0, data =      maildatapooled)     Residuals:      Min 1Q Median 3Q Max      -29.345 -9.325 -0.915 8.540 40.150      Coefficients:      Estimate Std. Error t value Pr(>|t|)      YearVar -7.089 1.252 -5.664 3.04e-07 ***     AgeGroupDummy1 142.292 9.211 15.447 < 2e-16 ***     AgeGroupDummy2 185.508 9.211 20.139 < 2e-16 ***     AgeGroupDummy3 218.170 9.211 23.685 < 2e-16 ***     AgeGroupDummy4 255.733 9.211 27.763 < 2e-16 ***     AgeGroupDummy5 278.180 9.211 30.200 < 2e-16 ***     AgeGroupDummy6 300.910 9.211 32.667 < 2e-16 ***     AgeGroupDummy7 282.325 9.211 30.650 < 2e-16 ***     YearVar:AgeGroupDummy2 -1.737 1.770 -0.981 0.3298      YearVar:AgeGroupDummy3 -2.401 1.770 -1.357 0.1792      YearVar:AgeGroupDummy4 -3.772 1.770 -2.131 0.0366 *      YearVar:AgeGroupDummy5 -2.915 1.770 -1.647 0.1040      YearVar:AgeGroupDummy6 -3.587 1.770 -2.026 0.0465 *      YearVar:AgeGroupDummy7 -2.372 1.770 -1.340 0.1846      ---     Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     Residual standard error: 14.97 on 70 degrees of freedom     Multiple R-squared: 0.9946, Adjusted R-squared: 0.9935      F-statistic: 917.9 on 14 and 70 DF, p-value: < 2.2e-16      My interpretation of the interaction term coefficients here is that they represent the difference in slope of YearVar between the corresponding AgeGroup in the interaction term and the YearVar result at the very top. That is, they are something like a period effect across the different age groups. Can anyone see a problem with what I've done here or with my interpretation? This second regression is the closest thing to period effects across distinct age groups that I've been able to muster. Any critiques/new ideas are welcome.


title: Predicting with GAM, using an offset

text: I came in a new post because I’m quite confused with the use of an offset in gam() function (mgcv package). I’m actually working on predicting marine top predator’s distribution using aerial surveys (predators’ observations) and secondary production model outputs. My data are shaped in a 0.25° resolution grid, over large areas such as south-west Indian Ocean. Surveys have been done in several sectors over that grid, that is, all the grid is not sampled. My secondary production covariates are from SEAPODYM model, developed by Lehodey et al. I have a set of about 12 covariates, and wish to find the best combination of them to _predict_ predators’ distribution. That said, because my goal is to predict predators’ distribution, I performed a model selection based on GCV score, after having excluded all combinations of correlated covariates. I am modelling observations, therefore I used a quasipoisson distribution (overdispersed data) with sampled area per pixel as an offset. I made my model selection with the model written as follows:               mod1<-gam(Y ~ offset(log(sampled area))+ covariate1 + covariate2 + covariate3 + covariate4, family=quasipoisson)      Until then, everything worked well, but then I started predicting predators’ observations over my full grid. Doing this with the previous model resulted in predicting only over sampled pixel. Following colleagues’ advices, I predicted using a different writing, thinking it would not change anything in the results:               mod2<-gam(Y ~ covariate1 + covariate2 + covariate3 + covariate4, offset(log(sampled area)), family=quasipoisson)      Indeed, it allows me to predict observation over the entire grid, but I am completely confused since I looked at the summaries of the models. The results are in fact quite different, with both GCV scores and Explained Deviances having higher values:               summary(mod1)     Family: quasipoisson      Link function: log           Formula:     species ~ offset(log(sampled area) + s(covariate1,k=4) + s(covariate2,k=4)          + s(covariate3,k=4) + s(covariate4,k=4)          Parametric coefficients:                 Estimate Std. Error t value Pr(>|t|)     (Intercept) -0.12268    0.08025  -1.529    0.127          Approximate significance of smooth terms:                      edf Ref.df     F  p-value         s(bathy_mnk_B) 2.722  2.951 15.89 4.83e-10 ***     s(meso_mnk_B)  2.995  3.000 17.52 3.52e-11 ***     s(bathy_mnk_P) 2.989  3.000 20.37 6.21e-13 ***     s(mmeso_mnk_P) 2.536  2.855 12.71 6.49e-08 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           R-sq.(adj) =  0.0838   Deviance explained = 16.4%     GCV score = 22.371  Scale est. = 22.188    n = 1491               summary(mod2)     Family: quasipoisson      Link function: log      Formula:     species ~ s(covariate1,k=4) + s(covariate2,k=4) + s(covariate3,k=4) +          s(covariate4,k=4)     Parametric coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  1.30093    0.07884    16.5   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Approximate significance of smooth terms:                      edf Ref.df     F  p-value         s(bathy_mnk_B) 2.802  2.974 17.37 5.18e-11 ***     s(meso_mnk_B)  3.000  3.000 16.35 1.86e-10 ***     s(bathy_mnk_P) 2.988  3.000 27.75  < 2e-16 ***     s(mmeso_mnk_P) 2.460  2.806 12.47 1.12e-07 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           R-sq.(adj) =  0.0938   Deviance explained = 17.5%     GCV score = 31.809  Scale est. = 31.547    n = 1491      Given that I based my model selection on GCV scores, the fact that it changed like that puzzled me, and I am wondering whether or not my model selection was correct. I re-performed a model selection with the second kind of model (Y~X1+X2,offset), and resulting selected models are not always the same than with my first procedure. Seeing this, I am wondering on which writing of the model should I base my modelling and predictions? Does anyone could explain me why those two models are giving such different results? Thanks in advance for your help and advices, Charlotte


title: Is there a way to correct standard errors and/or prediction intervals for multiple comparison after doing backwards selection?

text: It is well known that most model selection algorithms can easily fall into a multiple comparison trap. To quote Friedman: > Consider developing a regression model in a context where substantive theory > is weak. To focus on an extreme case, suppose that in fact there is no > relationship between the dependent variable and the explanatory variables. > Even so, if there are many explanatory variables, the R2 will be high. If > explanatory variables with small t statistics are dropped and the equation > refitted, the R2 will stay high and the overall F will become highly > significant. This is demonstrated by simulation and by asymptotic > calculation. Now, say you've got a multiple regression model based on multiple theories; you expect correlation between $y$ and many of your $X$'s, but you're not entirely sure how those correlations are affected by partial co-linearity between your $X$'s. Further, lets say that you suspect that there might be interactions between some of your $X$'s, but you're not sure. Because of collinearity, you want to estimate a model with all the $X$'s, rather than several separate models. Now, say you were fitting the model to the population, rather that a sample from the population. Suppose that some subset of your $X$'s and your interactions were in fact poor predictors. Since you're in the population, all coefficient estimates are unbiased, and those poor predictions represent real, but weak, correlations. But what if you're not sure about the model? Coefficient estimates suffer from omitted variables bias, or misspecification bias, or probably both. Now go from the population to the sample. You specify the same rich model based on your best guess of theory, expecting that some aspects of your theory are correct, and some are not. If you **do not** select out (say via AIC) some parameters, your model is overfit, and will generalize poorly. If you **do** select out certain terms (say via backwards selection), your model will suffer from a multiple comparison problem -- certain features will be deemed "significant" even if they do not represent real population correlations. So what is an applied guy to do? I've got a complicated dataset with a bunch of correlated variables, and a bunch of reasonably well-justified theories. I want to test whether there is any support for the theories, and I want to make predictions (with associated measures of uncertainty of predictions) based on the most robust model I can construct. I've been doing iterative backwards selection -- dropping terms that lower AIC the most when dropped. **Is there a way to correct standard errors and/or prediction intervals for multiple comparison after doing backwards selection**? (Note that I do NOT want to do ridge regression or lasso, because I am not willing to trade bias for variance. I know that mis-specification can lead to bias, but I am using semi-parametrics to avoid functional form mis- specification on features.)


title: In Excel, how do I plot two rows against each other?

text: So I have data like:               Cost    20  30  10   5     Rating   5   3   2   5      I want to make a chart of rating vs. cost, so the points would be               [(5,20), (3,30), (2,10), (5,5)]      I can't seem to get excel to do anything other than put the two rows as independent series. Am I missing something, or I do have to pivot the data somehow to make it do that? (Actually, I'm using an old-ish version of Numbers.app on OS X, but I'm hoping the concept will be the same. I have access to excel if need be.)


title: Interpreting Residual and Null Deviance in GLM R

text: How to interpret the Null and Residual Deviance in GLM in R? Like we say the for AIC that smaller AIC is better, is there any similar and quick interpretation for the deviances also? Null deviance: 1146.1 on 1077 degrees of freedom Residual deviance: 4589.4 on 1099 degrees of freedom AIC: 11089


title: Dunnett's test in R returning different values each time

text: I am using the R 'multcomp' library (http://cran.r-project.org/web/packages/multcomp/) to calculate Dunnett's test. I am using the script below:               Group <- factor(c("A","A","B","B","B","C","C","C","D","D","D","E","E","F","F","F"))     Value <- c(5,5.09901951359278,4.69041575982343,4.58257569495584,4.79583152331272,5,5.09901951359278,4.24264068711928,5.09901951359278,5.19615242270663,4.58257569495584,6.16441400296898,6.85565460040104,7.68114574786861,7.07106781186548,6.48074069840786)     data <- data.frame(Group, Value)     aov <- aov(Value ~ Group, data)     summary(glht(aov, linfct=mcp(Group="Dunnett")))      Now if I run this script through the R Console multiple times I get very slightly different results each time. Here's one example:                        Simultaneous Tests for General Linear Hypotheses          Multiple Comparisons of Means: Dunnett Contrasts               Fit: aov(formula = Value ~ Group, data = data)          Linear Hypotheses:                Estimate Std. Error t value Pr(>|t|)        B - A == 0 -0.35990    0.37009  -0.972  0.76545        C - A == 0 -0.26896    0.37009  -0.727  0.90019        D - A == 0 -0.09026    0.37009  -0.244  0.99894        E - A == 0  1.46052    0.40541   3.603  0.01710 *      F - A == 0  2.02814    0.37009   5.480  0.00104 **     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     (Adjusted p values reported -- single-step method)      And here's another:                        Simultaneous Tests for General Linear Hypotheses          Multiple Comparisons of Means: Dunnett Contrasts               Fit: aov(formula = Value ~ Group, data = data)          Linear Hypotheses:                Estimate Std. Error t value Pr(>|t|)         B - A == 0 -0.35990    0.37009  -0.972   0.7654         C - A == 0 -0.26896    0.37009  -0.727   0.9001         D - A == 0 -0.09026    0.37009  -0.244   0.9989         E - A == 0  1.46052    0.40541   3.603   0.0173 *       F - A == 0  2.02814    0.37009   5.480   <0.001 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     (Adjusted p values reported -- single-step method)      As you can see, the above two results differ very slightly, but it's enough to move the final group (F) from two stars to three stars, which I find worrying. I have several questions regarding this:   1. Why is this happening?! Surely if you put the same data in each time you should get the same data out.    2. Is there some kind of random number being used somewhere in the Dunnett's calculation?   3. Is this slight variation each time actually a problem?


title: Correlation and colinearity in nonlinear regression?

text: i have a new data set which is basically as bad as the last (same sort of data) and have been asked to try non linear regression on it, with the focus on partition (I will be using boosting and bagging). The dependent variable is continuous however. I require the input variables intact (as in not factorized) as i again require identification of what input variables affect the dependent variable of interest and how. Again the many X input variables are very different in terms of distribution, i have categorical inputs too. I do not see much information in general with regards to non linear regression and linear correlation. Leaving the data set as is causes multiple inclusion of what are strongly correlated variables particularly in bootstrapping (the reason seems obvious due to the nature of bootstrap). to what degree is this a problem? For example; if X1 and X2 are collinear and X1 is Out of Bag with X2 in bag you have effectively not taken that variable out. Is this even an issue?


title: Linear regression for time-series prediction

text: Say we have $N$ time series $X_t^i$ for $i=1...N$and we want to predict a separate time series $Y_t$. Let's consider the following model: $Y_t = \sum_{i} \beta_i X_{t-1}^i $ I am just trying to figure when such a model makes sense   * Does this model have a name? Is there an interpretation in terms of random processes? Are there any known applications?   * Does it make sense to solve the $\beta_i$ by using regression approaches (ridge, lasso, etc) using the training points $(<x_{t-1}^1... x_{t-1}^N>, y_t)$?


title: What is the impact of low predictor variance on logistic regression coefficient estimates?

text: Let's say I am using a logistic model to predict whether it rains (yes or no) based on the high temperature and have collected data for the past 100 days. Let's say that it rains 30/100 days. Furthermore, let's say that 70/100 days the temperature is actually 50 degrees. So there were only 30/100 days where the temperature was not 50 degrees. So in terms of the response variable (rain - yes or no), there is enough variation to develop the model. However, in regard to the predictor (temperature), with so many days where the temperature is 50 degrees, how does a recurring variable "element" (50 degrees) impact the calculations of a logistic model?


title: Estimating a joint probability distribution with antimode

text: I am looking for a certain multivariate probability distribution function to fit to my data, but the usual multivariate normal distribution is unfit, since my data has a dent (antimode) instead of a peak in certain cases. (In other cases the data does have a peak.) My data consists of data points with multiple variables, which are all constrained to be within the interval [0,1]. In the univariate case I would simply use the beta distribution, but the multivariate analogue, the Dirichlet distribution, seems unfit as well; it constrains the input variables to sum up to 1, instead of demanding all the input variables to lie within the interval [0,1]. Furthermore, it seems that the distribution doesn't depend on the covariance matrix - it seems to be unable to distinguish high covariance from low covariance given the same variances. What kind of distribution should I use?


title: Inference on original population based on sample - estimating sampling probability

text: In my research I am working on some data (let’s call it sales data), however we almost have no information on the original sample. My goal is to come up with possible explanations on how it could be sampled by consistency arguments. In more detail, the data provided is composed of sales transactions and it has 2 columns (buyer and receiver). One of the most reasonable ways to sample would be (we believe) to assign a probability (p) to include each transaction or not, and then include a subset of transactions according to random realizations ( sampling without replacement). This corresponds to binomial sampling. Our information on the original sample is limited (almost none), except that I have an estimate on total number of buyers there may be on overall (circa ~2 million). In addition, we assume that buyers are the customers that belong to a network of original data provider, however sellers may belong to other networks, thus our attention is limited to buyers. The consistency argument is based on the following: on the original network if the buyer is highly active (there are more transactions) then the probability that that specific buyer gets included in a sample is higher than a less active customer. Then by counting how many transactions each buyer makes, and grouping them by common number of transactions – we get the number of buyers who make $i$ number of transactions. Then, based on our estimate on the total number of buyers in the network (2 M), my goal to estimate the sampling rate (p). R code:                   set.seed(0);         ##generate random sales transactions              ss1= sample.int( 100, size = 100,replace = TRUE, prob = NULL)         ss2= sample.int(10000, size = 100, replace = TRUE, prob = NULL)              buyer=     sort(ss1,decreasing=TRUE);         receiver = sort(ss2,decreasing=TRUE);                   data= data.frame("buyer"= buyer, "receiver"= receiver);              ##compute the number of sales - each buyer makes               count= data.frame(table(buyer));              buyer.sales = count[order(count$Freq,decreasing=TRUE),];                   df = data.frame(table(buyer.sales$Freq)); ##number of buyers who make 1                 transaction, 2 transactions etc...                    ## We have a guess on total number of buyers - 2M. What would be the sampling rate(p)?      Here are the approaches I’ve tried so far: **Approach 1:** Pose Y_i the number of customers who make “i” calls in the sample (data), and similarly X_i the number of customers who make “i” calls in the original network. We know the Y’s and try to get X’s. Then if the transactions are sampled with a constant p rate, then $Y_i = sum_{k=i}^{k= max(i)} N_k * (choose(k,i))* p^k * (1-p)^{(n-k)}$; If we define the matrix A as follows:                A <- function(i,j,prob) ifelse(i > j, 0,(choose(j,i)*((1-prob)^(j-i))*(prob^(i))) );      The problem would be solve A*x= b, which has an easy solution since A is upper triangular. But for my current problem, the x values that I get do not make any sense – they highly range between plus-minus infinity and ideally I would think that x vector should have some kind of correlation with y’s. I solved the linear system of equalities by minimizing L1, L2 norms and also with the positivity constraints as well. Still, not correlated results- and the sum of x_i’s are highly different than 2M. **Approach 2:** This approach is on subsampling from the sample (data). I thought I could subsample from the sample with a given p, then count the number of buyers I get in that subsample. The probability value –p, which has generates the following ratio ( the number of buyers original / the number of buyers in sample), equal to (the number of buyers in sample / the number of buyers in subsample) would be the most reasonable choice. This approach gave me p~0.02% but I don’t think this is correct. I think approach 1 makes more sense but the results do not suggest that. I should note as well, that the data on hand is composed of monthly transaction data for a period of 12 months. Thus the amount of data I take as the sample is important as well. I tried poisson rate estimation (Estimate the size of a population being sampled by the number of repeat observations) but the amount of data I include, highly effects the estimate of total number of buyers – and I am really stuck. I would think that the original data would have some kind of a power law behavior, but that’s purely a guess without any justification. I am really stuck on which approach makes more sense. Any comments, ideas appreciated.


title: Can we approximate the distribution of S?

text: I want to understand how the sampling distribution of the whole covariance matrix behaves for large $n$. I am trying to use the delta method and multivariate CLT. I am trying to show that when the $X$'s are IID with nonzero finite fourth order moments then the operation (sample variance) takes place, as in this equation: $$\sqrt{n}(s_X^2 - \sigma^2) \underset{n\rightarrow \infty}{\overset{d}{\longrightarrow}} N(0, \mathbb{E}X_1^4 - (\mathbb{E}X_1^2)^2).$$ (For simplicity I have assumed $\mathbb{E}X_i = 0$). Help is solicited in understanding the steps I need to take. I started with studying the asymptotic behavior of the sample variance and then went on to use the delta method for the following function $g(a,b) = b - a^2$.


title: How to interpret the output of predict.coxph?

text: After fitting a coxmodel it is possible to make predictions and retrieve the relative risk of new data. What I don't understand is how the relative risk is computed for an individual and what is it relative to (i.e. the average of the population)? Any recommendations for resources to help understand (I not very advanced in survival analysis so the simpler the better)?


Should false discovery be controlled at the data acquisition level, or should this be at the data interpretation level? I have an experiment in which microarrays were used to quantify the expression of about 30,000 genes (variables) in two groups of biological tissues (group- sizes of 75 and 76). Raw array data was pre-processed to remove background signals and genes whose expression levels are undetectable, and to normalize values across arrays. The final data was then examined using the Mann-Whitney U test to compare gene expression between the two groups to identify differentially expressed genes, with false discovery rate (FDR) controlled with the Benjamini-Hochberg procedure. At FDR <5%, no gene is identified as differentially expressed, and I formally conclude that 'there is no differential expression of genes between the two tissues.' Now, suppose someone is interested in the expression of only one specific gene. Using my pre-processed gene expression data-set and the U test, they compare the two groups for the expression of only this gene and notice a P value <0.05, the commonly used significance threshold in my field of study. As this does not involve multiple testing, there is no false discovery control. Can this observer formally conclude that 'the gene is differentially expressed between the two tissues,' contradicting my conclusion?' Or should the observer have applied false discovery control because such a control has to be applied at the data acquisition level (as per which, data on multiple variables were collected) and not at the data interpretation level (as per which, data for only one variable was analyzed)?


title: How to fix weight in bayesan average

text: I am reading the article of wikipedia about bayesian average. This seems very interesting about on developing rating systems, but I have a doubt. How to fix the weight C of the average?


title: How to get a nice regress model?

text: I am doing a regression analysis which troubled me. My independent variable are 4 interplanetary condition components, and the dependent variable is the latitude of auroral oval boundary. So far, the specific relationship is still unknown in physical principle, what we want to do is to get a model (function expression) from the massive data which shows how these independent variables affect the dependent variable. I used the Matlab statistical toolbox to do the regression analysis, but the results were very bad. The p values of the F statistic and t statistic are very small, but the R² is also very low, about 20%. So how should I improve the R²? Are there good methods? I see that SVM (or LS- SVM) can do regression anaysis, is it a good way to manage the massive data, multiple independent variables regression anaysis? The following are the results:               mdl =      Linear regression model:         y ~ 1 + x1*x2 + x1*x3 + x1*x4 + x2*x3 + x2*x4 + x1^2 + x2^2 + x3^2 + x4^2          Number of observations: 18471, Error degrees of freedom: 18457         Root Mean Squared Error: 2.44       R-squared: 0.225,  Adjusted R-Squared 0.225       F-statistic vs. constant model: 413, p-value = 0        when we add another predictor, i.e., the independent variables become 5, the resuts of the regression analysis are: Linear regression model:               y ~ 1 + x1*x2 + x1*x3 + x1*x4 + x1*x5 + x2*x3 + x2*x4 + x2*x5 + x3*x5 + x2^2 + x3^2 + x4^2 + x5^2      Number of observations: 18457, Error degrees of freedom: 18439 Root Mean Squared Error: 2.21 R-squared: 0.366, Adjusted R-Squared 0.366 F-statistic vs. constant model: 627, p-value = 0


title: Confidence Interval for variance given one observation

text: This is a problem from the "7th Kolmogorov Student Olympiad in Probability Theory": > Given one observation $X$ from a $\operatorname{Normal}(\mu,\sigma^2)$ > distribution with both parameters unknown, give a confidence interval for > $\sigma^2$ with a confidence level of at least 99%. It seems to me that this should be impossible. I have the solution, but haven't read it yet. Any thoughts? I'll post the solution in a couple days. [Follow-up edit: Official solution posted below. Cardinal's solution is longer, but gives a better confidence interval. Thanks also to Max and Glen_b for their input.]


title: Rates affecting rates

text: Stats noob here. I've been trying to find out if there is a technique to perform a certain type of analysis but so far a lot of searching hasn't yielded anything. My question is, if a large peak or trough is identified in an organisation's overall rate, is there a way of identifying which of the organisation's production sites has contributed the most to this rate change? The data describes an organisation and it's performance against a rate based performance measure (let's say % wastage on a washer production line). The dataset contains 4 fields. These fields are:   1. Site (There are 100+ of these)   2. Month   3. Numerator (Number of washers lost to wastage)   4. Denominator (Total number of washers produced) The organisation monitors it's performance based on its overall rate (sum of all numerators / sum of all denominators) by month. This gives a trend line. If there is a sudden large increase/decrease in rate then the organisation wants to target the sites which have most contributed to this rate change. I'm not certain if there is an answer to this question, as one site could see a decrease in denominator and another could see an increase in its numerator value, meaning these sites both show a small change, whilst the organisation as a whole may see a larger scale change when the figures are combined and the overall rate is calculated. I was originally thinking about something like a logistic regression with weighting but wasn't sure if this was quite the way to go...


title: Web page change frequency

text: I want to know I often I need to crawl certain webpages, to make sure I have the newest content as fast as possible. The naïve approach to this, is to simply calculate the average over a certain period of time. However, different web pages are updated at different rates, so this wouldn't be a good fit. I did some googling, and apparently, the Poisson distribution is a good fit for this kind of estimation. My question is, how should I calculate lambda? And it should be constantly updated, correct? **Edit:** John's answer helped me out; for further reference, check this article: Estimating Frequency of Change Thank you


title: Is using the same data for feature selection and cross-validation biased or not?

text: We have a small dataset (about 250 samples * 100 features) on which we want to build a binary classifier after selecting the best feature subset. Lets say that we partition the data into: Training, Validation and Testing For feature selection, we apply a wrapper model based on selecting features optimizing performance of classifiers X, Y and Z, separately. In this pre- processing step, we use training data for training the classifiers and validation data for evaluating every candidate feature subset. At the end, we want to compare the different classifiers (X, Y and Z). Of course, we can use the testing part of the data to have a fair comparison and evaluation. However in my case, the testing data would be really small (around 10 to 20 samples) and thus, I want to apply cross-validation for evaluating the models. The distribution of the positive and negative examples is highly ill-balanced (about 8:2). So, a cross-validation could miss-lead us in evaluating the performance. To overcome this, we plan to have the testing portion (10-20 samples) as a second comparison method and to validate the cross-validation. In summary, we are partitioning data into training, validation and testing. Training and validation parts are to be used for feature selection. Then, cross-validation over the same data is to be applied to estimate the models. Finally, testing is used to validate the cross-validation given the imbalance of the data. **The question is:** If we use the same data (training+validation) used in selecting the features optimizing the performance of classifiers X, Y and Z, can we apply cross-validation over the same data (training+validation) used for feature selection to measure the final performance and compare the classifiers? I do not know if this setting could lead to a biased cross-validation measure and result in un-justified comparison or not.


title: Does the 'reference group' in a Cox proportional hazards model have to exist?

text: My understanding is that a hazard ratio from a Cox proportional hazards model compares the effect on the hazard rate of a given factor to a reference group. Does that actual group have to exist in the data? Say we enroll people in a study of how long before they buy a couch. We right- censor at 3 years. For this example we have two factors: age < 30 or >= 30, whether they own a cat. It turns out the hazard ratio of "owns cat" to the reference group (age < 30, "doesn't own cat") is 1.2, and significant (say p<0.05). Do I know that the data includes people with age < 30 who don't own a cat? For this example it seems a little crazy, but if there are many factors in the model, it becomes easier to imagine that there might be no data for some combination of factors, and that might end up as the reference group if chosen arbitrarily. In fact I think I've had that happen (in R, if it matters).


title: What is a trend line?

text: _Mathematics and Statistics are quite far behind in my educational curriculum, and yet now I seem to need them again, so this set of questions is very basic I suppose:_ ## What is a trend line? What various types of trend lines are there? And what are their respective use or relevance? ### Background: I'm working on a thesis and I am producing some visual rendering of a survey I ran, in which I tried to measure public knowledge of (colonial) history and heritage in a town of Kerala, India. I have data for instance, representing the result of a test, the average result of which is compared to various criteria of the population (age, education, etc.) With regard to criteria such as age-group, income or education, it is also interesting to find out whether their increase improves the results. I am working on this data with LibreOffice, which offers me the possibility to show a "trend line", and I suppose the trend-line is adequate to show the relation between _x_ and _y_ , i.e. criteria and result. However, there are 4 types of trend lines available: _linear, logarithmic, exponential_ and _power_. Which one should I chose? Which difference does it make? How would I figure out a trend line if LibreOffice would not do it for me? Wikipedia brings some light on those questions. Although I am not sure whether I should look at linear regression or technical trend line analysis. Could you please provide a short "Trend Line for Dummies" introduction? Thanks :)


title: As a reviewer, can I justify requesting data and code be made available even if the journal does not?

text: As science must be reproducible, by definition, there is increasing recognition that data and code are an essential component of the reproduciblity, as discussed by the Yale Roundtable for data and code sharing. In reviewing a manuscript for a journal that does not require data and code sharing, can I request that the data and code be made available   1. to me at the time of review   2. publicly at time of publication (the journal supports supplements) also, how might I phrase such a request? * * * _update_ : although I am interested in the general case, this particular case consists of a meta-analysis with all previously published data, and the code is simple linear models in SAS _side note_ the ability to make cross-study inference (as is the goal of meta- analysis) would be greatly enhanced if more studies provided raw data _update 2_ : I requested the data and code from the editor for purposes of review, the editor considered the request reasonable, and I have received the requested material (sufficient but with cryptic variable names, no metadata, and few inline comments) within a day.