title: Asymptotic normality of MLE in exponential with higher-power x

text: Given the distribution: $f(x;\theta) = \frac{3}{\theta}x^2e^{-x^3/\theta}$ if $x>0$ the MLE for $\theta$ is $\frac{1}{n}\sum_{i=1}^n x_i^3$. It's an unbiased estimator with variance $\theta^2/n$. The Fisher information number is $1/\theta^2$. Now, I am supposed to check whether or not this estimator shows asymptotic normality. The problem here is that my textbook has a condition which I can't find in any other source: "$\frac{d^2}{d\theta^2} \ln f(x;\theta)$ is a continuous function in $\theta$, uniformly in $x$." It's this uniformly continuity that I can't seem to find anywhere else. This seems to be the only condition not met for asymptotic normality. I've found that $\frac{d^2}{d\theta^2} ln f(x;\theta) = \frac{1}{\theta^2} - 2\frac{x^3}{\theta^3}$. Since uniform continuity is equivalent with the boundedness of the derivative, I derived this to $x$ to find $ - 6\frac{x^2}{\theta^3}$, which is not bounded for x. My textbook also mentions that "these conditions are in particular fulfilled for (identifiable) members of the exponential class". So, my question is, is this really a necessary condition, and am I right that this MLE is not asymptotically efficient? Thanks.


title: How to diagnose multicollinearity using the output of vif function in R?

text: I am running a logistic regression in R and am attempting to determine if multicollinearity is a problem with my model.   When I run `vif()` on my final model, I get `GVIF` and `GVIF^1/(2*Df)` columns. From what I have read `GVIF^1/(2*Df)` is what I should use to assess muticollinearity, but I have been unable to determine what values I should use as a cut-off point. Any help would be greatly appreciated.


title: Can I interpret the inclusion of a quadratic term in logistic regression as indicating a turning point?

text: In a Logistic Regression with linear and quadratic terms only, if I have a linear coefficient $\beta_1$ and quadratic coefficient $\beta_2$, can I say that that there is turning point of the probability at $-\beta_1 / (2\beta_2)$?


title: Worst-case error related to Cramer-Rao bound

text: Asked this previously on Math.Stackexchange, maybe this fits here. I would like to understand the relation (if any) between the Cramer-Rao Lower Bound of estimation theory and the following simple definition of "reconstruction accuracy" which doesn't use any statistical framework. Let $\theta\in\Theta\subset \mathbb{R}^n$ be an unknown parameter vector and let $x=f(\theta)\subset\mathbb{R}^m$ be a set of measurements, from which $\theta$ is to be reconstructed. Assume that the absolute error in each measurement is bounded from above by $\epsilon>0$, and otherwise nothing is known about the measurement error. Let $T_{\theta,\epsilon}=f^{-1}(x+B(0,\epsilon))$ (the full preimage, always well- defined). Obviously, $\theta\in T_{\theta,\epsilon}$. The "best possible" reconstruction accuracy is defined as $\delta_{\theta,\epsilon}=\sup_{t\in T_{\theta,\epsilon}} \|t-\theta\|$. (Informal justification for the definition: if the measurement error is bounded by $\epsilon$, the original parameter vector $\theta$ can only be recovered up to accuracy $\delta_{\theta,\epsilon}$ because any parameter vector $\hat\theta\in T_{\theta,\epsilon}$ could have produced a measurement vector $\hat x \in B(x,\epsilon)$.) In contrast, Cramer-Rao lower bound depends on some distribution $p(x;\theta)$ of the measurements, and gives a bound on the variance of any unbiased estimator $\hat \theta$ of $\theta$ satisfying some regularity conditions. Apparently, the two notions are different and in general not comparable. Still, are there cases in which the Cramer-Rao bound can be thought of as the $\delta_{\theta,\epsilon}$ defined above? ## Example Consider the simple case of a single parameter. Assuming $f'(\theta)\neq 0$, first-order approximation of the reconstruction accuracy as defined above is $$\delta_{\theta,\epsilon}\approx \frac{\epsilon}{f'(\theta)}. \qquad (1) $$ On the other hand, assuming that the measurements are modeled as $$ x = s(\theta) + w$$ where $w$ is Gaussian with zero mean and variance $\sigma^2$, we have (see e.g. Kay's book, section 3.5) $$ CRB\{\theta\} = \frac{\sigma^2}{\biggl(\frac{\partial s}{\partial \theta}\biggr)^2} \qquad (2)$$ Comparing $(1)$ and $(2)$ seems to suggest that in this case the notions are very closely related.


title: Differences between Baum-Welch and Viterbi training

text: I am currently using Viterbi training for an image segmentation problem. I wanted to know what the advantages/disadvantages are of using the Baum-Welch algorithm instead of Viterbi training.


title: Adding gamma and 3-parameter log normal distributions to L-moments ratio diagram lmrd()

text: How to adapt this piece of code but for: - gamma distribution - 3 parameter log normal More specifically, where can I find the specification of the parameter (lmom) for pelgam() and pelln3()? Lmom package info just gives: pelgam(lmom), lelln3(lmom), where lmom is a numeric vector containing the L-moments of the distribution or of a data sample.               # Draw an L-moment ratio diagram and add a line for the     # Weibull distribution     lmrd()     weimom <- sapply( seq(0, 0.9, by=0.01),     function(tau3) lmrwei(pelwei(c(0,1,tau3)), nmom=4) )     lmrdlines(t(weimom), col='darkgreen', lwd=2)      source: http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=lmom:lmrdpoints Thank you.


title: How to calculate ratings/rankings from Paired comparison / Pairwise comparison of large data-sets?

text: As part of my graduate thesis (area: psychology) I have gathered preference data. The data includes approximately 50000 heads-up comparison between elementX and elementX. I have a total of 15 elements. The participant would be shown 36 comparisons in random order, of which he/she would then choose the preferred element. I need to be able to calculate a meaningful "interval" scale-value for each element in order to see if the actual preferences of each element can predict theoretically generated scale values using regression analysis. The Data: 14 elements, forming a 14x14 matrix table 91 unique possible comparisons ((14*14-14)/2) = 50000 comparison (NOT evenly distributed across all comparisons) The problem: Since showing the participant all 91 possible comparison would be too time consuming I've just presented each participant 36 comparisons. Moreover, each participant was not eligible for each element, the eligible elements was determined if the participant owned the element. For example: ParticipantX owns element; 1,5,6,8 and 10. And would therefore only be shown comparisons that includes those elements. Furthermore, as owning some elements were considered more rare, I made a priority list for each element, so that element1 > element2 > .... > element14 For example: If the participant would be eligible for 50 possible comparisons (but restricted by the 36 comparison-cap) the comparisons shown to the participant would be based on the rarity order. I'm afraid I destroyed the data by making this rarity order, since it resulted in some some elements being compared unequally to some elements; For example: of 1000 comparisons that element 1 was included in, 30% was with element3, 10% with element4 etc... My questions: What method would you suggest would be the best for calculating "scale" values for each element? Currently I've calculated the basic probability that a given element has been choosen over another element. For example; element1 being chosen 95% of times to be preferred over any other element. However, I would prefer a more refined method to calculate the element values. Thanks for your time, Chris, Finland.


title: SVMPerf - only one vector in the model file?

text: Am I missing something? When I use SVMLight, there's always a lot of vectors in the model file. But when I create one with SVMPerf - there's only one vector in the file? (?) But still, it works fine. Am I missing something? Also the number of documents that it supposedly used is way smaller than I provided him with ( I supplied 400,000 examples and it says it used 165 training documents? )


title: Use of Wilcoxon test for non-normal data akin to Two One Sided T-test

text: I'm analysing paired data for equivalency and it's not normally distributed, i.e., the difference of the paired result is not normally distributed due to, amongst other things, outliers. If it were normally distributed I would use a Two One Sided T-test (TOST). My questions are:   1. Can two one-sided Wilcoxon tests with suitably defined limits ("mu" in R terminology) be used just like the Two One Sided T-test approach?   2. If so, is it possible to establish a single p-value? Or is the larger p-value for the limit specified used?   3. Are there any published references for doing this? (I have searched but couldn't find anything immediately available, although my search may not have been thorough enough.) Thanks in advance!


title: Measurement error in maximum counts

text: I'm familiar with the concept of a mean value of data and the variation around the mean. Is it possible to quantify variation around maximum values? For example, take the below data collected across 10 years. I want to present the maximum value for each month, but I also want to quantify the variation around the maximum value for each month across the 10 years:               counts <- data.frame(year = sort(rep(2000:2009, 12)), month = rep(month.abb,10), count = sample(1:500, 120, replace = T))      First 20 rows of data:               head(counts, 20)             year month count     1  2000   Jan    14     2  2000   Feb   182     3  2000   Mar   462     4  2000   Apr   395     5  2000   May   107     6  2000   Jun   127     7  2000   Jul   371     8  2000   Aug   158     9  2000   Sep   147     10 2000   Oct    41     11 2000   Nov   141     12 2000   Dec    27     13 2001   Jan    72     14 2001   Feb     7     15 2001   Mar    40     16 2001   Apr   351     17 2001   May   342     18 2001   Jun    81     19 2001   Jul   442     20 2001   Aug   389      What quantities can I use: Standard deviation? Interquartile range? Range of maximum values? Confidence intervals?


title: Why are distributions important?

text: This may as well go down as the silliest questions ever asked on this forum, but having received sound and meaningful answers to a previous question, i thought i will stretch my luck again. I have been very confused for some time on the importance of statistical distributions especially as they relate to asset returns and even more specifically in asset allocation. My question to be specific is this: Assume i have 20 years of S&P 500 monthly returns data, why should i need to assume a certain kind of distribution (i.e Normal/Johnson/Levy flight etc) for my asset allocation decision when i can simply just make my asset allocation decisions based on the historical data i have with me?


title: Gaussian Distribution Prediction

text: What are some common methods of making distribution predictions? I have a set of features $x_1,x_2,x_3$ which map to Gaussian distributions ($\mu,\sigma^2$). That is, the feature vector of a single example $X_1$ can be defined by a distribution with ($\mu_1,\sigma^2_1$). What are some ways to go about building a model which can predict ($\mu,\sigma^2$) given a feature vector $X$? I think my biggest issue here is that I don't know which terms to search for, so I'm having trouble finding relevant papers or software packages. Recommended R libraries would be a bonus!


title: Hypothesis Testing question

text: In a training center, a new strategy was implemented. After the implementation of this new strategy the training center claimed that on an average $80\%$ of the students passed. In order to find out whether this claim can be justified, we take a random sample of $30$ students and see that $60\%$ of the students passed. Do we have sufficient evidence to accept the $80\%$ claim?


title: Association rules vs decision tree vs rule learning

text: Is it ok to generate a classifier by collecting all association rules so that the conclusion part refers to the target variable? does it perform better than decision tree or rule learning?


title: Random variables for which the distribution of the sum of the RV with a Gaussian RV is known

text: For a counter example, I am searching for random variables $Y$ such that for a independent normal random variable $X$ the distribution of $Z=Y+X$ is known parametrically. Ideally, the Shannon entropy of $Z$ should also be known. The obvious case is when $Y$ is normal as well. Another example would be $Y$ being Bernoulli. In that case, the distribution of $Z$ is a mixture of two Gaussians. However, the entropy of this distribution is unknown. Are there other/better examples?


title: Cluster analysis on ordinal data

text: I'd like to perform a cluster analysis on ordinal data (Likert scale) by using SPSS. I have around 140 observations and 20 variables that are scaled from 1 to 5 (1: I strongly agree, 3: neutral, 5: I strongly disagree). As a result, I want to assign one cluster to each person, such as person 1 belongs to the group of technology-enthusiastic users... Since I'm new in that area, I want to ask you, if you might help me with regard to the procedure. I have already read something about cluster analysis in the Internet, but there are different opinions on how to treat ordinal data. Is it a common way to convert the ordinal data to interval data first and then perform the cluster analysis?


title: Linear Regression-Categorical variable with 93 levels

text: I am a beginner in SAS. I have a data set called cars with 25 variables and 93 observations. I have attached the image of the data set.![enter image description here](http://i.stack.imgur.com/YO4Fl.png) Here you can see that there is a variable `Model`: all values are different. The variable `Manufacturer` too has more than 20 levels. I am not sure what how to deal with these two variables, `Model` and `Manufacturer`. I want to build linear regression model for predicting the maximum price of a car.


title: Bayes' Rule - Law of Conditional Probability

text: A market research firms conducts studies regarding the success of new products. The company is not always perfect in predicting the success. Suppose that there is a 50% chance that any new product would be successful (and a 50% chance that it would fail). In the past, for all new products that ultimately were successful, 80% were predicted to be successful (and the other 20% were inaccurately predicted to be failures). Also, for all new products that were ultimately failures, 70% were predicted to be failures (and the other 30% were inaccurately predicted to be successes). What is the a priori probability that a new product would be a success? a) 0.50   b) 0.80   c) 0.95   d) 0.70   e) 0.60 I don't understand how to do this problem. This uses Bayes' Rule. Can someone show the steps to do this? This is a problem in my textbook. No its not homework, its for my own studying.


