title: How to get the points of a linear regression line?

text: I do a simple linear regression `mod <- lm(y ~ x)` and I plot its residuals, doing `plot(mod$residuals)`. Two questions:   1. I need to plot the line of the regression in the chart where I plotted the residuals, how can I do it?   2. I need to know the points of this line, because I should plot this line in another application, so I need those points (as a vector), how can I do? Thank you!


title: Extreme value simulation with Monte Carlo

text: I would like to seek your help with some questions to simulating extreme values. For example, I have written the following R code to generate QQplots for a normally distributed data, varying the size of the sample from 10 up to 1B. As seen in the graph as below, I had to go up to 1B samples to see extreme values @ 6 Sigma from the mean.               par(mfrow=c(2,3))     for(i in c(10, 100, 1e+3, 1e+4, 1e+5, 1e+6, 1e+7, 1e+8, 1e+9)){       data <- rnorm(i, mean = 0, sd = 1)       qqnorm(data, main=sprintf("Sample Size=%d", i)); qqline(data, col='red')     }      **Question1** : I am interested in extremes values only (tails). Though MC sampling generates most of the data close to the mean and I am wondering which method I can use to crate a few (up to a thousand maybe) samples and yet get some extreme values up to 6 sigma from mean. **Question2** : the QQplot seems to be significant only after ~ 1K points. Also, no matter the sample size, the QQplot always shows some tail-offs. this must be some sampling error but I'm not sure. Would you please point me to any literature that explains this behavior ? ![enter image description here](http://i.stack.imgur.com/mLTlR.png)


title: Correlated state-space models

text: I'm struggling with Reinsel's book "Elements of Multivariate Time Series Analysis," because I thought that it would be a good idea to switch from Vector ARMA to state-space representations; particularly after reading What are disadvantages of state-space models and Kalman Filter for time-series modelling?, and after struggling with parameter estimation of VARMA processes. The problem is that I can't see how correlation between two time series can be modeled in state-space representations. Is it that the correlation appears in the error covariance matrices? It sounds weird to me, particularly coming from VARMA. (Switching to state-space models is being harder than I thought!)


title: Expected value of a random variable

