title: CDF and confidence interval on non-parametric binomial data

text: I am having trouble with a series of system tests in which each test is assumed to be a normal approximations of binomial distribution. Each test sample, n, is a Bernoulli ("success", "failure") trial of the system. Each system test has approximately 550 IID observations. Each has a lifetime that may vary between 2-8 seconds, regardless of success or failure. During this lifetime a success or failure will occur. Therefore, for each system test an estimated probability of failure, $p$, is determined where $p = \rm{Pr}[X=0]$. However, since the test data in which the estimated probability of failure was determined is available, I wish to develop a CDF, F(t), that provides the probability of failure per unit time (0.1 seconds). More specifically, I wish to develop a graph to show when the failure occurred. I have a few questions that I am having problems with:   1. I would like to develop the CDF's, but I am not sure whether the lifetime differences of each sample varying from 2-8 seconds may bias the CDF curve. This curve will provide the cumulative estimated probability of failure in 0.1 second increments.    2. Do the number of failures within each time interval affect the significance of the estimated probability of failure within that particular time interval?   3. In addition I would like to develop two corresponding curves that illustrate the estimated probability of failure using a 95% confidence level. How do I do this, and what assumptions of the distribution should I consider? I realize that this may be a trivial problem, but I have doubt about some assumptions. Since my work has recently highly relied on these CDF responses to make important decisions, I would like to gain more confidence on how to complete these distributions and provide the resources/references that support these statistical conclusions.


title: Importance of normal distribution

text: Why did the normal distribution become such a popular (important) distribution? I know one reason is because of CLT. Can you please give more reasons?


title: Linear regression with group-dependent intercept and variance