title: Why is my R-squared so low when my t-statistics are so large?

text: I feel this may be a slightly dumb question but I ran a regression with 4 variables, and all are very statistically significant, with T values $\approx 7,9,26$ and $31$ (I say $\approx$ because it seems irrelevant to include the decimals) which are very high and clearly significant. But then the $R^2$ is only .2284. Am I misinterpreting the t values here to mean something they're not? My first reaction upon seeing the t values was that the $R^2$ would be quite high, but maybe that is a high $R^2$?


title: Using ANCOVA with data matched on propensity scores

text: I've read papers comparing them, but never seen a study that used them together. Is this done? Why or why not? Suppose you use ANCOVA to analyze a reduced sample of matched pairs generated using propensity scores (assuming the propensity scores are estimated from the complete set of confounders X that jointly effect the outcome variable Y and the treatment variable T). Suppose the model has T (considered as a grouping variable for ANCOVA) along with covariate Z predicting the continuous variable Y. Will estimates of the effect size for the treatment and the regression coefficient for Z be biased if Z is included in the set X or if Z and Y are not independent? I would love to find articles that discuss ideas along these lines.


title: Ensemble time series model

text: I need to automate time-series forecasting, and I don't know in advance the features of those series (seasonality, trend, noise, etc). My aim is not to get the best possible model for each series, but to avoid pretty bad models. In other words, to get small errors every time is not a problem, but to get big errors once in a while is. I thought I could achieve this by combining models calculated with different techniques. That is, although ARIMA would be the best approach for a specific series, it may not be the best for another series; the same for exponential smoothing. However, if I combine one model from each technique, even if one model isn't so good, the other will bring the estimate closer to the real value. It is well-known that ARIMA works better for long-term well-behaved series, while exponential smoothing stands out with short-term noisy series.   * **My idea is to combine models generated from both techniques in order to get more robust forecasts, does it make sense?** There might be many ways to combine those models.   * **If this is a good approach, how should I combine them?** A simple mean of forecasts is an option, but maybe I could get better predictions if I weight the mean according to some goodness measure of the model.   * **What would be the treatment of the variance when combining models?**


title: PACF manual calculation

text: I am trying to replicate the calculation that SAS and SPSS do for the partial autocorrelation function (PACF). In SAS it is produced through Proc Arima. The PACF values are the coefficients of an autoregression of the series of interest on lagged values of the series. My variable of interest is sales so I calculate lag1, lag2...lag12 and I run the following OLS regression: $$Y_t=a_0+a_1Y_{t-1}+a_2Y_{t-2}+a_3Y_{t-3}+\ldots+a_{12}Y_{t-12}.$$ Unfortunately the coefficients that I get are not even close to the PACF (lags 1 to 12) that SAS or SPSS provide. Any suggestions? Is there something wrong? What comes to my mind is that the least squares estimation of this model might not be appropriate and maybe another estimation technique should be used. Thanks in advance.


title: Power analysis ordinal probit regression

text: Suppose I have $n=100$ observations of ordinal data and get threshold coefficients $b_1, \dots, b_3$ and probit slope $b_4$. I want to test the hypothesis $H_{0}: \frac{b_{3}}{b_{4}} = \frac{1}{2}$ vs. $H_a: \frac{b_{3}}{b_{4}} \neq \frac{1}{2}$. So I want to count the proportion of times we fail to reject the null hypothesis. To do this, would we just keep sampling 100 observations similar to the previous ones and return the proportion of p-values greater than $0.05$?


What happens if you reject normality of residuals when estimating with least squares? Is it too important to have normality on the residuals?


title: How is the bayesian framework better in interpretation when we usually use uninformative or subjective priors?