text: Random variable $X$ has the probability density function \begin{equation*} f\left( x\right) =\left\\{ \begin{array}{ccc} n\left( \frac{x}{\theta }\right) ^{n-1} & , & 0<x\leqslant \theta \\\ n\left( \frac{1-x}{1-\theta }\right) ^{n-1} & , & \theta \leqslant x<1% \end{array}% \right. \end{equation*} show that if $k\in \mathbb{N}$ \begin{equation*} \mathrm{E}\left( X^{k}\right) =\frac{n\theta ^{k+1}}{n+k}+\sum% \limits_{i=0}^{k}\left( -1\right) ^{i}\binom{k}{k-i}\frac{n}{n+i}\left( 1-\theta \right) ^{i+1} \end{equation*} It is easy to find the first term of $\mathrm{E}\left( X^{k}\right) $ but i couldn't find the second one. i think i have to use beta distribution properties. i tried to simulate the integral \begin{equation*} \int\nolimits_{\theta }^{1}nx^{k}\left( \frac{1-x}{1-\theta }\right) ^{n-1}% \mathrm{d}x \end{equation*} to beta pdf using $u=\frac{x-1}{\theta -1}$ transformation but i couldn't get a reasonable result. After this transformation is applied i have to find \begin{equation*} \int\nolimits_{0}^{1}u^{n-1}\left[ 1+u\left( \theta -1\right) \right] ^{k}% \mathrm{d}x \end{equation*} but i couldn't. Also I tried to use the equality \begin{equation*} x^{k}=1+\left( x-1\right) \sum\limits_{n=0}^{k-1}x^{n} \end{equation*} but i couldn't get the result.


title: Is there an equivalent to Lower bound of Wilson score confidence interval for variables with more outcome

text: After reading "How Not to Sort by Average Rating" (http://www.evanmiller.org/how-not-to-sort-by-average-rating.html), I was curious to know if there was the same thing for variables with more than two outcomes (0,1) or even continuous variables. For example, how would you generalise the lower bound to the Amazon problem ? Clearly there are 5 outcomes (one for each number of star given by the user). What measure would you use to make the 4.5 stars with 2000 votes better rated that the 5 star with 2 votes ? Also, it seems to me that this kind of problem could have a bayesian interpretation. I mean using the formula in "How not to sort" is not far from setting a prior in the distribution, maybe a Bernoulli with parameter inferred on the whole dataset / category the item belongs to ? Does anyone know a reference for this particular problem ?


title: How low multiple R-squared value is enough to reject a model?

text: I'm doing a simple linear regression I tried:               > mod <- lm(rnorm(100,sd=2) ~ rnorm(100,sd=2.1))     > summary(mod)          Call:     lm(formula = rnorm(100, sd = 2) ~ rnorm(100, sd = 2.1))          Residuals:         Min      1Q  Median      3Q     Max      -4.0396 -1.0698  0.0803  0.9823  5.4893           Coefficients:                          Estimate Std. Error t value Pr(>|t|)       (Intercept)            0.0280     0.1868   0.150   0.8812       rnorm(100, sd = 2.1)  -0.1523     0.0900  -1.692   0.0939 .     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 1.868 on 98 degrees of freedom     Multiple R-squared: 0.02838,    Adjusted R-squared: 0.01846      F-statistic: 2.862 on 1 and 98 DF,  p-value: 0.09387       as you can see the multiple R-squared is very very low. I must to reject that model, but what is the level to understand if I can accept the model or not? Thank you


title: Confidence interval for the standard deviation of a Normal distribution with known mean

text: Suppose $\textbf{Y} = (Y_{1}, ... , Y_{n})$ is a random sample from the $N(\mu, \sigma_{0}^{2})$ distribution where $\mathrm{E}(Y_{i}) = \mu$ is unknown but $\mathrm{SD}(Y_{i}) = \sigma_{0}$ is known. It's possible to construct a confidence interval for the mean of a Normal distribution with known standard deviation using the pivotal quantity: $$Q = \frac{\bar{Y} - \mu}{\sigma_{0} / \sqrt{n}}$$ Then, suppose $\mathrm{E}(Y_{i}) = \mu_{0}$ is known but $\mathrm{SD}(Y_{i}) = \sigma$ is unknown, why can't we just use the following pivotal quantity to construct a confidence interval for the standard deviation with known mean? $$Q = \frac{\bar{Y} - \mu_{0}}{\sigma / \sqrt{n}}$$


title: How to speed up the training in neural network when mini-batch training is used?

text: Can anyone give me some ideas on possible techniques to speed up the training process of multilayer artificial neural network if the training involves mini- batch? So far, I understand that stochastic training probably leads to a faster convergence, but if we have to use mini-batch training, is there any way to make the convergence faster? (Some pointers to relevant papers will also help!) Thank you!


title: Differentiate two corelated variables

text: I am working on a prediction problem. I have a programming background and a small statistics one. If I have two dependent variables and one independent one, how do I work out if the two are co-related? For example, x1 has a correlation coefficient of .6 and x2 has one of .4. I hypothesize that x1 and x2 are not totally dependent (identity) and not totally independent to each other. Is there a way to analyze this situation mathematically to eventually decide weights of a final regression equation? I suppose I could put one x into bins and make a piecewise equation based on certain cutoffs. But, this seems not to be optimal and in reality there are many x's. I have not found any good literature on how to differentiate weights of dependent variables. I know machine learning is designed for this but I would like a 'old school slide ruler' approach where I can understand the differences of variance. I am in a exploratory phase where given two (and really many more) variables figure out the outcome. If the two variables are co-related then some of the variance is overlapping. I would like to know how much overlaps and how to weight the variables in an ending regression equation. (this my not be the right terminology so please bear with me.)


title: Unbiased estimator of variance of binomial variable

text: $Y_{1...n}\sim \operatorname{Bin}(1,p)$, iid, and I need to find an unbiased estimator for $\theta=\operatorname{var}(y_i)$. I did some calculations and I think that the answer is $p(1-p)-\frac{p(1-p)}{n}$   * Is this correct?   * If not, how can I find an unbiased estimator?


title: For data reduction, what is this technique called?

text: I have observations X along with their labels Y. I then create a histogram of Y. I then remove observations such that the histogram still retains the same distinct shape. Does anyone know what data reduction technique this is called?


title: tbats model for multiplicative seasonality

text: In De Livera, Hyndman & Snyder (2011), a TBATS model (exponential smoothing state space model with Box-Cox transformation, ARMA errors, trend and seasonal components) for additive seasonality has been discussed. Can TBATS be used with multiplicative seasonality? Are there any restraints or numerical difficulties in TBATS or BATS with multiplicative seasonality?


title: Generalised least squares: from regression coefficients to correlation coefficients?

text: For least squares with one predictor: $y = \beta x + \epsilon$ If $x$ and $y$ are standardised prior to fitting (i.e. $\sim N(0,1)$), then:   * $\beta$ is the same as the Pearson correlation coefficient, $r$.   * $\beta$ is the same in the reflected regression: $x = \beta y + \epsilon$ For generalised least squares (GLS), does the same apply? I.e. if I standardise my data, can I obtained correlation coefficients directly from the regression coefficients? From experimenting with data, the reflected GLS leads to different $\beta$ coefficients and also I'm not sure that I'm believing that the regression coefficients fit with my expected values for correlation. I know people quote GLS correlation coefficients, so I am wondering how they arrive at them and hence what they really mean?


title: Would two different pivotal quantities of the same parameter give the same confidence interval?

text: I couldn't think of an example but it would be great if someone could give one.


title: Sum of normal truncated random variables

text: Suppose I have $n$ independent normal random variables $$X_1 \sim \mathrm{N}(\mu_1, \sigma_1^2)\\\X_2 \sim \mathrm{N}(\mu_2, \sigma_2^2)\\\\\vdots\\\X_n \sim \mathrm{N}(\mu_n, \sigma_n^2)$$ and $Y=X_1+X_2+\dotsm+X_n$. How would I characterize the density of $Y$ if the distribution of each $X_i$ is each truncated to within $(\mu_i - 2\sigma_i, \mu_i + 2\sigma_i)$? In other words, I'm sampling from $n$ independent normal distributions, discarding samples not within $2\sigma_i$ of each mean, and summing them. Right now, I'm doing this with the R code below:               x_mu <- c(12, 18, 7)     x_sd <- c(1.5, 2, 0.8)     a <- x_mu - 2 * x_sd     b <- x_mu + 2 * x_sd          samples <- sapply(1:3, function(i) {       return(rtruncnorm(100000, a[i], b[i], x_mu[i], x_sd[i]))     })          y <- rowSums(samples)      Is there any method for generating the density of $Y$ directly?


title: Random Variable as a function of mean probability

text: Suppose that in a real communication system we make a scheduling decision $\omega$ at any time slot $t$ that affects the probability of loss of a packet (i.e. the probability of loss depends on our decision and we assume independence from slot to slot). We define as the mean packet loss probability as $P_{loss,t} = \frac{total\ packets\ lost\ until\ t}{total\ packets\ transmitted\ until\ t}$ Does it make any sense to define a random variable as $X_t=a\times P_{loss,t} + b$ with $a,b$ real constants?


title: Identifying formulas from MANOVA in a research paper

text: I need to identify the formulas for computing the group differences with Bonferroni corrections and Cohen's d effect size in the following paper (chart on the top of page 6): Research Paper What formula was used to do the group comparisons? Also, what were the degrees of freedom?


title: Efficiently finding approximate factors in expectation propagation

text: I've been trying to wrap my head around expectation propagation for a while now, but I'm struggling a bit with finding approximate factors in an efficient way and 'fully-factorized' approximations. I'm currently trying to write my own implementation (to check if I've fully understood everything) which is based on this paper of Minka and Infer.NET (as a reference for how it could be implemented). * * * **Finding approximate factors – computational complexity** My understanding is that factors in a factor graph represent the CPDs of my Bayesian network. For example – given random variables $x_1,x_2,x_3$ – I could define $x_3 = x_1 + x_2$. The corresponding factor $f(x_1,x_2,x_3)$ would then be defined as $p(x_3|x_1,x_2) = (pdf_{x_1} \cdot pdf_{x_2})(x_3)\rightarrow$ the convolution of $x_1$ and $x_2$. Depending on the distributions involved marginalizing out other variables might become intractable in a bigger network as the integral gets too complicated. Because of that, you try to approximate the factors (and the joint distribution) by distributions of the same exponential family (f.e. Gaussians). However, when refining the approximations and finding $\arg\min KL(f_a(\boldsymbol{x}) q^{\backslash a}(\boldsymbol{x}) || \hat{f}_a(\boldsymbol{x}) q^{\backslash a}(\boldsymbol{x})$ (equation 52 in the paper), I still have to take the integral over all random variables which still seems to be a computationally expensive task (especially because it has to be done several times). Even when I'm just matching the expectations of the features of the exponential family chosen, this results in a quite computational heavy integral. Am I misunderstanding something here or missing out some simplifications? * * * **Fully-factorized approximation** In the paper, many algorithms 'use a fully-factorized approximation'. After my understanding this means that all factors are split up into separate factors such that each factor only depends on a single variable. So $f(x_1,x_2,x_3) \approx q(x_1)q(x_2)q(x_3)$. Am I right when I assume that this is done 'automatically' by the algorithms? Do the messages from factors to variables represent the factors of a fully- factorized graph? * * * Thanks in advance for any replies. If I'm talking nonsense feel free to tell me – this whole thing is messing with my head a bit at the moment and I keep confusing myself ;)


title: What are good resources for discussing solutions to "Introduction to Statistical Learning in R"?

text: This question is more about learning statistics than about statistics _per se_ , so please suggest another forum if this one is not appropriate. I found this textbook (ISLR by James, Witten, Hastie, and Tibshirani) online and it seems like a great resource. I read a few chapters and then realized that I wasn't getting good comprehension. So now I've decided to answer the questions at the end of each chapter and write them up in LaTeX/knitr. I'm through chapter 3. It would be amazing to discuss the solutions with other people who are using the same book. Maybe I'm getting them all wrong! Can anyone suggest a resource/forum/website where such a discussion might take place?


title: Interpretation of NbClust result

text: The plots show the output of `NbClust()`. By looking at the plot, is that correct to say that `k=5` is the optimal number of clusters?               > nc <- NbClust(m3, distance="euclidean", min.nc=2, max.nc=20,      +               method="kmeans", index="dindex")     [1] "*** : The D index is a graphical method of determining the number of clusters.           In the plot of D index, we seek a significant knee (the significant peak in Dindex          second differences plot) that corresponds to a significant increase of the value of          the measure."            "All 10829 observations were used."      ![enter image description here](http://i.stack.imgur.com/aSX7T.jpg)


title: Variance of ARMA(2,1)

text: How do I find the variance of this ARMA(2,1) model? $$ X_t=0.5X_{t-2} + e_t + e_{t-1} $$ I know the formula for ARMA(1,1) but when trying to solve I just keep getting an endless path of higher auto-covariances to find.


title: Meta-Analysis of GEO data using GEO2R

text: I'm using data from the GEO and GEO2R online analysis tool provided by NCBI. The GEO2R analysis gives the following scores: adjusted P, P, and t scores. I have 5-10 studies using the GEO2R tool. Is there a way to combine all these scores and filter out only significant genes? Should I convert the P or T scores into Z scores? Then I could calculate the $z=\frac{E(zk)}{\sqrt{K}}$ to create a consensus z -score (union of genes).


title: error structure in Generalized linear models when y is continuous data and errors not normally distributed

text: Lets say I have continuous y and x variable and I run a linear regression:               mdl1<-lm(y ~ x)      A generalised linear model should also give me the same parameters value if I do not specify the error structure (i.e. by default it assumes that the error structure is Gaussian)               mdl2<-glm(y  ~ x)      Both the above model should give me the same results (since in the `mdl2`, by default the error structure is gaussian) My question is if the residuals are not normally distributed in the `mdl1` (i.e. I do a shapiro.wilk test on `mdl1` residuals, which gives me a p-value of 0.02),               shapiro.test(rstandard(mdl1)      then in the glm what error family do I specify considering both y and x are continuous. What my understanding was I can specify `family=poisson` or `family=binomial` if my response variable was either count or proportion.               mdl3<-glm(y  ~ x,family="poisson") # when y is count data     mdl4<-glm(y  ~ x,family="binomial") # when y is proportion data      But in case of response variable being continuous and errors not normally distributed what error structure do I need to give?               mdl5<-glm(y ~ x, family=?????) # y is continuous data


title: What is the z-value of glm model parameters?

text: I know what the z-value of a single observation is, that is explained in Wikipedia. But what is the z-value of a parameter in a glm model?


title: Anomaly detection using exponential weighted moving average

text: I would like to detect anomaly using exponential weighted moving average. I don't have series of data points. All I have is EMA(t-1) and the data point of the current time(t) DP(t). From these data, I can calculate the new EMA(t). The EWMA constant will be 0.85 (assuming). Now I have two EMAs={EMA(t-1) ,EMA(t)} and DP(t). Is it possible to determine whether DPt is anomaly or not? using any logic? * * * For Eg,                 Time(t) DataPoint(t)    EMA(t-1)    EMA(t)       1       120               0          102         .        .                .           .      .        .                .           .      .        .                .           .      .        .                .           .      .        .                .           .          10      300              150        277.5         I don't want to take series of data points into calculation. At any time instance say 10, I want to know whether the data point 300 (in this case) is an anomaly or not. I also have the EMA(9)=150 and EMA(10)=277.5 for calculation (if needed) Is there any logic to calculate this? Thanks in advance! * * * Update: I thought of the below logic but not sure whether it would work for sure ABS(DP-EMA(new))> 3 *SD {EMA(old),EMA(new)} (or) ABS(DP-EMA(new))> 3 *EMWSD {EMA(old),EMA(new)} DP - Data Point   EMA - Exponential Weighted Moving Average


title: Order statistics of equal correlated continuous random variables

text: Suppose that $X_1, \ldots, X_n$ are mutlivariate normal with equal correlation $\rho$ and each of them are marginally distributed as $N(0,1)$. Let $X_{(1)}, \ldots, X_{(n)}$ be the corresponding order statistics. The distribution of $X_{(1)}$ and $X_{(n)}$ are easily found. What about the distribution of the other order statistics? Can anyone give some reference on this? Thank you. Hanna


title: Simultaneous Equation System for logit/probit?

text: Is it possible have an SES where the component equations are probabilistic, say, logit or probit? I am evaluating a number of quality metrics of services provided by a number of providers. The metrics are binary (pass/fail quality certification). The obvious approach would appear to be to estimate the logit/probit equations, conditional upon provider characteristics, jointly in something like an SUR, but I can't find anything like that in the discrete choice literature.


title: Summing normal instead of beta distributions, consequences for the density function of the sum?

text: **Background:** I've modeled a project effort prediction as a Google Spreadsheet template. Details of the Model: http://sites.google.com/site/effortprediction/methodology .Google Spreadsheet does not implement beta distribution functions. In PERT a beta distribution is used to estimate the effort of a task in Ideal Person Days (i.e.: one worker 8h a day without distraction). **Problem:** What consequences has use of normal distributions instead of beta distributions when the distribution are summed up to a single distribution. **Details** A _project_ consists of a list of _tasks_. Each task is _beta distributed_ and _independed_ , the sum of the tasks therefore should obey to the _central limit theorem_. The _project_ therefore is normal distributed. How does Expected Value and SD change for the project when I use normal instead of beta distributions for each task? I assume that the tasks distribution is skewed more often to the right than to the left. **Questions**   * Are the tails fater?   * Is the Expected Value higher or lower with normal distributions?    * Is the dispersion higher or lower?


title: Fitting an exponential distribution to data and finding third quartile problem

text: I have a data-set with range of `0 to 1`. I'm fitting a distribution to it in `MATLAB` using `fitdist(Data_set, 'exponential')`. I want find third quartile with `log(4)/λ` . When I calculate this, the result is `4.28738`As i mentioned my data set range is between `0 and 1`. why this value is bigger than maximum of my data range? How can i solve this problem? I think this is from tail o exponential distribution.


title: Probability generating function for negative values of random variables?

text: What if we have negative integral values for a random variable?Then is it possible to write a probability generating function for it? All definitions I have seen so far is for non negative integer values. I hope someone could assist me. Thanks


title: How to treat sentences including "and so on" when analyzing open-ended responses?

text: I'm analyzing responses from open-ended questions, by item (specifically, sections of an entire answer). I were asked to mention features or actions about a subject. I have sections that I reach to categorize. But I'm not sure how I must work with responses where I can find actions or characteristics but where the sentence ended up with "and so on", "and other features", "and more", etc. I registered these "part" in some responses independently of the matter in question. How to cope with such a situation? Could we just ignore that "part" of the answer or do we need to create a new category for them like a "generalizations"?


title: Meta-analysis and homogeneity -- what did these guys do?

text: I appreciate any insight into this meta-analysis. This is a meta-analysis on alogliptin efficacy and safety. In the 2nd paragraph of the discussion that make this statement: > Although heterogeneity testing showed a statistically significant > dissimilarity in the results of the included studies, sensitivity analysis > has shown the stability of the overall odds ratios with the withdrawal of > any of the study from the analysis without a significant improvement of the > heterogeneity. Thus, the credibility of the results of this meta-analysis > did not seem compromised. This is because; when the number of included > studies is small and heterogeneity is large, the robustness of the results > is best assessed with a sensitivity analysis What in the world do they mean by "sensitivity analysis has shown the stability of the overall odds ratios with the withdrawal of any of the study from the analysis without a significant improvement of the heterogeneity" Do they mean when one outlier study is removed the ORs don't change? What are they trying to communicate? Please help! Citation: Efficacy of alogliptin in type 2 diabetes treatment: a meta-analysis of randomized double-blind controlled studies (PDF) - _BMC Endocrine Disorders_ **2013** , 13:9


title: Exponential distribution hypothesis testing

text: I want exponential distribution hypothesis testing with `MATLAB`. Is any test or method for that?


title: Alternative to LPM and logit regression in Stata

text: I'm dealing currently with complex dataset that was already tested with Linear probability model and logit regression. I'd like to find an alternative for original regressions. The sample uses a binary dependent variable explained by many regressors. Besides the sample is quite big (~50 000 observations). Any suggestions will be welcome!


title: Finding limit distribution $T_n$

text: Suppose $X_1,X_2,\ldots,X_{2n}$ is a random sample of $N(0,1)$. If $D_n=\displaystyle\sum_{i=1}^n X_{i}^2$ and $R_n=\displaystyle\sum_{i=1}^n \frac{X_{2i-1}}{X_{2i}}$ and $T_n=\displaystyle\frac{R_n}{D_n}$, what is the limit distribution of $T_n$?


title: Odds and odds ratios in logistic regression

text: I am having difficulties understanding one logistic regression explanation. The logistic regression is between temperature and fish which die or do not die. The slope of a logistic regression is 1.76. Then the odds that fish die increase by a factor of exp(1.76) = 5.8. In other words, the odds that fish die increase by a factor of 5.8 for each change of 1 degree Celsius in temperature.   1. Because 50% fish die in 2012, a 1 degree Celsius increase on 2012 temperature would raise the fish die occurrence to 82%.   2. A 2 degree Celsius increase on 2012 temperature would raise the fish die occurrence to 97%.   3. A 3 degree Celsius increase -> 100% fish die. How do we calculate 1, 2 and 3? (82%, 97% and 100%)


title: Type of input data (proportions or log-ratios) for compositional analysis

text: I am comparing tree species use vs. availability within bird territories. I am using the R `compana()` function which uses Aebischer et al. (1993) method (package adehabitatHS) to determine if there are significant differences in proportional availability and proportional use of different tree species. I am trying to determine if the input to the program needs to be proportions (of both use and availability) or log-ratios of use and availability. Additionally, I am also trying to figure out whether or not the `compana` code for compositional analysis automatically replaces missing values in the available matrix by using the weighted mean $\lambda$ when you specify `test="randomization"`, which means that randomisation tests are performed for both the habitat ranking and the test of habitat selection (check out CRAN website for more information). I just can’t figure out if `compana` does the replacement of zeros in the available matrix automatically or if you need to do something else with the code.


title: Separating the populations in a bimodal distribution

text: I have a data set which displays a bimodal distribution. This was determined by plotting a histogram of the frequency vs number. I now need to separate the two original populations and therefore find an intersection point of sorts. From the plot it looks like the point might be approx. -1.0 to -0.8. Is there a straight forward calculation or function that I can use to locate this point more accurately? ![Histogram](http://i.stack.imgur.com/2TNd2.jpg)


title: Finding a trend in variance

text: I have written a little function that samples from a normal distribution with linearly changing variance               rnormltv <- function(n,mu,rsd) {         vc <- vector(mode="numeric", length=n)         i <- 1         for (x in seq(from=rsd[1], to=rsd[2], length=n)) {             vc[i]=rnorm(1,mu,x)             i <- i+1         }         return(vc)     }      Suppose we obtain 1000 samples, with a decreasing variance from 9 to 1.               set.seed(1)     x <- rnormltv(1000,0,c(9,1))     plot(x)      The goal is to estimate the parameters 9 and 1 given this data, the mean of zero, assuming a normal distribution and a linear trend in variance. Conceptually I can see how this can be done with maximum likelihood, but how can it by done in practice, preferably in R ? Note: This is not homework.


title: Significance of initial transition probabilites in a hidden markov model

text: What are the benefits of giving certain initial values to transition probabilities in a Hidden Markov Model? Eventually system will learn them, so what is the point of giving values other than random ones? Does underlying algorithm make a difference such as Baum–Welch? If I know the transition probabilities at the beginning very accurately, and my main purpose is to predict output probabilities from hidden state to observations, what would you advise me?


title: Should the risk of Type II error associated with a n.s. factor be discussed if this factor significantly interacts with another factor?

text: When computing a 2(A)$\times$2(B)$\times$3(C) ANOVA, I observe that only factor A has a significant main effect. Factors B and C were not significant, yet display a moderate effect size and low power. Yet the factors B and C significantly interact with A and with eath other (A$\times$B, A$\times$C and B$\times$C). **Question:** When it comes to writing up the results, it is relevant to indicate that factors B and C may hide a risk of Type II error (i.e. due to their low power yet moderate effect size)? I have a sense that this is not appropriate since these factors significantly interact with A.


title: Difference Between Simultaneous Equation Model and Structural Equation Model

text: Can anybody please help me to understand what are the differences between Simultaneous Equation Model and Structural Equation Model (SEM)? It will be great if somebody can provide me some literature on it. Also, is there any literature where SEM has been used in timeseries context? The literatures I'm getting are mostly explained SEM in cross-sectional data context. Thank You!


title: How to perform Student's t-test having only sample size, sample average and population average are known?

text: Student's $t$-test requires the sample standard deviation $s$. However, how do I compute for $s$ when only the sample size and sample average are known? For example, if sample size is $49$ and sample average is $112$, I will then attempt to create a list of $49$ identical samples with values of $112$ each. Expectedly, the sample standard deviation is $0$. This will create a divide- by-zero problem in the $t$ test. **ADDITIONAL DATA:**   The average income of ACME North Factory workers is $\$200$. It is reported that a random sample of $49$ workers in ACME South Factory had an annual income of $\$112$. Is this difference statistically significant? Am I correct in saying that the population mean is $\$200$?


title: Effect of sample decimal points on power in rank-based tests

text: I perform different rank-based non-parametric tests (Friedman and Nemenyi) post-hoc and noticed that the decimal places of the samples have a strong effect on their power. Please consider the following toy example:               Samples 1: 0.12,0.23,0.22,0.17     Samples 2: 0.15,0.17,0.18,0.1     Samples 3: 0.19,0.29,0.27,0.19      Using all decimal places might show that samples 1&2 are not significantly different but samples 1&3 and 2&3 are. Nevertheless, omitting the last decimal place will make samples 1 equal to samples 3. Thus, omitting too many decimal places might make the samples too similar (they will share a rank). On the contrary, not omitting them makes them prone to small irregularities in sample computation or handling of floats (every sample gets a different rank, although being very close to each other). Questions:   1. What is the best way to deal with this, especially in a publication?   2. Is the difference between per samples described as "effect size"? In this case, rounding to a specific amount of decimal points would serve as a constraint to the minimum effect size to be considered. Thanks!


title: Outliers in proportional hazards model

text: How can one easily found an outlier in PH model and why are they so important? Can anybody give me some advice or suggest which R function may I use? Thanks for comments and help.


title: Question about the Bayesian Inference of a parameter

text: In order to understand the difference between the Frequentist and Bayesian inference, I was reading the presentation at: http://www.stat.ufl.edu/archived/casella/Talks/BayesRefresher.pdf . In order to explain the difference between the approaches, the author is using the following example(Page 6): If I understand correctly, there are three different cases of protocols (something related to toxicity, it seems) and there are different studies for each of the three protocols, which has examined the relation of the protocol to a specific cause, which is called something as "AAN". For example, for one of the three protocols, different studies showing the total number of cases (first) and relation to "AAN" (second) is as following : 1) 66,11 2) 1756,129 3) 272, 48 4) 151, 18 ... etc. Each of these number pairs belong to different studies. Now, the Bayesian model for these studies are given as such: Given $X_j$ is the AAN frequency of the $j$th study about the protocol $i$, $X_j$ is distributed as $X_j \sim Binomial(n_j,p_i)$. $n_j$ is the total number of incidences for the $j$th study and $p_i$ is the parameter of the distribution we want to infer. What I did not understand here is, the author says that, in this model $p_i$ can vary from study to study. In my understanding of the Bayesian approach, here, all of the studies constitute our data, $D=\left( n_1,X_1,n_2,X_2,...,n_K,X_K \right)$ where $K$ is the total number of studies. For $p_i$ we have a prior distribution of $P(p_i)$. So we try to find the posterior distribution $P(p_i|D)$. In my understanding $p_i$ cannot vary from study to study: $p_i$ is first generated from the prior distribution as $p_i \sim P(p_i)$ and then this generated value of $p_i$ is used to generate each $X_j \sim Binomial(n_j,p_i)$. Yes, $p_i$ is not fixed as in the frequentist approach, but it varies over different realizations of all $K$ studies, once a $p_i$ is generated from the prior, then all $K$ studies use the same $p_i$. So, it should not change from a single study to another. So, this slide has confused me at this particular point. Am I right with my thoughts here or did I misunderstand something? Thanks in advance.


title: Plottin lavaan SEM/CFA in APA-style

text: I like doing SEM's with lavaan. But I dont like reporting them because I haven't found a good way to plot them in APA-style. I usually rebuild the model in PowerPoint and add the parameters by hand. What a pain! To look at the results of the SEM I usually use semPaths from the semPlot package, which is great, but it's output is not in APA-style (or anything close). How do you deal with that issue when presenting an SEM? Has anyone found a good solution for this?


title: What can I do to if I have a nominal dependent variable and scale independent variables?

text: Our dependent variable is a nominal variable that asks respondents "would you share this video?" with the responses being either yes or no. The hypotheses we are testing are all relational directional hypotheses and all our independent variables are measure using a Likert scale from 1-7, therefore this is scaled data. The problem we are facing however, is that when we try run the hypotheses tests in SPSS, the output says that the data could not be calculated because our nominal independent variable is constant (i.e. it is not a scale dependent variable). Is there a solution to this problem? Advice would be greatly appreciated.


title: Is there a formula to find the MAPE for $Y_t$ if we know the MAPE for $\Delta$log($Y_t)$?

text: Is there a simple formula to find the MAPE for $Y_t$ if we know the MAPE for $\Delta$log($Y_t)$ ~ iid N($\mu$,$\sigma^2$)? Is there an algebraic relation between the two? What if I use RMSFE instead? Or some other measure of forecast accuracy.


title: What is the difference between Kaplan-Meier Model and Cox PH model in this case

text: Say you have a case where you want to perform survival analysis on a set of individuals with a single binary-valued covariate. Assume that the PH assumption is satisfied. Say that the goal is to analyse the difference between the two groups (corresponding to the 2 values of the covariate). In this case, what would be the difference between using Kaplan-Meier and Cox PH model? Also, intuitively, why would they be different? I know that Cox takes into account the value of the covariate, but in this case, since it is binary- valued, Kaplan-Meier would also do this in a sense.


What does the theoretical crosscorrelation function look like for two unrelated random walks? I know that for two random walks, they will be spuriously correlated. But what about the autocorrelations?


title: Statistical demand forecasting

text: How is batch demand forecasting done in retail like in Walmart where number of products to forecast are very large in number and products are short lived i.e have less than 36 months of historical data? I know there statistical methods like ARIMA, exponential smoothing, but in retail ( for example walmart) there are 100,000+ products ranging in various categories like detergents,grocery etc each having its own seasonality, geographical hierarchy, product hierarchy. Also some of them have very irregular demand i.e zero in most of the months and very high in some months for example spare parts


title: Include nesting factor as fixed effect in a GLMM

text: I have the following GLMM:               success ~ age + gender + group/task + (1 + group/task|school/subject), family = binomial      I want to know whether participants' probability to succeed in certain problem-solving tasks can be predicted by the type of task. I have 6 tasks which can be categorized into 3 groups (A, B, C) with 2 tasks in each group. Each participant received 3 tasks (one from A, one from B, one from C; combination and order counterbalanced). Cochran's Q- and post-hoc McNemar tests revealed that the three groups differ in their success rates: A is easier than B and C, and B and C are equally difficult. I used crosstabs to analyze whether the tasks within each group differ in their success rates and found that they are different for A and B. Now I would like to do a comparison of all individual tasks (not just those within one category). My question is: Is the equation above correct in terms of the fixed effects or is there any reason to include group as an extra fixed effect (e.g. to see whether task has an effect on top of the group effects found in the McNemar tests)? Would that be unnecessary? What does it mean to include both group/task and group?


title: How to model a count dependent variable with upper limitation

text: I have a dependent variable, which have 0, 1, 2, or 3 for its value. I asked participants to choose three items and coded 1 if it is in a certain category and 0 otherwise. I add the three binary variables to be my dependent variable. So 3 is the maximum value of the variable. I know that Poisson regression is used a lot for the count dependent variable and right-censored model for the dependent variable with upper limit. But checking examples available in the Internet, it does not seem appropriate to analyze my data.


title: VECM "goodness of fit" in R

text: I'm using ca.jo in R to perform the Johansen test on a given dataset. I obtain my VECM coefficients, cointegration rank, etc. However, it does not seem to give any notion of "strength of relationship" or "goodness of fit". So my questions are: 1 - What is a standard/reliable way to measure how strong the relationship is in a given VECM cointegration model? 2 - How do I obtain that measurement in R?


title: Measured moderator, which is influenced by an independent variable

text: I have a measured variable, which is found to be influenced by an independent variable (i.e., experimental conditions). I will report this in the first part of the result section. Also, as an additional analysis, I want to use the measured variable as a moderator, which I assume influences a causal relationship between the independent variable and another dependent variable. This kind model is not quite common, I guess, but it sounds not problematic to me. I might have the multicollinearity problem (the two independent variable are not independent to each other) but it wouldn't be a severe problem unless the IV and moderator are closely related. Am I right?


title: lm weights and the standard error

text: I am fairly new to R, but am experienced in other languages and also in data analysis in the physical sciences. I have a problem and will illustrate with a straight line fit using `lm()` In the physical sciences, if the dependent variable (y) has explicit uncertainties (also called errors) and those uncertainties all have the same value, then if we fit to a straight line then the values of the slope and intercept are unchanged but the uncertainties in the slope and intercept will be different than if the data do not have explicit uncertainties. I can not find an R fitting routine that does the correct thing. I've written fitters that do the right thing in C, Mathematica, Maple, and LabView, but hope to avoid doing so with R. Here's an example:               d.f = data.frame( x = c(1, 2, 3, 4), y = c(2, 3.9, 6.1, 7.9), u = c(.1, .1, .1, .1))      The variable `u` is the uncertainties in the values of `y`. Do a straight line fit ignoring the uncertainties:               fit1 <- lm(y ~ x, data = d.f)      This gives exactly what I think is the correct fit, including the "Std. Error" in the slope and intercept. The weights are conventionally $1 / (\text{uncertainty}^2)$, so I give that to `lm()`               fit2 <- lm(y ~ x, data = d.f, weights = 1/u^2)      The Residual standard error is different for these 2 fits, as it should be, but the Std. Error in the fitted parameters is not. For example, the intercept for both fits is $1.99$, which is correct, and the Std. error in the intercept for fit1 is $0.05196$, which is also correct. However, these values are exactly the same for fit2. The correct value for the uncertainty in the intercept for fit2 is actually $0.04472$. I've searched around a fair amount, and have failed to find a least-square regression fitter in R that does this correctly. In fact, the ideal such fitter would also accept uncertainties in both the `x` and `y` coordinates, and fit using an "effective variance" technique. Or perhaps I just don't understand how to use `lm()` or interpret its results.


title: Appropriate test for detecting a signal in normally distributed noise

text: I am doing some signal processing and I have a histogram which has a bell shape when there is only noise in the signal. (I have been advised that this is to be expected due to the central limit theorem) ![Histogram](http://i.stack.imgur.com/xwd6R.png) However when the signal that I am looking for is present The shape of the histogram changes to ![With Signal](http://i.stack.imgur.com/AIm6l.png) Note that I have removed the bars in the middle because they are so high it makes it hard to see the detail at the base of the curve. I would like a test to discriminate between these two cases. I am aware that there are some tests for normality however I am unsure of how appropriate they are in this case as explained here. One approach that I have considered is testing the smoothness > $\int_{-\infty}^{\infty}{[f''(x)]^2 dx}$ where f is the normal distribution **EDIT** As Glen_b sugested, I tried the Anderson Darling Test               float[] sortedYs;     sortedYs = Hist.OrderBy(a => a).ToArray<float>();     float[] cdf = new float[sortedYs.Length];          MathNet.Numerics.Distributions.LaplaceDistribution ld = new MathNet.Numerics.Distributions.LaplaceDistribution(U, (b));     for (int i = 0; i < sortedYs.Length; i++)     {         cdf[i] = (float) ld.CumulativeDistribution(sortedYs[i]);     }          float AD = 0;     for (int i = 0; i < cdf.Length; i++)     {         AD -= (float)((2 * (i + 1) - 1) * (Math.Log(cdf[i]) + Math.Log(1 - cdf[cdf.Length - 1 - i])));     }     AD /= cdf.Length;     AD -= cdf.Length;     AD *= (float)(1 + (0.75 + 2.25 / cdf.Length) / cdf.Length);      (I'm using C# and Math.Net) I'm still not quite sure if this is how anderson darling works. (about 8000 - 14000 for bubbly, about 4000- 12000 for non bubbly). Its also probably worth noting that I have quite a large data set, since my images have 1920*1080 pixels. The AD value calculated is much higher for moving video frames. The smoothnes calculation is done by estimating $\int_{-\infty}^{\infty}{[f''(x)]^2 dx}$ using the raw data for f and then by doing the same estimate but using a laplace distribution (using sample mean and b) for f I find the difference of these two estimates. I am getting about 0.01 and less for no motion, and higher for moving images. However there is overlap in these categories as well so its not as reliable as I would like yet. **EDIT** I'm going to post some more histograms No Motion ![no motion](http://i.stack.imgur.com/s1oZZ.png) Motion Present ![motion](http://i.stack.imgur.com/LZ2fS.png) I was looking at the noise in the tails all this time, but if I step back and look at the whole histogram, it can be seen that images with motion have a conspicuously lower variance in their histogram than images without motion. I think this might be partially due to the way pixels are categorised into bins in a histogram. (I'm using opencv and emgu so I will have to look at the source code for more clues here). Images with motion actually have a wider spread of gray scale values than images with no motion. However when calculating the histogram, the pixels have to be placed in bins from 0 to 255. So an image with widely varying grayscale values has lower resolution bins. And so, more pixels get placed in the same bin, even if they a bit different. This does kind of ruin my data, but at the same time, its an effect that I can take advantage of too.


title: What is a semantic map

text: Based on: J. P. Carvalho, "On the Semantics and the Use of Fuzzy Cognitive Maps in Social Sciences" (WCCI, 2010 -- PDF) and Richard Dagan's web page Cognitive Mapping. A cognitive map consists of nodes representing concepts and ties indicating relationship between the concepts. In Carvalho section IV, the semantics is used in perspective of the presented model. Dagan's definition of semantic map is that it is a tree like structure branching from a main concept. Also, the answer http://cs.stackexchange.com/questions/9800/are-semantic- networks-ontologies?rq=1 further adds to my confusion as to what is the significant different between cognitive map where nodes are concepts and semantic map where also the nodes are concepts. I am unable to find significant difference between concept map and semantic map. Can somebody please explain in simple language? Thank you


title: Forecasting model inputs that are both auto-correlated and are calibrated over time?

text: How does one account for model inputs that are both a) auto-correlated and b) calibrated over time? I'm interested in forecasting the outcomes of sporting events. Let's say that each team has a score that is assigned each week of the season that indicates how good or bad that team is. That score is updated each week based on how the team performed in that week's game. Each team's score at the beginning of the season is a "best guess" based on the previous year and their upcoming schedule. I could easily build a historical model from all of the previous weekly scores, but something seems off about this. Scores that are assigned early in each season have greater uncertainty about them than scores assigned later in the season. Obviously, each week's score is correlated with (at least) the previous week's scores, but the scores become increasingly precise as the season progresses. Treating each score as an equal input in the forecasting model seems wrong. What would be the best modeling approach for dealing with this? Or am I overthinking it?


title: Comparing pre and post test

text: I would like to compare the pre and post test scores of students from years 1-4 (scores divided by question and by year). What test should I use?


title: Making sense of principal component analysis, eigenvectors & eigenvalues

text: In today's pattern recognition class my professor talked about PCA, eigenvectors & eigenvalues. I got the mathematics of it. If I'm asked to find eigenvalues etc. I'll do it correctly like a machine. But I didn't **understand** it. I didn't get the purpose of it. I didn't get the feel of it. I strongly believe in > you do not really understand something unless you can explain it to your > grandmother -- Albert Einstein Well, I can't explain these concepts to a layman or grandma.   1. Why PCA, eigenvectors & eigenvalues? What was the _need_ for these concepts?   2. How would you explain these to a layman?


title: How to explain PCA in layman's terms

text: I would like to know how can I explain PCA to someone not familiar with statistics. I have categorical/binary data and when I apply PCA the representation actually makes sense. So, if I put PCA in context what is the intuition behind it? Say, you have variables x,y,z and row 1= 0 0 1 row 2= 0 1 0 (i.e. for 1st observation x=0 y=0 z=1 etc.). Let's say 1 indicates that a function of a server is on and 0 is off. x is function A, y function B and z function C. Assume I have checked what is the status of the functions n times. So, what covariance matrix will give me? What is the meaning of eigenvectors and eigenvalues in this context?


title: How are eigenvectors and principal components related?

text: > **Possible Duplicate:**   >  Making sense of principal component analysis, eigenvectors & eigenvalues I am currently going through a PCA tutorial. However, I am a bit confused. For PCA, we calculate the covariance matrix and then find the eigen vectors and eigen values of the covariance matrix. These vectors are the principal components.   * How are the eigenvectors and the principal components related?   * Why does finding the covariance matrix from the data and its eigenvectors are the principal components?


title: Principal Component Analysis (PCA) for binary data

text: First of all, I would like to note that I have read similar topics in CrossValidated but I am not fully satisfied. I have a dataset which consists of an $N\times M$ binary matrix. 1 means that an action is performed and 0 that it is not. I apply PCA to the dataset and surprisingly get very good results, especially when I reduce it to only two dimensions. I am looking for the intuition behind performing PCA on such a dataset (i.e. where each attribute contains categorical data; you can give whatever example you think is more understandable) and whether a more appropriate technique can be applied. I am working with MATLAB and I need the data in a clustering friendly form.


title: What to do when the predictors do not accord with common sense/literature, but the model is fine/best according to log likelihood and LRT?

text: I would try to clarify the problem and then ask the questions. * * * **The problem** (variable names are masked due to confidentiality): I ran a binary logistic regression, in which there were 5 independent variables (IVs): A, B, C, D, and E. A and B are my concern. C is also my concern and I would talk about it later. They have two factors each (A1 and A2, B1 and B2). When I ran the estimation, some significant findings appeared. Nevertheless, the direction of the coefficients (beta values) for A and B was the opposite of what I expected! According to the most of the literature (not all of it), I expected A and B to have a positive beta, while they both had a negative coefficient. First I checked the models for about three full days. There was no mistake in them, and the directions did not change whatever changes I made to the models (except that in none of those changes, I attempted to drop the interactions). The log likelihoods as well showed that I am in a good direction. Then I decided to put my subjective view against the strange results aside and trust the results of the regression analysis. Then I passed to discuss those strange results and tried to justify the controversial findings. While discussing, I came to this point that "those two variables are heavily interconnected. Firstly, they had a significant interaction. Secondly, the distribution of the predictor A was heavily affected on the B, and according to the literature, A and B could have opposite effects; it could be important in my sample which was not balanced in terms of B. This imbalance could confound the effect of A as well. So I though maybe these are causing some problems. Then since C was as well strange, I thought maybe the whole model is being affected in a bad way by problems such as multicollinearity. I asked myself "what will happen if I isolate only A and B in the model?" If the interactions between A, B, and C are some sources of bias, can reducing the number of IVs lead to different results? The answer was yes: when I excluded all the other IVs from the model and left only A, B, and A*B, one of the coefficients became favorable and more in line with the literature and common sense. Thus I might tell that some errors do exist in my model which disrupt the main model (such as multicollinearity maybe). The I decided to examine every strange predictor, in isolation. When I excluded the interactions and left only the five IVs, the results seemed Much more consistent. Apparently, the problem begins when some specific interactions (but not all of them) are added to the model. After adding them, the directions of betas for A and C get reversed. It is a little annoying since by adding those specific interactions, the log likelihood reduces considerably (from about -75 to -48), so I cannot easily ignore those interactions. * * * **Questions** ( **the main ones are 3 and 4** , but an answer to the rest is as well much appreciated):   1. When the model acts strangely, but LRT and log likelihood tell that it is fine, which one should we choose? The subjective common sense, or the objective statistical measures?   2. Do you think is there a "problem" in my case, in the first place? Maybe everything is fine. If you wished, I can provide the raw data too.   3. What would you do in my case? At least three choices can be made: A. Dropping the interactions. B. Not dropping them and reporting the strange model. C. Reporting both the models with and without interactions, and also models of limited numbers of IVs (for example only A and B), and then try to subjectively discuss that "it is the interactions that cause the main large model strange".   4. I am going to do the latter (3.C). But that would be so messy and not so good looking. I wonder if there is an elegant, objective way of finding the source of error in the main model (well if there is any errors, of course), so instead of subjective discussions, I can substantiate my claims on some objective statistical measures. For example, is there a way to highlight the problematic interactions according to some statistics?    5. Do you have any other valuable suggestions or ideas?


title: Convergence of multivariate ECDF

text: Gilvenko-Cantelli assures uniform a.s. convergence of univariate ECDF. My questions are:   1. Are there similar assurances for multivariate ECDF?    2. How is the rate of convergence dependent on the dimensions of the r.v.? Thank you in advance.


title: Convergence Rate for Multivariate Kernel Density Estimators

text: I know that univariate kernel density functions converge uniformly a.s. to the true distribution. Is this also true for multivariate kdf? Is there a theorem that gives the rate of convergence in the multivariate case? What about the univariate case? Thank you in advance.


title: How to predict demand from historical "continuous" event data (date, lat, lon)?

text: I am attempting to predict demand for our service, both quantity but maybe more important, location (hotspots). I am by no means an experienced statistician, so I need some help :) I have all the historic data for our service, date, latitude and longitude. As far as I understand, the first thing to do is not to deal with the latitude and longitude: somehow they need to be converted into a single dimension right? After that what type of analysis should be done to the data? I think dealing with the date directly might also be the wrong way to go. My idea here is to deal only with week days, so I can predict the demand for a type of day (any Tuesday) instead of a specific Tuesday. I am looking for some guidance as to how to achieve this. I am a good programmer, but I do need some help finding the right way.


title: Why don't statisticians use mutual information as a measure of association?

text: I've seen a couple talks by non-statisticians where they seem to reinvent correlation measures using mutual information rather than regression (or equivalent/closely-related statistical tests). I take it there's a good reason statisticians don't take this approach. My layman's understanding is that estimators of entropy / mutual information tend to be problematic and unstable. I assume power is also problematic as a result: they try to get around this by claiming that they're not using a parametric testing framework. Usually this kind of work doesn't bother with power calculations, or even confidence/credible intervals. But to take a devil's advocate position, is slow convergence that big of a deal when datasets are extremely large? Also, sometimes these methods seem to "work" in the sense that the associations are validated by follow-up studies. What's the best critique against using mutual information as a measure of association and why isn't it widely used in statistical practice? edit: Also, are there any good papers that cover these issues?


title: Is analysis of existing data always exploratory, or can it be used for hypothesis testing?

text: This is a question about the rhetoric of describing analysis done using a public data set or any other pre-existing data set. **Here is the hypothetical situation:** A researcher reports that they have a hypothesis. To test this, they take a sample, n, of individuals that meets study criteria from a database of N participants. They run a single, a priori, test of this hypothesis on the sample, and report that the hypothesized group difference, or correlation, or whatever, is greater than 0, p<.05. _Is this exploratory or is it a fine example of hypothesis testing?_ Because any single reported analysis of an existing data set might be one of many interrogations, any such report should be framed as exploratory rather than hypothesis testing. In other words, presenting an analysis of an existing data set using the rhetoric of a single sample hypothesis test seems to be misleading the reader to over-interpret the results. However, I can also see an argument being made for the researcher being given the benefit of the doubt as is similarly the case in all non-pre-registered studies. In other words, can we assume that the researcher did just test this one hypothesis on this single subset of the existing database?


title: Bound for Arithmetic Harmonic mean inequality for matrices?

text: **NOTE:** _This question has originallybeen posted in MSE, but it did not generate any interest. It was first posted there, because the question itself is a pure matrix-algebra question.   Nevertheless, since the motive has to do with statistics and econometrics, I am posting the question on Cross Validated also, in the hope that some statistics/matrix algebra savvy brain will have something to contribute._ The framework is as follows: We have a cross-sectional i.i.d. sample $\\{\mathbf y, \mathbf X\\}$, where $\mathbf y$ is a $N \times 1$ column vector, and $\mathbf X$ is a $N\times K$ matrix. We postulate a linear relationship between $\mathbf y$ and $\mathbf X$, $$\mathbf y = \mathbf X \beta + \mathbf u $$ where $\mathbf u$ is white-noise with variance $\sigma^2$, and exogenous to the regressors in the $\mathbf X$ matrix, and $\beta$ is a $K\times 1$ column vector of unknown constant coefficients. Under this assumption, the OLS estimator is unbiased and consistent. Now assume $N$ is "large", say $10.000$ or more (samples that large have started to appear in the econometrics field also). Then a researcher could conceivably entertain the following two options: $A$) Run _one_ OLS regression using the whole sample. This tactic can be thought of as appealing to the consistency property of the OLS estimator. Call this estimator $\hat \beta$. $B$) Divide the sample into $m$ _disjoint_ sub-samples (for simplicity, assumed of equal length, and note that their union equals the whole sample), run $m$ regressions, and calculate the average of the $m$ coefficient estimates she will thus obtain. This tactic can be thought of as appealing to the unbiasedness property of the OLS estimator. Call this averaging estimator $\bar b_m$. _(Note that tactic $B$ does not fall into any re-sampling approach, like bootstrap, subsampling, or jackknife -to be exact, it has been considered as a marginal case in applying jackknife in time series, but it is not trully a jackknife method)._ I have derived a nice-looking (to me) result that shows that the variance of the whole-sample estimator is always smaller than the variance of the averaging estimator: $$\text{Var}\left(\bar b_m\right) > \text{Var}\left(\hat \beta\right) $$ I say it is nice-looking because the result uses the Arithmetic-Harmonic mean inequality for PD matrices, proven here: Specifically, writing $Z_l= \left(X_l'X_l\right)^{-1}$ for the inverse moment matrix of the regressors from the _l_ -th sample, _l_ =$1,...,m$, denoting by $A_m$ the arithmetic mean and by $H_m$ the harmonic mean of these $Z$ matrices, it is not hard to arrive at the following: $$\text{Var}\left(\bar b_m\right) = \frac1m\sigma^2A_m > \frac1m\sigma^2H_m =\text{Var}\left(\hat \beta\right) $$ ...the inequality holding in the matrix sense. Note that $H_m$ is the harmonic mean of $\left(X_1'X_1\right)^{-1},...,\left(X_m'X_m\right)^{-1} $in the true matrix sense, not a matrix containing the harmonic means of the corresponding elements of the matrices it averages. So the averaging estimator $\bar b_m$ is always less efficient than the whole sample estimator $\hat \beta$. **My question: Are there any known bounds for the difference between Arithmetic-Harmonic means for matrices?** For real numbers they are (see the wiki article and the original resources here and here ). Why? because it will be helpful in order to move to the next step and compare estimators that may be neither unbiased nor consistent, and so one is left with a criterion such as Minimum Squared Error to compare them. Any suggestion, link or reference will be really appreciated.


title: Intuition on the definition of the covariance

text: I was trying to understand the Covariance of two random variables better and understand how the first person that thought of it, arrived at the definition that is routinely used in statistics. I went to wikipedia to understand it better. From the article, it seems that good candidate measure or quantity for $Cov(X,Y)$ should have the following properties:   1. It shoukd have a positive sign when two random variables are similar (i.e. when one increases the other one does to and when one decreases the other one does too).   2. We also want it to have a negative sign when two random variables are oppositely similar (i.e. when one increases the other random variable tends to decrease)   3. Lastly, we want the this covariance quantity to be zero (or extremely small probably?) when the two variables are independent of each other (i.e. they don't co-vary with respect to each other). From the above properties, we want to define $Cov(X,Y)$. My first question is, it is not entirely obvious to me why $Cov(X,Y) = E[(X-E[X])(Y-E[Y])]$ satisfies those properties. From the properties we have, I would have expected more of a "derivative"-like equation to be the ideal candidate. For example, something more like, "if the change in X positive, then the change in Y should also be positive". Also, why is taking the difference from the mean the "correct" thing to do? A more tangential, but still interesting question, is there a different definition that could have satisfied those properties and still would have been meaningful and useful? I am asking this because it seems no one is questioning why we are using this definition in the first place (it kind of feels like, its "always been this way", which in my opinion, is a terrible reason and it hinders scientific and mathematical curiosity and thinking). Is the accepted definition the "best" definition that we could have? * * * These are my thoughts on why the accepted definition makes sense (its only going to be an intuitive argument): Let $\Delta_X$ be some difference of for variable X (i.e. it changed from some value to some other value at some time). Similarly for define $\Delta_Y$. For one instance in time, we can calculate if they are related or not by doing: $$sign(\Delta_X \cdot \Delta_Y)$$ This is somewhat nice! For one instance in time, it satisfies the properties we want. If they both increase together, then most of the time, the above quantity should be be positive (and similarly when they are oppositely similar, it will be negative, because the $Delta$'s will have opposite signs). But that only gives us the quantity we want for one instance in time, and since they are r.v. we might overfit if we decide to base the relationship of two variables based on only 1 observation. Then why not take the expectation of this to see the "average" product of differences. $$sign(E[\Delta_X \cdot \Delta_Y])$$ Which should capture on average what the average relationship is as defined above! But the only problem this explanation has is, what do we measure this difference from? Which seems to be addressed by measuring this difference from the mean (which for some reason is the correct thing to do). **I guess the main issue I have with the definition is taking the difference form the mean**. I can't seem to justify that to myself yet. * * * The interpretation for the sign can be left for a different question, since it seems to be a more complicated topic.


title: how to do category regression

text: I have a set of data which contains students first term gpa, final gpa, their age, race, state and so on, I need to predict student's final gpa based on the others, my problem is how to convert race and state into numbers, I can't use dummy variables because there is 8 race and 20 states


title: What kind of model is used by 20 Questions?

text: Which kind of machine learning concept / model is used in 20 Questions?   Is this kind of thing best solved by a neural network?   Where I can read something about it?


title: First-differences vs. pooled OLS

text: Two-period panel regression model:   * $Y_{it}=\beta_0\delta_{t2} + \beta_1X_{it} +\alpha_i + \epsilon_{it}$ where   * $t=1, 2$    * $i=1,...,n$ indexes the cross sectional units drawn via random sampling.   * $\alpha_i$ is an individual-specific unobserved effect.   * $\delta_i=1$ if $t=2$; $0$ if $t=1$ In other words:   * $Y_{i1}=\beta_1X_{i1} +\alpha_i + \epsilon_{i1}$   * $Y_{i2}=\beta_0 + \beta_1X_{i2} +\alpha_i + \epsilon_{i2}$ Now for the first-difference estimator $\hat{\beta_1}^{FD}$ to be consistent, one must assume strict exogeneity :   * $E[\epsilon_{it}|X_{i2}, X_{i1}]=0$ for $t=1, 2$ For the pooled OLS estimator $\hat{\beta_1}^{POLS}$ to be consistent, one must assume weak exogeneity:   * $E[u_{it}|X_{it}]=0$ where $u_{it}=\alpha_i + \epsilon_{it}$ I know that $\hat{\beta_1}^{POLS}$ will be less efficient than $\hat{\beta_1}^{FD}$ because $u_{it}=\alpha_i + \epsilon_{it}$ is serially correlated, but **is it necessarily the case that first-difference estimators will always be superior to pooled OLS estimators?** Surely the assumption of strict exogeneity in the first-differences case rules out the possibility that $X$ might have a persistent effect on $Y$ thereby biasing $\hat{\beta_1}^{FD}$? In other words, although the weak exogeneity assumption makes $\hat{\beta_1}^{POLS}$ inefficient might it be the case that the strict exogeneity assumption biases $\hat{\beta_1}^{FD}$? **There is a trade-off between strict and weak exogeneity.**


title: How do I do a log-log regression in SPSS?

text: I am unable to find a way to do a log-log regression in SPSS? Could someone help me do this step by step? Thanks.


title: Margin of error for non-dichotomous data

text: I am trying to calculate the margin of error for specific survey questions. I know how to calculate this if the data are dichotomous (i.e., a 75-25 or 60-40 split). But how do you deal with questions that have more than one response? For example: > Which of the following age ranges best describes you?   >  () Younger than 18   >  () 18 to 30   >  () 30 to 50   >  () 50 to 70   >  () Over 70   * Is the margin of error calculated for this question by treating each response as it were a separate question? (That is, if 20 respondents out of 200 picked "Over 70", would I treat "Over 70" as a separate question that has a 90-10 split?)    * If not, how would margin or error be calculated?


title: Are very low ANOVA and Kruskal Wallis p-values reason to doubt the results?

text: I have four relatively small sets of data under five different conditions. For each set of data I want to see whether the differences between the conditions are statistically significant. The first two sets are by visual inspection (histograms) the least normally distributed, whereas the two latter sets look almost normally distributed, with a few peaks outside what would be the Gaussian function. All data sets fail Kolmogorov-Smirnov and the t-test, but succeed the Chi-square goodness of fit test. The quantile-quantile plot for the first and last set is plotted below, which also suggests almost normal distribution for the better sets (is linear enough, or must it be 45 degrees?): set1: ![qqplot1](http://i.stack.imgur.com/lm4Xi.png) set 4: ![qqplot2](http://i.stack.imgur.com/3wjD8.png) The p-values for an ANOVA and Kruskal-Wallis for the four data sets are:                       ANOVA        Kruskal Wallis     set 1   3.21e-07     5.31-05     set 2   1.23e-06     1.04e-04     set 3   3.39e-11     1.18e-16     set 4   1.31e-23     7.20e-17      The outcome of those tests indicate what I found by visual inspection, that the sets are more and more normally distributed from set 1 to set 4. However, in other research and examples I generally see p-values in e.g. ~0.00X. These values are tiny! It seems counter-intuitive given that the data sets failed tests for normality. Is that a reason for concern? Given these results, am I justified in running the ANOVA and plotting confidence intervals, regardless of the fact that the t-test and Kolmogorov-Smirnov test failed? I know that the ANOVA is relatively robust towards non-normal distributions, but is my data too far off? Edit: Results of t-tests for first and last set as requested:               t-test:     set 1     condition    1         2         3         4         5     H            1         1         1         1         1      P (1e-09 *)  0.0000    0.0007    0.5115    0.0000    0.0007          set 4     condition    1         2         3         4         5     H            1         1         1         1         1     P (1e-11 *)  0.0000    0.0000    0.0000    0.1729    0.0000


title: Generating values from copula using copula package in R

text: I have a bunch of questions concerning the use of the copula package in R. My overall aim is to generate synthetic values using copulas. I am analyzing a hydrological data: annual peak discharge [m³/s] and corresponding volume [m³]. I managed to apply tests on serial independence and dependence. Furthermore I identified and excluded ties and created pseudo-observations (transformation of copula values between [0,1]). Since I don’t know which copula is the best, I fitted the copula parameter first:               fg <- fitCopula(copula=gumbelCopula(), data=u) #  u is my data      I will do this for all available copulas in R. Afterwards I test the goodness of fit with the following function:               gofCopula(copula=gumbelCopula(fg@estimate), x=data[,2:3], N=1000, method="Sn",                estim.method="mpl", simulation="mult")      Using the “best” copula, I then want to create synthetic values. I found a function to create random samples, but I am not sure, if it does what I need.               random_samples <- rCopula(copula=gumbelCopula(fg@estimate), n=10000)      It seems to me that this function creates only random values, but is the dependency structure of my data set considered? There is also another function in the copula package `mvdc`, for the construction of multivariate distributions from copulas. What is actually the difference of `mvdc` and `rCopula`, both are generating synthetic values, aren’t they? One last question is: Once I am able to generate my synthetic values, how can I transform them back to their real units? From reading through the documentation I understood tat I have to multiply the values for (u,v) with the inverse of their particular cdf, is this true?


title: How to formulate a likelihood ratio test for an exponential distribution?

text: > Find the form of the likelihood ratio test of $H_0 : \lambda = \lambda_0$ > against $\lambda\neq\lambda_0$ when $X_1,X_2,...,X_n$ is a random sample > from $Ex(\lambda)$. Simplify it as much as possible. I've started as follows, where $\lambda_1$ is some alternative parameter $$ \frac{L_0}{L_1}=\frac{\prod f(x;\lambda_0)}{\prod f(x;\lambda_1)} $$ I think the p.d.f. for an exponential is $\lambda e^{-\lambda x}$ so $$ \frac{L_0}{L_1}=\frac{\prod \lambda_0 e^{-\lambda_0 x_i}}{\prod \lambda_1 e^{-\lambda_1 x_i}} $$ $$ =(\frac{\lambda_0}{\lambda_1})^n\frac{\prod e^{-\lambda_0 x_i}}{\prod e^{-\lambda_1 x_i}} $$ Which I think can be simplified like this by taking logs $$ =n log(\frac{\lambda_0}{\lambda_1})\frac{e^{-\sum\lambda_0 x_i}}{e^{-\sum\lambda_1 x_i}} $$ Not sure about this but, I think that becomes $$ =n log(\frac{\lambda_0}{\lambda_1}) e^{-\sum x_i ({\lambda_0}-{\lambda_1})} $$ Am I right so far? Not sure where to go from here. I think the critical region is where the above is less than some constant.


title: Setting up feature vectors

text: I am working on a classification project and I want to use SVM's and/or Clustering Algs. What I am having trouble with is deciding how to set up my feature vectors. I have already decided what my features will be and I have already taken into account things like correlation and overfitting. Some of my features are categorical (region, gender, occupation...) and some are numeric (age, income,...) I would like to use some if not all of these. But I am not sure how to quantify the features in such a way that is accurate and will not create noise. (One thing I thought to do was for each categorical variable, say gender, calculate something relevant that has to do with the variable (say probability of default given that you are of some gender). The idea is to find something highly correlated with the categorical variable and use that value instead.) Id appreciate any tips or references! Thanks


title: Probit regression in R giving singular Hessian matrix

text: I am trying to run a probit regression using panel data in R by first computing the log likelihood and then using the `optim` function to optimize.   1. Scale of predictor variables: The predictor variables that I have vary in scale significantly. Some of the variables are in tens while some are in '000s. When I am trying to run the `optim` function, it frequently gives the error message: `vmmin is not finite`. So I have to be extremely careful in choosing the initial values. And some of these initial values are as small as 0.0001. Is this standard in probit? Or do we need to normalize the predictor variables before using probit?   2. I have a few variables that are computed within the estimation code, i.e., some of the parameters that are being estimated are used to calculate these inbuilt covariates. When I run the full model with all such inbuilt parameters, the `optim` function gives me a singular Hessian matrix. Again, I am not sure what the singularity of the hessian implies. When I the add the inbuilt covariates one at a time, the Hessian is not singular.  Has anyone encountered such problems before? Please advise.


title: MPE and MAPE averages

text: I am a computer science graduate with a minor in mathematics, however it has been a long while since taking any statistical courses. I am looking at calculating the Mean Percent Error and Mean Absolute Percent Error between expected scores and actual scores. Lets use the school setting for a familiar analogy. If I calculate the MPE and MAPE for each school as a whole based on various test scores against expected test scores, could I then use those MPE and MAPE results (summed and averaged) to get the MPE and MAPE for larger subsets such as combining multiple School's MPE/MAPE results into a district level MPE/MAPE, or would that skew the results from a loss of proper weighting of data? Is there another piece of data I could retain alongside the school MPE/MAPE to address combining them into larger subsets, such as the number of students, or number of tests taken to maintain proper weighting?


title: Principal components using correlation matrix in R

text: My understanding is that `prcomp` and `princomp` work off the dataset itself (row of observations, across variables in the columns). Is there a function that will run a principal component analysis directly off a correlation or covariance matrix, without having the "raw" dataset?


title: Using Leave-One-Out Cross Validation with LARS

text: I have a kind of obscure question about using the Least Angle Regression (LARS) algorithm for variable selection. If I'm understanding it right, my professor formulates LARS as such: $$\mathbb{min}\ \hat{y}^T\hat{y}\ \mathbb{such}\ \mathbb{that}\ |t_j| \le t\ \forall\ j$$ where $t_j\ = (x_j^Tx_j)^{-1/2}x_j^T(y-\hat{y})/\sigma$. The dataset I'm working with deals with grades for a small class ($\mathbb{n = 22}$ students, $\mathbb{p = 16}$ graded assignments) and the response is the students' grade on the final exam. I stumbled across a result that's either mildly interesting or a coincidence, and was wondering if anybody knows about something mathematically similar to this that's related to variable selection? I standardized the predictors and centered the response and ran a Leave-One- Out Cross Validation for each of the $\mathbb{p}$ steps in the LARS algorithm. I wanted to figure out a way to tell which step I should stop at. One thing I thought of was to create a $\hat{y_i}^T\hat{y_i}$ vector with $\mathbb{n}$ elements where $\hat{y_i}$ is the predicted model with the $\mathbb{i^{th}}$ observation missing. $\hat{y}^T\hat{y}$ is the value of the full predicted model at each step. From there I made $\mathbb{p}$ "z-scores" calculated so that: $$z_j= \dfrac{\hat{y}^T\hat{y}-mean(\hat{y_i}^T\hat{y_i})} {sd(\hat{y_i}^T\hat{y_i})}$$ What I found was that $\mathbb{z_j}$ decreased until it reached a minimum at $\mathbb{z_8}$ and then increased. It so happened that the model with 8 variables was the exact same model we reached when we decided to only include positive $\mathbb{\beta_j} \mathbb{'s}$ ($\beta_9$ was negative) since that made realistic sense. I also ran this data through normal stepwise regression and got something like a minimum AIC at 3 variables and minimum BIC at 13 variables. I'm trying to make my code for this technically correct and run this with different data, but my long winded question is essentially if anyone thinks that this something mathematically related to an idea from regression/variable selection or just a strange uninteresting coincidence?


title: Meaning of latent features?

text: I'm trying to understand matrix factorization models for recommender systems and I always read 'latent features', but what does that mean? I know what a feature means for a training dataset but I'm not able to understand the idea of latent features. Every paper on the topic I can find is just too shallow. Edit: if you at least can point me to some papers that explain the idea.


title: Calculating posterior probabilities in R?

text: I asked a question earlier on here. Essentially, I am trying to evaluate a warning system that consists of a light bulb (of a specific color) being switched on, to indicate a predicted warning threat level. Srikant came up with a simple solution which involved calculating the posterior probabilities of the warning system, using bayes theorem. I am now contemplating implementing this in R. However, I am not quite sure how to proceed, since I am an R newbie. The data looks like this:               warning light, event level occured?     red, yes     none, yes     green, no          etc ...      I would like to write a simple R script that will help me to calculate: a). the posterior probablity of an event occuring, for a given state of the light bulb b). a C.I to attach to the posterior probability obtained in (a) above Since I am new to R, I would be grateful for the steps (and commands) required to do the above. My data will be in a simple csv file in the format described, so I can simply scan() it into R. How to proceed from there, is what I need help with..