text: This is my first question on CrossValidated and I'm not a professional statistician (although I am trained in theoretical probability) so please be indulgent. I have data of the form $(X_i, Y_i, K_i)_{1 \leq i \leq n}$ where $X_i,Y_i$ are 1-dimensional continuous variables and $K_i$ is categorical (ordinal in fact), representing the group in which the $i$-th observation belongs. I want to fit a linear regression model in order to explain the dependent variable $Y_i$ by a multiple of $\alpha X_i$, plus an intercept that depends on the group, so $$ Y_i = \alpha X_i + \beta(K_i) + \varepsilon_i$$ I have been able to estimate $\alpha$ and the $\beta(k)$s by introducing dummy regressors $X^k_i = 1_{K_i=k}$, and running a simple OLS estimator on the formula `Y ~ X + Z1 + ... + Z8 - 1`. But I have reasons to think that the residual variance is different from one group to another : $\mathrm{Var}(\varepsilon_i|X) = \sigma^2(K_i)$ and obviously unknown. I think this will have a material effect on the estimates but I don't know how to implement this (say in R). My googling has led me to the package "RFGLS" but I'm not sure how to apply it to my situation since I'm not familiar with the whole "fixed effects" / "random effects" terminology. I would very much appreciate if anyone could provide me with tips as to how to estimate this model, or with a formulation of my problem in the canonical statistical language. Thanks. PS : If it helps, the function $k \mapsto \beta(k)$ is expected to be decreasing, and this is true of the OLS estimates $\hat{\beta}_k$. I expect it to be also true for the variance function $k \mapsto \sigma^2(k)$. PS2 : Actually my ultimate hope is to prove that $\alpha$ is "reasonably" close to zero, for instance by not being able to reject $H_0 : \alpha=0$ at the $5\%$ level (I'm aware using hypothesis testing to "prove" the null is not exactly good practice). Edit : Sorry I didn't make this clear in the first place : my working hypothesis is that the $\varepsilon_i$s are centered and mutually independent given $X$.


title: Connecting the dots in a graph

text: Given a graph which is the product of an equation, we can meaningfully calculate any points that solve the equation, and consequently, also a line that goes through the points. The line, at any of its points, is the answer. But what about measures like "number of cars per hour / parking lot." Does it make sense to connect the dots of each hourly measure? Wouldn't this lead us to false results, or at least misleading results, since we do not know how many cars were there in the meantime?


title: not concave iterations on a zinb model

text: my name is vincenzo and i have that type of problem with zinb that you intorduce in this discussion (ZIP converges but ZINB does not. Should I drop this model?): the iterations continues to be not concave. If I take out one variable from my model, the second iterationis ot concave, the third is "backed up", then the model run normally. What is going on? can you help me please? thank you in advance! Regards, Vincenzo


title: Phylogenetic GLS with categorical response variable

text: I want to perform phylogenetic comparative analysis on a multilevel (3 level) categorical response variable, preferentially in R. Does anyone know how to do that or if there is a way to do at least? Thanks very much in advance! Best wishes, Jasmine


title: Chi-square test problem (and some more maybe)

text: A colleague of mine has asked me for help on this problem, but I am confused with it as well, the main question he asked me was that the chi-square result he was getting was very high. He conducted a survey, and among other questions, one of them is a question which has 8 sub-questions about certain processes in company of examinees and the examinees have to give their grade of the quality of those processes by choosing values from 1 (bad) to 5 (great) (classical Likert scale). There is one other value (sixth) that can be selected which is causing me trouble aside from the main problem. That sixth option is titled "N/A" and means that the examinees don't have that process in their company. Enough for intro (hopefully). So the data table he made in Excel (for one subquestion) looks like this (there were 195 examinees in total): ![Image 1.](http://i.stack.imgur.com/Y6HA4.jpg) I wonder why did he assume that theoretical frequencies (f$_t$) are fixed to 32,5 (shown as rounded to 35 in data table) for all values? My opinions is that he is wrong, and that the theoretical frequencies should be these (and if not these, for sure not the same for every mark): ![Image 2.](http://i.imgur.com/hJdBA.jpg) The chi-square value is even greater now. I am not even sure why does he calculate it in this fashion in Excel. My opinion that for this data set chi-square value should be calculated with CHITEST function of the following syntax (cell addresses are as in Image 2.):               =CHITEST(C9:C14;D9:D14)      That way i get a chi-square value of 1,26177E-44, which I think is ok. Please help me, all I thought I know or knew on topic seems to be vanished now. Sorry if any terms are not understandable, English is not my primary language. If any clarification is required, I'll do my best. Thank you in advance for your time and knowledge and sorry that images are not automatically displayed, I need to gain more reputation in order to enable them.


title: A/B testing with multiple conversion goals

text: Suppose I run A/B testing with two variants on my website and track multiple conversion goals (e.g. goal a: post a comment, goal b: click facebook like button). To make it simple for now, a conversion goal is a binominial variable that can take a value of either 0 or 1. The system I use for capturing the results lets me only see the totals, e.g.:               _Variant A_       Total samples: 10000       Total conversions goal a: 400       Total conversions goal b: 600                 _Variant B_       Total samples: 10000       Total conversions goal a: 750       Total conversions goal b: 500      In the case of having only one conversion goal, I could use a z-test or chi- square to determine if both variants differ significantly. The number of samples and the number conversions is all the input I need. In the case of having multiple conversion goals, what method can I use to determine if two variants are significantly different?


title: What's the name of this discrete distribution (recursive difference equation) I derived?

text: I came across this distribution in a computer game and wanted to learn more about its behaviour. It comes from the decision as to whether a certain event should occur after a given number of player actions. The details beyond this aren't relevant. It seems applicable to other situations, and I found it interesting because it's easy to calculate and creates a long tail. Every step $n$, the game generates a uniform random number $0 \leq X < 1$. If $X < p(n)$, then the event is triggered. After the event once occurs, the game resets $n = 0 $ and runs through the sequence again. I'm only interested in one occurrence of the event for this problem, because that represents the distribution that the game is using. (Also, any questions regarding multiple occurrences can be answered with a single occurrence model.) The main "abnormality" here is that the probability parameter in this distribution increases over time, or put another way, the threshold rises over time. In the example it changes linearly but I suppose other rules could apply. After $n$ steps, or actions by the user, $$ p(n) = kn $$ for some constant $0 < k < 1$. At a certain point $n_{\max} $, we get $p(n_{\max}) \geq 1 $. The event is simply guaranteed to occur at that step. I was able to determine that $$ f(n) = p(n)\left[1 - F(n - 1)\right] $$ and $$ F(n) = p(n) + F(n-1)\left[1 - p(n)\right] $$ for PMF $f(n)$ and CDF $F(n)$. In brief, the probability that the event will on the $n$th step is equal to the probability $p(n)$, less the probability that it has already happened on any preceding step. Here's a plot from our friend Monte Carlo, for fun, with $k \approx 0.003$. The median works out to 21 and average to 22. ![enter image description here](http://i.stack.imgur.com/R1s2A.png) This is broadly equivalent to a first-order difference equation from digital signal processing, which is my background, and so I found that quite novel. I'm also intrigued by the notion that $p(n)$ could vary according to any arbitrary formula. My questions:   1. What's the name of this distribution, if it has one?   2. Is there any way to derive an expression for $f(n)$ without reference to $F(n)$?   3. Are there other examples of discrete recursive distributions like this? **Edits** Clarified process about random number generation.


title: Finding suitable rules for new data using arules

text: I am using R (and the arules package) to mining transactions for association rules. What I wish to do is construct the rules and then apply them to new data. For example, say I have many rules, one of which is the canonical `{Beer=YES} -> {Diapers=YES}`. Then I have new transactional data where one of the records has purchased beer but not diapers. How can I identify a rule where the LHS is met, but not yet the RHS? R example:               install.packages("arules")     library(arules)          data("Groceries")     **#generate Rules omitting second record**          rules <- apriori(Groceries[-2],parameter = list(supp = 0.05, conf = 0.2,target = "rules"))      Rules generated are:               > inspect(rules)       lhs                   rhs                   support confidence     lift     1 {}                 => {whole milk}       0.25554200  0.2555420 1.000000     2 {yogurt}           => {whole milk}       0.05603010  0.4018964 1.572722     3 {whole milk}       => {yogurt}           0.05603010  0.2192598 1.572722     4 {rolls/buns}       => {whole milk}       0.05664023  0.3079049 1.204909     5 {whole milk}       => {rolls/buns}       0.05664023  0.2216474 1.204909     6 {other vegetables} => {whole milk}       0.07484238  0.3867578 1.513480     7 {whole milk}       => {other vegetables} 0.07484238  0.2928770 1.513480      The second transaction shows this customer, since they have yogurt but not whole milk perhaps should be sent a coupon for milk. How can any applicable rules in "rules" be located for new transactions?               > LIST(Groceries[2])     [[1]]     [1] "tropical fruit" "yogurt"         "coffee"


title: Total area under any probability density function

text: What's the name of the theorem that tells us that the total area under any probability density function, discrete or continuous, equals 1? My stats book actually _defines_ a PDF by requiring that $$\sum_{x}f(x)=1\quad\text{or}\quad\int_{-\infty}^{\infty}f(x)=1$$ In other words, $f(x)$ is a PDF **only if** the above is true (along with a few other requirements). Rather than using this as a definition, is it possible to _prove_ that it's true for any PDF? If so, has it been done before? I'd like the name of the theorem/proof so I can reference it.


title: How can I generate data from a SVM?

text: I would like to generate samples from a support vector machine, with a Gaussian kernel and a fixed `C` and `sigma`? For examples `(x1,x2)` and corresponding class label `y` Effectively I would like to generate a dataset such that I know a `C` and `sigma` that work well on it. Then knowing these values I can access various method of hyper-parameter optimisation methods? I've seen how it can be done with Gaussian processes but they are probabilistic so I'm not really sure how or if it can be done with SVMs? Or is there a different way to construct a dataset to assess hyper-parameter optimisation routines?


title: Identifying feature values that influence an outcome

text: I have a data set which has data about 1 million people. Data about each person consists of a 'Score' and about 100 features (each of which refers to some characteristic of the person - example - age, location, ethnicity, level of education, access to internet, etc. etc. ) I want to find the clusters/pockets within this population that have a very low 'Score. I also want to get the set of feature values that are representative of these clusters. An example output would be - (1) People with Location = X,Y,Z and Level of Education < A have low score AND (2) People with Level of Education = B and Location = C have low score AND (3).. so on. What statistical analysis/ML method is best suited for this?


title: What's the difference between a component and a factor in parallel analysis?

text: The `psych` package in R has a `fa.parallel` function to help determine the number of factors or components. From the documentation: > One way to determine the number of factors or components in a data matrix or > a correlation matrix is to examine the “scree" plot of the successive > eigenvalues. Sharp breaks in the plot suggest the appropriate number of > components or factors to extract. “Parallel" analyis is an alternative > technique that compares the scree of factors of the observed data with that > of a random data matrix of the same size as the original. fa.parallel.poly > does this for tetrachoric or polychoric analyses. When I run the function I get the following output:               Parallel analysis suggests that the number of factors =  7                                   and the number of components =  4      What is the difference between a _factor_ and a _component_?


title: Effective variance for two independent variables

text: When you're fitting a model $y=f(x)$ to data (${x_i, y_i}$) with errorbars on both the independent ($x$) and response ($y$) variables, it's standard that you can define an 'effective variance' when calculating the chi-squared $$(\delta y_i)^2 + \left(\delta x_i \left(\frac{df}{dx}\right)_{x=x_i} \right)^2$$ Here the differential is evaluated at the value of $x_i$ of interest. My question is if this can be extended to the case where there are two independent variables $(x, z)$, so that you're fitting a function $y=g(x, z)$. Something like: $$(\delta y_i)^2 + \left(\delta x_i \left(\frac{dg}{dx}\right)_{x=x_i} \right)^2 + \left(\delta z_i \left(\frac{dg}{dz}\right)_{z=z_i} \right)^2$$


title: How robust is estimation using sum of powers?

text: Minimizing sum of squares for scale and location estimation is optimum wrt efficiency for normal distributions, and doing so but using the sum of absolute values is maximum likelihood estimation scheme for the Laplace distribution, i.e. for having stronger exponential tails. The Huber estimator is switching smoothly from quadratic to abs deviations and is optimum for something in between, but it comes with the need for a scale estimate as additional parameter to deal with :-(. To avoid that but still having a potential even higher robustness (like for distributions with power law tails) we can also look to the sum of sqrt(|x-xo|). This way outliers have even less impact, but I wonder how efficient is this compared to MLE or M-estimates? Or is there a distribution where this sum of sqrt is already an ML estimator? And we can extend the sum of sqrts method, i.e. avoiding still the need for a scale parameter plus improving/tailoring efficiency? I found no real information on this in the literature :-(. Here are some of my results: I made a little comparision of sum of squares versus sum of abs value plus also sqrt and x^4 and x^0.25 on different distributions. For normal data of course x² is best, but of course for a uniform x^4 is much better, and e.g. sqrt is better than x² already on a Student-4, for a Cauchy x^0.25 is best amoung this little benchmark, but still well below ML efficency!! So I wonder if the scheme of sum of powers can be generalized and improved further? Best regards Stephan


title: Too many controls?

text: I am purposely keeping this question vague since this has happened to me with multiple datasets. I have a DV and a few IV then several controls, and all looks fine. Then I add one more controls and the output of my OLS looks the same BUT I see the following               F(X, Y)  = .     Prob > F = .      This is happening in Stata but, I am hoping it is a theory/application problem and not a program problem. What I do not understand is that the control I am adding applies to 25% of the sample, so it is not a fine grained control. What does this mean? Can I use the point estimates? Is there any documentation about this happening?


title: Reducing the number of variables by clustering or pca

text: Assume you have much more variables than samples and you want to reduce the number of variables (because you think some of them are redundant). The natural approach to do this is PCA, but is it possible to apply clustering algorithms where you consider the features as samples and try to regroup them? If so: is there a 'better' approach between pca and clustering of the variables?


title: Empirical Bayesian Approach in Normal Normal Model with Unknown Variance

text: For $i=1, \ldots, K$ and $j=1, \ldots,n$, assume the following model. \begin{align}\label{model} X_{ij} \mid \mu_i & \sim N(\mu_i, \sigma^2) \nonumber \\ \mu_i \mid \sigma^2 & \sim N(\mu_0, \sigma^2/\kappa_0) \nonumber \\ \sigma^2 & \sim Inv-\chi^2(\nu_0, \sigma^2_0) \end{align} We assume everything is independent. Then the posterior is as below \begin{align*} \pi(\mu_i \mid \mathbf{X_i}) \sim t_{\nu_n}(\mu_i \mid \mu_n, \sigma^2_n/\kappa_n) \end{align*} where \begin{align*} \mu_n &= \dfrac{\kappa_0}{\kappa_0+n}\mu_0 + \dfrac{n}{\kappa_0+n} \bar X_i\\ \kappa_n & = \kappa_0 + n\\ \nu_n & = \nu_0 + n\\ \nu_n \sigma^2_n & = \nu_0\sigma^2_0+(n-1)s^2+\dfrac{\kappa_0 n}{\kappa_0+n}(\bar X_i - \mu_0)^2 \end{align*} I would like to make inference using empirical Bayes approach for which I need to estimate $\mu_n, \sigma^2_n, \kappa_n$ and $\nu_n$. Is there any reference on this or could anyone give some suggestions? Thank you!


title: Which test is useful to assess if a source is the same in two samples of binary data?

text: I am looking for a statistical test that assesses if two samples which have binary values come from the same source or from different ones. Example: I have one set that has 10 0s and 6 1s and another set which has 50 0s and 29 1s. I would like to test the two sets from belonging from the same source or are significantly different. What is the best one? Is t-test appropriate for this specific case?


title: pymc random seed doesn't guarantee the same posterior samples?

text: If pymc.numpy.random.seed(0) guarantee the same random number sequence to initialize a stochastic variable (say a Uniform distribution), why does its posterior samples (from trace plot) don't have the same values for multiple runs with the same seed=0 ? Is there any internal random seed encoded in the pymc module ? Or, this is cause by the assigned "probability α" ? ( See http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm#Intuition ) Have anyone encountered this problem ? Is this related to the older version of pymc that I'm using ? Currently, I'm using pymc version 2.2 because failed to install version 3 to my computer with windows 7 platform.


title: Nonlinear model wirh effects over linear term

text: I have been working in `R` with nonlinear models such us: $Y = \alpha_{0}\text{(varia)} + \alpha_{1}\text{Time}\text{(varia)} + \alpha_{2}\sin(\frac{2\pi\text{Time}}{\alpha_{3}}) + \alpha_{4}(-1^{\text{Time}})$ And I would line whether a Bernoulli variable affects to the $\alpha_{1}\text{Time}$ and intercept or not. Such Bernoulli variable could be:               varia<-rep("BEFORE","AFTER"),each=30)      Using:               Y<- -2.5+rnorm(60)+2*sin(2*pi*time/8)+2.5*(-1)^time     time<-seq(1,60)      I found that `nls` can fit this model but without the effect of this categorical variable, I also found that the `nlme` package can estimate effects over variables not over terms of the nonlinear model. My question is: Which of these package could help me? and, how could I add this Bernoulli variable in the code? Regards and thanks.


title: Are slopes of calibration experiment different?

text: I'm running a calibration experiment for an instrument that measures soil moisture content. I have 6 buckets filled with different soils that I weigh daily and at the same time take a measurement with the soil moisture probe. I now have 6 different XY linear models all with different slopes/intercepts/R2 etc. and want to test to see if the slopes for the different buckets are different. To be clear, and I think the reason other answers aren't appropriate, is that n = 1, as each bucket represents a different soil type (sand, clay, silt, mixes of those for the other threee); however, the calibration curves have ~18 points each (18 different moisture contents and their corresponding instrument measurement). I've been using R, so my data are in 12 columns AX, AY, BX, BY, ... FX, FY and would appreciate the code to run the test as well. (A = sample A (sand), B = clay etc. with X and Y being moisture content and probe reading, respectively, taken at that same time (e.g., AX is moisture content for sand)). Moisture contents are different between A to F (i.e., not the same value at each weighing time). Thoughts?


title: Multiple hypothesis tests, "comparing the expected values"

text: If I have two sets of random samples and I'd like to test whether there is some difference or not between them (or if one is better than the other) I'd test the null hypothesis $H_0: \mu_1 = \mu_2$ with a one-tailed test. So if it turns out that the null hypothesis is rejected on some significance level I could say the one has a significantly higher mean than the other (It was some time ago I studied this so please correct me if I'm wrong). But what if I have more than two sets of random samples, is it possible to do a similar (I would assume: set of) test(s) so that I could draw the conclusion that one or some of has a significantly higher mean than the others. If possible you someone please direct me to some keywords, papers and/or briefly explain how it's done. Thanks in advance!


title: Latent Dirichlet Allocation

text: Nowadays I am trying to understand LDA. I have read a lot but something is missing. For instance, when I have a vector with three different variables (buy order, sell order, volume) respectively: [14,34,23] And have over 12 Million records per each customer. I want to profile them by using LDA. How should I organize my data? What will be the analysing process? Furthermore, how will I interpret the results? I could not find any numerical (not topic modeling) examples. Thanks in advance.. Teitia


title: Normal approximation to binomial

text: What do I do when the normal approximation is not valid? Here's the question I'm trying to answer: A student guesses all 15 answers on a multiple-choice test. There are 5 choices for each of the questions. Find the probability (correct to 6 decimal places) that the student passes the test. thanks


title: Maximum Likelihood

text: Let $x$ have an exponential density $p(x|\theta) = \theta e^{-\theta x} \text{ if }x \ge 0;\quad 0 \text{ otherwise.}$ a) Plot $p(x|\theta)$ _versus_ $x$ for $\theta = 1$. Plot $p(x|\theta)$ _versus_ $\theta$, $(0\le \theta \le 5)$ for $x = 2.$ How would I do this?


title: Random forest on grouped data

text: I am using random forest on high-dimensional grouped data (50 numeric input variables) which have a hierachical structure, since the data have been collected with 6 replications at 30 positions of 70 different objects resulting in 12600 data points which are not independent. It seems random forest is over-fitting the data, since the oob error is much smaller than the error which we get when leaving data from one object out during training and then predicting the outcome of the leaved out object with the trained random forest. Moreover we have correlated residuals. I think the overfitting is caused since random forest is expecting independent data. Is it possible to tell the random forest about the hierarchical structure of the data? Or is there another powerful ensemble or shrinkage methode which can handle high-dim grouped data with strong intereaction structure? Any hint how I can do better?


title: Time series digitization

text: I am looking for a viable method to create an alphabet out of multidimensional time series data. My intention is to use pattern/string matching algorithms (potentially genetic) to find profitable patterns. Extra points if I can do this in R. Any ideas?


title: How to generate user-friendly summaries of cluster analysis?

text: I used BIRCH and HAC to do clustering on my data. I want to now what type of information I can include in reports that my users can generate to get more insights on the clusters. I would have to dumb down the statistical terms and represent them as much as possible visually in these reports for my users.


title: Spatial econometrics -- computing residuals

text: Consider the model in spatial econometrics denoted SAR by James P LeSage: $y = \rho W y + X \beta+ \epsilon$ I use the `R` package `spdep` and the econometric toolbox by LeSage from www.spatial-econometrics.com. However, when comparing the way these two compute the residuals there are differences. The `R` package by Roger Bivand computes them by $\dot{y}=y - \rho W y$ and then linearly regressing $\dot{y}$ on the $X$ matrix. Obtaining residuals from this fit. (`lm(y-rWy~X-1`). Yet LeSage computes the residuals the way (I originally thought would be the right way): $y-(\rho W y)^{-1}Xb$ So here're my questions:   1. Why do the two authors compute the residuals differently?   2. Is there an objection against one or the other way of computing them? Thank you for kind help!


title: How can you implement a two-way ANOVA with nesting in R or SPSS?

text: I currently am conducting a study where I have three variables: one that is binary, and two numerical ratio variables. Each of the subjects in my study has values for each of the three variables. Variables:   * Condition (binary): Values 0 and 1   * Pre (ratio)   * Post (ratio) I want to test if there is a significant difference between the the pre and post variables of the `0` control group and the `1` experimental group. Both groups have 103 subjects. The data meet all typical ANOVA assumptions such as normality and the like. I was thinking of nesting the variables as follows and then running a two-way ANOVA. Variable 1 is `exposed` / `not` and variable 2 is `pre` / `post`. People are nested in Variable 1 (meaning that each person gives both pre and post information for either exposed or not exposed conditions) Would this be the correct way to approach this problem? Also how would I implement this statistical analysis, preferably in R or SPSS?


title: Multivariate outlier analysis of data with multimodal distribution

text: I have a data on methylated gene expression which follows a bimodal distribution. Can I use a multivariate outlier analysis to detect differentially methylated genes? What are the various methods of doing this and what assumptions come with them? I know similar analyses have been used in analyzing gene expression data, but they claim the data are assumed to follow a multivariate normal distribution.


title: Meaning of “post hoc” multiple comparisons

text: My original thought was “post hoc” in a post hoc multiple comparison procedure means that the procedure is taken after an omnibus test (such as ANOVA F test) gives a significant result. But after reading Howell's _Statistical Methods for Psychology_ , “post hoc” in a post hoc multiple comparison procedure means the procedure is planned after data collection, and “prior” in a prior comparison means the procedure is planned before data collection. So was my original thought wrong? Is Tukey–Kramer method (Tukey's HSD) post hoc? How about Scheffe method? * * * A table from Howell's book (here is link to Google Books, although not having the table): ![enter image description here](http://i.stack.imgur.com/uyS4J.png)


title: Computing confidence intervals for prevalence for several types of infection

text: I have a dataset about a population in hospital, and what type of infections patients have. Let say the number of patients is 100, 10 of them have pneumonia `(group A)`, 20 of them have urinary tract infection `(group B)`; keep in mind that group A and B can be overlapping, that is a patient suffering from pneumonia can also have urinary tract infection. I need to estimate the prevalence of different infection type in this population (i.e., prevalence of pneumonia, prevalence of urinary tract infection). I'm not sure whether assuming binomial distribution, like the one below from here, is appropriate: $$\text{SE}=\sqrt{\frac{p\times (1-p)}{n}}\times\sqrt{1-f}$$ Using this formula, I will compute multiple "binomial" estimates (i.e., one for each type of infection). I would feel comfortable to use this if only I need to describe prevalence of one type of infection, but in this case, I need to describe several ones from the same population. I am not sure if using the formula is appropriate or not in this case. Can anyone here enlighten me? Thanks!


title: How to do ANOVA on data which is still not normal after transformations?

text: I'm looking at the effect defeat and entrapment inducing conditions have on subjective ratings of defeat and entrapment at three different time points (among other things). However the subjective ratings are not normally distributed. I've done several transformations and the squareroot transformation seems to work best. However there are still some aspects of the data that have not normalized. This non- normality manifests itself in negative skewness in High entrapment high defeat conditions at the time point I expected there to be the highest defeat and entrapment ratings. Consequently I think it could be argued that this skew is due to the experimental manipulation. Would it be acceptable to run ANOVAs on this data despite the lack of normality, given the manipulations? Or would non-parametric tests be more appropriate? If so is there a non parametric equivalent of a 4x3 mixed ANOVA?


title: Model for continuous response and a mix of continuous and categorical predictors

text: Which statistical model is appropriate when the response is continuous and the predictors are a mix of continuous and categorical? What is the disadvantage in using GLM combined with gaussian family? Here is my dataset and model in `R`:               df <- structure(list(as.factor.pred. = structure(c(1L, 1L, 5L, 3L,      2L, 8L, 3L, 5L, 2L, 2L, 3L, 2L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,      1L, 2L, 4L, 1L, 1L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 1L, 7L, 3L, 3L,      6L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 5L, 2L, 2L, 1L, 1L, 2L, 1L, 1L,      1L, 1L, 1L, 1L, 2L), .Label = c("A", "B", "C", "D", "E", "F",      "G", "H"), class = "factor"), res = c(33, 33, 37, 32, 32, 26,      33, 28, 25, 34, 29, 35, 26, 20, 27, 19, 30, 33, 27, 24, 26, 28,      27, 23, 26, 25, 24, 26, 24, 25, 21, 21, 23, 24, 23, 27, 23, 20,      21, 22, 22, 22, 22, 23, 23, 21, 22, 21, 21, 23, 23, 18, 20, 18,      18, 18, 19)), .Names = c("as.factor.pred.", "res"), row.names = c(NA,      -57L), class = "data.frame")          names(df)[1] <- "pred" ## fix up the names to match formula          model <- glm(res ~ pred, data = df)          summary(model)          Call:     glm(formula = res ~ pred, data = df)          Deviance Residuals:         Min      1Q  Median      3Q     Max       -6.625  -3.000  -0.625   2.000  10.000            Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept)  2.300e+01  9.080e-01  25.331   <2e-16 ***     predB        2.625e+00  1.471e+00   1.784   0.0806 .       predC        4.714e+00  1.971e+00   2.391   0.0207 *       predD        3.500e+00  3.397e+00   1.030   0.3080         predE        6.333e+00  2.823e+00   2.243   0.0294 *       predF       -8.283e-15  4.718e+00   0.000   1.0000         predG        1.000e+00  4.718e+00   0.212   0.8330         predH        3.000e+00  4.718e+00   0.636   0.5278         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for gaussian family taken to be 21.43562)              Null deviance: 1277.6  on 56  degrees of freedom     Residual deviance: 1050.3  on 49  degrees of freedom     AIC: 345.85


title: Automated Normality Test on large samples

text: I have a automated measurement system which measures the volume of particles. In our case it is known that the distribution of the volume particles follow the lognormal distribution. But under some circumstances this assumption can be violated (errors during the image recording etc). These cases should be recorded. My idea is to do an automated normality test on the log-transformed sample values. When the p-value indicates a rejection of the null hypothesis, the user of the system should be informed. The problem is, that the sample size is usually very big (several thousand values) and since the tests are known to be very sensitive for large samples ( i use the Shapiro-Francia test), i came up with the following idea. Instead of using the whole sample, i create x sub samples by randomly picking n values for each sub sample. The test is than applied on the sub samples and mean p-value for all sub samples x is used to test the null hypothesis. Is this from the statistical point of view a possible solution for my problem? Greetings Kevin


title: A good book with equal stress on theory and math

text: I have had enough courses on statistics during my school years and at the university. I have a fair understanding of the concepts, such as, CI, p-values, interpreting statistical significance, multiple testing, correlation, simple linear regression (with least squares) (general linear models), and all tests of hypothesis. I had been introduced to it much of the earlier days mostly mathematically. And lately, with the help of the book Intuitive Biostatistics I have grasped and unprecedented understanding towards the actual conceptual theory, I believe. Now, what I find I lack is understanding of fitting models (estimating parameters to model) and the like. In particular, concepts such as maximum likelihood estimation, **generalized** linear models, bayesian approaches to inferential statistics always seem foreign to me. There aren't enough examples or tutorials or conceptually sound ones, as one would find on simple probabilistic models or on other (basic) topics on the internet. I am a bioinformatician and I work on RNA-Seq data which deals with raw read counts towards finding, lets say, gene expression (or differential gene expression). From my background, even if I am not familiar with statistical models, I am able to grasp the reason for a poisson distribution assumption and negative binomials and so on.. But some papers deal with generalized linear models and estimate a MLE etc.. which I believe I have the necessary background to understand. I guess what I am asking for is an approach some experts among you deem useful and (a) book(s) which helps me grasp these concepts in a more intuitive way (not just rigorous math, but theory backed with math). As I am mostly going to apply them, I would be satisfied (at the moment) with understanding what is what and later, I can go back to rigorous mathematical proofs... Does anyone have any recommendations? I don't mind buying more than 1 book if the topics I asked for are indeed scattered to be covered in a book. Thank you very much!


title: Generalized RBF Kernels

text: There is the notion of Generalized RBF Kernels, for example in "Towards Optimal Bag-of-Features for Object Categorization and Semantic Video Retrieval" from Jiang (1) or in formula (2.72) in http://agbs.kyb.tuebingen.mpg.de/lwk/sections/section23.pdf (2). They define: Gen. RBF Kernel: $e^{-pd(x,y)}$, where (1) says d can be any distance function and (2) requires it to be a metric. I now have been wondering for a long time: Is a function defined by $e^{-pd(x,y)}$ where d is a metric always a Mercer's Kernel? I think in (2) it sounds like it, but if it were true, not many people know of it, since they are still looking to prove that $e^{-EMD(x,y)}$ is a kernel, where EMD = Earth Mover's Distance, which is a metric. If its not a Mercer's Kernel - do you know of any counter example?


title: Correlation based on one object

text: **Background:** I have two random variables - height and weight. We often estimate correlation between them based on random sample of individuals. However, what would happen if I took only one individual and measured his weight and height in many points in time. **Questions:** What is the limitation of such approach? Is correlation based on one individual still correlation between two random variables? Thank you.


title: How to fit distribution having expected value and 0.95 quantile?

text: I would like to estimate parameters of lognormal distribution having only two values: expected value and 0.95 quantile. What do you recommend? Thanks!


title: Testing endogeneity against decision probability

text: I am working with a model where the literature suggests there is potential endogeneity between the dependent variable and the primary independent variable of interest (a binary treatment). However, there is a logical argument that selection may be made into the treatment across the spectrum of the continuous dependent variable. I have a dataset (N=513) where I am trying to determine if selection to treatment is essentially random with respect to the dependent variable. Is there a way to test the probability of selection to treatment to demonstrate that there is no endogeneity? I have already done the regular tests for endogeneity and the results have been inconclusive. I want a "positive" result that will support the lack of endogeneity, or that will demonstrate endogeneity does in fact exist.


title: P-values and null hypothesis

text: I am making a study guide for a statistics test at the end of the summer semester. My professor went through some stuff that will be on the test. There is one type of question that I just don't understand. It gives you a p-value and asks you to retain or reject the null hypothesis, for example:               For p = 0.1500, what is the decision if α = 0.05?      A. Retain the null hypothesis.      B. Reject the null hypothesis.           For p = 0.1500, what is the decision if α = 0.01?      A. Retain the null hypothesis.      B. Reject the null hypothesis.       I don't know how to tell what to do with the null hypothesis. Could someone explain this example to me because there are about 7 or 8 of these on the study guide?


title: Estimating a p-value when you can't compute it for the whole set

text: I have a two large lists and I want to compute the p-value that they are similar. The similarity algorithm is a black-box, but for this application we will trust it gives accurate p-values. My problem is my list is too large and the algorithm won't give answers for them. However, I think that if I take a smaller random sample of each list the algorithm will give a p-value. Is it permissible to randomly sample both lists a large number of times and average the p-values? Should I randomly sample with or without replacement?


title: Is there a good reason for the name "Missing at Random"?

text: Non-statisticians can be easily misled by the phrase "missing at random". Their natural interpretation tends to be closer to MCAR than to MAR. I had to explain the concept to an ecologist today, and his reaction was something along the lines of "But why on Earth would you call it that?". I must admit, I couldn't come up with a better answer than "that's just the name we've always used". Does anyone here have a better response?


title: What are some interesting and well-written applied statistics papers?

text: What are some good papers describing _applications_ of statistics that would be fun and informative to read? Just to be clear, I'm not really looking for papers describing new statistical methods (e.g., a paper on least angle regression), but rather papers describing how to solve real-world problems. For example, one paper that would fit what I'm looking is the climate paper from the second Cross-Validated Journal Club. I'm kind of looking for more statistics-ish papers, rather than machine learning papers, but I guess it's kind of a fuzzy distinction (I'd classify the Netflix Prize papers as a bit borderline, and a paper on sentiment analysis as something I'm _not_ looking for). I'm asking because most of the applications of statistics I've seen are either the little snippets you seen in textbooks, or things related to my own work, so I'd like to branch out a bit.


title: Exploratory Analysis - finding the most important factor

text: I have a dataset of 113 variables. In exploratory analysis the first thing I want to know is what are the most important factors on a single variable (revenue). I learned that naive Bayes would work, but I do not know how. I am using `e1071` package and coded according to the example:               library(e1071)     data(iris)     classifier <- naiveBayes(iris[,1:4], iris[,5])      table(predict(classifier, iris[,-5]), iris[,5])      If you know any other algorithm that is not computationally expensive and suitable for exploratory analysis please let me know.


title: Remeasuring "bad" values

text: We have a large number of samples whose concentration we measure twice, averaging the two values. Typically, the coefficient of variation (cv) for each sample is < 5%, but for a few samples the cv is high. We assume that in these cases something went wrong with one or both concentration measurements. We can afford one more concentration measurement for the samples with high cv's. My question is, how to use the three measurements to achieve the "best" estimate of the true concentration? Average all three measurements? Pick the two with the lowest cv? Or...? Many thanks for any insights or pointers to literature.


title: 3 variables and conditional probability

text: Following a previous question, lets say we now have 3 variables: $L$, $B$, $S$:                    S        /   \       L     B      So $L$ depends on $S$ $B$ depends on $S$ $P(S) =$ 0.5 $P( \lnot S) =$ 0.5 $L$ that depends on $S$: $P(L|S) =$ 0.10   $P( \lnot L|S) =$ 0.90   $P(L| \lnot S) =$ 0.01   $P( \lnot L| \lnot S) =$ 0.99 $B$ that depends on $S$: $P(B|S) =$ 0.60   $P( \lnot B|S) =$ 0.40   $P(B| \lnot S) =$ 0.30   $P( \lnot B| \lnot S) =$ 0.70 I have gotten $P(L) =$ 0.055 $P(\lnot L) =$ 0.945 $P(B) =$ 0.45 $P(\lnot B) =$ 0.55 As the previos question to get $P(S)$ after I observe that $P(L)=1$   $P(S∣L)= P(L∣S)P(S) / P(L) =$ (0.10)(0.5) / (0.055) = 0.9091   so $P(S)=$ 0.9091 Similarly when we have $P(B)=$1   $P(S∣B)= P(B∣S)P(S) / P(B) =$ (0.40)(0.5) / (0.45) = 0.3636   so $P(S)=$ 0.3636 But What do you do when you observe both events: $P(L)=$1 and   $P(B)=$1 How do you modify the above formula to get $P(S)$?


title: How to compare mean distance traveled by two groups?

text: I have two independent groups who travel with different modes. I want to compare the distance traveled (which is a continuous variable) by different modes for these two groups. How do I compare them in SPSS? My sample size is small: each group has only 32 observations. * * * Sorry, I am a newbie in stats. Here is the actual question. I collected travel distance and the choice of travel mode in two different locations. There are 32 people in each location with 6 different modes. Each person can choose one mode and give the travel distance. My query is, how should I compare to see if the travel distances of these two groups are significantly different? If it were categorical data, I could use chi square, but does the t-test work here? I believe the t-test is to compare two groups with one mean in each group, but I have 6 $\times$ 2 separate mean values (one for the each mode and two locations). I hope I have made it clear now. Thanks


title: Difference between sampling distribution and data distribution

text: I went to a stats course refresher last week and the instructor talked about data distribution and sampling distribution. I am just practicing with the exercises shown in class. Based on my dataset below, For the data in Week2, what can I explain about the ‘data distribution’. and the ‘sampling distribution’ of the sample mean? I would appreciate any explanation of these two terms. Thanks               fishy <- structure(list(week1 = c(2.011, 1.994, 10.332, 7.056, 4.926,                                       6.12, 2.039, 5.948, 7.731, 8.41,                                        11.055, 7.157, 5.855, 25.243,                                       36.553, 76.281, 38.902, 64.689, 42.934,                                        80.373, 115.858, 145.981,84.735,                                        163.084, 190.472, 295.79, 254.71,                                        273.446, 582.495),                            week2 = c(535.013, 513.534, 824.442, 1130.764,                                       1396.367, 1122.016, 1263.061, 1449.587,                                         1588.527, 680.99, 1861.677, 1432.656,                                       2921.025, 3595.931, 2071.98, 1666.726,                                       1594.989, 1522.255, 2496.464, 2169.722,                                       1870.255, 1039.203, 54.847, 0.266, 60.603,                                      601.822, 244.916, 124.749, 74.059)),                            .Names = c("week1", "week2"), class = "data.frame",                            row.names = c(NA, -29L))     boxplot(fishy,week1~week2)


title: Evaluation of a classification model in Rapidminer

text: I have some trouble with the evaluation of a classification model in Rapidminer. I train a naive bayes and svm model by using a split validation (first training of naive bayes, then application of the model and evaluation by the performance operator). The result is one performance window with the results of the testing data. But, as far as I know the software should return two performance windows. One for the results of the training and one for the testingdata. How can I see if my model is overfitting if I just see the results of the testingdata? Thank you


title: Classification with 3 groups, repeated measurements, missing values, more predictors than subjects

text: I am working on a classification problem with the following characteristics:   * Individuals belong to one of three groups. The groups are "somewhat ordinal": controls, subclinical and clinical group.   * Each individual is measured five times on each predictor. These are technical replicates taken at the same time. All predictors are numerical.   * Many predictor measurements are missing. For some individuals and some predictors, all five measurements are missing.   * I have more predictors than subjects. What I already did or am thinking about:   * I already tried an approach by Yeung & Bumgarner (2003), which addressed all four of my issues (e.g., it downweights predictors with high average intra-individual variation). Unfortunately, there does not seem to be any code for it, so I coded it myself in R. When I optimized the shrinkage parameter $D$ via cross-validation, the optimal value was very high, so that in effect, the algorithm discards all predictors and assigns everything to the most common group. This may just be what is in my data, or I may have made an error in coding, but I would like to try at least one alternative approach before giving up.   * Partial least squares and similar approaches that essentially rely on matrix algebra will - as far as I understand - have problems with missing data. And I guess they would only deal with my replicates by averaging the observations from each ID/predictor combination, which sounds wasteful to me.    * Given that my groups are somewhat ordinal, I am tempted to use a polytomous Rasch model next, where the link functions could be analogous to mixed models (to account for the replications), with some sort of regularization/shrinkage/lasso/ridge regression (to account for $p>n$). I haven't found an implementation of this anywhere, so it looks like I would need to "roll my own" here. Any thoughts or pointers? Anything that is already implemented in R (or somewhere else) would be especially welcome. I paste some example training and test data in R below.               training.data <- structure(list(ID = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L,        3L, 4L, 4L, 4L, 5L, 5L, 5L), .Label = c("a", "b", "c", "d", "e"       ), class = "factor"), group = structure(c(1L, 1L, 1L, 1L, 1L,        1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L), .Label = c("A", "B",        "C"), class = "factor"), predictor.1 = c(0.537, 0.42, 0.15, 0.245,        0.738, 0.513, 0.612, 0.633, 0.685, 0.29, 0.358, 0.624, NA, NA,        NA), predictor.2 = c(NA, NA, NA, NA, 0.846, 0.489, 0.954, 0.627,        0.709, 0.296, 0.975, 0.858, 0.904, 0.104, 0.242), predictor.3 = c(NA,        0.636, 0.007, NA, 0.396, NA, NA, 0.5, 0.776, NA, 0.222, NA, 0.398,        NA, 0.197), predictor.4 = c(0.364, 0.406, 0.761, 0.195, 0.695,        0.205, 0.215, 0.336, 0.311, 0.907, 0.91, 0.267, 0.815, 0.703,        0.486), predictor.5 = c(0.601, 0.568, 0.038, 0.246, 0.118, 0.992,        0.61, 0.486, 0.893, 0.004, 0.123, 0.373, 0.436, 0.595, 0.289),        predictor.6 = c(0.261, 0.75, 0.603, 0.921, 0.398, 0.362,        0.384, 0.603, 0.563, 0.591, 0.387, 0.918, 0.904, 0.965, 0.076       ), predictor.7 = c(0.251, 0.724, 0.169, 0.658, 0.702, 0.687,        0.474, 0.769, 0.081, 0.19, 0.798, 0.717, 0.514, 0.672, 0.911       ), predictor.8 = c(0.502, 0.868, 0.261, 0.424, 0.884, 0.569,        0.213, 0.781, 0.384, 0.771, 0.657, 0.501, 0.833, 0.357, 0.924       ), predictor.9 = c(0.379, 0.375, 0.01, 0.316, 0.747, 0.585,        0.316, 0.054, 0.463, 0.731, 0.239, 0.377, 0.87, 0.895, 0.373       ), predictor.10 = c(0.737, 0.397, 0.491, 0.906, 0.111, 0.61,        0.902, 0.77, 0.321, 0.137, 0.49, 0.561, 0.131, 0.202, 0.1       )), .Names = c("ID", "group", "predictor.1", "predictor.2",        "predictor.3", "predictor.4", "predictor.5", "predictor.6", "predictor.7",        "predictor.8", "predictor.9", "predictor.10"), row.names = c(NA,        -15L), class = "data.frame")          test.data <- structure(list(ID = structure(c(1L, 1L, 1L), .Label = "f",       class = "factor"),        group = c(NA, NA, NA), predictor.1 = c(0.203, 0.568, 0.458       ), predictor.2 = c(0.729, 0.531, 0.156), predictor.3 = c(0.584,        0.995, 0.079), predictor.4 = c(0.307, 0.085, 0.152), predictor.5 = c(0.966,        NA, 0.108), predictor.6 = c(0.637, 0.97, 0.699), predictor.7 = c(0.597,        0.672, 0.806), predictor.8 = c(NA_real_, NA_real_, NA_real_       ), predictor.9 = c(0.943, 0.106, 0.71), predictor.10 = c(0.804,        0.198, 0.408)), .Names = c("ID", "group", "predictor.1",        "predictor.2", "predictor.3", "predictor.4", "predictor.5", "predictor.6",        "predictor.7", "predictor.8", "predictor.9", "predictor.10"), row.names = c(NA,        -3L), class = "data.frame")


title: Multiplying two entropy values

text: I have seen formulas that either sum or subtract entropy values (eg Information Gain). However, I have not seen a "multiplication/division of entropy values" anywhere. I would like to know if there is any formula/theorem that multiplies/divides entropy values. Additionally, I would like to know if there is any mathematical reasoning/meaning behind such operation. Thanks.


title: Intraclass correlation coefficient, non parametric data

text: Does the ICC require that my data are normally distributed? If more information about my data is needed, just ask.


title: How to generate predicted survivor curves from frailty models (using R coxph)?

text: I want to compute predicted survivor function for a Cox proportional hazards model with frailty terms [using survival package]. It appears that when frailty terms are in the model, the predicted survivor function cannot be computed.               ## Example      require(survival)     data(rats)          ## Create fake weight     set.seed(90989)     rats$weight<-runif(nrow(rats),0.2,0.9)          ## Cox model with gamma frailty on litter     fit <- coxph(Surv(time, status) ~ rx+weight+frailty(litter,dist="gamma"),     data = rats)           ## Compute survival curve from the cox model for rx=0 and weight=0.5 kg     plot(survfit(fit, newdata=data.frame(rx=0,weight=0.5)),xlab = "time",     ylab="Survival")           ## Running this line, I get following error message:     Error in survfit.coxph(fit, newdata = data.frame(rx = 0, weight = 0.5)) :      Newdata cannot be used when a model has sparse frailty terms      I tried using both sparse and non-sparse computation methods by using `sparse=TRUE, Sparse =FALSE, sparse =0, sparse=5` options. However, none worked. How do I compute predicted survivor curve based on my frailty model?


title: AR(1) selection using sample ACF-PACF

text: The following graph shows the ACF (sample autocorrelation function) and PACF (partial autocorrelation function) of the residuals in a linear regression. There is a sinusoidal decay in the ACF and two spikes at lag 1 and 4 in the PACF. Considering the spike at lag 4, can we still assume AR(1) for these residuals? I should add that, the best model based on AIC is ARIMA(2,0,1) and based on BIC is AR(1) using the `auto.arima` function in package forecast. Here are the codes:               auto.arima(resid(lm.3),ic=c("aic"))     Series: resid(lm.3)      ARIMA(2,0,1) with zero mean               Coefficients:          ar1     ar2     ma1       0.0579  0.5125  0.8500     s.e.  0.2228  0.1981  0.1753          sigma^2 estimated as 0.0008776:  log likelihood=207.36     AIC=-406.71   AICc=-406.29   BIC=-396.33               auto.arima(resid(lm.3),ic=c("bic"))     Series: resid(lm.3)      ARIMA(1,0,0) with zero mean               Coefficients:          ar1       0.7771     s.e.  0.0616          sigma^2 estimated as 0.0009202:  log likelihood=205.11     AIC=-406.23   AICc=-406.1   BIC=-401.04      ![enter image description here](http://i.stack.imgur.com/HEQvS.png) The time series plot of the original residuals is ![enter image description here](http://i.stack.imgur.com/gNNwc.png) And the ACF and PACF plot for the first difference is ![enter image description here](http://i.stack.imgur.com/5f8zY.png)


Suppose all the simple correlations between $x_i$ and $x_j$ are $r$ for all $i,j=1,2,\dots,p, i\neq j.p>8$. Find the partial correlation coefficient $r_{1p.2468}.$ By definition, $$r_{1p.2468}=\frac{cov(e_{1.2468},e_{p.2468})}{\sqrt{var(e_{1.2468})}\sqrt{var(e_{p.2468})}}$$ how can I find $cov(e_{1.2468},e_{p.2468}),var(e_{1.2468}),var(e_{p.2468})$ in terms of $r$?


title: Is a sequence of random variables indexed by a homogeneous Poisson process process strictly stationary?

text: I'm revising for an exam and have no idea how to approach this question: Let $\{N_t\}_{t\geq 0}$ be a homogeneous Poisson process of parameter $\lambda > 0$. Let $\{X_k\}_{k\geq 0}$ be a sequence of random variables, identically independently distributed, with $E[X_k] = 0$ and $Var[X_k] = \sigma^2 < +\infty$. Is the process defined by $Y_t = X_{N_t}$ strictly stationary? I'm guessing the fact that the Poisson process has independent increments and we have: $P(N_{t_2} - N_{t_1} = k) = \dfrac{\lambda (t_2 - t_1))^k}{k!} e^{-\lambda (t_2 - t_1)}$ is important, but not sure how it relates, any hints appreciated!


title: Reliable factors

text: Briefly, my study ($n=188$) consists of using PCA on a 111-item questionnaire claiming to have 21 different disorders (factors). I did a MAP test (Velicer's Minimum Average Partial test) to decide on the number of factors to extract (it said 16), followed by PCA with varimax rotation. Using .4 as minimum loading. The tenth and eleventh components of the rotated component matrix have only 3 items (.698, .616, .542 and .572, .568, .483). Are these components reliable?


title: What are the different types of averages?

text: 1. I know there are three types of averaging methods in statistics: Mean, Mode and Median. Are there any other type of averaging methods that statisticians use?    2. How do I know which method is best for a particular set of data?


title: How to do an instrumental variables regression with an instrumented interaction term in Stata?

text: I'm having a bit of a problem with Stata syntax. I need to do the following regression: $$y = ax + bz + c(xz) + e$$ where both $x$ and $z$ are instrumented and also the interaction term $xz$ uses the instrumented values of $x$ and $z$. Just generating the predicted values for $x$ and $z$ and using them as regressors yields incorrect standard errors. Edit: I also need to do a similar regression with only one of the variables instrumented and with this one instrumented variable being in the interaction term.


title: Univariate feature ranking in classification

text: Scikit-learn has function to evaluate the F-statistics for univariate feature importance feature selection. According to the web page they are calculating ANOVA F value. If I understood correctly, the univariate feature importance checks how important a specific feature is to the target without taking into consideration the other features. If I have the features $f_1,f_2,f_3,...f_n$, the univariate feature importance will check how well $f_1$ fits into the target, $f_2$ fits with the target,... $f_n$ fits into the target? Why do we use ANOVA F when there are only two variables that are checked against each other? I think I need some basic background here to understand the use of ANOVA F when doing univariate feature selection.


title: State Space formulation of Hodrick-Prescott ﬁlter

text: I would like to apply the Kalman filter in order to get a causal Hodrick- Prescott ﬁlter. The Hodrick-Prescott ﬁlter models a time series $(y_t)_{t=0}^T$ as $$ y_t = \tau_t + c_t $$ where $\tau_t$ is a trend component and $c_t$ is a cyclical component. This reference defines a state space formulation of the form $$ y_t = \tau_t + c_t $$ as the measurement equation and $$ \tau_t = 2 \tau_{t−1} − \tau_{t−2} + \epsilon_t $$ for the unobservable trend. I have three questions on this: A) $c_t$ is assumed to be a random error here, right? Normally distributed with constant variance. This seems a difficult assumption to me. B) What's the logic behind the equation for the trend? C) Does anybody know a different state space formulation for this problem? Or a nother reference? Thanks!


title: Var(X) is known, how to calculate Var(1/X)?

text: If I have only $\mathrm{Var}(X)$, how can I calculate $\mathrm{Var}(\frac{1}{X})$? I do not have any information about the distribution of $X$, so I cannot use transformation, or any other methods which use the probability distribution of $X$.


title: Expected value from a regression table

text: My question is related to the discussion that take place in this post. The question is: how can I compute the expected values of conditions GroupB/Cond2 and GroupB/Cond3 in that case? I know that GroupA/Cond1, 2 and 3 are calculated from coef A (6.1372) + coefs from conditions 1, 2 and 3) and I know that GroupB/Cond1 is calculated from coef B (6.0758). The problem is that the coefficients from the R output to GroupB/Cond2 is -0.1055. So, the number in the table should be 6.1913 (6.0758 + 0.1055), but it is 6.0853. Similarly with GroupB/Cond3: it should be 6.1758 (6.0758 + 0.1000), but it is 6.1149. I don't understand that difference. Is there an error in the table or is it that I really don't understand how to calculate the values? In other words, my question is about how to calculate the coefficients for the interactions terms in the lmer function output from R.               ---------------------------------     Group          A - 6.1372     B - (6.1372 - 0.0614) = 6.0758     ---------------------------------          Condition          1 - (6.1372)     2 - (6.1372 + 0.1150) = 6.2522     3 - (6.1372 + 0.1000) = 6.2372          ---------------------------------     Interactions          A1 (6.1372)          = 6.1372     A2 (6.1372 + 0.1150) = 6.2522     A3 (6.1372 + 0.1000) = 6.2372          ---------------------------------          B1 (6.0758) = 6.0758     B2 (6.0758 + ??????) = 6.0853     B3 (6.0758 + ??????) = 6.1149          ---------------------------------


