title: Way to correct sample selection bias with unknown selection?

text: I would greatly appreciate some advise on a statistical problem that haunts me. Suppose you wish to estimate the effect of $x$ on $y$, but the probability to observe $\\{y_i, x_i\\}$ also depends on the effect of $x$. Say, if $x_i > z_i$, then you don't observe $\\{y_i, x_i\\}$. If $x_i < z_i$, then you do. If this was a standard textbook example of sample selection models, we would observe $z$, and could adjust, but in my case I don't. I have no information on how many individuals did not make it into my sample. To be concrete: I don't know they ever existed. That is, they died before the sample was collected, and thus do not appear in my sample. My sample consists of survivors, so to say. Following the comments below, I augment my question to make it a bit more precise. Suppose a group of $m$ people was exposed to a disaster, say. 20 years later I want to see how much this experience shines through today, but those that were most severely affected died along the way. I nevertheless collect data on all remaining $n$ survivors. Here, $n < m$. Unfortunately, I have no information whatsoever about the exposed population and I do not observe their characteristics $\\{x_i, y_i\\}$. I do not know $m$, and I do not know how many of the initial population survived. All I have is my group of survivors. If I want to estimate the effect of the disaster on some outcome, that estimate will be an biased because the whole sample is conditioned on survival. Does anyone know any solution to this problem? Are there any methods out there to, say, compute bounds maybe? Any hints are highly appreciated.


title: Can you ever have known parameters?

text: Maybe a bit of a philosophical question - but can you ever truly have known parameters in data? I have a set of data for which the dataset is complete, but the parameters will still be estimates i believe?


title: Suggestions for tools or strategies to visually compare results

text: This is not exactly a technical question related to any method, but this seemed like the best forum to ask this. The question I am asking is more along the lines of experimental data management. I often run into situations where I have to compare several sets of similar results VISUALLY, for example experiments run with different combinations of methods/datasets/parameters etc. The quantities that I am comparing are numbers (e.g. accuracy, AUC ROC scores etc.). An example of this kind of data would be               method1/dataset1/param1 0.01     method1/dataset1/param2 0.20     ....     method9/dataset5/param4 0.40      If I want to compare the results, the performance of method1 vs method2 (on dataset5). I would like to filter out everything else and keep only the relevant information, or something along those lines, so that the clutter is minimized when visually inspecting the results. Now it is possible to manually do this for every comparison or put the results in a database and run the relevant queries or by graphing those values. What I am asking is, is there an approach that you use to streamline the process, so that the overhead is minimal?


title: 25 percent chance of winning, but only 3 wins in the past 20 events

text: If anyone is interested, this question is prompted by the recent NBA Lottery results. In the past 20 years, the team with the worst record in the NBA has won the lottery exactly 3 times. The team with the worst record every year has a 25% chance of winning the number one overall pick. Of course, many people will talk about how the lottery is rigged and cite this statistic as an example of why. However, it doesn't seem like this is a statistically significant deviation from the expected value of 5 wins. Does anyone know which test I should use to find out whether the worst record team winning the overall pick only 3 or fewer times in the past 20 years is statistically significantly unlikely? Mostly interested in the process here, since I should be familiar with this. My first instinct tells me to use the binomial probability function and thus, $\sum_{x=0}^3 {n \choose x} (.25)^3 (.75)^{(20-x)}$, which yields around 22.5%, which is greater than an $\alpha$ of say, 0.05, so this isn't statistically significant. However, I feel like there is a better way to do this involving the expected value of 5. Does anyone have a clue on what I'm looking for here?


title: What's the difference between these ordered factors?

text:
str(test)     'data.frame':   767 obs. of  2 variables:     $ datefield: Ord.factor w/ 59 levels "1984-04-01"<"1984-07-01"<..: 1 2 3 4 5 6 7 8 9 10
         $ somevar : num  43.7 55.6 43.5 54.1 42.8 ...          > str(Italy)     'data.frame':   1008 obs. of  2 variables:     $ year: Ord.factor w/ 48 levels "1951"<"1952"<..: 1 1 1 1 1 1 1 1 1 1 ...
         $ gdp : num  8.56 12.26 9.59 8.12 5.54 ...      I turned the datefield into an ordered factor because I try to reproduce an example with own data. Now I wonder what the difference between is, or put it differently, what does 1 2 3 4 5 6 7 8 9 10 mean opposed to 1 1 1 1 1 1 1 1 1 1 ... ?


title: Approaches for extracting single-year age estimates from age-buckets

text: I'm trying to estimate a binary function of $Y$~$Bernouli(p_{i})$ where $p_{i}=f(age)$ where f is an unknown continuous function. There are obviously a lot of methods to try to reconstruct f (splines, Gaussian processes, etc), but what if my age observations are collapsed into buckets (18-29, 30-44,etc)? Obviously, it's impossible to reconstruct f even with infinite data. And reconstruction of f is going to depend fairly strongly on assumptions about the nature of f. But this seems like it should be a common problem and I'm curious what approaches people have used/would use. edit 7/14/2014:I should clarify that in this case, I have a pretty good distribution of the overall distribution of age in the population.


title: variable depending on 2 variables conditional probability

text: Previosuly I have asked how to compute some conditional probabilities, but I am missing this particular case: lets say we now have 3 variables: $T$, $L$, $E$:                 T     L        \   /          E      So $E$ depends on $T$ and $E$ depends on $L$ I have these probabilities: $P(T) =$ 0.0104   $P(L) =$ 0.055   $P(E) =$ 0.0648 $P(E|L,T) =$ 1.0   $P(E|L,\lnot T) =$ 1.0   $P(E|\lnot L,T) =$ 1.0   $P(E|\lnot L,\lnot T) =$ 0.0 How is modified the probability of $T$ and $L$ when we observe that $P(E)$ = 1?


title: Estimating the number of registered interaction designers in a given territory

text: I plan to conduct a research study about Interaction Designers. Prior to starting on sampling methods I am finding it very hard to get my hands on some demographics and statistics, e.g. How many people are registered to be employed as Interaction Designers, User Experience Designers etc. in a particular territory (UK in my case). Any help regarding this would be highly appreciated.


title: Estimating sample size required for optimal performance of latent semantic indexing?

text: Latent semantic indexing seems to work well; e.g. it is independent of language, etc. However, it appears to use the similarity of frequencies of terms in the corpus to categorize them. If this understanding is correct, is there a way to measure the size of the dataset that will give optimal performance?


title: Choosing between Latent Dirichlet Allocation and Association Rule Mining

text: Data Mining and Statistical Analysis has a general discussion on stats vs data mining. If I may, narrow down the question a bit - are there any general demarcations that allows you to decide which approach is more suited for using - LDA or ARM? (other than appreciation for the maths in the former)


title: p-value reported as 3.2531E-129 in excel does not work in CHIINV function

text: When I insert the p-value `3.2531E-129` into excel's `CHIINV` function I get `#NOMBRE` as a result.   How do I find or report the chi-square value if this doesn't work?   Someone suggested I wrote out the 0 when the value was `6.1203E-05` to `0,0000061203`.   So, how can I change `3.2531E-129` to `0.000000......32531`?   Thank you for your help.


title: Finding minimax estimator $\theta$

text: Suppose $X_1,\ldots,X_n$ is a random sample from $N(\theta,1)$ so that $\theta\geq a$ . How can I calculate minimax estimator $\theta$ under loss function $L(\delta,\theta)=(\delta-\theta)^2$


title: Hierarchical variable: checking whether two subgroups have to be distinguished or can be pooled together

text: I am doing a study to analyze the effects of knocking out a particular gene on a certain behaviour in flies. After a log transformation, my observations have a normal distribution so I'm going to do a three-way ANOVA, with age as the first, sex as the second and genotype (mutant vs wild-type) as the third variable, also looking at interactions. The problem is, my mutants come from two slightly different strains that have never been shown to differ behaviourally but I should still check first to determine whether I need to distinguish between the two in the subsequent analysis or not. I'm wondering whether this could be done with a nested ANOVA (all the examples of nested ANOVAs I've seen look quite different from what I have, which is why I haven't tried it, but I'm probably just not getting it). So far, what I've done is run a three-way ANOVA (type III ss) on just the mutants, with age as the first, sex as the second and genotype (mutant strain 1 vs mutant strain 2) as the third variable. Age, sex and the interaction of age and sex are both significant while mutant strain isn't, nor are any of the interaction terms containing it. My question is, can I just report the fact that the mutant strain has no significant effect and then only worry about the age and the sex when I'm running the ANOVA on all my flies (this time without distinguishing between the two strains of mutants) or do I have to report it and take it into account now too? And what if mutant strain had come out significant? Would that change anything as to whether or not I should report the significance of age and sex at this stage? I would have just done a t-test to compare the two mutant strains but I was afraid there might be an effect that was masked by either age or sex as both seem to have an effect on the behaviour I'm looking at. Thank you so much in advance, I know this is a newbie question (it's my very first real statistical analysis but I'm enjoying every second of it) but if any of you had a second to answer I'd be very grateful.


title: Reproducing table 18.1 from "Elements of Statistical Learning"

text: Table 18.1 in the Elements of Statistical Learning summarizes the performance of several classifiers on a 14 class data set. I am comparing a new algorithm with the lasso and elastic net for such multiclass classification problems. Using `glmnet` version 1.5.3 (R 2.13.0) I am not able to reproduce point 7. (the $L_1$-penalized multinomial) in the table, where the number of genes used is reported to be 269 and the test error is 13 out of 54. The data used is this 14-cancer microarray data set. Whatever I have tried, I get a best performing model using in the neighborhood of 170-180 genes with a test error of 16 out of 54. Note that in the beginning of Section 18.3, on page 654, some preprocessing of the data is described. I have contacted the authors -- so far without response -- and I ask if anybody can either confirm that there is a problem in reproducing the table or provide a solution on how to reproduce the table.


title: Correct variance estimation for random sample post-stratification weights

text: I'm dealing with a random sample and I have a question about its poststratification weights. Let's assume that:   * this is a random sample of addresses, but post-stratification weights have been calculated (using Census estimates) because the sample doesn't match the target population with respect to some important variables.   * there are no major non-response problems.  Should the post-stratification weights be used:   1. as probability weights (for example using [pweight] in Stata or Lumley's survey package in R) when estimating variance, SE, etc... **OR**   2. should they simply be used to re-adjust the sample (changing the distribution but not the total sample size) which can then be analysed with random sample variance estimation procedures? I'd opt for the first option (probability weights, but I may be wrong). Most discussions I find are about clustered samples or more generally samples that were not designed as random in the first place. Other discussions show how post-stratification weights are combined with other types of weights (e.g., to correct for missing data) to produce the final probability weights, but they don't discuss post-stratification in itself. In a survey with good response rates, I see post-stratification as a way to correct non-response bias when this is not very clearly driven by some characteristics and to make the population representative at the same time, two things which are I think often fairly hard to distinguish in some surveys. In other words even when they are just issued as just post-stratification weights and no major non-response problem is mentioned they should be used as probability weights when estimating standard errors etc....