title: Is it possible to do time-series clustering based on curve shape?

text: I have sales data for a series of outlets, and want to categorise them based on the shape of their curves over time. The data looks roughly like this (but obviously isn't random, and has some missing data):               n.quarters <- 100     n.stores <- 20     if (exists("test.data")){       rm(test.data)     }     for (i in 1:n.stores){       interval <- runif(1, 1, 200)       new.df <- data.frame(                       var0 = interval + c(0, cumsum(runif(49, -5, 5))),         date = seq.Date(as.Date("1990-03-30"), by="3 month", length.out=n.quarters),         store = rep(paste("Store", i, sep=""), n.quarters))       if (exists("test.data")){         test.data <- rbind(test.data, new.df)           } else {         test.data <- new.df       }     }     test.data$store <- factor(test.data$store)      I would like to know how I can cluster based on the _shape_ of the curves in R. I had considered the following approach:   1. Create a new column by linearly transforming each store's var0 to a value between 0.0 and 1.0 for the entire time series.   2. Cluster these transformed curves using the `kml` package in R. I have two questions:   1. Is this a reasonable exploratory approach?   2. How can I transform my data into the longitudinal data format that `kml` will understand? Any R snippets would be much appreciated!


title: Ordering of time series for machine learning

text: After reading one of the "Research tips" of R.J. Hyndman about cross- validation and time series, I came back to an old question of mine that I'll try to formulate here. The idea is that in classification or regression problems, the ordering of the data is not important, and hence _k_ -fold cross-validation can be used. On the other hand, in time series, the ordering of the data is obviously of a great importance. However, when using a machine learning model to forecast time series, a common strategy is to reshape the series {y_1, ..., y_T} into a set of "input-output vectors" which, for a time t, have the form (y_{t-n+1}, ..., y_{t-1}, y_{t}; y_{t+1}), sorry for the latex-sy notation. Now, once this reshaping has been done, can we consider that the resulting set of "input-output vectors" do not need to be ordered? If we use, for example, a feed-forward neural network with n inputs to "learn" these data, we would arrive at the same results no matter the order in which we show the vectors to the model. And therefore, could we use k-fold cross-validation the standard way, without the need to re-fit the model each time? Sorry if this seems as a silly question or if it is already answered in the site (I haven't been able to find it). Thanks!