title: Comparing noisy data sequences to estimate the likelihood of them being produced by different instances of an identical Markov process

text: _(Prompted to some extent by the answers already given by Shane and Srikant, I've rewritten this to try to clarify what I'm getting at, if only to myself.)_ Suppose we have several similar systems, each with behaviour that approximates a continuous time Markov process. That is, there are some number of discrete states the system can be in and associated probabilities of transitioning from one state to another at any instant, depending solely on the current state. For now, consider the processes to be _stationary_ , ie the transition probabilities do not change over time, and unaffected by seasonality or other external considerations. Unfortunately, we cannot measure the state of any system directly, but have instead to measure a proxy quantity, which varies with state but is _not_ discrete and is subject to various sources of noise and error. The principal question is this: **Q1: Given two data sequences produced independently from two such systems, how can we decide whether the underlying Markov processes are the same?** Now, it _may_ be that the best way to approach this is as two separate problems:   1. Convert the imperfect proxy sequence into an idealised time series of (categorical) states   2. Determine whether the state sequences correspond On the other hand, such a separation _might_ involve discarding some information from the data in step 1 (eg, about its variability) that would be useful in step 2. Which leads to: **Q2: Does it make sense to decompose the problem in this way or is it better to compare the proxy data directly?** If such a decomposition does make sense, that opens up a whole other issue about how to do the idealisation, but that's definitely a question for another day. Shane, below, mentions goodness-of-fit and distributional tests such as Anderson-Darling, and that seems like a promising approach. But I'd like to check I'm understanding the idea correctly. Given sufficient samples in a sequence, we would expect the proportion of time spent in each state to tend to the stationary distribution. So one could test the distributions of occupancies in the two sequences for similarity. (I have the vague sense a two-sample Kolmogorov-Smirnov might suit for this, but please set me right about that.) The thing is, I'm not sure how good this can be as evidence. If the distributions are very different, that seems like a reasonable strike against the underlying processes being the same, but what if they're very similar? Can we draw conclusions in that direction? **Q3: Does a good fit of occupancy distributions tell us anything useful?** It seems like there could be an infinite number of processes that will tend to the same stationary distribution. I _think_ this is very unlikely in practice, and that different systems will tend to have distinctly different behaviour, but still it's worth considering. Finally, we will often have a model of the underlying process that we are looking for, although it may not be perfect. So we could compare each sequence with the expected behaviour of the model, instead of with each other. We may also have more than two sequences to test. **Q4: Is it better to compare multiple sequences to a single model, even an approximate one, or are we better off comparing data directly?**