title: Iterative PCA R

text: I have a largish data set (400,000 variables of 1000 samples). I would like to identify what is the best set of these variables for capturing most of the variance between samples. What's the best way to perform an iterative principal component analysis that will perhaps go through X rounds of between Y and Z variables each time and result in a list of which round leads to the fewest number of PCs to capture a certain percentage of variance?


title: When can the variance be written as a function of expectation?

text: 1. What kinds of distribution have variance as a function of expectation, i.e. let $X$ be a random variable of a distribution, $Var(X) = f(E(X))$ for some function $f$ ? Sufficient conditions including typical examples of distribution (such as exponential family distributions, or even natural exponential family distributions), and/or necessary conditions will be fine. I saw this is assumed in Generalized Estimating Equations. In particular, in the natural dispersed exponential family with pdf $f (y; θ, φ) = \exp\left(\frac{yθ − b(θ)}{a(φ)}+ c(y, φ)\right)$, is there always a function $V( · )$ such that $Var(X) = V(E(X))a(φ)$? See this note.   2. Does "variance structure or variance function" refer to that the variance can be written as a function of the mean?   3. If $E(X) = f(\theta)$, and $Var(X) = f'(\theta)$, is $Var(X)$ a function of $E(X)$? One example that satisfies the condition is a canonical exponential family distribution.


title: How to test for "parallelness" of two distributions?