title: Generalized Linear Model estimating multiple parameters concurrently

text: Is there an extension of the generalized linear model that is able to estimate (using maximum likelihood) more than one parameter of a distribution at once. For example, I know that for a gamma distribution with shape parameter $k$ and scale parameter $\theta$, a GLM can estimate the coefficient matrix $\beta$ for predicting $\theta$ from given covariates, responses, and $k$. Is there a way to estimate both parameters as a linear combination of predictor variables ? In chapter 10 of Generalized Linear Models (McCullagh & Nelder, 1989), Chapter 10 deals with joint modelling of mean and dispersion from covariates. However, they do not go into detail about a general algorithm to do so.


title: Modelling a variable that is both categorical and continious

text: I am attempting to design a study in which the impact of a variable on patient decision making is evaluated having adjusted for some other clinical variables. To give more detail, a clinical continuous variable (call it X) that is known to be associated with disease recurrence is used along some other clinical variables. Suppose, I have 100 patients and over the course of the study, the patients are randomized in such a way that half of them receive the variable X along other variables and choose to undergo treatment or not. While, the other half do not receive the info on variable X and based on other variables revealed to them choose to undergo treatment or not. Now, I'm interested to evaluate the impact of X on the decision (treatment vs no treatment), but as a continuous variable and the problem is half the patients did not receive the variable X. So, what is the appropriate approach here? Can I use a logistic regression model where I adjust for 1)other clinical variables, 2)whether someone received X or not (binary), and then also 3) X as continuous variable (for those who did not receive X, this value would be 0). I would appreciate your thoughts on this. thanks