title: Why do categorical predictor variables in regression need to be recoded as multiple predictors?

text: I'm learning about machine learning using Python's library scikit learn, and in their tutorial here they mentioned about a categorical variable `color` which can have values `purple`, `blue` and `red`. What is the reason behind using 3 boolean variables `color#purple`, `color#blue` and `color#red` instead of having the single variable `color` but mapping the values `purple`, `blue`, `red` to `1`, `2`, `3`? Will doing either way have any effects on the regression fitting/prediction?


title: PCA and exploratory Factor analysis on the same data set

text: I would like to know if it makes any logical sense to perform principal component analysis (PCA) and exploratory factor analysis (EFA) on the same data set. I have heard professionals expressly recommend:   1. Understand what the goal of analysis is and choose PCA or EFA for the data analysis   2. Having done one analysis there is no need to do the other analysis. I understand the motivational differences between the two, but I was just wondering if there is anything wrong in interpreting the results provided by PCA and EFA at the same time?


title: Interpreting plot.lm()

text: I had a question about interpreting the graphs generated by plot(lm) in R. I was wondering if you guys could tell me how to interpret the scale-location and leverage-residual plots? Any comments would be appreciated. Assume basic knowledge of statistics, regression and econometrics.


title: What is the distribution of maximum of a pair of iid draws, where the minimum is an order statistic of other minima?