text: 23 patients suffering from bruxism (grinding of teeth) have been given a treatment. Their baseline bruxism has been measured for three nights, then they were treated for 23 nights, then again the "post" baseline was measured for three nights. The raw data looks like this:               patient night   score   duration    minutes of sleep     1       1       20      25          480     1       2       28      67          480      The values are:   * patient = patient number   * night = night (1 to 26)   * score = total number of bruxism attacks during that night   * duration = total seconds of bruxism (sum over all attacks) for that night   * minutes of sleep = the observed timespan (rounded to a quarter hour) In all, there are 598 datasets, 26 consecutive nights for 23 different patients. The interesting variables have been defined as: nph = number of bruxism attacks per hour (score * 60 / minutes of sleep) sph = duration of bruxism in seconds per hour (duration * 60 / minutes of sleep) The nightly means over all patients are:               nph <- c(8.808858, 8.364923, 11.932373, 9.108704, 8.135258, 6.886013, 6.379688, 6.034062, 5.731728, 5.823831, 5.128421, 4.526408, 5.145101, 5.860569, 7.691637, 5.645240, 5.151750, 7.496109, 5.375605, 5.994595, 5.650951, 6.459269, 4.490204, 7.891916, 7.684742, 7.583042)     sph <- c(23.528665, 16.382689, 29.492815, 23.752084, 19.966185, 16.104159, 15.135596, 13.218658, 11.626881, 11.321739, 13.418337, 8.391212, 9.062977, 10.665424, 9.535756, 9.925929, 9.040313, 14.635182, 8.549451, 10.346057, 9.742318, 16.437902, 7.443140, 14.204343, 12.648306, 13.401150)      These are mean number of bruxism attacks per hour (nph) and total duration of bruxism in seconds per hour (sph), over 26 nights, with treatment from night 4 to 23. A plot of these two variables shows that their progress over the 26 nights is roughly "parallel", indicating that both react in roughly the same way to the treatment. ![enter image description here](http://i.stack.imgur.com/M2OJG.png) The filled circles are the three nights pre and three nights post treatment. Due to the scaling the maximum value for both variables fall on the same point (third night). My question is: How can I analyse how "parallel" these two "curves" are?


title: 'systemfit' package and systems of non-linear equations/regression

text: I am trying to estimate a system of non-linear equations using 'systemfit' package in R. I have had issues with it. The two equations share the same parameters i.e. "sigma", "al" and "ae". I expect them to take the same values across the two equations. By theory, I expect 'sigma's to be the same. Though, I am getting different values for the coefficients. I am aware that 'systemfit' allows you to impose linear restrictions on coefficients, but it seems impossible to do this with 'nlsystemfit'.               gammae=exp(-0.7)      gamma0=exp(-0.7)      #beta=0.6      quale=q*e           mod1=lsesv ~ log(gammae)+((sigma-1)/sigma)*(ae*trend1620-0.6*al*trend1620+log(quale)-0.6*log(labour)-(0.4)*log(kapital))      mod2=loutput~log(gamma0)+(sigma/(sigma-1))*log(((exp(al*trend1620)*labour^0.6)*(kapital^(0.4)))^((sigma-1)/sigma)+gammae*((exp(ae*trend1620)*quale)^((sigma-1)/sigma)))      eqSys=list(mod1,mod2)           start.values=c(ae=0.0016,al=0.0016, sigma=0.1)      model.ols=nlsystemfit("OLS",eqSys,start.values,data=uk,maxiter=1000)      # Estimated coefficients               >model.ols$b*1                 ae         al      sigma          0.01647120 0.04431699 0.50860670      >model.ols$b*2             ae         al      sigma           0.03294239 0.08863398 1.01721341       In the above result, the coefficient estimates of equation 1 are half those of equation 2. I was wondering if there is any procedure that allows me to equate the respective coefficeint estimates via restrictions or alternative procedures. Thank you, mrrox


title: Seasonal dummies significance issues

text: I have a question concerning the significance of the seasonal dummies in my ARIMA-model (I do not use seasonal differencing or seasonal AR/MA as I have quite regular seasonality and I get better performance results): What to do, if say 1 or 2 are significant but the rest are not? This happens to me particularly when I use monthly dummies. Shall I exclude the insignificant dummy variables? If I do so (say I exclude from 11 monthly dummies the dummies for February, March), this will lead to a change again in the significance and coefficients of the remaining monthly dummy variables significantly! I could not find any information on that in my standard forecasting textbooks. * * * Finally I found a reference which I wanted to share to conclude the post: > The exclusion of some seasonal dummies because their estimated coefficients > have low t-scores is not recommended. Preferably, testing seasonal dummy > coefficients should be done with the F-Test instead of with the t-test > because seasonality is usually a single compound hypothesis rather than 3 > (or 11 with monthly data) individual hypothesis having to do with each > quarter (or month). To the extent that a hypothesis is a joint one, it > should be tested with the F-Test. If the hypothesis of seasonal variation > can be summarized into a single dummy variable, then the use of the t-test > will cause no problems. Often, where seasonal dummies are unambiguously > called for, no hypothesis testing at all is undertaking. Studenmund, A.H. and Cassidy, H.J. (1997). Using econometrics: A practical guide. The Addison-Wesley series in economics, 3rd ed., p. 257, Addison- Wesley, Reading, Mass.


title: "Interestingness" function for StackExchange questions

text: I am trying to put together a data-mining package for StackExchange sites and in particular, I am stuck in trying to determine the "most interesting" questions. I would like to use the question score, but remove the bias due to the number of views, but I don't know how to approach this rigorously. In the ideal world, I could sort the questions by calculating $\frac{v}{n}$, where $v$ is the votes total and $n$ is the number of views. After all it would measure the percentage of people that upvote the question, minus the percentage of people that downvote the question. Unfortunately, the voting pattern is much more complicated. Votes tend to "plateau" to a certain level and this has the effect of drastically underestimating wildly popular questions. In practice, a question with 1 view and 1 upvote would certainly score and be sorted higher than any other question with 10,000 views, but less than 10,000 votes. I am currently using $\frac{v}{\log{n}+1}$ as an empirical formula, but I would like to be precise. How can I approach this problem with mathematical rigorousness? In order to address some of the comments, I'll try to restate the problem in a better way: Let's say I have a question with $v_0$ votes total and $n_0$ views. I would like to be able to estimate what votes total $v_1$ is most likely when the views reach $n_1$. In this way I could simply choose a nominal value for $n_1$ and order all the question according to the expected $v_1$ total. * * * I've created two queries on the SO datadump to show better the effect I am talking about: Average Views by Score Result: ![Views by Score](http://i.stack.imgur.com/AjwgO.png) Average Score by Views (100-views buckets) Result: ![Score by Views](http://i.stack.imgur.com/2dVzr.png) * * * The two formulas compared Results, not sure if straighter is better: ($\frac{v}{n}$ in blue, $\frac{v}{log{n}+1}$ in red) ![Formulas](http://i.stack.imgur.com/v8RVc.png)


title: Optimize a boolean function

text: I have some data that a downstream system needs an optimized function of boolean logic for. Essentially I have data similar to:               User cat1 cat2 cat3 cat4       1    0    0    0    1       2    1    0    0    1       3    0    1    1    0        I must optimize cat4 as a function like this: `(cat1 or cat2 or cat3)`. "Or" is the only operation that I can use. Does anyone know a particular technique / strategy for optimizing a union of categories in this way? (Note, I will have many categories.)


title: Regression Model - forcing a coefficient/effect to be 1

text: Currently we have a linear model which includes 3 independent variables and the dependent variable "Y" which is the predicted values (range from -0.5 to 2000) from the model. Now we want to leverage the same model by just removing one of the predictor (lets call that as X1). One of my colleague suggested to run a model by have a new dependent as (Y-X1) against the two predictors. Is this a correct approach or is there any better way like "offset" technique for continuous dependent variable? I really appreciate your help in advance.


title: Code for evaluating perplexity of an LDA topic model

text: I'm working on a research project where we are trying to assess model fit improvement of our LDA model (the basic Blie et al., 2003 LDA approach). I'm looking for code that allows us to assess the perplexity of our LDA model as we tune it. I know the R Package "topicmodls" has a perplexity function but my PI wants other options outside of R. Any recommendations? thanks all!


title: Test for treatment change

text: I have a study where the goal is to see if a new diagnostic tool provides any addition to the old procedure. A set of patients have been diagnosed and classified into 7 categories, each corresponding to a given treatment. The new diagnostic tool is then used, and the patients are diagnosed and classified again. Patients are then classified into two categories:   1. the suggested treatment remains the same; or    2. the suggested treatment is changed. It seems as if this should be a simple question, but still I can't figure out what kind of test I could use to see if the diagnostic tool provides any additional information leading to revision of treatment. Given that there is a test, I would also like to do some sample size calculations for the next part of the study.


title: How does ggplot compute confidence intervals for regressions?

text: The R plotting package ggplot2 has an awesome function called stat_smooth for plotting a regression line (or curve) with the associated confidence band. However I am having a hard time figuring out exactly how this confidence band is generated, for every time of regression line (or "method"). How can I find this information?


title: Influential residual vs. outlier

text: First, I should state that I have searched on this site for the answer. I either didn't find a question that answered my question or my knowledge level is so low I didn't realize I already read the answer. I am studying for the AP Statistics Exam. I have to learn linear regression and one of the topics is residuals. I have a copy of _Introduction to Statistics and Data Analysis_ on page 253 it states. > Unusual points in a bivariate data set are those that fall away from most of > the other points in the scatterplot in either the $x$ direction or the $y$ > direction > > An observation is potentially an **influential observation** if it has an > $x$ value that is far away from the rest of the data (separated from the > rest of the data in the $x$ direction). To determine if the observation is > in fact influential, we assess whether removal of this observation has a > large impact on the value of the slope or intercept of the least-square > line. > > An observation is an **outlier** if it has a large residual. Outlier > observation fall far away from the least-square line in the $y$ direction. Stattreck.com states four methods to determine an outlier from residuals: > Data points that diverge in a big way from the overall pattern are called > outliers. There are four ways that a data point might be considered an > outlier. > >   1. It could have an extreme X value compared to other data points. >   2. It could have an extreme Y value compared to other data points. >   3. It could have extreme X and Y values. >   4. It might be distant from the rest of the data, even without extreme X > or Y values. > These two sources seem to conflict each other. Could anyone help clear up my confusion. Also, how does one define extreme. The AP Statistics uses the rule if the data point is outside of (Q1-1.5IQR,Q3+1.5IQR) the it is an outlier. I don't know how to apply that from just a graph off the residuals.


title: Bootstrapping hierarchical/multilevel data (resampling clusters)

text: I am producing a script for creating bootstrap samples from the `cats` dataset (from the `-MASS-` package). Following the Davidson and Hinkley textbook [1] I ran a simple linear regression and adopted a fundamental non-parametric procedure for bootstrapping from iid observations, namely **pairs resampling**. The original sample is in the form:               Bwt   Hwt          2.0   7.0     2.1   7.2          ...          1.9    6.8      Through an univariate linear model we want to explain cats hearth weight through their brain weight. The code is:               library(MASS)     library(boot)               ##################     #   CATS MODEL   #     ##################          cats.lm <- glm(Hwt ~ Bwt, data=cats)     cats.diag <- glm.diag.plots(cats.lm, ret=T)               #######################     #   CASE resampling   #     #######################          cats.fit <- function(data) coef(glm(data$Hwt ~ data$Bwt))      statistic.coef <- function(data, i) cats.fit(data[i,])           bootl <- boot(data=cats, statistic=statistic.coef, R=999)      Suppose now that there exists a clustering variable `cluster = 1, 2,..., 24` (for instance, each cat belongs to a given litter). For simplicity, suppose that data are balanced: we have 6 observations for each cluster. Hence, each of the 24 litters is made up of 6 cats (i.e. `n_cluster = 6` and `n = 144`). It is possible to create a fake `cluster` variable through:               q <- rep(1:24, times=6)     cluster <- sample(q)     c.data <- cbind(cats, cluster)      I have two related questions: How to simulate samples in accordance with the (clustered) dataset strucure? That is, **how to resample at the cluster level?** I would like to sample the clusters with replacement and to set the observations within each selected cluster as in the original dataset (i.e. sampling with replacenment the clusters and without replacement the observations within each cluster). This is the strategy proposed by Davidson (p. 100). Suppose we draw `B = 100` samples. Each of them should be composed by 24 possibly recurrent clusters (e.g. `cluster = 3, 3, 1, 4, 12, 11, 12, 5, 6, 8, 17, 19, 10, 9, 7, 7, 16, 18, 24, 23, 11, 15, 20, 1`), and each cluster should contain the same 6 observations of the original dataset. How to do that in `R`? (either with or without the `-boot-` package.) Do you have alternative suggestions for proceeding? The second question concerns the initial regression model. Suppose I adopt a **fixed-effects model** , with cluster-level intercepts. **Does it change the resampling procedure** adopted? [1] Davidson, A. C., Hinkley, D. V. (1997). _Bootstrap methods and their applications_. Cambridge University press.


title: In a multi-level model, what are the practical implications of estimating versus not-estimating random effect correlation parameters?

text: In a multi-level model, what are the practical and interpretation-related implications of estimating versus not-estimating random effect correlation parameters? The practical reason for asking this is that in the lmer framework in R, there is no implemented method for estimating p-values via MCMC techniques when estimates are made in the model of the correlations between parameters. For example, looking at this example (portions quoted below), what are the practical implications of M2 versus M3. Obviously, in one case P5 will not be estimated and in the other it will. Questions   1. For practical reasons (the desire to get a p-value through MCMC techniques) one might want to fit a model without correlations between random effects even if P5 is substantially non-zero. If one does this, and then estimates p-values via the MCMC technique, are the results interpretable? (I know @Ben Bolker has previously mentioned that "combining significance testing with MCMC is a little bit incoherent, statistically, although I understand the urge to do so (getting confidence intervals is more supportable)", so if it will make you sleep better at night pretend I said confidence intervals.)   2. If one fails to estimate P5, is that the same as asserting that it is 0?   3. If P5 really is non-zero, then in what way are the estimated values of P1-P4 affected?   4. If P5 really is non-zero, then in what way are the estimates of error for P1-P4 affected?   5. If P5 really is non-zero, then in what ways are interpretations of a model failing to include P5 flawed? Borrowing from @Mike Lawrence's answer (those more knowledgeable than I are free to replace this with full model notation, I'm not entirely confident I can do so with reasonable fidelity): M2: `V1 ~ (1|V2) + V3 + (0+V3|V2)` (Estimates P1 - P4) M3: `V1 ~ (1+V3|V2) + V3` (Estimates P1-P5) _Parameters that might be estimated:_ **P1** : A global intercept **P2** : Random effect intercepts for V2 (i.e. for each level of V2, that level's intercept's deviation from the global intercept) **P3** : A single global estimate for the effect (slope) of V3 **P4** : The effect of V3 within each level of V2 (more specifically, the degree to which the V3 effect within a given level deviates from the global effect of V3), while enforcing a zero correlation between the intercept deviations and V3 effect deviations across levels of V2. **P5** : The correlation between intercept deviations and V3 deviations across levels of V2 Answers derived from a sufficiently large and broad simulation along with accompanying code in R using lmer would be acceptable.


title: F test and t test in linear regression model

text: F test and t test are performed in regression models. In linear model output in R, we get fitted values and expected values of response variable. Suppose I have height as explanatory variable and body weight as response variable for 100 data points. Each variable (explanatory or independent variable, if we have multiple regression model) coefficient in linear model is associated with a t-value (along with its p value)? How is this t-value computed? Also there is one F test at the end; again I am curious to know about its computation? Also in ANOVA after linear model, I have seen a F-test. Although I am new statistics learner and not from statistical background, I have gone through with lots of tutorials on this. Please do not suggest for going me with basic tutorials as i have already done that. I am only curious to know about the T and F test computation using some basic example. Thanks !!


title: How do we usually select the best combination of parameters of a machine learning model (for a given dataset)?

text: Am I wrong, or the standard way of optimizing a machine learning model is by evaluating the algorithm over the (initial) dataset for all possible combinations of parameters, and then pick up the one that gives the 'best' evaluation score? If we still don't get an acceptable score, we should then slice and dice the (initial) dataset and start the process from the beginning, but this time applying it on the new dataset? Does this sounds right to you or I am missing something here?


title: Compare the slopes and intercepts between 2 variables in SPSS

text: There are 4 groups, each one having two 2 variables, X and Y. I am trying to understand the regression between them. For X variable I have 2 categories, positive and negative. I am not understanding how to compare the slopes and intercepts between X (positive and negative) and Y. Do I have to invert the negative values?


title: Time-delayed neural networks - Learning the "max" layer

text: In Unified Architecture for NLP paper time-delayed neural network proposed as a way to deal with variable length input. Input window slides over sequence and label each output with "time", then next layer takes max over time. But no explanation given how to back-propagate error through this max layer. We able to compute error only for specific "time", potentialy different for each max layer output node. Reading original paper for TDNN didn't clarify anything. For me it looks like completly different NN architecture without any mention of max layer. Can someone explain or point to explanation? My own googling skills failed.


title: case down-weighted to emulate equal presence and absences

text: In GAM-based logistic regression pseudo-absence data is necessary to fit the model with corresponent presence data. In many cases, pseudo-absence data are generated 10 times more than presence records in order to ensure sufficient landscape variety in the pseudo-absences. To ensure that the presences and pseudo-absences contributed equally to the model and avoid bias, I would like to know how that case-weighting option can be get with R. For example, I want weight each pseudo-absence as n/1000, where n is the number of presence sites and 10000 are pseudo-absences sites generated randomly.


title: When would maximum likelihood estimates equal least squares estimates?

text: Under what situations would MLE (Maximum Likelihood Estimates) equal LSE (Least Squares Estimates)? I got an impression that under norm 2 ($L_{2}$), MLE and LSE are equal. For example, the process of solving $\mathrm{min}||y−Ax||_{2}$ is actually the MLE estimation of parameter $A$ for random variable $y=Ax+\epsilon$ where $x$ and $\epsilon$ is normal. However, is that generally true that the minimization problem under $L_{2}$ norm is the same as maximum likelihood estimation? For example, consider a quadratic function $f(X) = XAX$, can minimizing the distance between $f(X)$ and some value $Y$ under $L_{2}$ can be solved by MLE?


title: Free internet or downloadable resources for sample size calculations

text: Today I noticed this question, and I thought it would be helpful if we had a thread that listed resources that people could conveniently access for power analysis / sample size calculations, perhaps analogous to this thread: Resources for learning R.


title: Variables are often adjusted (e.g. standardised) before making a model - when is this a good idea, and when is it a bad one?

text: In what circumstances would you want to, or not want to scale or standardize a variable prior to model fitting? And what are the advantages / disadvantages of scaling a variable?


title: How to determine correlation between stationary and non-stationary time series

text: I have three time series of economic data based on quarterly observations; A, B and C, and I would like to ascertain the correlation (or not) between A and C as well as the correlation between B and S. The 1st order differenced time series of B passes the ADF and KPSS tests for stationarity. Series A & C fail the stationarity tests even after trying detrending (based on linear regression and centred moving averages of order 4) as well as 1st and 2nd order differencing. What is the best way to ascertain the extent of correlation between A & C and between B & C under these circumstances? Below are the time series data (tab delimited):               A   13,603  15,062  (22,984)    14,704  14,285  15,585  (17,460)    21,145  20,926  28,117  (6,524) 31,190  25,610  34,311  (21,376)    14,140  3,416   26,526  (21,159)    30,874  51,579  50,426  (19,874)    52,980  30,338     B   (92,345)    19,415  13,045  104,693 214,196 (39,180)    (29,979)    112,499 5,914   60,787  92,253  124,716 23,638  362,566 (66,896)    209,127 103,986 (13,418)    389,962 (161,400)   177,945 (36,645)    148,722 189,477     C   819,019 716,641 767,830 1,177,339   1,254,122   1,254,122   985,382 716,641 614,264 806,222 819,019 844,613 1,018,655   1,108,235   1,261,801   1,474,234   1,412,807   1,678,988   1,638,037   1,228,528   1,023,773   1,279,717   1,279,717   1,023,773


title: Gibbs sampling to produce posterior pdf

text: Suppose we have the following classical normal linear regression model: $$y_i = \beta_1 x_{1i} + \beta_2x_{2i} + \beta_3x_{3i} + e_i$$ where $e_{i} \sim iid.N(0, \sigma^2)$ for all $i = 1, 2, \cdots, n$ and $x_{1i} = 1$ for all $i = 1, 2, \cdots, n$. Assume that we have known data values for both $x_{2i}$ and $x_{3i}$ for all $i = 1, 2, \cdots, n$. Defining $\boldsymbol{\beta} = (\beta_1, \beta_2, \beta_3)'$ and assuming a non-informative prior of the form $p(\boldsymbol{\beta}, \sigma) \propto \frac{1}{\sigma}$, then we can show that the conditional posterior pdf for $\boldsymbol{\beta}$ and $\sigma$, that is, $p(\boldsymbol{\beta}|\sigma, \mathbf{y})$ and $p(\sigma|\boldsymbol{\beta},\mathbf{y})$ are normal and inverted gamma, respectively. The question is: Use a Gibbs Sampler and estimate the posterior pdf of the parameter function: $\displaystyle \psi = \frac{\beta_2 + \beta_3}{\sigma^2}$. Now I have run a Gibbs sampler (in R) and after a burn in period of 100 draws, I have obtained 1000 draws of $(\boldsymbol{\beta}^{(i)}, \sigma^{(i)})$, that is, I have a sample of $(\beta_1^{(i)}, \beta_2^{(i)}, \beta_3^{(i)}, \sigma^{(i)})$ for $i = 1, 2, \cdots, 1000$, how can I use these draws to produce an estimate of the posterior pdf of $\psi$? In other words, how can I estimate $p(\psi|\mathbf{y})$? * * * EDIT: I'm new to MCMC logarithms. I do understand what you mean but I am still not too sure how to use it in the context of this question. From what I've learnt so far, say we have $\boldsymbol{\theta} = (\theta_1, \theta_2)'$ and $p(\boldsymbol{\theta}|\mathbf{y}) = p(\theta_1, \theta_2|\boldsymbol{y})$ is the joint posterior, then the marginal posterior of $\theta_1$ is given by $p(\theta_1|\mathbf{y}) = \int_{\theta_1} p(\theta_1|\theta_2,\mathbf{y})p(\theta_2|\mathbf{y})d\theta_2$, now say I have a sample of $M$ draws of $(\theta_1^{(i)}, \theta_2^{(i)})$ from $p(\theta_1, \theta_2|\boldsymbol{y})$, then to estimate $p(\theta_1|\mathbf{y})$, we use the following sample mean: $\widehat{p(\theta_1|\mathbf{y})} = \frac{1}{M} \sum_{i=1}^M p(\theta_1|\theta_2^{(i)}, \mathbf{y}) $, that is, we estimate the marginal density by averaging over the conditional densities. How can I implement that here? * * * EDIT2: After the help of @Tomas and @BabakP, I have tried to code the problem into R myself (I'm quite new to R, only learnt the language in the last couple of weeks). The following is my code: (Note: just for practice, I purposely went the other route of drawing from each conditional rather than @BabakP's more efficient decomposition)               ###########################################################     ########       Generation of data = y               #######     ###########################################################               true_beta = matrix(c(2,0.5,0.7),3,1)      ## True value of beta (vector) used in the data generating process (dgp)          true_sig = 0.3                            ## True value of sig used in the data generating process (dgp)               ## Setting the random number seed          set.seed(123456, kind = NULL, normal.kind = NULL)               ## Number of observations          nobs=20               ## Generating the values for x1, x2, x3 and y. Using Gaussian distribution here for convenience.          x1 = matrix(1,nobs,1)     x2 = matrix(c(rnorm(nobs,0,1)),nobs,1)     x3 = matrix(c(rnorm(nobs,0,1)),nobs,1)     u  = matrix(c(rnorm(nobs,0,1)),nobs,1)%*%true_sig          x = cbind(x1,x2,x3)          y = x%*%true_beta + u               ###########################################################     ### Specification of sample statistics     ###########################################################               invxx = solve(t(x)%*%x)     betahat = invxx%*%t(x)%*%y          df = nobs - 3     sighat = sqrt((t(y-x%*%betahat)%*%(y-x%*%betahat))/df)               ###############################################################     ###     ###  Generation of draws of beta/sigma/psi via Gibbs sampling     ###     ###############################################################               B = 100+1        ## Extra +1 since the Gibbs sampler coded below starts at j=2, hence burn-in period is from j=2 to j=101     M = 1000         ## Number of draws after burn in period.     repl = B+M       ## Total number of Gibbs iterations including burn in and number of draws after burn in period.               ## Dimensioning the matrix of beta/sigma/psi draws ##                betav = matrix(0,repl,3)           sigv = matrix(0,repl,1)     psiv = matrix(0,repl,1)               ## Generating underlying normal random draws to be used in Gibbs algorithm               stnbeta = matrix(c(rnorm((repl*3),0,1)),repl,3)                ## Specifying starting value for Gibbs chain          sigv[1] = true_sig               ## Gibbs sampling          for(j in 2:repl){            ## Generation of beta (vector) via its known multivariate normal conditional            meanbeta = betahat       varbeta = (sigv[j-1])^2*invxx            betagen = chol(varbeta)%*%stnbeta[j,] + meanbeta        betav[j,] = betagen                 ## Generation of sigma via its known inverted gamma conditional            a = t(y-x%*%betav[j,])%*%(y-x%*%betav[j,])            zvec = matrix(c(rnorm(nobs,0,1)),nobs,1)            chisq = t(zvec)%*%zvec            sigv[j] = sqrt(a/chisq)                 ## Draws of psi from its marginal posterior pdf            psiv[j] = (betav[j,2]+betav[j,3])/(sigv[j])^2                    }          #################################################################################     ## Estimation of the marginal posterior pdf of psi using kernel density smoothing     #################################################################################          ## Remove burn in period of psiv which is from j=1 to j=101 to create the vector psi which conists of 1000 draws after the burn in period           psi = psiv[102:repl]               ## Histogram of draws of psi          hist(psi,30,ylab="Percentage frequencies",xlab=expression(psi),          main=expression(paste("Histogram of draws of ",psi)))               ## kernel density plot of the draws of psi          d=density(psi)     plot(d,main=expression(paste("Kernel density estimate of marginal of ",psi)),          xlab=expression(psi),ylab=expression(paste("p(",psi,"|y)")))      My question is now how can I use the 1000 draws of psi from the above code to estimate the 5% quantile of $p(\psi|\mathbf{y})$?


title: Confusion related to linear regression

text: I learned a linear regression model with some input xs and output ys. My ys are always positive. However, now when I test the model on my test data, I get negative values for ys. I know I can either take the absolute value of ys or make them 0. But that's not the right way to do it. What is the correct approach to deal with this situation?


title: Contingency table with 2 different time series

text: I am trying to put the two time series of my Excel spreadsheet (see link below), into a contingency table. https://www.dropbox.com/s/2f96oylxj97fuih/example.xls The first series is the number of subscriptions per hour for 24 hours. The second series consists of the hours when ads for different radios stations have been displayed (e.g. ad for Channel 1 at 2:00PM). I have doubts about the feasibility of putting these 2 time series into a contingency table. Is this reasonable?


title: Random forest is overfitting

text: I am trying to use Random Forest Regression in scikits-learn. The problem is I am getting a really high test error: `train MSE, 4.64, test MSE: 252.25.` This is how my data looks: (blue:real data, green:predicted): ![Forest regression cleaned](http://i.stack.imgur.com/qmhbs.png) I am using 90% for training and 10% for test. This is the code I am using after trying several parameter combinations:               rf = rf = RandomForestRegressor(n_estimators=10, max_features=2, max_depth=1000, min_samples_leaf=1, min_samples_split=2, n_jobs=-1)      test_mse = mean_squared_error(y_test, rf.predict(X_test))     train_mse = mean_squared_error(y_train, rf.predict(X_train))          print("train MSE, %.4f, test MSE: %.4f" % (train_mse, test_mse))     plot(rf.predict(X))     plot(y)      What are possible strategies to improve my fitting? Is there something else I can do to extract the underlying model? It seems incredible to me that after so many repetitions of the same pattern the model behaves so badly with new data. Do I have any hope at all trying to fit this data?


title: Data Mining / ML applications in hydrodynamics?

text: I have a question about Scientific Data Mining.   1. Do you know successful case studies of applying Data Mining / Machine Learning techniques in hydrodynamics?   2. In general, does it make actually sense to try to apply DM/ML techniques to such deterministic systems as gas/fluid flows which are described by Navier-Stokes equations? I guess that answer is 'no' in the case of very simple flows. But maybe it makes sense if we have a deal with complex data from turbulent/multiphase/... flows?   1. What could be the problem formulation in this area? I will appreciate to your opinion, links to the papers and web-pages.


title: Process for Standardising and Normalising data

text: First: I'm not well versed in statistics terminology so please forgive me - I'll try to be as verbose as possible with my problem. This is a problem which I've previously solved very naively. I'm looking to apply more standard statistical theory to this problem in an attempt to get more accurate results. Any help with validating my approach or (rejecting it) and methods I could use would be appreciated (even the correct terms). The problem is a commerce related one: I have a collection of products (say 100) which have a variety of attributes. For this example I'll pick 3 of them:   1. Discount Percentage (Current Price/RRP) * 100   2. Views (Number of times viewed)   3. Stock Level Using these attributes I wish to be able to configure various sorting rules in order to place these products into an order. The sorting rules are where the 'configuration' happens - but the statistical part of my problem comes before that. In order to compare these attributes of products properly I need to: ## Standardize & Normalize. So the values and distribution of these attributes:   1. Discount Percentage: 0 to 100 with a distribution where most products have values of 0(50% of products), 20(25%), 30(10%), 50 (15%)   2. Views: 0 to 30000 with a fairly even distribution   3. Stock Level: 0 to 40 with a fairly even distribution (but different from the distribution for views) I take the entire population for each attribute, and using the max and min values I 'reset' each value to place in in a scale between 0 and 10 - where the maximum value is 0 and the minimum value is 10. For each member I work out if it is in the 0-10% percentile, 10-20 percentile etc. and give it a score of 1 to 10. However I feel this is a very naive solution as it doesn't take into account the distribution. For some attributes I get 99 members of the population with a score of 10 and 1 member with a score of 0. I'd like a more even distribution where outliers cannot change the end result so drastically (is this called smoothing?). I _could_ do this after the step above - but too much data is destroyed by this point (imo) to do a good job here. Is there a technique for doing both these steps at the same time (does this have a name?). Am I right in thinking that the normalization method for each attribute might end up being different? Should I just use trial and error on my data until I get to something that I deem satisfactory or can I consistently achieve a more equatable distribution without having to manually tweak this process? I'd like to have a range of scores between 0 and 1 (rather than 0 and 10). I with probably end up using R to help me with this as by all accounts it seems to be the best tool for the job. I'm not sure if it's overkill for what I see as being a fairly simple problem - but the complexity involved in the 'correct' Normalisation process may be beyond what I can easily do in my scripting language. Before I jump in with that though I'd like to validate my assumptions/theory. ## Weighting Once the above steps have taken place I allow domain knowledge experts to configure the weighting of these various attributes in order to determine the sort order (for example they may want to say that the number of views is 3 times more important than then stock level, in which case I'll multiple the normalised 'views' figure by 3, and then add it to the others before sorting the entire population). I've read that sometimes this knowledge could be incorporated into the standardisation process - so that the figures themselves take into account their importance relative to other attributes. I don't want to do this. In the example above I've talked about a population of 100. My real population is between 1000 and 50000, but I still think this is small enough to be able to work with the whole set of data rather than taking a sample. Any pointers would be appreciated as I've spent hours reading and I feel now like I'm going round in circles. -Thanks


title: Appropriate Application of Survival Analysis

text: I have an experiment that will produce observations of the time until an event occurs. Some basic properties are that   1. We count the number of events that have occurred at some point $t_1,...,t_n$.   2. Event times are interval censored, between $(t-1,t]$,   3. Individuals will not leave the trial between $t_1,...t_n$, that is an individual either experiences the event by $t_1$ or does not, at which point they are censored,   4. A large proportion of the individuals will not receive the event by $t_n$, when we terminate the experiment, and   5. I cannot assume any underlying parametric decay models. It appears to be a natural application for survival analysis. However, it is complicated by the fact that it is trivial to repeat the experiment from the same initial configuration multiple times. In effect, we will have a set of $m_{1,...,s}$ event counts (where $s$ is the number of samples) for each observation time $t_i$. I'm relatively new to statistics, and I'm struggling to see how to apply survival analysis to this data (if it is even applicable and there aren't more appropriate methods to measure this kind of time-to- event data). My inclination is to construct the survival function around the mean number of observed events at each interval (i.e. $\bar{m}_{1,...,n}$), which should better approximate the expected number of events at each interval in the population, however I have no idea if this is appropriate or of its implications. I've searched to no avail on Google Scholar, if anybody could point me to more material (or give me the correct nomenclature for what I'm trying to do), it would be appreciated. ### Edit Given that the $(t-1, t]$ intervals are uniform across all samples, say that I had the following matrix describing the cumulative number of individuals for which an event has occurred in each interval $M = \left(\begin{array}{ccc} 0 & 24 & 35 & 52 & 60 & 71 \\\ 0 & 22 & 38 & 57 & 64 & 75 \\\ 0 & 26 & 34 & 55 & 62 & 72 \\\ 0 & 21 & 32 & 52 & 61 & 73 \end{array}\right)$ where each row gives the event count for the same set of individuals at risk at $t=0$ across all samples (i.e. multiple instances of the experiment), and each column is an observation interval. I assume that by taking the mean number of events for each interval, I can get a better estimation of the expected population survival, so let $n$ denote the number of time intervals, $s$ denote the number of samples (experiment instances), then the vector $\bar{M} = \left[ {{\sum_{i=1}^{s}M_{it}}\over{s}} \right]_{t=1...n}$ will be the mean number of observed events for each time interval. My goal, then, is to use this as the input to the survival estimation. Let $f$ be the number of individuals at risk when $t=0$. Using the naive estimator (for now, given that the event intervals are uniform across all samples and there is no censorship until $t_n$), the survivor function could be estimated as: $S(t) = {{f - \bar{M}_t}\over{f}}$ Which would (hopefully) be a better estimate of the population survival than any individual sample (a single row from $M$). To reformulate my question:   1. Is $\bar{M}$ an appropriate input to a survival function estimation? I haven't seen this approach in any of the material I've read.   2. As I'm really, painfully a novice at statistics, can somebody point me to some material (academic papers, textbooks, wikis etc. would be fine) on estimating the confidence interval and variance for this survival function estimation? I presume it will not be identical to standard formulations. Apologies if my original question was confusing, I probably didn't include sufficient information.


title: fisher.test() and Chi-square bug in R

text: I'm getting this error when conducting fisher.test() in R:               freqTable =          1  2  3  4  5  6  7      1 14  0  2  0  0  0  0      2  0  0  0  9  0  0  0      3  0  6  0  0  0  0  0      4  0  0  0  0  6  0  0      5  0  0  0  0  0 10  0      6  0  0 10  0  0  0  0      7  2  4  9  1  0  0 30          > fisher.test(freqTable, workspace=2e+07,hybrid=TRUE)$p.value     Error in fisher.test(freqTable, workspace = 2e+07, hybrid = TRUE) :        FEXACT error 30.     Stack length exceeded in f3xact.     This problem should not occur.      Chi-square is also problematic:               chisq.test(freqTable)          Pearson's Chi-squared test          data:  freqTable     X-squared = 466.81, df = 36, p-value < 2.2e-16          Warning message:     In chisq.test(freqTable) : Chi-squared approximation may be incorrect      Looks like there's an open bug report: https://bugs.r-project.org/bugzilla/show_bug.cgi?id=1662 I'm wondering if anyone has found an appropriate workaround for these cases? Thanks


title: How to efficiently manage a statistical analysis project?

text: We often hear of project management and design patterns in computer science, but less frequently in statistical analysis. However, it seems that a decisive step toward designing an effective and durable statistical project is to keep things organized. I often advocate the use of R and a consistent organization of files in separate folders (raw data file, transformed data file, R scripts, figures, notes, etc.). The main reason for this approach is that it may be easier to run your analysis later (when you forgot how you happened to produce a given plot, for instance). What are the **best practices for statistical project management** , or the recommendations you would like to give from your own experience? Of course, this applies to any statistical software. ( _one answer per post, please_ )


title: "Ties" in Bayesian Survival Analysis

text: I have been reading about Bayesian methods for performing survival analysis. One of the things that strikes me is that in the three or so books that I have read, there has been no mention of handing ties, e.g. Efron, Breslow, or exact. Is this complication limited to frequentist statistical treatments?


title: How to detect when a regression model is over-fit?

text: When you are the one doing the work, being aware of what you are doing you develop a sense of when you have over-fit the model. For one thing, you can track the trend or deterioration in the Adjusted R Square of the model. You can also track a similar deterioration in the p values of the regression coefficients of the main variables. But, when you just read someone else study and you have no insight as to their own internal model development process how can you clearly detect if a model is over-fit or not.


title: What are some good resources to understand Hierarchical Statistical Models like Mixed Models?

text: The problem I am working on not only requires me to use several levels but also identify, from a set of attributes, which ones could be treated as levels. What is a good resource to understand Mixed Models? Also, what are the concepts that one must understand before putting such models to use? And what are other models that could be useful to in such hierarchical setting.


title: Deciding about dimensionality reduction with PCA

text: I have 2D data (I have a zero mean normalized data). I know the covariance matrix, eigenvalues and eigenvectors of it. I want to decide whether to reduce the dimension to 1 or not (I use principal component analysis, PCA). How can I decide? Is there any methodology for it? I am looking sth. like if you look at this ratio and if this ratio is high than it is logical to go on with dimensionality reduction. **PS 1:** Does PoV (Proportion of variation) stands for it? **PS 2:** Here is an answer: PCA and proportion of variance explained does it a criteria to test it? **PS 3:** Does Kaiser criterion is for that? **PS 4:** If you want to involve to this question, you don't read the comments below just read question it is ok.


title: Is Lorenz curve the same as QQ-plot?

text: Are Lorenz curves and QQ-plots the same? If not, where are the differences? I read about both of them and they appear to be two terms for the same type of plot / statistical technique to compare distributions. I was not able to find any confirmatory source for this. Perhaps you know?


title: Descriptive statistics in a repeated measures study

text: I'm doing a retrospective study looking at how the mortality of piglets differs with increasing age of pigs. Since my data is for each sow over seven years and each sow may only have mortality in one or two years, I'm running into a lot of trouble with descriptive stats. My Friedman's test showed a significance over the seven years but my standard deviations are all over the place! I've read a bit about Cousineau's confidence intervals in within-subject designs, but I can't find any reference to it being used in anything but psychology studies. Would it still apply here and, if not, does anyone have any suggestions on what to do? I can't figure out how to get error bars that don't make my data seem insignificant!


title: Distance Metrics For Binary Vectors

text: I have vectors of same length consisting of 1 and 0. I am trying to find out how similar they are. So far I am using hamming distance that I calculate sum of one vector then sum of second vector and the difference between this is the difference of the days. With 1 and 0 it works pretty well. My problem is that it doesn't reflect in any way where is the difference in the vectors and what is the variance of the error. I have thought of counting of how many 1 been misplaced to 1 of the next vector and how many 0 have been misplaced. It gives little bit more of information but still doesn't tell anything about the variance of the error. The vectors are used to represent occupancy of house in time, with every 1 indicating that house is occupied and 0 that it is not. From this I am trying to predict how next day will look.


title: Calculate mutual information of twitter's users interactions

text: I have a research project about finding out user's interaction on Twitter. If I want to calculate the mutual information between two users A,B: I(A,B), to know the probability of A and B appearing in the same discussion. So how can I calculate this? What can be the p(A), p(B) and p(A,B) in this case? As my thought, p(A) is (number of discussions have A)/(total number of discussions) p(A,B) is (number if discussions have A and B)/(total number of discussions) Is it correct? I'm not very good at statistical learning so I hope you can explain for me a little bit clearly and in details. Thanks.


title: One inflated beta regression with random effects using GAMLSS

text: I am new to modelling percentage data, and I would be greatfull for some advice. I have proportion data (0,1] on a percentage of money sent by Player B to Player A. Participants received an amount of money, and could decide what percentage they will send back. I have two categorical predictors (1st with 3 factors, 2nd with 2); one continuous predictor; and one nesting factor (class). Since the data are bound between 0 and 1, I figured out that the best option would be Beta regression. I tried to use `hglm` package which fitted well, however, since the data are one-inflated (many people chose to send back the full amount), I am looking for other options. As most appropriate seem to be `gamlss` package, which can use BEOI (Beta One Inflated) distribution. I used this code:               m1 <- gamlss(percent~cat1+cat2+continous, random(class), family=BEOI, data=dat,                   mixture="gq", K=1)      From what I understand from package help files, this should be the simplest option. However, it produces very different results from the `hglm`command. Especially the standard errors are higher than beta coefficients, leading to non-significant results. I tried to specify other other functions in the model (e.g., `K, sigma.formula, nu.formula, mixture` etc.), but these are beyond my understanding, and I am not really sure what I did there. I would very much appreciate any suggestions regarding either how to better specify the model, or simple explanations of `gamlss` function.


title: Failure rate for repairable units from usage and maintenance data

text: What are the approaches for characterizing the failure rate from usage and maintenance data of a fleet of repairable units? The usage data contains risk indicators (typically vibration levels, or operation limit exceedances), operational time durations, operating conditions, and recording times. The maintenance data contains maintenance events, start / end times. Plan to evaluate three models for repairable systems: the constant repair rate (HPP/exponential) model, the power law (Duane) model and the exponential law model (Ref Section 8.4.5. How do you fit system repair rate models?, http://www.itl.nist.gov/div898/handbook/toolaids/pff/apr.pdf


title: Five-point or seven-point Likert scale?

text: I am designing a questionnaire by using five-point Likert scale. What is the advantage of five-point Likert scale over seven-point Likert scale? Thank you.


title: CI for monthly production based on daily production data

text: Say I have 16 work days in a month, and I have 16 days of production data so that I have a mean and variance for this daily data (e.g. mean is 3750 and variance is 500). I can construct a confidence interval for the daily data, but how do I inflate the interval width to construct a confidence interval of the same level of confidence for the mean monthly production? That is, if I make a 95% confidence interval for the daily mean I also want a 95% confidence interval for the monthly mean. [You could created a very large CI around the one monthly production datum, but that is not what I am interested in.] Thanks!


title: Help with Anova of categorical and continuous variable in R and SPSS output

text: I am having some trouble running an Anova on categorical variables in R and matching SPSS output. What I need to do is run an anova on the dataset below (its a made up data set). But, I need to know if the mean of each category is significantly from the total mean of all races.               Satisfaction    Race     3   Asian     4   Cacasion     5   African American     2   Other      5   African American     3   African American     4   African American     5   African American     2   Asian     3   African American     1   Cacasion     1   Cacasion     1   Cacasion     5   Other      5   Other      5   Other      5   African American     5   Asian     4   Asian     5   Other      5   Other      5   Other      1   Cacasion     4   Cacasion      For example, the mean of all races is 3.5 :               > mean(test$Satisfaction)     [1] 3.5       What I would like to know is if the mean score for each race is significantly different from the total mean of 3.5 and the p-value. I ran an Anova in R with the following model, but R will set one catagory as the refernce and test is against the others :               > lm.test <- lm(test$Satisfaction ~ test$Race)     > summary(lm.test)          Call:     lm(formula = test$Satisfaction ~ test$Race)          Residuals:         Min      1Q  Median      3Q     Max      -2.5714 -1.0000  0.4286  0.8482  2.0000           Coefficients:                       Estimate Std. Error t value Pr(>|t|)         (Intercept)         3.8571     0.5023   7.679 2.18e-07 ***     test$RaceAsian     -0.6071     0.8330  -0.729   0.4745             test$RaceCacasion  -1.8571     0.7394  -2.512   0.0207 *       test$RaceOther      0.7143     0.7103   1.006   0.3266         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 1.329 on 20 degrees of freedom     Multiple R-squared: 0.391,  Adjusted R-squared: 0.2997      F-statistic:  4.28 on 3 and 20 DF,  p-value: 0.01732       The output is telling me that the mean for African American is 3.8571 and is significantly different from the mean of the caucasian group. It is not different from the mean of group Asian and Other. Is there a way for me to set the intercept to 3.5 in R and get significant compared to the mean and not the reference group. Or should I be using another tests altogether? My stats isn't that great so if its another tests a brief explain on which test and how to run it in R would be great.


title: How to resample in R without repeating permutations?

text: In R, if I set.seed(), and then use the sample function to randomize a list, can I guarantee I won't generate the same permutation? ie...               set.seed(25)     limit <- 3     myindex <- seq(0,limit)     for (x in seq(1,factorial(limit))) {         permutations <- sample(myindex)         print(permutations)     }      This produces               [1] 1 2 0 3     [1] 0 2 1 3     [1] 0 3 2 1     [1] 3 1 2 0     [1] 2 3 0 1     [1] 0 1 3 2      will all permutations printed be unique permutations? Or is there some chance, based on the way this is implemented, that I could get some repeats? I want to be able to do this without repeats, guaranteed. How would I do that? (I also want to avoid having to use a function like permn(), which has a very mechanistic method for generating all permutations---it doesn't look random.) Also, sidenote---it looks like this problem is O((n!)!), if I'm not mistaken.


title: Betweenness centrality applied to Amazon books graph

text: I made a visualization of Amazon related products. Every link in visualization means two products are often bought together. Now I'm applying various graph analysis techniques and am fascinated by the results. The biggest problem though is translating scientific terms to human language (please pardon me if it sounds snobberish). For example, I calculated indegree centrality of nodes, and I called it "What's popular here". After all the most connected node means it's most often bought with all other products in a given graph. Currently I'm working on _betweenness centrality_ , and it yields quite interesting results, often contradicting _degree centrality_. But how would you interpret these results? Most important products? For example, books graph of The art of R programming has the following top 3 nodes: **Indegree Centrality:**   1. The art of R programming \- 18 outgoing edges   2. R Cookbook (O'Reilly Cookbooks) \- 14 outgoing edges   3. Doing Bayesian Data Analysis: A Tutorial with R and BUGS \- 10 outgoing edges **Betweenness Centrality:**   1. The art of R programming \- centrality value of 1210   2. What is a p-value anyway? \- centrality value of 896   3. Visualize This \- centrality value of 784 The graph itself looks like this: ![enter image description here](http://i.stack.imgur.com/m3X1J.png)


title: R versus R Studio output differences

text: I ran the same code in R and in R-studio, but got two different results. Does anybody know why this is occurring, and if there is a fix for this? If there is no remedy, which is the "correct" program to use? Some information about the code I am running: I am using smooth.spline() as the only built in r-function (and the rest is my own code, which was copy/pasted from one program to the other). I am not generating random numbers (unless smooth.spline() does somewhere in the way it functions), and when I run the code multiple times in the same program I get the same results. Also, the results are the same on different computers in the same program. So much for reproducible results! **EDIT:** The only difference I can find is RStudio is 64-bit and R is 32-bit. Here is the code that I am running. More background, this is an assignment to code our own ACE algorithm for maximal correlation. I know there is already an r-function that exists (and I am sure it is more efficient than this).               my.ace <- function(x, y, tol=1e-6, maxiter=500) {            # inputs:       #       #   x       #    - numeric vector of length n       #   y       #    - numeric vector of length n       #   tol       #    - maximum tolerance for convergence (>0)       #   maxiter       #    - maxmimum number of iterations (>1)       #       # outputs:       #       #   fx       #    - optimal transformation of x       #   gy       #    - optimal transformation of y       #   maxcor       #    - correlation between fx and gy       #   iter       #    - number of iterations performed            # check inputs            if(!is.numeric(x) || !is.numeric(y) || length(x) != length(y)) {         stop('x and y must be numeric vectors of the same length')       }       if(tol <= 0 || maxiter < 2) {         stop('tol must be greater than 0 and maxiter must be greater than 1')       }            # setup variables            fx <- x-mean(x)       fx <- fx/sqrt(sum(fx^2))       gy <- y       corOld <- -100       corNew <- sum(fx*gy)       iter <- 0            # run ace algorithm            while(abs(corNew-corOld) >= tol && iter < maxiter) {         iter <- iter + 1         corOld <- corNew         gy <- smoother(fx,gy)         gy <- gy-mean(gy)         gy <- gy/sqrt(sum(gy^2))         fx <- smoother(gy,fx)         fx <- fx-mean(fx)         fx <- fx/sqrt(sum(fx^2))         corNew <- sum(fx*gy)       }            # return results            return(list(fx=fx,gy=gy,maxcor=corNew,iter=iter))     }


title: Removing seasonal and non-seasonal oscillations with least-square method

text: I have data on sea temperatures at different depths. With these data I need to remove seasonal and non-seasonal oscillations by fitting a function that consists of two sinusoids with periods of 12 and 6 hours on temperatue data (for every depth). The function needs to be put in seasonal diagram. And this adjustment needs to be done with least- square method.


title: computing confidence interval for ranges, confidence level 0.95

text: I have some ranges of values and frequencies of each range. Participants choosed an entire range when asked, and not a particular value. Example: Range [0,20] : frequency = 20 Range (20,60], frequency = 10 Range (60, 100], frequency = 15 I'm trying to figure out how to get a confidence interval, something like get an average interval of the values that i have, having a confidence level of 95%. Assuming that ranges represent for example how many minutes a person waited in line at the supermarket, I would like to say "Hey...given the data, you will probably wait 30-40 minutes in line if you go to the supermarket" Can you help or point me to some links where I can read how to achieve this? Thanks!


title: Best modelling packages in R?

text: In R I use a lot the packages plyr, stringr, Hmisc and ggplot2. Each of these packages take the base code and make functions that are more intuitive and easier to work with. Each of these packages tend to focus on the data munging and exploration side of data analysis. Are there any packages that make working with linear models easier and more intuitive?


title: Coefficients from a penalized Cox PH

text: I'm using the R package Penalized (0.9-42) on a Cox PH model. I'm using L2 (Ridge) on the grounds that I don't want to shrink my coefficients to 0. I don't understand why when I ask for:               fit <- penalized(Surv(ttocomplete2, event)~ x1 + x2 + as.factor(x3), lamda2=15.96))     print(coefficients(fit, "all"))      I get coefficients for the baseline/ref level of my categorical variables which, of course, I don't get when using `coxph{survival}` or `cph{rms}`.


title: Comparison of Kaplan-Meier curves across ordered groups

text: I am familiar with the log-rank test for comparing multiple Kaplan-Meier curves, but I am looking for a test that will compare across ordered groups (an ordinal variable). A significant result from the log rank test indicates that at least one of the groups has a survival curve different than the others. I want a test where a significant test indicates that there is a monotonic trend in the survival curves over the levels of the ordinal grouping variable. As a more concrete example, consider survival curves for cancer of different stages, and I want to show not just that the curves differ, but that the survival worsens as stage increases. As a bonus, pointers to implementations of the algorithm in R would be appreciated. * * * Edit to address `@dardisco` comment: When I tried to find a suitable test, I did find a reference so something that seemed appropriate in the documentation for MedCalc software at http://www.medcalc.org/manual/kaplan-meier.php. Two quotes that sounded like what I wanted from that are > Options: Linear trend for factor levels: Allows testing for a linear trend > across levels of the factor. It is appropriate if factor levels have a > natural ordering (for example, factor codes represent doses applied to > different groups). Kaplan-Meier assumes that the factor levels are equally > spaced. and > Logrank test for trend: If more than two survival curves are compared, and > there is a natural ordering of the groups, then MedCalc can also perform the > logrank test for trend. This tests the probability that there is a trend in > survival scores across the groups. But no specific algorithm for this test is stated.


title: Test if two samples of binomial distributions comply with the same p

text: Suppose, I have done:   * $n_1$ independent trials with an unknown success rate $p_1$ and observed $k_1$ successes.   * $n_2$ independent trials with an unknown success rate $p_2$ and observed $k_2$ successes. If, now $p_1 = p_2 =: p$ but still unknown, the probability $p(k_2)$ to observe $k_2$ for a given $k_1$ (or vice versa) is proportional to $\int_0^1 B(n_1,p,k_1) B(n_2, p, k_2) \text{d}p = \frac{1}{n_1+n_2+1}\binom{n_1}{k_1}\binom{n_2}{k_2}\binom{n_1+n_2}{k_1+k_2}^{-1}$, so if I want to test for $p_1 \neq p_2$, I only need to look in which quantile of the corresponding distribution my observations are. So far for reinventing the wheel. Now my problem is that I fail to find this in literature, and thus I wish to know: What is the technical term for this test or something similar?


title: How to combine the results of several binary tests?

text: First off let me say that I had one stats course in engineering school 38 years ago. So I'm flying blind here. I've got the results of what are essentially 18 separate diagnostic tests for a disease. Each test is binary -- yes/no, with no threshold that can be adjusted to "tune" the test. For each test I have what is ostensibly valid data on true/false positives/negatives when compared to the "gold standard", yielding specificity and sensitivity numbers (and anything else you can derive from that data). Of course, no single test has sufficient specificity/sensitivity to be used alone, and when you "eyeball" the results of all tests there's frequently no obvious trend. I'm wondering what is the best way to combine these numbers in a way that will yield a final score that is (hopefully) more reliable than any single test. I've so far come up with the technique of combining of the specificities of TRUE tests using               spec_combined = 1 - (1 - spec_1) * (1 - spec_2) * ... (1 - spec_N)      and combining sensitivities of the FALSE tests the same way. The ratio               (1 - sens_combined) / (1 - spec_combined)       then seems to yield a reasonably good "final score", with a value over 10 or so being a reliable TRUE and a value under 0.1 or so being a reliable FALSE. But this scheme lacks any true rigor, and for some combinations of test results it seems to produce an answer that is counter-intuitive. Is there a better way to combine the test results of multiple tests, given their specificities and sensitivities? (Some tests have a specificity of 85 and sensitivity of 15, other tests are just the opposite.) ## OK, my head hurts! Let's say I've got tests 1-4 with sensitivities/specificities (in %):   1. 65/50   2. 25/70   3. 30/60   4. 85/35 Tests 1 and 2 are positive, 3 and 4 negative. The putative probability that 1 is a false positive would be (1 - 0.5), and for 2 (1 - 0.7), so the probability that both are false positives would be 0.5 x 0.3 = 0.15. The putative probability that 3 and 4 are false negatives would be (1 - 0.3) and (1 - 0.85) or 0.7 x 0.15 = 0.105. (We'll ignore for the moment the fact that the numbers don't add up.) But the presumed probabilities that 1 and 2 are true positives are 0.65 and 0.25 = 0.1625, while the presumed probabilities that 3 and 4 are true negatives are 0.6 and 0.35 = 0.21. Now we can ask two questions:   1. Why **_don't_** the numbers add up (or even come close). (The sens/spec numbers I used are from "real life".)   2. How should I decide which hypothesis is (most likely) true (in this example it seems to be "negative" for both calcs, but I'm not sure that's always the case), and what can I use for a "figure of merit" to decide if the result is "significant"? ## More info This is an attempt to refine and extend an existing "weighting" scheme that is entirely "artistic" in nature (ie, just pulled out of someone's a**). The current scheme is basically on the lines of "If any two of the first three are positive, and if two of the next four, and either of the next two, then assume positive." (That's a somewhat simplified example, of course.) The available statistics don't support that weighting scheme -- even with a crude weighting algorithm based on the measured stats I come up with significantly different answers. But, absent a rigorous way of evaluating the stats I have no credibility. Also, the current scheme only decides positive/negative, and I need to create a (statistically valid) "ambiguous" case in the middle, so some figure of merit is needed. ## Latest I've implemented a more-or-less "pure" Bayesian inference algorithm, and, after going round and round on several side issues, it seems to be working pretty well. Rather than working from specificities and sensitivities I derive the formula inputs directly from the true positive/false positive numbers. Unfortunately, this means that I can't use some of the better quality data that isn't presented in a way that allows these numbers to be extracted, but the algorithm is much cleaner, allows modification of the inputs with much less hand calculation, and it seems pretty stable and the results match "intuition" fairly well. I've also come up with an "algorithm" (in the purely programming sense) to handle the interactions between interdependent observations. Basically, rather that looking for a sweeping formula, instead I keep for each observation a marginal probability multiplier that is modified as earlier observations are processed, based on a simple table -- "If observation A is true then modify observation B's marginal probability by a factor of 1.2", eg. Not elegant, by any means, but serviceable, and it seems to be reasonably stable across a range of inputs. (I'll award the bounty to what I deem to have been the most helpful post in a few hours, so if anyone wants to get a few licks in, have at it.)


title: When was the cross-correlation function invented? And what was the motivation behind it?

text: In a previous question I asked about the history of auto-correlation In this question I'm interested in the discovery process behind cross- correlation.   * When was it invented?    * Was it independently invented multiple times?   * What motivated its invention?


title: Effect size for rectangular frequency table

text: I have run several Chi-square tests and want a measure of effect size. I have some 1x2 tables, 1x4 tables, 4x9, and 3x4 tables of frequency data. What measure of effect size can I use on these rectangular frequency tables, if each cell has a different expected frequency? (I would like to use something like Cramer's V, but it only applies to tables with the same expected frequencies.) Sample data: two groups, 'false' observes 535 (expected: 3/4), 'true' observes 222 (expected: 1/4).


title: When was the autocorrelation function invented? And what was the motivation for it?

text: I'm just very curious about the discovery process behind the autocorrelation function.   * When was it invented?    * Was it independently invented multiple times, for example?   * What was the motivation for it?


title: Is there an R optimization package that can handle integer constraints and non-linear objective functions?

text: I am looking for an optimization routine that can optimize a non-linear objective function with integer constraints. NuOPT for S-Plus, CPLEX, or Matlab include powerful optimization packages for these kinds of optimizations? Is there any similar kind of package for R? Or is there an optimization procedure (perhaps genetic algorithms) that can solve this problem? For what it's worth I'm performing a portfolio mean-variance optimization. Thanks!


title: Video lectures about data mining?

text: I want to learn data mining. Are there any free video lectures out there which explain the process of data mining in depth?


title: How to fit a longitudinal model with binary outcomes

text: I'd like to fit a **longitudinal model** for where _multiple subjects experience binary outcomes over time_. To accomplish that, I'd like to use an additive random effect for each subject and an autoregressive error process to model the temporal stability in the "success probability." To illustrate, I'm simulating data from the following model: $\text{logit}(\mu_{it}) = \alpha_i + \beta x_i + \epsilon_{it}$ where $\alpha_i \sim \text{N}(0, \sigma_{\alpha}^2), \quad \epsilon_{it} \sim \rho \, \epsilon_{i, t-1} + \nu_t, \quad \nu_t \sim \text{N}(0, \sigma_{\nu}^2)$ and we observe: $Y_{it} \sim \text{Bernoulli}(\mu_{it}),$ for subjects $i = 1, \ldots, N$ and time periods $t = 1, \ldots, T$. It's kind of a lame longitudinal model because there aren't any fixed growth terms (e.g., a function of time), but I think it will suffice for the purposes of this discussion. The lme4 package in R does a really nice job with everything in this problem except for the temporal stability issue. So I went to **nlme** , but I don't think it ever work, at least the way I'm trying to use it. I'd like to know what went wrong, and any strategies to make this a reality. _Data Generating Process_               N <- 500     n <- 15     ar = 0.8     dat.list <- list()     ids.list <- list()     fixed.list <- list()     ran.effs <- .75 * rnorm(N)     x <- rnorm(N)     for (i in 1:N) {       lin.pred <- 0 + ran.effs[i] + .7 * x[i] + arima.sim(n = n, list(ar = ar, sd = .5))       dat.list[[i]] <- as.numeric(runif(n) <= plogis(lin.pred))       ids.list[[i]] <- rep(i, n)       fixed.list[[i]] <- rep(x[i], n)     }          longit.df <- data.frame(id = unlist(ids.list),                             x = unlist(fixed.list),                             y = unlist(dat.list))      _Model Fitting_ (tons of errors)               library(nlme)     my.nlme <- lme(y ~ x, data = longit.df,                    random = ~ 1 | id,                    weight = varPower(form = ~ fitted(.)*(1-fitted(.)), fixed = 1),                    cor = corAR1(0.6, form = ~ 1 | id))


title: Can I test hypothesis based on estimated parameters?

text: I estimated the mean and variance of two latent variables through two groups of data. I can't use the data to do hypothesis testing, because I am interested in the latent variable. Is there a way to test the whether the two latent variables are significantly different?


title: Classification of data with incomplete label sets

text: There is such a problem: we have to process multi-label classification (assignmet of tags) of text articles, using some pre-labeled training set. But for many texts in the training set, should be assigned more tags, that it was done by their authors. For example, there is a text about war medicine, and there is assigned only the text "war", and tag "medicine" is omitted. Are there any common, conventional methods to deal with such data?


title: Dimension Reduction to a Relative Single Value After PCA in R?

text: After reducing the dimensions of a matrix using PCA, I want to transform the "big" components (i.e. they cumulatively account for 90% of the variance) of each row into a single value. The value does not have to be unique, only relative so that they can be compared. Is there a statistical technique that will allow me to do this (especially in R)? Pros and cons would be helpful. I apologize if my wording is off--I am new to statistics.


title: How to find correlation among dependent variables?

text: If I want to find how strongly a dependent variable is related to another dependent variables in a study, do I make use of multiple regression? The reason I am asking is because the book mentions that for instance to conduct a standard regression, I will need to input all the independent variables and the corresponding dependent variable into SPSS. Not sure the procedure if I need to find a relationship between two or more dependent variables. Kindly clarify. Thanks


title: Using ANOSIM or PERMANOVA to abiotics data?

text: I'm new here and have some doubts. I'm working with a large data set of water quality and want to test for differences between the hydrological stations and between stretches of the river. There is a large internal variability into class factors. For this, I intend to do the following analyses: Exploratory PCA: to select variables that explain the variability in dataset. DA: to sort and visualize the a priori groups. Manova (transformed data) or ANOSIM: to detect statistical significance between a priori groups. In this analysis lies my main doubt, because I frequently use ANOSIM or PERMANOVA for species abundance data. Could I use these analyses for abiotic data? Another doubt: When I use PERMANOVA or ANOSIM, I do SIMPER analysis with a complement, to track what variables (species) that caused the difference detected. Could I use this analysis (SIMPER) too? Or should I test each variable separately by univariate analyses as Kruskal-Wallis or ANOVA?


title: How would you explain statistical significance to people with no statistical background?

text: Background:   I had to perform a data analysis for a client (some kind of lawyer) who was an absolute beginner in statistics. He asked me what the term "statistical significance" means and I really tried to explain it... but since I'm not good at explaining things I failed ;)


In simple linear regression, what is the covariance between the error term and the residual? Model: $y_i = \beta_0 +\beta_1 x_i + \varepsilon_i$ What will be the $\rm {cov}(\varepsilon_i,\ e_i)$, where $e_i = y_i - \hat y_i$?


title: Distribution of a logit transform of a normal variable

text: If I have a normal distribution with a given mean and variance and apply a logistic transform to it, what is the mean and variance of my transformed variable? This seems like it has to be a well known problem, but I haven't been able to find a quick reference anywhere.


title: Interpreting intercepts in mixed effect model with categorical predictors

text: Trying to fit a linear mixed effects model with 2 categorical predictors (group & worker) where worker is a random effect and group a fixed effect. I'm trying to figure out 1) whether I should specify intercept=0 and 2) why these 2 model results seem to give different conclusions about the effect of group. Model1: tps ~ group + (1 | worker) Model2: tps ~ group + (1 | worker) + 0               summary(Model1):     Linear mixed model fit by REML ['merModLmerTest']     Formula: tps ~ group + (1 | worker)         Data: mydata           REML criterion at convergence: 3489.872           Random effects:      Groups   Name        Variance Std.Dev.      worker   (Intercept) 1866     43.20         Residual             3165     56.26        Number of obs: 318, groups: worker, 18          Fixed effects:                 Estimate Std. Error     df t value Pr(>|t|)         (Intercept)    70.15      15.59  11.27   4.501 0.000848 ***     group phone   -20.85      21.75  10.83  -0.959 0.358586         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Correlation of Fixed Effects:                 (Intr)     group phone -0.717          summary(Model2):     Linear mixed model fit by REML ['merModLmerTest']     Formula: tps ~ group + (1 | worker) + 0         Data: mydata           REML criterion at convergence: 3489.872           Random effects:      Groups   Name        Variance Std.Dev.      worker   (Intercept) 1866     43.20         Residual             3165     56.26        Number of obs: 318, groups: worker, 18          Fixed effects:                    Estimate Std. Error    df t value Pr(>|t|)         group computer    70.15      15.59 11.27   4.501 0.000848 ***     group phone       49.30      15.17 10.40   3.251 0.008291 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Correlation of Fixed Effects:                 grpcmp     group phone 0.000       In the first model the 'phone' effect is the same as the difference between the two groups' effects in model 2 (this makes sense because in model1 the 'computer' group is the baseline). In model2, both groups' effects are significant, whereas in model1 only the intercept is significant. Which is the "right" model for a situation where the group predictor is binary? It must be only one or the other (seems to indicate that model1 is correct, because there the intercept "is the same as" the computer group, right? Model2 allows a "zero" value for group which doesn't make sense). Am I right about this? And how to interpret the fact that in model1 the intercept is significant but 'phone' is not?


title: How to test stability/variability of statistic over time, using R

text: I have a long time series data (running several years). I have split the data into several series, each of 1 month duration. I then perform some custom calculations on each of the 'shortened' time series and collate some robust summary stats (mode, MAD etc on my calculated data). I then want to see if there are any statistically significant differences between the stats calculated for the shortened series - so I can answer questions like:   1. Does the mode change over time   2. Does the MAD change over time ... etc? I am a little rusty, but I think this is what is referred to as Panel data? My questions then is which statistical test would be most appropriate to use to check if the summary stats change over time. I would be grateful if someone could post a little snippet/pseudocode in R to show how to do this.


title: Calculate standard error of difference between groups with different N

text: I have two datasets $X$ and $Y$ that have paired data, however one of the sets is larger than the other. I am interested in the standard error of the difference $Z = X - Y$. I guess I could use the approach described here, but I would also like to weight the standard deviation according to the sample sizes. So I came across this formula: $$ s^2_{pooled} = \frac{(n_1-1)s^2_1 + (n_2-1)s^2_2}{n_1+n_2-2} $$ which gives the common standard deviation of $Z$. Now the question: in the formula above, can I just substitute $s^2$ with the respective standard error of X and Y and take the square root in the end? Or do I have to calculate $s^2_{pooled}$ first, and the divide by $\sqrt{n}$? And if the latter, which $n$ do I take for that? Is it ok to derive an "average N", e.g. $(n_1 + n_2) / 2$? Thanks.


title: Additional parameter adds entropy to the centre of the density function

text: When an extra shape parameter is added in the distribution than a statement in favour of this parameter is written "extra parameter can control both tail weights and adding entropy to the centre of the density function". I am quite clear about tail weight but cannot understand adding entropy to the centre and how to prove this additional shape parameter add entropy to the centre. Also what advantages are obtained by adding entropy.


title: Interpretation of the constant in random-effects estimation with AR(1) process

text: I'm estimating a Random Effects model with an AR(1) disturbance for the errors (xtregar in Stata). The reported results show that the constant (in Stata terminology) is positive and significant. The question is, in general terms, how do you interpret this constant term?


title: Initiator follower analysis with time series data sets

text: I am a newbie to this forum. I searched different white papers and codes on google but couldn't find a solution, that's when I registered on this forum.. Please share in case you guys have a idea as to how to go about this problem.. I have a data set which consists of 30 prop. traders, and their historical trades on each symbol with trade dates..I want to understand the relationships that exist among each of these traders.. For eg: I want to understand who first trades on a stock and who follows next and so on.. basically derive the initiator-follower relationships.. And predict the same for the future... Any ideas as to how to go about? I started by taking a difference between trade dates for each ticker for any two traders and then taking a correlation of those dates/or difference in trade dates... I do not think this is a right approach..Any suggestions?


title: How to standardize data in Excel (mean = 0, SD = 1)

text: I have a large set of variables measured on different scales/units, and want to standardize them to the same scale, with a mean of 0 and standard deviation of 1, so that I can run a PCA on them. I have heard that the `NORM.S.DIST` function in Excel is appropriate - can anyone verify this?


title: Are large data sets inappropriate for hypothesis testing?

text: In a recent article of _Amstat News_ , the authors (Mark van der Laan and Sherri Rose) stated that "We know that for large enough sample sizes, every study—including ones in which the null hypothesis of no effect is true — will declare a statistically significant effect.". Well, I for one didn't know that. Is this true? Does it mean that hypothesis testing is worthless for large data sets?


title: Is it possible to have too much data? When should I subsample?

text: Social science PhD student here looking for a little help. I set out to collect around 300 online surveys for a study, but ended up with over 4,000. At that sample size should I sub-sample my results because it's more likely that significant results will be found because of the sample size? Thanks.


title: Question regarding significance level

text: Say, for instance, I'm estimating a model with about 5 million observations using linear regression or MLE. Given that the estimates are consistent, using the standard rule of rejecting the null on a 5% significance level will pretty much never make me entitled to reject the null (of the parameter being equal to zero), i.e. everything is significant. I should probably use another significance level, which is smaller than 5%, and my question is if there is any literature on picking the right significance level, or how would you guys handle this problem? Should I simply stick with the 5% significance level? Thanks.