title: Testing significance of correlation between two autocorrelated series

text: Say I collected the shin bones of _N_ different skeletons; they are all around 30cm long, and I measured different properties P1, P2, P3, P4 and P5 along these bones every 3mm (so I have 100 data points for each P) to see how well they measure bone density. I also measured the amount _F_ of fungus along these bones (every 3mm as well). I am especially interested in P5, which seems a pretty good indicator of bone density, and I want to see how well P1, P2, P3 and P4 are correlated with P5. I would like to have 1) a global measure of correlation between each P and P5 taking into account the measures along the whole bone (with significance test). Since the bone density varies in a particular way (non-linear, non cyclic, ...), and the amount of fungus can influence the other P, I would also be very interested in calculating 2) the local correlation (say along 3cm segment) between each P and P5 along the bone. The problem for 1) is that all the series are autocorrelated (since all the properties vary slowly in space). I would have done a simple cross- correlation, but then the p-value would need correction for the fact that the series are autocorrelated. Then what can I do? And should I average the values of all skeletons at each point along the bone before computing the correlation? For 2) I was thinking of calculating a "running/moving correlation" on a window, say 10-points long (3cm) so as to have a measure of correlation at each point. Then I would not only have to take into account the autocorrelation; I would also have to correct for the multiple testing caused by a "running/moving correlation". Any idea on how I could do this (Bonferroni corrections would be too conservative)? One of the things I would like to show for the running/moving correlation (but not only) is that low values of correlation between P1, P2, ..., P4 and P5 are obtained when the fungus measure _F_ is high (just saying). I use R, but I am a beginner in statistics, so I would be grateful if you could take that into account in your answer :) Thanks. P.S. I am not studying bones, but that seems a simpler story.