text: Consider $n\cdot m$ independent draws from cdf $F(x)$, which is defined over 0-1, where $n$ and $m$ are integers. Arbitrarily group the draws into $n$ groups with m values in each group. Look at the minimum value in each group. Take the group that has the greatest of these minima. Now, what is the distribution that defines the maximum value in that group? More generally, what is the distribution for the $j$-th order statistic of $m$ draws of $F(x)$, where the kth order of those m draws is also the pth order of the n draws of that kth order statistic? All of that is at the most abstract, so here is a more concrete example. Consider 8 draws of $F(x)$. Group them into 4 pairs of 2. Compare the minimum value in each pair. Select the pair with the highest of these 4 minima. Label that draw "a". Label the other value in that same pair as "b". What is the distribution $F_b(b)$? We know $b>a$. We know a is the maximum of 4 minimums of $F(x)$, of $F_a(a) = (1-(1-F(x))^2)^4$. What is $F_b(b)$?


title: Can statistics really be used in education to determine test scores?

text: I am a middle school teacher. My district prides itself on having the kids take 3 tests during the year to predict their score, and my teaching, on "the real test" at the end of the year. My question: Is this practice even real? Or is the district fooling themselves in thinking these numbers mean more than what they do?


title: how to add instances data into test arff file with the same attribute and relation with training file