text: It is often argued that the bayesian framework has a big advantage in interpretation (over frequentist), because it computes the probability of a parameter given the data - $p(\theta|x)$ instead of $p(x|\theta)$ as in the frequentist framework. So far so good. But, the whole equation it is based on: $p(\theta|x) = {p(x|\theta) . p(\theta) \over p(x)}$ looks to me little suspicious for 2 reasons:   1. In many papers, usualy uninformative priors (uniform distributions) are used and then just $p(\theta|x) = p(x|\theta)$, so bayesians get the same result as frequentists get - so how is then bayesian framework better in interpretation, when bayesian posterior and frequentists likelihood are the same distributions? It just yields the same result.   2. When using informative priors, you get different results, but the bayesian is affected by the subjective prior, so the whole $p(\theta|x)$ has the subjective tinge too. In other words, the whole argument of $p(\theta|x)$ being better in interpretation than $p(x|\theta)$ builds on a presumption that $p(\theta)$ is kind of "real", which normally is not, it is just a starting point we somehow choose to make the MCMC run, a presumption, but it is not a description of reality (it can't be defined I think). So how can we argue that bayesian is better in interpretation?


title: Intuition for consequences of multicollinearity

text: So we have a regression equation with one explained variable and 10 explanatory variables. What I have read so far:   1. Multicollinearity doesnt affect the regression of the model as a whole.   2. But if we start looking at the effect of individual predictor variable Xs on the explained variable, then we are going to have inaccurate estimates.  I have tried to start thinking intuitively about it like follows:   * High Multicollinearity means that in a matrix, two or more rows/columns are linearly dependent on each other. In other words, in a 3-dimensional space there are 2 vectors which have (almost) the same direction only different magnitudes (is this right?) I'd appreciate it if someone could explain how this translates into "multicollinearity not affecting the regression as a whole but only individual variable's coefficient estimates". Also, could someone explain the statement in bold? I cant make sense out of it: > One other thing to keep in mind is that the tests on the individual > coefficients each assume that all of the other predictors are in the model. > **In other words each predictor is not significant as long as all of the > other predictors are in the model. There must be some interaction or > interdependence between two or more of your predictors.** which was an answer to this question: How can a regression be significant but all predictors insignificant?


title: Convergence proof for perceptron algorithm with margin

text: I was looking for the proof of convergence for perceptron algorithm with margin. I was not able to find it in any pattern classification text book or over the internet. Can anyone here please point to a text or reference for the above?


title: Bootstrap confidence interval for a minimization problem

text: Based on a simulated empirical distribution of $Z$, find $\mu^*_X$, the smallest $\mu_X$ such that $P(Z > c) = p$, where $Z$ is given by $Z = \prod^5_{i=1}(aX_i + bY_i)$ where the $X_i$'s are i.i.d $LN(\mu_X,\sigma_X^2)$, the $Y_i$'s are i.i.d $LN(\mu_Y,\sigma_Y^2)$ and $LN(,)$ is the Log-Normal distribution. (Meaning that $X_i = e^{\mu_X + \sigma_XW}$ where $W$ is standard normal, and similarly for $Y_i$.) Furthermore, the $X_i$'s and $Y_i$'s are independent. If I simulate a large number, $10M$ say, of i.i.d $N(0,1)$ outcomes, I can compute $M$ values of $Z$, $Z^1, \ldots Z^M$. I can do this for different values of $\mu_X$, with the same simulated data set since all that changes is the value of $\mu_X$ Then, I can estimate $P(Z > c)$ by $\#\{i:Z^i > c\}/M$ for different $\mu_X$. We know that the probability that a simulated sample $Z^i$ is greater than $c$ is by definition $p$. Hence, for every $\mu_X$ I use, the estimator $\#\{i:Z_i > c\}/M$ is a binomial proportion, for which we can compute confidence intervals with the Normal approximation, but this is not satisfactory as we cannot make inferences about $\mu^*_X$ that way (or can we?). Each time I simulate $10M$ standard normal variates, I expectedly get different $\mu^*_X$ , but what kind of method allows me to quantify that uncertainty? How would it work in this case? Thanks.


title: IV and Mediator have reciprocal relation

text: I am trying to create a mediation model for: X1=employee capability; X2 (Mediator)=Company capability; Y= employee earnings Although, I am not sure if I can still say that this will be a mediation model when I believe and support that X1 and X2 have reciprocal relation. Where, employee capability (X1) explains the company capability (X2) and vice-versa then both affect employee earnings(Y). First, is this considered a usual mediation model? Second, whatever model it will be, how do I proceed with the analysis/process? Thanks.


title: How to decide hyperparameters using empirical Bayes, with prior distribution a mixture of 2 normals?

text: I have a batch of(~200) velocity data of celestial bodies that I know they are either type A or B, supposedly with respective normal distribution of velocity. I want to use these observed data to decide the mean and variance of the two, as well as ratio between 2 types of stars. Any suggestion? PS I'm quite unfamiliar with either empirical Bayes or mixture models, so I really hope for some step-by-step guidance instead of tons of rigorous derivations or advanced techniques...


title: EMD-SVR method for time series

text: Has anyone attempted to use Empirical Mode Decomposition(EMD)-Support Vector Regression(SVR) in nonstationary time series forecasting? It seems quite interesting. As I observed it has high performance compared with single SVR. But I didn't find it as a real approach for time series.


title: Categorical predictor variables in multiple regression?

text: Can I run multiple regression with mostly categorical independent variables? I'd appreciate some examples to cite this approach as one of my reviewers is insisting that I justify the use of multiple regression for categorical IVs.               2 DVs (test scores) interval          IVs -     1 score (1-5 scale)     Gender dummy coded     Age dummy coded 3 levels     Member of ethnic minority dummy coded     Time spent out of country dummy coded 4 levels     Academic Discipline dummy coded 4 levels      Since the independent variables are mixed (mostly categorical) is multiple regression still acceptable? I want to explore the relationship of the IVs independently and in combinations on the test scores (DVs).


title: Account for both between and within variance in panel data

text: I have a general question on fixed-effects and mixed-effects models for panel data. I am doing a logistic regression on panel data, with data measured on the individual level. I understand that by using fixed-effects in my model (`xtlogit...,fe` in Stata) I basically ignore or do not look at the difference **between** individuals, but rather use the **within** variance of the IVs to calculate the coefficients. **Question:** But what if there's also considerable difference between the individuals that I would like to capture? Is it even possible to account for both the unobserved heterogeneity (through fixed-effects) and the variance between individuals? I read a lot about Random Effects, Population-averaged models, between- estimator and mixed-effects but I can't bridge my intellectual gap on the question "when to use which model?". I feel that the mixed-effects logistic regression (`xtmelogit` in Stata) might apply, but the term _random intercept_ got me thinking if that really captures the variance between individuals in panel data. Also my data is not _nested_ as in most of the examples on mixed-effects models which speak of clusters or group levels (e.g. students in schools).


title: How can I calculate an average speed rate that ignores occasional anomalies?

text: Horrible title, so let me explain - I think this must be a very basic statistical problem, but I don't remember what little statistics I once knew. I'm writing a low-level emulator of an old IBM CPU. Each machine instruction of the IBM1620 was accomplished by a sequence of various "machine cycles". On the original machine, each cycle took exactly 20 microseconds. I'm emulating the various cycles with subroutines; they are much faster than 20 microseconds, and they do not all take the same amount of time to execute. The 1620 had a magnificent light console that constantly showed the state of over 100 internal bits (data, registers, control switches) as the machine ran. I'm emulating that, too. Unlike the 1620, I can't manage to update my Java GUI 50000 times per second, as repainting after every 20 microsecond "cycle" would require. But I get a fine visual approximation simply by repainting the GUI every **N** cycles. [We can assume that repainting the GUI takes constant time.] Without repainting the GUI at all, the emulator runs much too fast; repainting after every cycle makes it run much too slow. The graphical repainting is the ideal "friction" to apply to the emulator, to get it to approximate the true speed of the machine. The question is, what is **N**? I have probes to count cycles and sample the system clock. The problem is that, due to several factors, the most important of which is that _my emulations of the various "cycles" do not take equal time_ , the number of cycles in a measurement interval can vary a good deal, depending on the nature of the instructions the CPU was carrying out at that time. I'd like, over time, to be able to converge on a "good" value of **N** , using some sort of method that can somewhat overlook anomalous readings. "Good"? "Anomalous"? I don't know how to define them. But if anyone can give guidance, even along the lines of "oh! you should check the Wikipedia article for **The Backwards Freight Train Problem** " [I'm making that up, of course], I'd appreciate it. EDIT: to illustrate the problem, I'm appending some measurements. Here are a series of measurement intervals, each lasting 10000 cycles, showing:   * interval time in microseconds (usecs),    * the number of cycles executed between redrawing the updated GUI (cycs/redraw=51), currently constant,   * the average GUI redraws per second (redraws/sec),   * the cycles per second, expressed as a percentage of a real IBM1620, which was always exactly 50000 cycles per second. And the question is: how might one vary **cycs/redraw** to bring **Rate** closer to 100%? Especially given that, as is obvious below, emulated cycles vary in duration. * * *               CPU: Interval: Cycles=10000; usecs=174595; cycs/redraw= 51; redraws/sec:1123.  Rate:114%     CPU: Interval: Cycles=10000; usecs=197842; cycs/redraw= 51; redraws/sec: 991.  Rate:101%     CPU: Interval: Cycles=10000; usecs=179968; cycs/redraw= 51; redraws/sec:1089.  Rate:111%     CPU: Interval: Cycles=10000; usecs=157880; cycs/redraw= 51; redraws/sec:1241.  Rate:126%     CPU: Interval: Cycles=10000; usecs=151724; cycs/redraw= 51; redraws/sec:1292.  Rate:131%     CPU: Interval: Cycles=10000; usecs=150892; cycs/redraw= 51; redraws/sec:1299.  Rate:132%     CPU: Interval: Cycles=10000; usecs=151908; cycs/redraw= 51; redraws/sec:1290.  Rate:131%     CPU: Interval: Cycles=10000; usecs=172101; cycs/redraw= 51; redraws/sec:1139.  Rate:116%     CPU: Interval: Cycles=10000; usecs=159090; cycs/redraw= 51; redraws/sec:1232.  Rate:125%     CPU: Interval: Cycles=10000; usecs=165604; cycs/redraw= 51; redraws/sec:1184.  Rate:120%     CPU: Interval: Cycles=10000; usecs=135410; cycs/redraw= 51; redraws/sec:1448.  Rate:147%     CPU: Interval: Cycles=10000; usecs=130933; cycs/redraw= 51; redraws/sec:1497.  Rate:152%      * * *


title: Using a deterministic trend in OLS

text: I have a time series data, see jpg. It goes straight up and then straight down, later tailing off. This seems to fit my textbook’s description of a deterministic trend being an almost exact function of time on the way up and on the way down. Using the D-F test I found that I could reject the existence of a unit root at the 5% level but not at the 1% level, and if I split the series into two parts, up and down, I cannot reject the unit root at even the 10% level. Is it valid to use this time series, as it is (ie not differenced), in an OLS regression on the basis that it is a deterministic trend and therefore ‘trend stationary’? I want to run it as an explanatory variable against another series which is difference stationary. ![nbl full](http://i.stack.imgur.com/LepDH.jpg) When I did it, I found the answer to be reasonable (ie something similar to what I was expecting), the residuals did not show autocorrelation in a D-W test and I can reject a unit root in a D-F test at the 1% level.


title: how the uncertainty of the coefficient in the model affect the uncertainty of the fit

text: I am wondering how to estimate the uncertainty of the fit if the coefficients in the model have uncertainties as well. for example: I try to use model **y=exp(-(a+b)*x)+c** to fit a data set **x** and **y** to find out the value of **b** and **sigma_b** , where **a** and **c** are predetermined coefficients and have uncertainty **sigma_a** and **sigma_c** , how should I do it?? Any help is greatly appreciated.


title: Checking assumptions lmer/lme mixed models in R

text: I ran a repeated design whereby I tested 30 males and 30 females across three different tasks. I want to understand how the behaviour of males and females is different and how that depends on the task. I used both the lmer and lme4 package to investigate this, however, I am stuck with trying to check assumptions for either method. The code I run is               lm.full <- lmer(behaviour ~ task*sex + (1|rat/task), REML=FALSE, data=dat)     lm.full2 <-lme(behaviour ~ task*sex, random = ~ 1|rat/task), method="ML", data=dat)      I checked if the interaction was the best model by comparing it with the simpler model without the interaction and running an anova:               lm.base1 <- lmer(behaviour ~ task+sex+(1|rat/task), REML=FALSE, data=dat)     lm.base2 <- lme(behaviour ~ task+sex, random= ~1|rat/task), method="ML", data=dat)     anova(lm.base1, lm.full)     anova(lm.base2, lm.full2)      Q1: Is it ok to use these categorical predictors in a linear mixed model?   Q2: Do I understand correctly it is fine the outcome variable ("behaviour") does not need to be normally distributed itself (across sex/tasks)?   Q3: How can I check homogeneity of variance? For a simple linear model I use `plot(LM$fitted.values,rstandard(LM))`. Is using `plot(reside(lm.base1))` sufficient?   Q4: To check for normality is using the following code ok?               hist((resid(lm.base1) - mean(resid(lm.base1))) / sd(resid(lm.base1)), freq = FALSE); curve(dnorm, add = TRUE)


title: What sampling method is this?

text: Please provide your input on what kind of sampling methodology is followed here.   * Data being collected for a questionnaire.   * the sample population consists of people in employment and IT professionals   * they were targeted by sending mails and posting messages on professional groups on social media sites. _**what sampling method would this be ?_** Its not random sampling because chose only to invite people on professional groups and discussion forums for IT professionals.


title: Pooled time series

text: Can somebody tell me   1. what are and how to do **pooled time series**?   2. and in what way it can help me to **build forecast models from clusters of time series?** Maharaj and Brett do this in their article http://ideas.repec.org/p/msh/ebswps/1999-9.html Thanks.


title: Is machine learning less useful for understanding causality, thus less interesting for social science?

text: My understanding of the difference between machine learning/other statistical predictive techniques vs. the kind of statistics that social scientists (e.g., economists) use is that economists seem very interested in understanding the effect of a single or several variables -- both in terms of magnitude and detecting whether the relationship is causal. For this, you end up concerning yourself with experimental and quasi-experimental methods, etc. Machine learning or statistical modeling that is predictive often entirely neglects this aspect and in many cases doesn't give you a specific degree to which one variable affects the outcome (logit and probit do seem to do both). A related question is to what extent do theoretically inspired economic or behavioral models have an advantage over atheoretical models when predicting to new domains? What would a machine learning or prediction-oriented statistician say to the criticism that without an economic model, you wouldn't be able to correctly predict new samples where covariates were very different. I'd be really happy to hear people's take on this from all perspectives.