title: Method to compare ratings from multiple different sources with missing data

text: I want a method to compare ratings from multiple sources and find a single measure that best reflects all the ratings. To give a specific example, let's call it "The fellowship review committee problem" (but please don't hate me, I am not part of a fellowship review committee, this is just an example): A college fellowship review committee tries to decide which of their students best deserves a fellowship. It reviews 1000 student applications and wishes to focus on the grades of each student in up to 5 specific courses (let's say the courses are A, B, C, D, E). Not all students are required to have taken all 5 courses, some have taken 3, some 4, and some 5 (but all applicants have taken at least 3 out of A, B, C, D, E). What is the best single number that measures which student is better? My initial idea to solve this problem was to use principal components analysis (PCA). If all applicants had taken all 5 courses, that would work well. Now, I see in this post that PCA does not work well with missing data. I have two questions:   1. If there is no missing data (all students have taken all courses), is there a better way than PCA to find a single measure?    2. If I do have missing data, is there a better solution than the so-called DINEOF procedure described in the link I posted above?


title: How to report most important predictors using glmnet?

text: I want to find the most important predictors for a binomial dependent variable out of a set of more than 43,000 independent variables (These form the columns of my input dataset). The number of observations is more than 45,000 (these form the rows of my input dataset). Most of the independent variables are unigrams, bigrams and trigrams of words, so there is high degree of collinearity among them. There is a lot of sparsity in my dataset as well. I am using the logistic regression from the glmnet package, which works for the kind of dataset I have. Here is some code:               library('glmnet')     data <- read.csv('datafile.csv', header=T)     mat = as.matrix(data)     X = mat[,1:ncol(mat)-1]      y = mat[,ncol(mat)]     fit <- cv.glmnet(X,y, family="binomial", type.measure = "class")     betacoeff = as.matrix(fit$glmnet.fit$beta[,ncol(fit$glmnet.fit$beta)])      betacoeff returns the betas for all the independent variables. I am thinking of showing the predictors corresponding to the top 50 betas as the most important predictors. My questions are:   1. glmnet picks one good predictor out of a bunch of highly correlated good predictors. So I am not sure how much I can rely on the betas returned by the above model run.   2. Should I manually sample the data (say 10 times) and each time run the above model, get the list of predictors with the top betas and then find those which are present in all 10 repetitions? Is there any standard way of doing this? What is the standard way of sampling in this case?   3. My other question is about cvm (cross validation error) returned by the above model. Since I use type.measure = "class", cvm gives the misclassification error for different values of lambda. How do I report the misclassification error for the entire model? Is it the cvm corresponding to lambda.min?


title: Rating how closely one graph models another

text: I'm using a genetic algorithm to generate a string that produces certain results I map into a line/bar graph. I'm trying to rate how closely the results produced by the genetic string compares to a graph I have produce manually by hand. Currently I am comparing the graphs be finding the area of difference between the two and getting a percentage representing where it sits between 0 difference and a worst case difference being where they don't overlap at all. Though I don't feel this rating really represents how closely the shape matches the other. As this is what concerns me the most, as two different results could have the same difference in terms of area but one may be a closer shape then the other, and I would like the one with the closer shape to have a higher rating! My question is, is there a standard way of rating graphs where i can get a percentage of how closely one matches the other and can you provide the way of going about this?


title: Can I incude an effect size as an independent variable in a meta-regression?

text: My question is whether I can use an effect size $X$ as a dependent variable and another effect size $Y$ as the independent variable in a meta-regression. For example, I conducted a meta-analysis for the effects of exercise in drinking problems and I found significant results and high heterogeneity. I want to do a meta-regression and to use the effect size of those interventions in anxiety as an independent variable and the effect size of drinking problems as a dependent variable (assuming that each study assessed both anxiety and drinking problems and I calculated the effect sizes as Hedges's $g$). Does this make sense to you?


title: MCMC model from CDF values?

text: I am looking for a working example of an MCMC Bayesian model (for JAGS, etc.) when random variates, from a given distribution with unknown parameters, are observed only via their CDF values. For instance, given x ~ Gamma(shape, scale), we observe only {x, CDF(x)} pairs. Could anyone point me to an example of this, with pseudocode? I have found nothing to date except some references to quantile regression which is not the sort of model I'm after. I am looking for some valid way to use the CDF in the likelihood instead of the PDF. Clarification: This is not a case where, e.g., two observed CDF values completely determine the two unknown parameters. It is assumed that the observed CDF values are uncertain and the parameters of the model unknown. If possible, I would like the set of observed CDF values to find the posterior for the model parameters in the same way that a set of random variates finds it when the likelihood is a product of PDF values. I thought that there might be a way to use the CDF in the likelihood factor instead of the PDF.


title: How do I identify the "Long Tail" portion of my distribution?

text: I have a number of series that would typically be described as normal skewed or Gamma distributed. For example, say I have a group of customers and have calculated their spend over a fixed length of time. I then create a histogram to see distribution of spend and find an extremely long tail for the small group of high-spenders. My question is, since I want to identify these high spenders, are there methods to empirically inspect a distribution of values and approximate the point at which the distribution becomes "long-tailed" to create a cut in the data? I am not looking at inspecting a histogram to find the long tail. Just looking for a consistent method to systematically cut the data.


title: How to process triaxial accelerometer data for finger taps in real time?

text: I have a triaxial accelerometer mounted on my index finger and I want to differentiate in real time when it taps on a touch screen vs when the middle finger taps. For this I take the accelerometer data in a time window when the tap is detected on the touch screen. The problem is that the data is not clearly distinguishing the two cases, because even when the middle finger does not tap, it is following the downward motion when the index finger goes to tap. There are also dependencies on whether the 1st and 2nd finger go alternately or same finger again and again. I am new to this and am unable to figure out how do I go about analysing. First of all, how exactly should I mount the accelerometer so that it gives me the most clear data? - for instance if I mount it loosely, it might jiggle more when tapped but it might give increase noise as well. Secondly, do I need to try and integrate the data to distance values, or get it into Fourier transform? Thirdly, if I need to use machine learning (decision tree/svm etc.) which features will help me the most? There are just too many things to try out and I can't code every possibility in time - some guidance will help.


title: How to get true odds from bookmaker's prices? Removing overround

text: I've developed a model for predicting the probability that each horse will win a race. The output of the model is the predicted probability that each horse will win, the sum of all the probabilities will be 1. If I were to look at bookmaker prices for a horse race they would look something like this: `3.9, 5.4, 3.95, 6.7, 9, 14` If we took the sum of the probabilities that each horse would win, we'd get something like this: `1/3.9 + 1/5.4 + 1/3.95 + 1/6.7 + 1/9 + 1/14 = 1.027` The odds come out to be over than 1 due to the bookmaker overround. This is problematic for me because I need to compare the odds offered at the bookmaker and the probability from my model to determine if there is any value placing a bet. In order to compare them accurately I will need to either   1. Remove the overround from the bookmaker prices (deflate the probabilities)   2. Inflate the output of my model to include the bookmaker's overround Can anyone suggest how to do this accurately from a theoretical point of view? I believe that doing this linearly - applying the same factor to each horse - is incorrect. I believe that each horse should take a percentage of the over round based on their odds of winning, the horses that are more likely to win should take a bigger proportion of the overround. Does this sound more correct? Thanks for your help.