text: I have training arff file, and also I have instances data at hand. Now I want to add the instances data onto another file, called test file, and this test file has the same relation header and attribute with previous training file. How to realize that?


title: Fitting GARCH Model

text: I'am getting more and more familiar with this kind of model (and others models too). I'm now used to fit this model with my data (rmgarch package in R). How is it done ? What is the theory behind parameters estimators ? I think about a maximum likelihood process, but I would like a confirmation and precision.


title: GRM or mixed effect models

text: I need your help. In experiment I have measured the grow of one plant on 40 locations. At one location different number of plants were measured, but always the same species. The distribution of response variable is log-normal. I wanted to show the distribution in the post but it is not possible because I'm new user of stats.stackexchange. I have additional explanatory variables like altitude, soil type, precipitation ect. and would like to prepare the model of variables that explain the grow of this plant. Because the distribution of response variable is not normal I thought to use the GLM (`glm` function in R with `family=gamma`). Should I include also the grouping variable "location" and because of this the use of Mixed Effect Models (`nlme` function in R package `lme`) would be more appropriate? How would the code look like? Could you recommend me some literature? I hope that the question is well described. Thank you in advance.


title: Which is the better method to do forecast..1-step or h-step ahead?

text: I am using forecast() package in R to predict future values. I have a time series data for approx 6-7 years. First, I split the data into training set and test set. Test set contained values of the last 12 months of original data. Then, applied ets() and auto.arima() on the training set. This can be done by iterating 12 times using 1-step ahead OR by running once to get 12 points at the same time, which is my question no. 1 Then, the output received is compared with the test set to check which of the two (ets or auto.arima) gives minimal error in comparison with test set. The selected model class is then applied on the original data set to predict future 12 values. So, here arises my question no. 2 In Short, I have two things to do   1. Compare accuracy measures for test set of ets and auto.arima to select one   2. Predict future values using the derived model from the first step **Questions**   1. Should I use one-step ahead or 12-step ahead forecasts to generate data for test set?   2. Should I use one-step ahead or 12-step ahead forecasts to predict future values?   3. Should the method be same in both the steps mentioned above? Kindly help.