title: How does the Goodman-Kruskal gamma test and the Kendall tau or Spearman rho test compare?

text: In my work, we are comparing predicted rankings versus true rankings for some sets of data. Up until recently, we've been using Kendall-Tau alone. A group working on a similar project suggested we try to use the Goodman-Kruskal Gamma test instead, and that they preferred it. I was wondering what the differences between the different rank correlation algorithms were. The best I've found was this answer, which claims Spearman is used in place of usual linear correlations, and that Kendall-Tau is less direct and more closely resembles Goodman-Kruskal Gamma. The data I'm working with doesn't seem to have any obvious linear correlations, and the data is heavily skewed and non-normal. Also, Spearman generally reports higher correlation than Kendall-Tau for our data, and I was wondering what that says about the data specifically. I'm not a statistician, so some of the papers I'm reading on these things just seem like jargon to me, sorry.


title: How do you compare OR between different groups using SAS?

text: How do you get a confidence interval for an odds ratio between two different groups (both of which are not the reference group) in SAS? Suppose $x_1$ is the reference level and $x_2,x_3,x_4$ are the other variables. Odds ratios in SAS give them for each level compared to the reference (e.g. $x_2$ vs $x_1$, $x_3$ vs $x_1$ etc..). How do you get the odds ratio between $x_2$ and $x_3$ in SAS including confidence intervals?


title: Goodness-of-fit test without analytical PDF and CDF

text: I have closed form moment-generating function and characteristic function of a distribution, which describes waiting time of a continuous univariate random process. However, I cannot analytically invert the characteristic function to get PDF. Now, I have run simulations of the original process in order to verify that my model adequately describes the process. Because data in this case are simulated, I exactly know model parameters. Therefore, I am not interested in parameter estimation for now. By eyeballing histogram of simulated data and numerically computed PDF, my model seems to fit the data very well. However, I want to say something more concrete and quantitative. Is there any goodness-of-fit test with arbitrary distribution, whose PDF and CDF are not analytically available? If not, what other approaches can I take?


title: Odds ratio interpretation

text: Does it make sense to say that the odds of some risk for a person in group A is $0.4$ times lower that that of someone in group B? Or would it be better to phrase it in "higher than" language?


title: How do you predict a response category given an ordinal logistic regression model?