title: Using Cleveland dot plots to visualize time-series data

text: When I'm dealing with time-series data I'm generally thinking about visualizing that data with a bar graph (small n) or a line plot (large n). For example, I might create something like the below: ![enter image description here](http://i.stack.imgur.com/au0Li.png) However, is there any instance in which a dot plot can be used to visualize time-series data? ( a clevland dot-plot is what I mean, with dates on y-axis and values on the x-axis) Here's the data and code I used for the previous graph.               conv = c(10, 4.76, 17.14, 25, 26.47, 37.5, 20.83, 25.53, 32.5, 16.7, 27.33)     click = c(20, 42, 35, 28, 34, 48, 48, 47, 40, 30, 30)          dat <- data.frame(date=c("July 7", "July 8", "July 9", "July 10", "July 11", "July 12", "July 13",                       "July 14", "July 15", "July 16", "July 17"), click=c(click), conv=c(conv),                       stringsAsFactors = FALSE)          dat          ggplot(dat, aes(as.character(date), conv)) +  geom_bar(fill="#336699", colour="black") + ylim(c(0,50)) +             opts(title="Conversion Rate") +             opts(axis.text.y=theme_text(family="sans", face="bold", size=10)) +             opts(axis.text.x=theme_text(family="sans", face="bold", size=8)) +             opts(plot.title = theme_text(size=15, face="bold")) +             xlab("") + ylab("")      EDIT: My question may not have been clear. I'm NOT asking how to generate a cleveland dot plot. I'm asking whether it's all right to use a cleveland dot plot to visualize time series data. According to the 'statistical visualization rulebook', are cleveland dot plots a good way to represent time series data?


title: Inferring parameters for a regression with features of both multivariate probit and ordinal regression?

text: I am dealing with data which is generated by a complex process, which I elaborate below; I am trying to answer one or more of the following questions- a) what is the right literature to look for previous work addressing my problem b) has someone solved my problem already c) if this problem is novel (which I personally doubt), what is the right approach to solve it? For $N$ individuals, with $M$ associated covariates, I observe outcomes in $J$ possible classes; thus my design matrix $X$ is of dimension $M$ x $N$ and my response matrix $Y$ is of dimension $N \times J$. For a given class $j$ there is a latent variable $U_{i,j}$ s.t. when $U_{i,j} < 0$, $Y_{i,j}=0$ and $Y_{i,j}=1$, otherwise . Up to this point, if $U \sim MVN(\beta X, \Sigma)$, I have described a standard multivariate probit regression scenario. However, if a given $U_{i,j}$ is greater than all other classes $U_{i,-j}$ _and_ greater than some constant $c_1$, then $Y_{i,j} =t$ if $c_{t-1} < U_{i,j} < c_t$ where $t\in \{1,2,...,T-1\}$; $Y_{i,j}=T$ if $U_{i,j} > c_{T-1}$. Here the cutpoints $c_t$ are unknown. On its own, the latter half corresponds to an ordinal regression problem. However, together this problem seems to be more complex than either probit or ordinal regression. Below is R code to generate data from this process, where I am generating the values of $U$ directly without specifying a linear relationship with $X$:               library(MASS) # for mvrnorm          N <- 1e3 # num. of observations     J <- 12  # num. of classes     cuts <- c(1.5, 2, 2.25, 2.5, 2.7, 3) # T=6          S <- rWishart(n=1, Sigma=diag(J), df=20) # sample cov. mat.     U <- mvrnorm(n=N, mu=rep(0, J), Sigma=S[,,1]) # generate U     Y <- apply(U, 1, function(u) { # generate Y from U          u <- as.numeric(scale(u)) # s.d. is 1 for identifiability          y <- ifelse(u < 0, 0, 1)  # U's extreme enough are observed          cond <- u > cuts[1] & u == max(u) # most extreme U's are ranked          if(any(cond))             y[which(cond)] <- findInterval(u[which(cond)], cuts)          return(y)       }     )      In my case, $N$ is approximately $3500$, $J=12$ and $M=30$, and $T=6$. The outcomes of the classes are correlated. Furthermore, the individuals are hierarchically related. Initially, I thought I could use JAGS to solve this problem. Before implementing the full model, I tried a hierarhcial multivariate probit regression. However, JAGS seemed to take an inordinately long time to sample from the posterior. Thus, I have since been searching through the simulated likelihood literature e.g. GHK, but have found no work reproducing the model I have described above. I am investigating approaches such as approximate Bayesian computation (ABC), indirect inference, and simulated method of moments. Am I missing something obvious? Are there pre-existing solutions to my problem? Are there other clever approaches I could use instead (perhaps something like PLS or CCA)? Thank you for your help CrossValidated folks!


title: Moving average of irregular time series data using R

text: I'm attempting to dig out some metrics that look at how reliably clients connect to a service. The raw data is in the form of "client A, came online|offline at time X". The connection is highly unreliable, and I want some type of moving average to show whether the connection is improving or not over time. Clients are not always connected, so simply going offline does not mean it's a fault. So far, I've taken then data and applied some assumptions to help simplify it, I assume that if a client reconnects within a minute of disconnecting then that is a fault. These I've modelled as a simple impluses, ie. "client A had fault at time X". The part I'm struggling with is how to turn this plot into a moving average (I'm playing with R to crunch the numbers). I believe I should be able to do this with a low pass filter, or use the zoo package and rollmean. However, I don't know how to handle the cases where the client simply didn't want to be online. Any suggestions?


title: Question on cumulative probability

text: In an article I'm reading the probability of an observed point falling in bin i is written as: $P_i=\frac{m_i}{\sum\limits_j m_j}$ where $m_i$ is the number of **model** points in bin i. Then, the **cumulative probability** of obtaining the entire data set $n_i$ is written as: $P=\prod\limits_i \left( \frac{m_i}{\sum\limits_j m_j} \right)^{n_i}$ where $n_i$ if the observed value of bin i. My question is: how was equation 2 obtained? Edit: this is the article in question. I didn't add much information in the question because I wanted to keep it as simple as possible. Equations in the article are (14) and (15).


title: What is the relationship between two points on probability density function?

text: The Wikipedia entry for Probability Density Function states that the PDF "describes the relative likelihood for this random variable to take on a given value." Two questions:   1. Does that mean that the ratio of two points reflects the difference in probability? For instance, on a standard normal distribution f(0.0)=0.3989, while f(1.0)=0.2420. Does this mean that an outcome of 0.0 is 1.65 times more likely than an outcome of 1.0? (I know that the actual probability of any exact value is zero.)   2. If one can compare points in this way, can you then compare points from two different densities? For instance, on a standard normal distribution f(1.0)=0.2420, whereas f(1.0) for a Student's-T distribution with 1 degree of freedom is 0.1592. Does this mean an outcome of 1.0 is 1.52 times more likely on a standard normal distribution than a T distribution?


title: Correcting experiment results

text: We have performed a microarray screening of about 200 samples. In each sample we measure about 100 different variables. For technical reasons the screening of these 200 samples was divided into two batches with a couple of weeks interval between them. When all the data has been collected, I have performed principle component analysis (PCA) on the 200 x 100 table. When we look on linear projection of first 4 components (responsible for ~70% of the variability), we see a clear division between the two experiment batches. An illustration to what I have can be seen here: http://img594.imageshack.us/img594/3687/pca.png What are the accepted techniques to address this issue?


title: Few machine learning problems

text: In a particular application I was in need of machine learning (I know the things I studied in my undergraduate course). I used Support Vector Machines and got the problem solved. Its working fine. Now I need to improve the system. Problems here are   1. I get additional training examples every week. Right now the system starts training freshly with updated examples (old examples + new examples). I want to make it incremental learning. Using previous knowledge (instead of previous examples) with new examples to get new model (knowledge)   2. Right my training examples has 3 classes. So, every training example is fitted into one of these 3 classes. I want functionality of "Unknown" class. Anything that doesn't fit these 3 classes must be marked as "unknown". But I can't treat "Unknown" as a new class and provide examples for this too.   3. Assuming, the "unknown" class is implemented. When class is "unknown" the user of the application inputs the what he thinks the class might be. Now, I need to incorporate the user input into the learning. I've no idea about how to do this too. Would it make any difference if the user inputs a new class (i.e.. a class that is not already in the training set)? Do I need to choose a new algorithm or Support Vector Machines can do this? PS: I'm using libsvm implementation for SVM.


title: How to interpret residual plots from time series regression

text: I am doing a time series regression between 2 variables. I used the dynlm library in R. I'm trying to understand how to interpret the results. Could you please point out where I am getting it wrong: 1) The R squared seems very low -- this indicates a weak linear relationship. Or too many outliers. 2) Normal QQ plot shows that there are a significant number of outliers -- at a later date in the time series? 3) Does the 'Residuals vs Fitted' plot show that the model is a pretty good fit for most of the data, except for those outliers? Any suggested readings (esp open-sourced materials available online) also appreciated. * * *               Call:     dynlm(formula = P ~ L(B))     Residuals:         Min      1Q  Median      3Q     Max      -63.711 -27.687 -14.907   2.364 146.157      Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept) 17.485051  13.500833   1.295    0.197          L(B)         0.019422   0.002384   8.146 5.84e-14     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1      Residual standard error: 48.17 on 182 degrees of freedom     Multiple R-squared:  0.2672,    Adjusted R-squared:  0.2632      F-statistic: 66.36 on 1 and 182 DF,  p-value: 5.838e-14      ![enter image description here](http://i.stack.imgur.com/xJbHF.png) ![enter image description here](http://i.stack.imgur.com/6cpT2.png) ![enter image description here](http://i.stack.imgur.com/JZbmX.png)


title: lasso and cross-validation (theoretical results)

text: is there any theoretical result which says that use the minimum of the cross- validation as value for the lasso penalty is a good choice? I would like something like $P(S_0 \subset \hat S_{lasso}(\lambda_{cv}))\rightarrow 1$ where $S_0$ is the set of true variable. Where can I find it?


title: Bayesian inference and degrees of freedom

text: While learning frequentist linear regressions, one thing the professors always talked about was about the _number of degrees of freedom_ , I never saw this expression in a bayesian book though. Perhaps because bayesian methods don't need this number to infer things like variance and such? My question is: is the number of degrees of freedom equals the number of parameters in an hierarchical bayesian model and if it's not, is there something equivalent one can calculate? In particular, I'm interested in when a model is overidentified in a hierarchical framework. For example, if I have 1000 observations and about 10 possible competing models with about 100 parameters each, if mix them all in an hierarchical model using, for example, trans-dimensional MCMC/Bayes factor, the will I have an overidentified model? My intuition says that it's possible that it won't, although the total number of parameters is greater than the number of observed parameters.


title: How to get percentage agreement between a group of factor columns?

text: I have a data.frame with 18 columns all of which are categorical variables containing the same factors e.g. "a", "b", "c". What I would like to produce is a matrix which is 18x18 and contains the proportion/percentage agreement between each column and the rest of the others. For example:- Input:               yoda <- data.frame(one=c("a","a","a","a"), two=c("a","b","a","a"), three=c("a","a","b","b"))      Output:               structure(c(1, 0.75, 0.5, 0.75, 1, 0.25, 0.5, 0.25, 1), .Dim = c(3L,      3L), .Dimnames = list(c("one", "two", "three"), c("one", "two",      "three")))      I could loop over the column names twice but I am sure that there must be a better way. I have looked into creating a distance matrix and correlation but none of the methods seem to fit the bill. -- EDIT -- Here is what I am doing currently:               library(plyr)     result <- aaply(data, 2, function(col) {       aaply(data, 2, function(col2) {         sum(col==col2)/length(col)       })     })


title: Using R for GLM with Gamma distribution

text: I currently have a problem understanding the syntax for R for fitting a GLM using the Gamma distribution. I have a set of data, where each row contains 3 co-variates ($X_1, X_2, X_3$), a response variable ($Y$), and a shape parameter ($K$). I want to model the scale of the Gamma distribution as a linear function of the 3 covariates, but I don't understand how to set the shape of the distribution to $K$ for each row of data. A situation that I think is analogous is that for a binomial distribution, the GLM requires that the number of trials ($N$) is known for each data entry.


title: Detecting pattern in panel Data R

text: I am struggling on panel data in R. I detected some serial correlation within my data, and now I would like to detect some pattern within it. There are 166 observations and 34 variables and the data looks like this:               id  1975  1976  1977  1978  1979  1980  1981  ... 2007      5   10    10    10    10    11    12    NA      29   1     2     2     2     2     2     2        NA      12  NA    NA    NA    NA    NA     5     6         5      15   1     1     4     3     2    NA    NA        NA      This panel is rather unbalanced and I reshaped it in order to analyze the specific years in order to detect some patterns. Is there any way to calculate or plot the data? Before I reshaped it, I tried some ggplot2 oder googleVis with the full data, but there are just too many groups (id) for a clear visualization. Is there a smart way to detect patterns of serial correlation within this kind of data frame?


title: Help with Theil-Sen regression function (written in R)

text: I have written a function (in R-language) to perform a Theil-Sen regression with large sample approximation (please see below). It seems to work well for the one-sided case as illustrated in Hollander & Wolfe (2nd Edition Ch. 9)               NP.lm <-function(dat, X, Y, alpha) {            dat <- dat[order(dat[,X]),]       n <- nrow(dat)           combos <- combn(n, 2)       i.s <- combos[1,]       j.s <- combos[2,]       Y.num <- dat[j.s,Y] - dat[i.s,Y]       X.dom <- dat[j.s,X] - dat[i.s,X]             Z.p <-  qnorm(alpha/2, lower.tail=F)       N <- (n*(n-1))/2       s <- (Y.num/X.dom)[X.dom != 0]       C.stat <- sum(sign(s))              slope <- median(s, na.rm=TRUE)                        intercept <- median((dat[,Y] - slope*dat[,X]), na.rm=TRUE)                 C.star <- C.stat / sqrt(n*(n-1)*((2*n)+5)/18)       p.val <- pnorm(-abs(C.star))            out <- data.frame(colnames(dat[X]), colnames(dat[Y]), slope, intercept, p.val)            out          }       # Example usage               set.seed(123) # to make the example fully reproducible     sigma <- as.matrix(data.frame(c(1, .7), c(.7, 1)))     N <- 50     dat <- data.frame( matrix(rnorm(N*nrow(sigma)), N, nrow(sigma)) %*%          chol(sigma) ) # cor structure via Choleski Decomp.     colnames(dat) <- c("X", "Y")     summary(lm(Y~X, data=dat))     NP.lm(dat=dat, X=1, Y=2, alpha=0.05)            colnames.dat.X.. colnames.dat.Y..     slope  intercept        p.val     1                X                Y 0.6719583 0.08337517 1.227402e-07          # This p-value can't be right      Could someone please tell me what I am doing wrong, and how to make the test 2-sided. Thank you so much for your help!


title: Simulate data for power analysis of logistic regression model - include covariance variance of variables?

text: I've tried to simulate data for a power analysis of a logistic regression. The results of the power analysis look reasonable: power=90% for a sample of 6000 persons. But I feel that the analysis lacks something. So, my question is: when generating the data should I include something about how the variables are correlated, or their covariance, other than just defining their linear relationship as I have done in the example below, and if so where do I write that into the code? I know other questions look like this but I'm not confident that they answer this question.               library(plyr) # functions     ## Define Function     simfunktion <- function() {        # Number in each sample       antal <- 6000       beta0 <- log(0.16) # logit in reference group       beta1 <- log(1.1)  # logit given smoking       beta2 <- log(1.1)  # logit given SNP(genevariation)       beta3 <- log(1.2)  # logit for interactioncoefficient for SNP*rygning        ## Smoking variable, with probabilities defined according to empirical studies.       smoking  <- sample(x = 0:2, size = antal, replace = TRUE, prob = c(40,25,40))        ## SNP variables with probabilities defined according to empirical studies       SNP      <- sample(x = 0:2, size = antal, replace = TRUE, prob = c(40,40,20))        ## calculated probabilites given the model:       pi.x     <- exp( beta0 + beta1*smoking + beta2*SNP + beta3*smoking*SNP) /                    ( 1 + exp(beta0 + beta1*smoking + beta2*SNP + beta3*smoking*SNP) )        ## binoial events given the probabilities:       sim.y    <- rbinom( n = antal, size = 1, prob = pi.x)         sim.data <- data.frame(sim.y, smoking, SNP)        #################### p-value of the interaction is extracted:        ## the model is run:       glm1     <- glm( data = sim.data, formula = sim.y ~ smoking + SNP + smoking:SNP,                        family=binomial )        ## p-value of the interactionterm is extracted:       summary(    glm( data = sim.data, formula = sim.y ~ smoking + SNP + smoking:SNP,                        family=binomial ))$coef[4,4]     }     pvalue     <- as.vector(replicate( 100 , simfunktion()))     mean(pvalue < 0.05)


title: What specification when dependent variable is a fraction?

text: Basically what the title says; I'm carrying out an exploratory regression on the determinants of board independence. Is using a GML regression the way to go? Thanks in advance


title: Closed form solution for the expectation of the square root of a hypergeometric variate

text: As the title states, is there a closed form formula for the expectation of the square root of a hypergeometric variable. Edit: Closed form approximate solutions based on related distributions or expansions are also welcome.


title: Using components from PCA as index / measurement?

text: I have used Principal Component Analysis to find main components that desribe the perceptions of the environment. The goal is to use those main components as index / measurement for another study. My problem is if I find a "parsimonious description of the variables" and use that as an index for the next study, how can I be sure that people what that term actually means (which variables it included)? So, I found that e.g. pleasant, comfortable, relaxing, welcoming etc. to be in the same component and decide to call them "positive perception" as use it in the next study as "How do you find those environments in terms of "positive perception"? How are people going to know that "positive perception" means pleasant, comfortable, relaxing and welcoming? Do I need another validation study to confirm that the new term I came up with is appropriately describing the variables? Would it be enough to ask a number of people what word they think would describe those variables best? Or what they think "positive perception" means? Many thanks in advance!!!


title: How one defines the expected value of two parameter distribution?

text: I have been given two discrete distributions $X$ and $Y$ and values $\mathbb{P}(x_i,y_j)$ for every $i,j$. What is the formula to compute for example $\mathbb{E}(X+Y)$ and $\mathbb{E}(XY)$? My lecture notes do not explain this but I have to count those. In the one parameter case I know the formula is $\mathbb{E}X=\sum_{x\in\mathbb{R}} xf_X(x)$ and I found from the notes that $f_{XY}(x,y)=\mathbb{P}(X=x\cap Y=y)$ but what is $f_{X+Y}(x,y)$. And I guess that $xf_X(x)$ should be of the form $g(x,y)f_{X+Y}(x,y)$ but what is $g$?


title: Can EFA be conducted on dependent variables?

text: I'm working on my dissertation which is a multiple case study examining the use of performance information. My dissertation chair states that I need to use confirmatory factor analysis; however, after months of learning SPSS and completing other analyses through the system, I am just learning the CFA cannot be completed in SPSS without AMOS which the school does not support AMOS. I've been working with a committee member and have been examining 2 and 3 factors related to my dependent variables which would help examine the underlying dimensions and associations with each factor but I'm concerned that I'll still be required to use CFA and have to learn a new statistical system. So my questions are: 1) can EFA be conducted on dependent variables and 2) is there a strong argument to use EFA over CFA? Thanks!


title: Software needed to scrape data from graph

text: Anybody have any experience with software (preferably free, preferably open source) that will take an image of data plotted on cartesian coordinates (a standard, everyday plot) and extract the coordinates of the points plotted on the graph? Essentially, this is a data-mining problem and a **reverse** data- visualization problem.


title: Tool for reading values from the chart

text: I am given a paper with Kaplan-Meier curve. I'd like "reverse-engineer" the chart to get the original data (without the censoring events, as they are difficult to discern in my low quality copy). I've seen the tool for this task on the Net, but it was long time ago and forgot its name. Can someone please remind me, how to find it? Thanks.