text: I want to predict a health problem. I have 3 outcome categories that are ordered: 'normal', 'mild', and 'severe'. I wish to predict this from two predictor variables, a test result (a continuous, interval covariate) and family history with this problem (yes or no). In my sample, the probabilities are 55% (normal), 35% (mild), and 10% (severe). In this sense, I could always just predict 'normal' and be right 55% of the time, although this would give me no information about individual patients. I fit the following model: the cut point for (y >= 1) = -2.18 the cut point for (y >= 2) = -4.27 $\beta$(test) = .60 $\beta$(family history) = 1.05 Assume there is no interaction and everything is fine with the model. The concordance, c, is 60.5%, which I understand to be the maximum predictive accuracy the model affords. I come across two new patients with the following data: 1. test = 3.26, family = 0; 2. test = 2.85, family = 1. I want to predict their prognosis. Using the formula: exp(-$X\beta$ -cutPoint) / (1+exp(-$X\beta$ -cutPoint)) (and then taking the differences amongst the cumulative probabilities), I can calculate the probability distribution over the response categories conditional on the model. R code (n.b., due to rounding issues, the output doesn't match perfectly):               cut1 <- -2.18     cut2 <- -4.27     beta <- c(0.6, 1.05)     X <- rbind(c(3.26, 0), c(2.85, 1))     pred_cat1 <- exp(-1*(X%*%beta)-cut1)/(1+exp(-1*(X%*%beta)-cut1))     pred_cat2.temp <- exp(-1*(X%*%beta)-cut2)/(1+exp(-1*(X%*%beta)-cut2))     pred_cat3 <- 1-pred_cat2.temp     pred_cat2 <- pred_cat2.temp-pred_cat1     predicted_distribution <- cbind(pred_cat1, pred_cat2, pred_cat3)     predicted_distribution      Namely: 1. 0 = 55.1%, 1 = 35.8%, 2 = 9.1%; and 2. 0 = 35.6%, 1 = 46.2%, 2 = 18.2%. My question is, how do I go from the probability distribution to a predicted response category? I have tried several possibilities using the sample data, where the outcome is known. If I just pick max(probabilities), accuracy is 57%, a slight improvement over the null, but below the concordance. Moreover, in the sample, this approach never picks 'severe', which is what I really want to know. I tried a Bayesian approach by converting null and model probabilities into odds and then picking the max(odds ratio). This does pick 'severe' occasionally, but yields worse accuracy 49.5%. I also tried a sum of the categories weighted by the probabilities and rounding. This, again, never picks 'severe', and has low accuracy 51.5%. What is the equation that takes the information above and yields optimal accuracy (60.5%)?


title: Fit an arbitrary PDF

text: I'm looking for a way to fit several arbitrary discrete distributions to another discrete distribution using maximum likelihood, but I'm having trouble finding a function to do this in `numpy` or `scipy`. I'm ok with using another library, so long as   * it's in C++ or Python   * it is designed specifically for fitting (trying to avoid GUIs and massive packages)  The fitted distributions come from toy Monte Carlo, so they are represented by a series of ordered pairs.


title: Multiple comparison with genotype data

text: I am looking at the effect of genotype on a continuous variable. In total 16 comparisons were done involving 9 different SNPs. The number of significant events was 3 out of 16 tests (Mann-Whitney U-test). Is it necessary to correct the p-values to account for multiple comparisons?


title: Gaussian Like distribution with higher order moments

text: For the Gaussian distribution with unknown mean and variance, the sufficient statistics in the standard exponential family form is $T(x)=(x,x^2)$. I have a distribution that has $T(x)=(x,x^2,...,x^{2N})$, where N is kind of like a design parameter. Is there a corresponding known distribution for this kind of sufficient statistics vector? I need samples from this distribution so it is kind of crucial for me to get exact samples from the distribution. Thanks a lot.


title: What methods can be used to specify priors from data?

text: ## Background I am generally interested in learning appropriate methods of using data to specify priors. A previous question asks how to elicit priors from experts and received some good recommendations. Here, I am interested in learning how to specify a prior using data. I plan to use these priors in a meta-analysis to synthesize additional data that I collect. _update_ Although John provides a 'correct' answer, in my case, it would require substantial modification of the original model to implement, so I would prefer to find a way to estimate the prior as a discrete step. ## Questions   1. What is the best way to specify such a data informed prior?   2. If I am working with parameters for a particular species (monkeys), and this species belongs to a group of organisms(primates), and data are available for primates but not for monkeys, would it be appropriate to fit a distribution based on the primate data? ## Example cases, first with proposed solution   1. I have 100 observations from 100 primate species of primate thumb length:              set.seed(0)     thumb <- rgamma(100, 4, 0.1)     library(MASS)     fitdistr(thumb, 'gamma')      Indeed, when there is no apriori reason to select a particular distribution, the distribution can be chosen by maximum likelihood:              for(dist in c('gamma', 'lognormal', 'weibull') {     logLik(fitdistr(thumb, dist))     }        2. I have collected 50 means, standard errors, and sample sizes from 50 different primate species, and 50 independent observations from another 50 species of eye diameter:              eye <- data.frame( diameter = rgamma(100, 4, 0.1),                        se       = c(rlnorm(50, 0.5,1), rep(NA, 50)),                        n        = c(rep(1:5, 10), rep(1, 50)))           eye <- signif(eye, 3)      How can I incorporate the sample statistics into my calculation of a prior?


title: For which distributions are the parameterizations in BUGS and R different?

text: I have found some distributions for which BUGS and R have different parameterizations: Normal, log-Normal, and Weibull. For each of these, I gather that the second parameter used by R needs to be inverse transformed (1/parameter) before being used in BUGS (or JAGS in my case). Does anyone know of a comprehensive list of these transformations that currently exists? The closest I can find would be comparing the distributions in table 7 of the JAGS 2.2.0 user manual with the results of `?rnorm` etc. and perhaps a few probability texts. This approach appears to require that the transformations will need to be deduced from the pdfs separately. I would prefer to avoid this task (and possible errors) if it has already been done, or else start the list here. **Update** Based on Ben's suggestions, I have written the following function to transform a dataframe of parameters from R to BUGS parameterizations.               ##' convert R parameterizations to BUGS paramaterizations     ##'      ##' R and BUGS have different parameterizations for some distributions.      ##' This function transforms the distributions from R defaults to BUGS      ##' defaults. BUGS is an implementation of the BUGS language, and these      ##' transformations are expected to work for bugs.     ##' @param priors data.frame with colnames c('distn', 'parama', 'paramb')     ##' @return priors with jags parameterizations     ##' @author David LeBauer          r2bugs.distributions <- function(priors) {            norm   <- priors$distn %in% 'norm'       lnorm  <- priors$distn %in% 'lnorm'       weib   <- priors$distn %in% 'weibull'       bin    <- priors$distn %in% 'binom'            ## Convert sd to precision for norm & lnorm       priors$paramb[norm | lnorm] <-  1/priors$paramb[norm | lnorm]^2       ## Convert R parameter b to JAGS parameter lambda by l = (1/b)^a       priors$paramb[weib] <-   1 / priors$paramb[weib]^priors$parama[weib]       ## Reverse parameter order for binomial       priors[bin, c('parama', 'paramb')] <-  priors[bin, c('parama', 'paramb')]            ## Translate distribution names       priors$distn <- gsub('weibull', 'weib',                            gsub('binom', 'bin',                                 gsub('chisq', 'chisqr',                                      gsub('nbinom', 'negbin',                                           as.vector(priors$distn)))))       return(priors)     }          ##' @examples     ##' priors <- data.frame(distn = c('weibull', 'lnorm', 'norm', 'gamma'),     ##'                     parama = c(1, 1, 1, 1),     ##'                     paramb = c(2, 2, 2, 2))     ##' r2bugs.distributions(priors)


title: Quantile regression power analysis

text: New to the site and to stats here! This may be a silly question, but I haven't been able to find a satisfactory answer on the procedure for a power analysis (or general guidelines about sample size) with a quantile regression. I plan on looking at 3 predictors along a single outcome. What would a recommended sample size be, or how can I conduct a power analysis? I appreciate any help, thanks!


title: Repeated coin flip experiments: what counts as a sample?

text: Apologies if this is a little confusing. If anything needs clarifying, I'll edit the question as best I can. Thanks. * * * If we take a sample of, say, 100 people, we could measure their heights, and find the mean and standard deviation. We can calculate an estimate of the standard error from this single sample to find where the population mean is likely to lie. If we repeated this sampling and measured further samples of 100 people from the population, then finding the means for each sample would give a sampling distributions of means. The mean of this sampling distribution estimates the population mean, while the standard deviation of this sampling distribution also estimates the standard error. This all seems fairly clear. * * * Now, consider a large number of trials, e.g. 20000 trials of 100 coin flips, where the number of heads in each trial is recorded and the fraction of heads obtained in each trial found. To me, it seems less obvious as to what constitutes a sample here. Would each one of these 20000 trials be considered as a single sample, so that we have 20000 samples? **If that's right, then when we average the fractions of heads found over all 20000 trials, then — _just as for the height example_ — the mean is an estimate of the population mean, and the standard deviation would be an estimate of the standard error.** The alternative I'm thinking of would be to treat each set of 100 coin flips as a single data point, which would give a single large sample of 20000 measurements. Treating this data as 20000 individual samples does seem like it should be correct since it's analogous to the height example. However, it doesn't seem immediately intuitive to me. Maybe this is because each trial results in a single fraction value, instead of dealing with the 100 individual flips of each trial. So, thinking of the result of each trial as a individual sample mean seems less obvious here, as the trial effectively yields a single value. On the other hand, I can see that the fraction of heads _is_ effectively a mean value, where we code heads as 1, tails as 0 and then divide by the total number of flips in a trial. What confuses me further is if we then decide to work with the raw number of heads found in each trial instead of the fraction of heads. The total number of heads in a trial is a single value, rather than a mean. However, my thinking is that the same would be true: the mean of the 20000 raw numbers of heads would be an estimate of the population mean, and that the standard deviation would estimate the standard error of the mean. I'm not sure if or why this would be the case though. * * * Edit: Really, I think what I'm trying to understand about this:   1. If I take repeated trials of a fixed number of flips and compile the fraction of heads, do these give a distribution of sample means? I think from Michael Chernick's answer that this is correct.   2. Next, if we look at the distribution of the raw scores of heads obtained in the repeated trials of a fixed number of flips, can this be again considered as a distribution of sample means? This is exactly the same as in 1, except the measured values are multiplied by the sample size. Because of this, I assume the answer is yes, but it just seems a little dissonant. When you reframe the data like that, you're using a total value, rather than the mean number of heads (the mean being equivalent to the fraction of heads).


title: Vector autoregression - number of lags

text: I am constructing a Vector autoregression model and I have used AIC to find how many lags I should use. Does 7 lags seem unreasonable? I am trying to find the impact the property market has had on the economy and my variables are GDP, inflation, employment rates, housing prices, number of housing complete, long-term interest rates and mortgage rates. How many lags makes economic sense?


title: Identifying the time lag between cause and effect

text: What approaches exist to observe the time lag between two variables? I need to analyze the relationship between blood pressure and some other factor, such as exercise. The data set I am drawing from has around 1800 individuals, with an average of 100 entries a piece. It is generally known that there is a strong relationship between exercise level and blood pressure. However, if a person increases their steps to 8000+ a day, how long will it take for their blood pressure to drop? I am new to this type of analysis, and this is a challenge I have been thinking about for weeks. I don't know if anyone wants to comment on possible approaches to this challenge or any issues surrounding it. Some issues I have been dealing with:   1. Is it better to treat this as a times series analysis or longitudinal data analysis? My understanding is that time series usually focuses on one variable with no missing data and is observed at consistent intervals, where as longitudinal is over a longer period and has inconsistent time intervals, dropouts, and missing data. The data I have seems to fit the longitudinal description more, but it also seems like time series could be used if I averaged the values by week so there would be no missing entries. I'm not sure about the pros and cons of each approach.   2. Should I be fitting a causal model, or would some other method like regression be more helpful? I've been looking at various possible causal models, for example Marginal Structural Models (MSM) or Structural Nested Models (SNM), but there seem to be very little information on their application. I did find one R package that applied inverse probability weights and then used Cox proportional hazards regression model on a survival object (MSM), but that seemed to be focus on weighting for confounding and right censoring. Its result was a correlation coefficient, which I don't think helps me. So I'm not sure if fitting a causal model is what I want, because that seems to be more focused on the making intellectually satisfying assumptions about relationships within the data and then determining the degree of causality, rather than providing information about time lag. If anyone knows about MSM, SNM, their use in R, or how they might relate to this problem, that would be awesome to hear.   3. What about survival analysis or SEM? I haven't explored these options very in-depth yet but they sound potentially relevant. I've kind of stalled, so any hints about what direction I might want to go would be really appreciated. Thanks in advance.


title: Best regression model, given coefficient of variation $R^2$ and mean-squared error

text: If you have 3 separate models in a multiple regression problem (and 3 ANOVA tables), which would be best given that you have the coefficient of determination, $R^2$, and mean-squared error values? So you have 1 table, with just $X_1$, another with just $X_2$ and a third with $X_1$, $X_2$ combined. Which is best given specific $R^2$ values and MSE values? I assume you're looking for the highest $R^2$ values, but how does MSE play into this?


title: Statistical measure for if an image consists of spatially connected separate regions

text: Consider these two grayscale images: ![river](http://i.stack.imgur.com/2ndLs.png) ![random](http://i.stack.imgur.com/3nWSz.png) The first image shows a meandering river pattern. The second image shows random noise. **I am looking for a statistical measure that I can use to determine if it is likely that an image shows a river pattern.** The river image has two areas: river = high value and everywhere else = low value. The result is that the histogram is bimodal: ![enter image description here](http://i.stack.imgur.com/3L9Lv.png) Therefore an image with a river pattern should have a high variance. However so does the random image above:               River_var = 0.0269, Random_var = 0.0310      On the other hand the random image has low spatial continuity, whereas the river image has high spatial continuity, which is clearly shown in the experimental variogram: ![enter image description here](http://i.stack.imgur.com/7Ytbd.png) In the same way that the variance "summarizes" the histogram in one number, I am looking for a measure of spatial contiuity that "summarizes" the experimental variogram. I want this measure to "punish" high semivariance at small lags harder than at large lags, so I have come up with: $\ svar = \sum_{h=1}^n \gamma(h)/h^2 $ If I only add up from lag = 1 to 15 I get:               River_svar = 0.0228, Random_svar = 0.0488      I think that a river image should have high variance, but low spatial variance so I introduce a variance ratio: $\ ratio = var/svar $ The result is:               River_ratio = 1.1816, Random_ratio = 0.6337      My idea is to use this ratio as a decision criteria for if an image is a river image or not; high ratio (e.g. > 1) = river. **Any ideas on how I can improve things?** Thanks in advance for any answers! EDIT: Following the advice of whuber and Gschneider here are the Morans I of the two images calculated with a 15x15 inverse distance weight matrix using Felix Hebeler's Matlab function: ![River_M](http://i.stack.imgur.com/FJFQl.png) ![Random_M](http://i.stack.imgur.com/8qeOd.png) I need to summarize the results into one number for each image. According to wikipedia: "Values range from −1 (indicating perfect dispersion) to +1 (perfect correlation). A zero value indicates a random spatial pattern." If I sum up the square of the Morans I for all pixels I get:               River_sumSqM = 654.9283, Random_sumSqM = 50.0785       There is a huge difference here so Morans I seem to be a very good measure of spatial continuity :-). And here is a histogram of this value for 20 000 permutations of the river image: ![histogram of permuations](http://i.stack.imgur.com/A7NcR.png) Clearly the River_sumSqM value (654.9283) is unlikely and the River image is therefore not spatially random.


title: Weird denominator degrees of freedom in SAS glimmix

text: I´m trying to learn generalized mixed models in SAS and I have now bumbed into a situation that I cannot figure out by myself. So I was thinking that maybe some of you could help me forward. I´m running a model where I have normally distributed measurements of traits like tarsus etc in birds from different populations. I have few fixed factors (sampling area, sex and year) and their interactions in the model and random term "sampling site clustered within the sampling area" (2 sampling areas with both 2 sampling sites). I´m using SATTERTH as a method for computing the denominator degrees of freedom. My problem is that for one particular trait I get really low degrees of freedom always when the sampling site is in the model, and for that particular factor. I noticed that playing around with the methods of computing the degrees of freedom, I can get "better" results, but I cannot understand why in this model, one fixed factor eats up my degrees of freedom. If anyone can and would like to explain to me why this happens and how I should select the method for computing the DDF´s in "dummies" way, it would make my day. I'm happy to provide you with more information, I was not certain how I should present my problem. Thank you in advance.


title: Estimate the parameters of beta exponential distribution via L-Moments

text: Estimate 3 parameters of beta exponential distribution in the case of censored type 1 samples via L-moments


title: Hyperparameter estimation in Gaussian process

text: I am trying to optimize the hyperparameters for a Gaussian process. As a starter I choose the squared exponential function for covariance where iI have to optimize 3 parameters $\sigma_f$, $\sigma_n$ and the length parameter $l$. $$k_y(x_p,x_q) = \sigma^2_f \exp\left(-\frac{1}{2l^2}(x_p-x_q)^2\right) + \sigma^2_n\delta_{pq}$$ As described in [1] I try to maximize the marginal likelihood using the parameters. To achieve this I use Commons Math optimizers and compute the gradients of the marginal likelihood function. Everything works fine except for a tiny example if I only use $ \sigma_f $ and $ l$ for the optimization. As soon as I try to optimize $ \sigma_n$ as well, I run into numerical problems and the the covariance matrix of the process is NOT positive-definite anymore after a few iterations of the optimizer and thus produces `INFINITY` or `NAN` errors within the optimizer. Can anyone explain that behavior and the connection to the hyperparameter $ \sigma_n$? [1] Gaussian Processes for Machine Learning Carl Edward Rasmussen and Christopher K. I. Williams The MIT Press, 2006. ISBN 0-262-18253-X.


title: Structural Equation Models (SEMs) versus Bayesian Networks (BNs)

text: The terminology here is a mess. "Structural equation" is about as vague as "architectural bridge" and "Bayesian network" is not intrinsically Bayesian. Even better, God-of-causality Judea Pearl says that the two schools of models are almost identical. So, what are the important differences? (Amazing to me, the Wikipedia page for SEMs doesn't even include the word "network" as of this writing.)


title: R correlation between two time series analysis

text: This is my very first post on CV so comment if I can improve my post. I have two webs that sell very similar products and also have a very similar group of customers. I was trying to prove that concept that "Fast Moving Inventory From Web A Will Also Move Fast on Web B" I have the daily snapshots of the inventory of both WebA and WebB. For example, for the product A from Web A, data looks like this:               TimeStamp  Product  Stock  UnitPrice     Oct 1st    A        100    1.2     Oct 2nd    A        90     1.2     Oct 3rd    A        40     1.2     Oct 4th    A        240    1.2     ..      For those two Webs, some of their products are exactly the same. For example, WebA and WebB both have product A. I am wondering is there a way to use R to do some time series analysis so I can prove that for those products that they have in common. There is a very high correlation (WebA sell product A crazy last weekend, which is also a good sale reflected on Web B's data)? Is there some R functions to do what I want? Then I can tell Web B that you need to carry those top selllers of Web A but Web B doesn't carry now.


title: Censored data when censoring time is uneven

text: I have a data set relating to credit defaults. The data set contains around 100 predictors (some categorical and some numeric). I am interested in predicting the time to default, if the person ever defaults. The term of the loan is 36 months, so a person could default at months 1-36, or never. The problem is for one segment of the data, I have 36 months of information, another segment of the data I have 24 months of information and another segment I have only 12 months of information. I realize that this is a censored data set, but my specific questions are:   1. How does the fact that the time of censoring is different for different segments of the data (12 months, 24 months, 36 months) affect the analysis? In the examples of censored data analysis I have read, the time of censoring is always the same. For example, a medical study that ended at month X, so all the data is right censored at this time.    2. If I am solely interested in maximizing my predictive power (according to some metric; I am aware of the dangers of data dredging) what are the best techniques for evaluating a dataset like this? Thanks for the help!


title: Why is this not a geometric distribution?

text: > > A recruiting firm finds that 20% of the applications are fluent in both > English and French. Applicants are selected randomly from a pool and > interviewed sequentially. Find the probability that at least five applicants > are interviewed before finding the first applicant who is fluent in both > English and French The answer for this question is $(0.8)^5$. However, since the questions asks for the probability that at least 5 applicants are interviewed _before_ finding the first applicant, wouldn't this be a geometric distribution, so $(0.8)^5 (0.2)$? Any help would be appreciated, thanks.


title: In an experiment w/binomial responses, some subjects gave the same answer for all trials. [How] can a Mixed Effects model (R's lmer) deal with this?

text: I recently ran a pilot of an experiment on Amazon's Mechanical Turk. In the experiment, participants read 5 items, and answered a yes/no question about each one. A between-participants factor was manipulated, which I'll call "Intention." After completing the pilot of the experiment, I attempted to analyze the results in R using the 'lme4' package.               lmer_max    = glmer(formula= Resp ~ 1+Intention + (1|PartCode)+ (1|Item), data=df, family="binomial")     lmer_nomain = glmer(formula= Resp ~ 1           + (1|PartCode)+ (1|Item), data=df, family="binomial")      Unfortunately, lme4 seems to dislike the data: using an older version of the package, the model fails to converge, and using a more recent version of the package, the model fit is worse than expected. Looking at the results, it's apparent that many subject's responded "yes" across all 5 vignettes. There's nothing wrong with them doing this within the context of the experiment, but I suspect it's one major source of problems in analyzing the data. One piece of confirmation for this possibility is plotting the lmer model using the `effects` package (Resp=0: No ; Resp=1: Yes): ![a ceiling effect](http://i.imgur.com/srhI2zY.png) This ceiling effect is unique to how the model is dealing with the random effects: the actual split between Yes's and No's is about Yes=2/3, No=1/3. My question is basically: how should I deal with this? I have three main options in mind, but am open to others.   1. **Re-design experiment so that there's more variation within each subject.** Obviously this is ideal, but it's also the riskiest and most difficult to do: the yes/no question participants are being asked is sort of philosophical in nature, and people come in with strong pre-conceptions. The manipulation in question seems to push their responses around, but there's still a lot of across-subject variability, and ultimately it may be difficult to guarantee that no subjects gives the same response on all trials.   2. **Re-design experiment so that there's FORCED variation within each subject.** Currently, subjects are asked this yes/no question for five items, and yes or no is a perfectly valid response for each. However, I could add a "catch" item in which "no" is the only plausible response, then remove any subjects who answered "yes" on all 6 items. This would guarantee that subjects, at maximum, answered 'yes' only 5 out of 6 times. However, this feels like a slightly odd strategy: instead of having some participants who gave the same response across all items, there will now be an item that has the same response across all participants. Is this a problem? I suppose the difference would be that, right now, about half of subjects had the same response on every item, whereas under this strategy, only 1/6 items would have the same response for every participant.   3. **Re-do analysis so that there's forced variation within each subject.** This is the laziest possibility, but I thought I'd ask. Is it in any way kosher to introduce a dummy item/trial into the data frame, in which all participants responded "no," simply to remove the difficulties that lmer is having? That is, instead of doing possibility 2, I would simulate possibility 2 with the current data. Or perhaps there's a way of telling lmer not to freak out over subjects who gave the same response across items?  Any help would be greatly appreciated.


title: Parameter estimation for the sum of two Independent (not necessarily i.d.) Gamma RVs

text: I'm a bit of a stats newbie so take it easy on me if this ends up being somehow trivial. I'm working on a problem that involves parameter estimation for the sum of two independent gamma distributions (which do not necessarily have the same scale or shape parameters.) I'm interested in what would be the most convenient way to numerically estimate the parameters of the model. I'm primarily interested in approximating MLEs (though I would be interested in hearing alternative methods.) The convolution expression for the density seems a bit messy, and while I found a paper giving an expression for the density in terms of an infinite series with a nice truncation error estimate, the objective function is still pretty messy. Particularly, I'm curious to know where to start looking in terms of a) asymptotic approximation for the MLE and b) software that is well-suited for this type of application. Any links or references that you find useful will be most appreciated. Thanks!


title: How to test whether the trend line between two time periods significantly differ?

text: I have chemical A concentration measured at various locations from 1993-2010. I have broken the time into 2 periods 1993-2004 and 2005-2010. **How can I test whether the trend line between the two periods significantly differ?** Someone told me to use the t-test. However, I thought that the t-test is for testing differences in means between 2 population and that it does not involve trend lines. But I might be wrong.


title: What is an appropriate classification method for analysis of down-hole geophysical data?

text: I am looking for a method that may be used to classify and 'cluster' some multi-dimensional data gathered from drill holes, to impute some (only partially measured) physical properties of the rocks within the drill holes. The data are (1) categorical qualitative: broadly classified lithology and stratigraphy, and (2) continuous or semi-continuous quantitative: where a set of instruments records a set of measurements of things like natural gamma response, electrical resistivity, acoustic velocity, and so on. A brief snippet of a set of data would be like this:               DEPTH   LITH    K   CALI    DT  GR  LAT     474.1173    12  1.290   240.1594    392.5612    54.2610 9.0883     474.2697    12  1.290   241.1248    389.7305    53.9297 9.2322     474.4221    12  1.290   242.9032    382.1031    53.5055 9.6901     474.5745    12  1.290   254.1457    388.4041    53.0907 9.7129     474.7269    12  1.290   248.7685    395.1500    53.8861 9.5200     474.8793    12  1.290   241.3833    393.8504    53.7665 9.8039     476.0985    12  1.194   236.8224    397.4059    55.1024 9.6624     476.2509    12  1.194   237.8029    397.2021    53.5209 9.5717      These geophysical well-logs have a response based on the physical properties of the rocks, ultimately on the mineralogy and physical state. I want to find the best set of well-log data that may be used to impute another physical property that is much harder and more expensive to obtain. For this problem it was suggested to me to use a SOM method implemented in some software. The primary benefit of the suggested implementation is its ease of use, and that it has been developed to cater to geoscientists. In my child-like naivete the SOM method seems a fine one to use. However I'm sure there are pros and cons I'm unaware of, and there might be better techniques. Its primary drawback seems to me to be that it is a little 'black- box'. Are there other methods suitable and appropriate that you can suggest?


title: Two sample t-test with mean centralized data

text: I have two sets of equal sized and normally distributed data, which are output from two different experiments. I wish to do a two-sample t-test with these two data sets to see if the results from the two experiments agree. **Is it appropriate to mean (or median) centralize the data sets before doing the t-test?** The variance of the data are equal and passes the var.test. My understanding is that if I centralize the data and the mean will be the same for the two data sets and the result of a t-test will always favor the null hypothesis of equal means. Can the same thing be said about centralizing data or not for a one-way ANOVA if I have more than two datasets?


title: Best method to validate a multiply imputed Cox model with R?

text: This question is with regards to using a test data set to validate an imputed Cox model using R. With a non-imputed data set I would use `val.surv()` from `rms`, but I'm not sure how/if I can use it with my multiply imputed data set. Further explanation: I created a predictive Cox PH model for 5-year RFS, and also used the `mice` package to multiply impute some missing data in the training data set. I then used `fit.mult.impute()` in Dr. Frank Harrell's excellent `Hmisc` package to obtain the pooled model. I have a data set that I would like to test the model on, but I am not sure how best to validate the pooled model. Multiple imputation is a common procedure that many researchers utilize, so there must be a way that R users are validating their models with imputed data. I would like to know what functions/avenues are available for me to test this pooled model with my validation data set? Here is some sample code to work with:               library(rms)     library(survival)     library(mice)          remove(veteran)     data(veteran)     veteran$trt=factor(veteran$trt,levels=c(1,2))     veteran$prior=factor(veteran$prior,levels=c(0,10))          #Set random data to NA      veteran[sample(137,4),1]=NA     veteran[sample(137,4),2]=NA     veteran[sample(137,4),7]=NA          impvet=mice(veteran)     survmod=with(veteran,Surv(time,status))          #make a CPH for each imputation     for(i in seq(5)){         assign(paste("mod_",i,sep=""),cph(survmod~celltype+karno,             data=complete(impvet,i),x=T,y=T))     }          #Now there is a CPH model for mod_1, mod_2, mod_3, mod_4, and mod_5.          pooled_mod=fit.mult.impute(survmod~celltype+karno,cph,impvet,data=veteran,surv=T)          #Here is a test data set.     remove(veteran)     test_dat=data.frame(trt=replicate(500,NA), celltype=replicate(500,NA), time=replicate(500,NA), status=replicate(500,NA), karno=replicate(500,NA), diagtime=replicate(500,NA), age=replicate(500,NA), prior=replicate(500,NA))     for(i in seq(8)){     test_dat[,i]=sample(veteran[,i],500,replace=T)     }          #Now there is a pooled model, "pooled_mod", and a test data set, "test_dat".      I'm looking forward to hearing about the R methods that can help in this situation.


title: How much information can you mine out of a name?

text: A name: first, possibly a middle, and surname. I'm curious about how much information you can mine out of a name, using publicly available datasets. I know that you can get the following with anywhere between a low-high probability (depending on the input) using US census data: 1) Gender. 2) Race. Facebook for instance, used exactly that to find out, with a decent level of accuracy, the racial distribution of users of their site (https://www.facebook.com/note.php?note_id=205925658858). What else can be mined? I'm not looking for anything specific, this is a very open-ended question to assuage my curiousity. My examples are US specific, so we'll assume that the name is the name of someone located in the US; but, if someone knows of publicly available datasets for other countries, I'm more than open to them too. I'm not quite sure if this is the right place for this, if it's not, I'd appreciate if someone could point me to a more appropriate place. I hope this is an interesting question, and this is the appropriate place!


title: Number of tried for a favorable result

text: I am new to statistics so I apologize for the simplicity of my question. I have a bag with 27 numbers inside; I will pull one number out until I get number 17. THe numbers will be returned to the bag after every try so the prob of success will always be 1/27. I am trying to determine on average how many pulls I will need.


title: How to find a meaningul threshold for relevant elements?

text: Is there any basic rule on which elements are over a certain threshold (and thus being important for further inspection) and which are below? Or can I set a "random" threshold just as I find it suitable for my current project? _Case example: I have a list of 94 different elements. The most frequent element occurs 2497 times, the least frequent elements just 1 time. Overall there are 19518 occurences. The elements have the following distribution (if helpful):_ ![distribution of elements](http://i.stack.imgur.com/RBI2I.jpg)


title: Applications of principal component analysis versus principal component regression?

text: I am investigating many different kinds of PCA versions, I am trying to find out whether PCR will apply to my analysis thus the question on use of PCR.


title: Testing difference between two means with pairwise data and absence of normality

text: I am trying to compare the difference between two means with two pairwise samples. Unfortunately, my data are very far of being normal. What test would you recommend to use in this situation? Should I revert to a nonparametric test?


title: How to simulate data based on a linear mixed model fit object in R?

text: Assume, I have a data set, which is similar to               require(nlme)      ?Orthodont       and my model is               fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)       **How can I use the model fit object`fm2` to generate several datasets, which have sample sizes 300, 400, 500, ... ?** I read this great answer on r-sig-mixed-models help but it seems incomplete.


title: How much is moderate violation to normality for one sample t-test?

text: I want to test if my data have a mean equal to zero (Ho: mu=0). The sample size in not large (n=21). My variable is numerical with average -0.10 and sample_sd 0.05. I don't know the population variance. The only test I know to do this, is t-test. I checked for normality using Shapiro Wilks and it gets rejected (therefore no evidence of normality). I understand that t is relatively robust to moderate violations of the normality assumption; what is not clear to me is, how much is moderate? I do need to run other test like this for other samples. May someone be so kind to provide me with examples of how much can I deviate in skewness and kurtosis for sample sizes between 12 and 25? What else can I do to test against zero if I can't use t-test? Thank you in advance ![enter image description here](http://i.stack.imgur.com/KBgTg.gif)


title: QQ plot of sigma2 from an MCMC regression?

text: I've been using R's `lm` to do some linear regression, but decided to give `MCMCregress` a try to get a feel for how it works. As expected, I got basically the same coefficients, but the extra `sigma2` value puzzles me. When I do a `qqmath` plot of the coefficients, I get the following graph, and I'm puzzled by the sigma2 plot. It's obviously not linear, but I'm not sure if that's meaningful in this context. I assume it's sigma _squared_ , and when I took the square root and plotted it, the line was straighter, but still curved. I guess my question boils down to: what is sigma2 telling me about the MCMC regression fit, and is a graph of it useful or should I ignore the graph and focus on something else? (All of the diagnostics and graphs I've done on my original `lm` fit seem to indicate that the fit is good, so I'm also wondering if the MCMC regression gives me more information or not.) ![qqmath plot of MCMCregress results](http://i.stack.imgur.com/wayi4.png) (If I need to provide the actual data, I can. I'm hoping that an answer depends more on what sigma2 is rather than on specific values.)


title: Multivariate logistic distribution

text: The normal distribution can be generalized into the multivariate normal distribution. Can the logistic distribution also be generalized into a similar multivariate distribution? Is there a multivariate generalization of the logistic distribution which depends on the covariance matrix $\Sigma$, similar to the multivariate normal distribution? The multivariate distribution should be such that its marginals are univariate logistic distributions.


title: Test Significance Between Two Groups (not mean?)

text: I'm conducting my bachelor thesis and I need to do research. This morning I received mail from my supervisor and I'm not quite getting wat he wants me to do. I hope someone can help me. PS: I will try to keep it as simple as possible. ![enter image description here](http://i.stack.imgur.com/Y0Sru.png) My supervisor send me the following: "The other is to compare the local versus foreign analyst recommendation statistically. For example, for each stock, calculate the difference between CAR_local analyst and CAR_adr analyst. Then report whether the differences are significant" My expectation is that the Foreign will be significantly higher than the local. What do you suspect I need to perform for kind of test. PS: I have SPSS and Eviews on my pc. Every help is much appreciated! EDIT: My N will be around 300 (so 300 stocks)


title: Maximum Likelihood estimation with eps ~ t5 while eps ~ N(0,1)

text: If one estimates an ML model with eps following a t5 distribution, while the eps are actually standard normally distributed, intuitively I would say the estimated standard errors would be too large. Is my intuition correct or am I missing something? The reason for my question is that if standard errors are systematically too large and you still find significant results it would be defendable to say that the results can actually only be "more" significant.


title: Time series degree of slope: calculating what I see

text: I want to calculate the degree of slope at each point in a time series. Different time series have different scales. The final number should be normalized in the range of +/-90 degrees. Basically, when I plot my time series in Excel, I can see the degree of slope up or down, 0=flat, 70=very steep up, -20=gradual slope down. I want to calculate the "number" for what I am seeing. I thought using the arctangent(P-P1), P=current point, P1=previous point would work. Not at all. For example on one time series: atan(1.166031374-1.168266667) yields -0.00224. On another times series, atan(11373.92-11342.05) = 1.539431. Certainly not normalized across different value scales nor producing values between +/-90. Visually, it’s so easy to see the degree of slope in my chart! Yet, over the last year I’ve tried more than a hundred work arounds, mostly complex. They approximate what I want but seem very convoluted and inelegant. I’d appreciate any insights into solving this problem.


title: Single variable changing over time - compare trajectories

text: I am running an experiment in the lab where a variable is measured over time (e.g: 60 seconds) I run several times the experiment without any disturbance. Then I run the experiment again by introducing a distrubance in the starting conditions. The trajectory of the variable I am measuring will start in a similar way and then will behave slightly different. My goal is to be able to identify that a disturbance was made just by looking at the trajectory. What's the best method to make this comparison? I want to be able to say "this experiment is significantly different from the non disturbed one, so something must be wrong". I was thinking of using some mvda techniques, such as PCA, and just use each time point as a variable, and then i can see clusters of observations. Is this a good approach?


title: Bias of Panel Generalization of Durbin-Watson

text: I'm working with an unbalanced panel dataset. (Country-Time) of approximate dimensions H=100 individuals i and average time length over individuals $mean(T_i)\approx7.5$. And about n= 8 regressors (appart from the fixed effects dummies). I want to use the Durbin-Watson panel generalization from _Bhargava, Alok, Luisa Franzini, and Wiji Narendranathan. "Serial correlation and the fixed effects model." The Review of Economic Studies 49.4 (1982) 533-549_ (http://restud.oxfordjournals.org/content/49/4/533.short) $d_P = \frac{ \sum_{i=1}^H \sum_{t=2}^T \left( \tilde{u}_{it} - \tilde{u}_{it-1} \right)^2 }{ \sum_{i=1}^H \sum_{t=1}^T \tilde{u}_{it}^2 }$ Critical values are given in the paper and are dependent on (T, H, n) with n the number of regressors. Which i linearly interpolate. $cv(7.5, 100, 8) \to d_{PL}=1.8561,\; d_{PU}=1.9039$ Now if i simulate data _without_ autocorrelation I tend to find very low values for $d_P$ that reject $H0: \rho=0$               #[R]     require(data.table)     set.seed(1)     DT <- data.table(i=c(rep(1:50, each=7), rep(51:100, each=8)),                       t=c(rep(1:7, 50), rep(1:8, 50)), u=rnorm(100*7.5))     DT[, ':='(du2=c(NA, diff(u))^2, u2=u^2), by=i] # difference by individual          #       i t          u          du2         u2     #  1:   1 1 -0.6264538           NA 0.39244438     #  2:   1 2  0.1836433 0.6562573681 0.03372487     #  3:   1 3 -0.8356286 1.0389152808 0.69827518     #  4:   1 4  1.5952808 5.9093205817 2.54492084     # ---                                              #746: 100 4 -1.3457937 0.1434371834 1.81116064     #747: 100 5  1.0336654 5.6618254864 1.06846414     #748: 100 6 -0.8117765 3.4056556180 0.65898102     #749: 100 7  1.8017255 6.8303923814 3.24621470     #750: 100 8  1.7715420 0.0009110448 3.13836092          DT.sumT <- DT[, list(sumdu2=sum(du2, na.rm=T), sumu2=sum(u2)), by=i] #sum over T_i     #       i    sumdu2     sumu2     #  1:   1 12.239715  4.688696     #  2:   2  8.925133  8.698205     #  3:   3  2.463443  4.030220     # ---          d_P <- DT.sumT[, sum(sumdu2)/sum(sumu2)] #sum over i     d_P # Durbin-Watson Statistic     # 1.708489     # next in seed gives --> 1.617848, 1.735762, 1.614137      Which are all < 2 so positive autocorrelation and all simulated $d_P< d_{PL}$ so 5% significant. But the data is IID... If a negative correlation of 0.2 is forced the test does flip side but still looks downward biased:               DT[, u:=u-0.2*c(0, u[1:(length(u)-1)])]     # d_P --> 2.038816      Am I missing something here? The tabulated 5%c.v. are for T=6 and T=10 so the short+wide property should not be a problem. Also balancing the simulation to e.g. $T_i=T=7$ gives low $d_P$


title: Calculating power for a survey to compare different metrics

text: I am currently planning a survey. The data will be used to test which of three metrics best captures some "true" underlying characteristic of the respondents. Getting data on the "true" underlying characteristic is costly, so I want a simplified metric that best approximates the "true" measure. For the current survey, I am gathering four data points for each observation. 1) metric 1 2) metric 2 3) metric 3 4) the "true" underlying characteristic In the future, I won't have data on the "true" characteristic, so I will rely on one metric. For future work we will collect just the data for the "best" metric. I am not clear about what effect size (aka difference in prevalence of indicator) to use to calculate sample size. Specifically, the hypothesized effect that I seek to identify, I believe is the difference between the metrics in their correlation to the true measure. For instance, if metric 1 has correlation 0.7 with the true characteristic, and metric 2 has a correlation 0.9, then the "effect size" of switching from metric 1 to metric 2 is 0.2. Am I on track?


title: Standardizing before/after/at all when using multi-class LDA for pre-processing step

text: If a multi-class Linear Discriminant Analysis (or I also read Multiple Discriminant Analysis sometimes) is used for dimensionality reduction (or transformation after dimensionality reduction via PCA), I understand that in general a "Z-score normalization" (or standardization) of features won't be necessary, even if they are measured on completely different scales, correct? Since LDA contains a term similar to the Mahalanobis distance which is already implying normalized Euclidean distances? So it would not only be not necessary, but the results after an LDA on standardized and non-standardized features should be exactly the same!?


title: How do we know the LMM model has to add random slope with random intercept?

text: When I read the literature about random intercept, random slope, and random intercept and slope model, I usually see the graph that they have lines for each participant. I do not know how to plot such graph and what this graph is called. Can anyone give me a link? or could you tell me what function in R can do this kind of graph to test if my model needs only random intercept, random slope, or both?


title: R - Test for homogeneity of regression slopes results in singular model

text: I am trying to check the assumptions of a two-way ANCOVA. So in my model I have   * two factors (F1, F2)   * one dummy coded two level covariate (C)   * one dependent variable (D) In order to check the assumption of homogeneity of regression slopes I tried to perform an ANOVA with type 3 sums for the model D ~ F1*F2*C to see whether any interactions with the covariate might be significant. Using the Anova function from the car package this corresponds to               modd<-aov(D~F1*F2*C)     Anova(modd,type=3)      However, I encounter the following Error message:               Error in Anova.III.lm(mod, error, singular.ok = singular.ok, ...) :       there are aliased coefficients in the model      My question is, whether it makes sense for the homogeneity test to force R to compute the ANOVA anyway by supplying the singular.ok=T option or what else I should do in order to check the assumption.


title: Sample sizes for differences between three groups

text: I am doing a study of two schools, and at each school I am sampling three groups. I am trying to determine why the person chose to go to that school and not another. I will be asking each group various questions about the school and options they may have had for other schools, etc. I want to be able to detect:   1. a difference of 20% in responses between the schools (e.g., proportion of students for whom current school was first school of preference, or proportion for whom there was no other choice), and   2. _within_ each school, 20% difference in responses between the three different groups (similar questions, but for students who enrolled in different eras). In both instances, I would like power of 80% for 0.05 significance. So, if I do what I think is the appropriate calculation, I effectively come up with: \begin{align} n &= \frac{0.5(0.84 + 1.96)^2}{0.2^2} \\ &= 98 \end{align} Can I just assume then that I need three groups of 99 within each school? I guess I'm confused because nobody ever seems to talk about calculating sample sizes when comparing more than two groups. Furthermore, is there anything wrong with sampling the three groups of 98 at each school and assuming that the total sample of 294 at each school will be sufficient to detect the 20% difference between the two schools?


title: Confusion between caret randomForest predict() results and reported model performance

text: This question seems related, but the consensus was that the issue had to do scaling the data, which I do prior to training, so I don't think that's the issue:   * Issue on prediction with FinalModel of RandomForest in R using the CARET package I've uploaded a sample data set, and here is how I generated my model:               library(randomForest)     library(caret)     library(ggplot2)          data <- read.csv("http://pastebin.com/raw.php?i=mE5JL1dm")          data_pred <- data[, 1:(ncol(data) - 1)]     data_resp <- as.factor(data$y)          data_trans <- preProcess(data_pred, method = c("center", "scale"))     data_pred_scale <- predict(data_trans, data_pred)          trControl <- trainControl(method = "LGOCV", p = 0.9, savePredictions = T)          set.seed(123)     model <- train(x = data_pred_scale, y = data_resp,                    method = "rf", scale = F,                    trControl = trControl)      Here's what `caret()` reports as the model performance:               > model          Random Forest           516 samples      11 predictors       5 classes: '0', '0.5', '1', '1.5', '2'           No pre-processing     Resampling: Repeated Train/Test Splits Estimated (25 reps, 0.9%)           Summary of sample sizes: 468, 468, 468, 468, 468, 468, ...           Resampling results across tuning parameters:            mtry  Accuracy  Kappa  Accuracy SD  Kappa SD       2     0.747     0.663  0.0643       0.0853         6     0.76      0.68   0.0507       0.068          11    0.758     0.678  0.0574       0.0763            Accuracy was used to select the optimal model using  the largest value.     The final value used for the model was mtry = 6.       In my "real" model, I have a training/hold-out set, and creating two sets of plots showing model predictions for the training/hold-out sets vs. the corresponding true observations. That's when I noticed something that seemed odd to me.               # data set of model predictions on training data vs. actual observations     results <- data.frame(pred = predict(model, data_pred_scale),                           obs = data_resp)          table(results)               obs     pred    0 0.5   1 1.5   2       0   148   0   0   0   0       0.5   0 132   0   0   0       1     0   0 139   0   0       1.5   0   0   0  38   0       2     0   0   0   0  59      And here's a plot confirming 100% accuracy on the training set:               p <- ggplot(results, aes(x = pred, y = obs))     p <- p + geom_jitter(position = position_jitter(width = 0.25, height = 0.25))     p      ![training predictions](http://i.stack.imgur.com/aFtX4.png) But if I look at the saved tuning predictions, subsetting only those where `mtry = 6` (what `caret` reports as the final model), I don't get anywhere near that performance:               model_resamples <- model$pred[model$pred$mtry == 6, c("pred", "obs")          table(model_resamples)                  0 0.5   1 1.5   2       0   296  69   5   0   0       0.5  51 228  48   0   0       1     3  28 255  24   9       1.5   0   0  16  32  15       2     0   0   1  19 101      And the same sort of plot:               p <- ggplot(model_resamples, aes(x = pred, y = obs))     p <- p + geom_jitter(position = position_jitter(width = 0.25, height = 0.25))     p      ![resampled predictions](http://i.stack.imgur.com/Z4gbm.png) Is this just a case of over-fitting where holding out 10% of the data can create a 25% decrease in performance, yet the final model trained with the same parameters but all rows can yield 100% performance? It seemed unlikely, but that's the only thing coming to mind at the moment. I just want to make sure there's nothing wrong in my training or predicting methods where I'm creating a problem where there shouldn't be one. * * * _Note:_ I created the tables/plots prior to adding `set.seed()` in the model training code above. The exact table and plot may differ slightly, but in re- running, they general result is the same (perfect re-prediction vs. ~77% reported by the model). It didn't seem to warrant re-doing the results/plots above, so I left them.


title: What is the difference between ex ante and a priori, if any?

text: In the context of estimating geophysical quantities from remotely sensed data (inverse theory), what do the terms _ex ante_ and _ex post_ mean? For context, see for example this paper by T. von Clarmann (emphasis mine): > The actual value of $χ^2$ allows to conclude if the differences > $\mathbf{\hat{x}}_\textrm{val} - \mathbf{\hat{x}}_\textrm{ref}$ are > consistent with the **ex ante** estimates of the uncertainties of the > difference, represented by its covariance matrix $\mathbf{S}_\textrm{diff}$ > , or if there is a significant inconsistency. and: > Instead, we use independent reference measurements, which allow to infer > **ex post** estimates of bias, precision, and accuracy. Is there any difference between the terms _ex ante_ and _a priori_?   Likewise, is there a difference between _ex post_ and _a posteriori_?   If yes, what are those differences?


title: Power Spectral Density

text: I'm currently writing an spectrogram to be able to plot Audio signals. So far, I have done the following:   * Split the data into blocks of `256` overlap of `128`   * Applied a Hanning window on each block    * DFT of each of the blocks    * Taken the first half of the vector (128) This has given me a 2D vector of size: `146x128` in order to check the results I have been getting, I used `matplotlib` in Python, the results were not the same and I spent weeks trying to find out why. It become obvious to me that `matplotlib` returns the Power Spectral Density of the signal, and, not the magnitude. I would like to therefore find the Power spectral density of the result. I've looked at the following: Here And to me, the algorithm looks like the following: I have 146 blocks containing 128 samples (complex):               block1 = [0, 5, 6, 8, ..... n] *where `n` is 128          block2 = [1, 4, 7, 8, ..... n]           .....     .....     .....,          block[146] = [10, 5, 6, 7]       If I therefore take the squared magnitude (average) of each these blocks: `average = sum(re*re+im*im)/128` This will give me a 1D vector, containing all of the averages for the blocks that I have in my 2D vector. I don't know where to go from here, basically. I've thought about, if I take the magnitude of each actual signal block `re*re+im*im` and then divided each of the numbers by the `average` Essentially, I'm using this to plot a spectrogram, so it will need to be in a 2D vector in order to do this. Any help would be greatly appreciated


title: Uncertainty in Genetic Algorithm output

text: I'm applying a simple genetic algorithm to an optimization problem (I need to find a 4-parameter function's global minimum) I start the generations run with say `100` chromosomes which I then evaluate, select, breed (crossover + mutate) and repeat for `N_gen` generations. When the algorithm is done (the `N_gen` generations were processed) I'm left with the best fitted solution: the one which evaluates to the minimum value from the last generation's `100` chromosomes/solutions. But this also leaves me with `99` other well fitted solutions that I have _a priori_ no use for. Is there a way to extract some information from these other solutions? Perhaps an uncertainty estimate for the best fitted solution? Intuitively I'd calculate the standard deviation for the parameters values obtained from the last generation and use that as an uncertainty estimate for the best fitted solution. Kind of like using those other solutions as a _bootstrap_ for my best solution.