title: Comparing two slopes

text: I have a dataset structured as follows:               Variable1  Variable2  Variable3  Variable4  Variable5  Variable6  Variable7  Variable8       60         70         75          73        57         69         85         90       55         66         47          79        63         76         82         97      etc. Each variable represents the detection performance for a certain position and condition. Variables 1 to 4 are meant to represent condition 1, Variables 5 to 8 are meant to represent condition 2. Detection performance for the same locations (but different conditions) are represented by variables 1 & 5, variables 2 & 6, variables 3 & 7, variables 4 & 8. I now have plotted the detection performance for both conditions as a line chart. What I want to do, is to compare the slope of both curves using SPSS. What would be the best way to accomplish this? I would like to test if one slope is steeper than the other. Both curves are representing the detection performance of a target at different locations. I have the directional hypothesis, that the slope of the one curve would be flatter than the other. I am aware that there are different slopes within the curve, so what I am looking for is the "overall slope" for each curve. I suppose this could be accomplished by calculating a regression line. However I am not aware about how to do this with my dataset. Usually predictor and predicted variable are positioned in two fields of your SPSS file. But now they are in the same field: the location of the target should predict the detection performance. And variable 1 to x represents the _detection performance_ on _location x_. Do I have to split those two pieces of information _detection performance_ and _location_ , or can this problem be solved with the existing structure of my dataset? This is how the plot looks like: ![enter image description here](http://i.stack.imgur.com/IlDvP.jpg) And I want to test if the curve of one condition is steeper / flatter than the curve of the other condition. Detection performance for condition 1 (4 possible target locations) is represented by variables 1 to 4. Detection performance for condition 2 (4 possible target locations) is represented by variables 5 to 8. The explanatory variables would be the position of the target, the response variable would be the detection performance. Both pieces of information are represented by the same field: Variable1 stands for _detection performance_ on _position a_ , Variable2 stands for _detection performance_ on _position b_ , and so on. My study design is a within subjects design, each participant was tested on both conditions. In the sample data I provided, each row stands for one participant. Each number stands for the average detection performance of the respective participant for the respective condition (condition 1: variable 1 - 4; condition 2: variable 5 - 8) and the respective location (the same location is represented in Variable 1 & 5, Variable 2 & 6, Variable 3 & 7, Variable 4 & 8).


title: How to compute lower edge of a credible interval for a beta distribution in python?

text: Suppose I flip a coin 584,369 times and observe 5,257 heads and 579,112 tails. I wish to know the lower edge of the 90% credible interval for the probability that the next flip will come up heads. I believe:   * I could use a beta distribution with alpha=5257 and beta=579112.   * I want the value of the cumulative distribution at 0.05.   * I could calculate that in Python with SciPy using the betainc and beta special functions. I installed SciPy, and wrote:               from scipy.special import betainc, beta          numer = betainc(5257, 579112, 0.05)     denom = beta(5257, 579112)          print numer     print denom     print numer / denom      I get               1.0     0.0     ./foo.py:6: RuntimeWarning: divide by zero encountered in double_scalars       print numer / denom     inf      That is not what I want. I believe this web page calculator says I want 0.008793824176376785502, which is plausible (close to 5257/(5257+579112)=0.00899602819451). So: are my beliefs correct? How do I calculate this lower edge in Python?


title: Bayesian inference with the wrong distribution

text: When an observation $x$ is generated by $P(x|\theta)$ for a parameter $\theta$ the Bayesian optimal estimator for the value of $\theta$ is $\hat\theta_{BEST}=\mathbb{E}[\theta|x]=\frac{1}{P(x)}\int d\theta P(\theta)P(x|\theta)$. Now assume that instead of using $P(x|\theta)$ in the above formula we have $Q(x|\theta)$ which is different from $P$ but have the same mean and variance for all $\theta$; that is $\hat\theta_{EST}=\frac{1}{Q(x)}\int d\theta P(\theta)Q(x|\theta)$ and $\mathbb{E}[Q(x|\theta)]=\mathbb{E}[P(x|\theta)]$ and $\mathbb{V}[Q(x|\theta)]=\mathbb{V}[P(x|\theta)]$ for all $\theta$. Are there any known boundaries on $MSE[\hat\theta_{EST}]-MSE[\hat\theta_{BEST}]$? What if we further assume Gaussianity of $Q$? Thanks!


title: Definition of dynamic Bayesian system, and its relation to HMM?

text: From Wikipedia > A Dynamic Bayesian Network (DBN) is a Bayesian Network which relates > variables to each other over adjacent time steps. This is often called a > Two-Timeslice BN because it says that **at any point in time T, the value of > a variable can be calculated from the internal regressors and the immediate > prior value (time T-1)**. DBNs are common in robotics, and have shown > potential for a wide range of data mining applications. For example, they > have been used in speech recognition, protein sequencing, and > bioinformatics. DBN have shown to produce equivalent solutions to Hidden > Markov Models and Kalman Filters.   1. I was wondering if "the immediate prior value (time T-1)" means the time index in a DBN is always discrete?   2. Does "at any point in time T, the value of a variable can be calculated from the internal regressors and the immediate prior value (time T-1)" mean a DBN is a discrete-time Markov process?   3. If I understand correctly, a HMM is a discrete-time Markov process too, if ignoring the output from state at the same time. So I wonder if HMM and DBN are the same concept? But another Wikipedia article says  > hidden Markov model (HMM) is a statistical Markov model in which the system > being modeled is assumed to be a Markov process with unobserved (hidden) > states. An HMM can be considered as the simplest dynamic Bayesian network. and there is another quote from the first article: > DBN have shown to produce equivalent solutions to Hidden Markov Models and > Kalman Filters. Thanks!


title: Handling forecast mean with complex seasonality

text: I'm using a time-series model to do weekly forecasts on the number of **incoming calls** to a company. This variable has a weekly 'in-month' pattern and a monthly 'in year' pattern, and i have data from 2008-01-01 to 2012-08-31. To explain it better: Mondays show peaks along weeks, and rainy months show peaks along the year. My idea is to take the **week-forecast point** and 'distribute' it along the week days, using the weekly-pattern information. What is a good technique to do this? Empirical density distributions maybe? My tools: R language. Edit: Below, weekly forecasts using 2 different methods. Results from `tbats` using the `forecast` package: ![enter image description here](http://i.stack.imgur.com/8ngG8.png) Results from `tslm` using the `forecast` package: ![enter image description here](http://i.stack.imgur.com/oN7JJ.png)


title: Formula for probability-probability plot

text: I am programmer developing a Software that features the display of a probability plot. As far as I know there are 2 different ways to display a probability plot (please correct me if I say something wrong):   * A Percent / Percent Chart   * A Quantile / Quantile Chart My customers want a Percent / Percent Chart. That brings some problems. There is no charting control I know that features a probability plot chart. ![enter image description here](http://i.stack.imgur.com/uxYS4.jpg) That's why I have to calculate the positions of the data points by hand. **Are there any formulas or calculations to place the data points and Y-axis values in a normal coordinate system?**


title: Extreme values for f-statistic

text: I have a question regarding my dat. It appears that i have extremely large values of the F statistics, ranging from 20 to 20,000. How is this possible? I should mention that it is panel data, controlled for fixed effects as well as clustering. thanx in advance, Luisa


title: How to assign degrees of freedom for two-way ANOVA with two within-subjects factors?

text: I have an experiment with 44 subjects and two within-subject factors, condition (3 levels) and round (6 levels). It is fully crossed such that each person completes one of each task condition in each round. I have done the analysis (and ANOVA) in R using the lmer function. My lmer model is:               lmer(dist ~ condition*round + (1|id) + (1|condition:id) + (1|round:id),data=data_dist)      For various reasons that are described around the web, lmer doesn't provide denominator degrees of freedom and thus it doesn't compute p-values. Because my design is simple I would like to assign the degrees of freedom just as I would for any other repeated-measures ANOVA. However, I have looked all over the web and can not find any resources for how to construct a repeated-measures ANOVA when both independent variables are within-subject. My current thinking for degrees of freedom is as follows: > condition: 2 (# of memory conditions - 1)   >  round: 5 (# of rounds - 1)   >  interaction: 10   >  subject error: 43 (# of subjects - 1)   >  error: 628   >  total: 688 (# of observations - 1) I think I also need degrees of freedom for the two interactions from the model (condition:id and round:id). Also after I calculate these degrees of freedom, I am unsure which one should be used in the denominator for my p-values. Can somebody please point me in the right direction? Thank you! **EDIT:** My views have evolved. Here is my current idea: > subj: 43   >  cond: 2   >  subj X cond: 86?   >  rnd: 5   >  subj X rnd: 215?   >  cond X rnd: 10   >  err: 688 - all of above   >  tot: 688 What's confusing me now is if the DoF should equal the number of non-empty cells for those subj interactions. Keep in mind I'm using lmer from R and thus being unbalanced is not a problem. This means that some subjects did not experience some rounds, for example, so those cells are blank. Should the subj interaction DoF match the actual number of observations, or just be the product of the relevant DoFs?


title: How to compute p-value from a permutation test?

text: I'm running a series of permutation tests inside 2 for loops (2 for loops will calculate some new data and then this data will be shuffled and then i apply a linear regression based on the shuffled data). Here is my code:               for(i in seq(1,12000,200))     {         for j in seq(1, 12000, 200))            {                 for(ind in 1:nrow(df))                 {                 ###calculate something (called X) depending on i,j, ind and df                 }           for (i in 1:100) # 100 sample             {                     sample(df$X,nrow(df), replace=FALSE)                       fit=lm(X ~ Y, df)   #linear regression calculation                      regerror=sum(residuals(fit)^2)                  res <- rbind(res, data.frame(shufflederr = regerror))                      }             }     }      The p-value is calculated based on the sum squared error of the linear regression in R which is:               > sum(res$shufflederr < observed$regerror)    #observed df contain the same error value but without shuffling      I have 3 questions:   1. Can a p-value of a pair `(i,j`) be equal to 0 ??is this statistically possible ?   2. Is there any optimal way to create a good cutoff of the `p-value` ?   3. Is there any better idea to do the permutation(or any other significance) test in my case ?


title: How to analyze a 2x3 design in which one level doesn't differ on one of the factors?

text: I have a 2x3 factorial design for my experiment: 3 levels of information given to participants (None, Moderate, Extreme), and 2 levels of time that the information focuses on (2050 or 2100), for those who received information. I'm curious about how to deal with the fact that my control group, receiving no information, basically does not differ on the Time factor. I'm a little worried that it reduces the power of my design to detect differences in the groups who did receive different (2050 or 2100) information. What's the best way to approach this? Should I instead do a 2 (Moderate, Extreme) x 2 (2050, 2100) design, essentially excluding the control group? If so, what test should I do to then compare the control group to those who received different levels of information?


title: Multiple regression with dummy variables

text: I need to calculate MR but my data consists of many dummy variables. I had a question in my survey that asked for the budget/month. But instead of giving the option to fill in a value, I offered ranges (steps of 250, so I have like 0-250; 250.01-500 and so on) I created dummy variables for these different steps so now I have a total of 6 dummy's and the problem is that I can't calculate the ln-values of these dummy's because all values are either 0 or 1. So any ideas how I can go around this problem?


title: Can Agresti-Coull binomial confidence intervals be negative?

text: According to http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm Agresti-Coull intervals cannot be negative; however using the formula from Wikipedia as well as the binom.confint function in R with x = 0 and n = 5 produces an interval of (-0.05457239, 0.4890549). Are there conflicting definitions of Agresti-Coull? Edit - Using binom.confint:               binom.confint(0, 5, .95, 'agresti-coull')              method x n mean       lower     upper     1 agresti-coull 0 5    0 -0.05457239 0.4890549      Manually calculating the lower bound (using the formula from Wikipedia):               > n <- 5     > x <- 0     > n2 <- n+qnorm(0.025)^2     > p <- 1/n2 * (0 + .5 * zn^2)      > zn <- qnorm(0.025)     1/n2 * (0 + .5 * zn^2) + zn * sqrt(1/n2 * p * (1 - p))     [1] -0.05457239


title: Assessing Significance in Neural Networks

text: I'm trying to determine statistical significance of inputs on the dependent variable using neural networks. I'm trying to test the null hypothesis that the partial derivatives are equal to zero. I found a method called Bootstraping but it requires to retrain the neural network a 1000 times. Is there any alternative less time consuming solution to this problem? Thanks,


title: Selection of lme models using AIC & appropriate random effects & variance structure

text: I am using three categorical predictor variables $X_1$, $X_2$, $X_3$ and one continuous dependent variable $Y$, and I want to treat $X_3$ as a random effect. The simplest model I could come with:                 M1<-lme(Y~X1*X2, random=~1|X3,method="ML")      Here are some more complex models:                 M2<-lme(Y~X1*X2, random=~1+X1|X3,method="ML")     M3<-lme(Y~X1*X2, random=~1|X3, weights = varIdent(form =~ 1 | X1*X2),method="ML")     M4<-lme(Y~X1*X2, random=~1+X1|X3, weights = varIdent(form =~ 1 | X1*X2),method="ML")             Model df       AIC       BIC    logLik   Test  L.Ratio p-value     M1   1   20  4211.434  4343.094 -2085.717                             M2   2   40 -4278.814 -4015.495  2179.407 1 vs 2 8530.249  <.0001     M3   3   37  2737.569  2981.139 -1331.784 2 vs 3 7022.383  <.0001     M4   4   57 -4863.952 -4488.722  2488.976 3 vs 4 7641.521  <.0001        * Please note that I had to increase iterations using lmeControl, to achieve convergence for M2 and M4.     1. Is it actually possible to use a categorical variable in a random slope model such as M2 & M4?     2. Should I be worry of the drop in AIC from high positives to low negatives?     3. Given that I have no missing values, do I have other possibilities than lme/lmer with equivalent structures of for random effects & variances?


title: How to plot a "differential word cloud" in R?

text: I have a questionnaire with one "open" text box, and 20,000 responses. I'm not planning to read all of them, but I'm thinking of interesting things to do with them. A general tag cloud would not be very interesting (I can predict which words would be the biggest), but it would be fun to split it along another categorical variable and see a tag cloud, or even just a word cloud, of the words that are the most different or most distinct for each categorical variable... Let's say this is the evaluation for a psych course... Of course "psychology", "lectures", "videos" etc will be the biggest words in general. But if we have a second question called "user category", and say "which words are the most common among "lifelong learners", and the least common among "students""... That would automatically eliminate the super-common words I mentioned above, and might show something interesting... If there are other interesting ways of visualizing words used compared with different categorical variables (maybe some kind of 2D graph, or network graph?), I'd love to hear about them as well. Thanks


title: Standardising the removal of outliers from a small data set

text: I have been wading through the many discussions on outliers on this site but I am still unfortunately having difficulty determining what to do with my data set. My study consists of a simple pre-post test setup whereby i conduct 6 tests prior to and 6 tests following my treatment. The purpose of the repeated tests is to decrease the effect of anomalous performance on the results and increase the sensitivity. We therefore expect that generally 1 or 2 of these 6 tests to be a bit different to the other 4. At present I've been removing abnormal scores by identifying those results that lie outside +/- 1 standard deviation. This practice has proven somewhat satisfactory (in that it detects most of those scores that upon visible inspection appear anomalous), however I read more and more that the use of standard deviation is not appropriate for determining outliers, and sometimes I've found that the scores deemed as outliers are not always those that appear the most anomalous. Therefore I'm wondering whether you can offer any suggestions as to how i can improve/better standardise my selection process? Is my current method acceptable? If so, do you know of any published material that has validated this approach (or does this approach have a name that i could google?) Alternatively, would the use of an interquartile range approach using the median value prove more satisfactory?


title: What does "covariance risk budget" mean in fPortfolio package?

text: I have been using the fPortfolio package and have found some inconsistency in the usage of terminology in the package. What does a "covariance risk budget" in the result returned by mean variance portfolio optimization- mean?


title: Estimating the variance of poker win rates

text: Suppose you have a casino with _n_ poker players. Each player has a win rate - the amount of money he wins or loses per hand. We assume that these win rates are normally distributed with a mean of 0. (We also assume that the players don't pay the casino any money.) Our goal is to estimate the variance _V_ of the distribution. For each player _x_ , we have observed a number of hands; we know how much money _x_ has won or lost on each of these hands. How would you go about estimating _V_? Can we get a better estimate if we add some empirical assumptions (in the vein of "in the long run, no-one can sustain a win rate of more than 1$/hand")? EDIT: Let me try to clarify what I mean by "win rate". If a player wins 500.000$ by playing a million hands then his _observed_ win rate is 0.5$/hand. With a million hands it's also likely that his _actual_ win rate is close to 0.5$/hand. The idea is that a player has an _actual_ win rate which cannot be observed directly, but which is a function of the player's skill. For example, if all players are equally skilled they will all have an actual win rate of 0; in this case, we also have _V_ =0. The question above is concerned with actual win rates. EDIT: My motivation for asking this question was to estimate how many players have an actual win rate of, say, more than 0.3$/hand. If you disagree with the assumptions made above, feel free to base your estimate on other assumptions.


title: Help with data analysis of small datasets

text: I have the following dataset (triplicate values of 5 independent measurements and duplicate values of a control):               Sample 1    Sample 2  Sample 3    Sample 4    Sample 5    C     181.8       58.2       288.9       273.2       290.9      53.9     120.3       116.8      108.9       281.3       446        39.6     86.1        148.5      52.9        126         150.3      The six conditions are independent production of a chemical in six different microorganisms. My aim was to find if the production of the chemical is significantly higher in Samples 1 - 5 than in C (control). At first, I carried out mean, SD and t-test (one-tail). Although SD error bars are large, two of the samples have mean values that are significantly higher than that of C. I decided to carry out median and Wilcoxon-Mann-Whitney test because of my worry about the size of my data and the high variations in the replicates. I was surprised however, that median and Wilcoxon-Mann- Whitney tests did not reveal statistical significant results. I will be happy if anyone could advise me on the best way to analyse this small dataset.


title: ANOVA and post-hoc analysis using R when we only have means and SD of three populations

text: I want to perform an ANOVA using R. I have three populations, represented by their respective means and SD: Pop.1: 5.5 +- 0.4 (n=100) Pop.2: 5.9 +- 0.3 (n=150) Pop.3: 6.2 +- 0.5 (n=200) Which is the exact code using R to perform the ANOVA using exclusively these data? Moreover, how can I perfom subsequently a post-hoc analysis? Thanks in advance


title: How does the RMS package's nomogram calculate points for continuous variables?

text: I have been reading a number of papers where researchers have created risk scores based on logistic regression models. Often they refer to "Sullivan's method" but I have no access to this paper and the explanations provided are far from clear. I noticed that Dr. Harrell's excellent RMS package provides a nomogram function which is in a way similar to creating a risk score (albeit with a very pretty graphical output). It seems after tinkering around with it that the way it works is by dividing the beta coefficients by the smallest beta coefficient and then multiplying a constant to create points for categorical variables. However I cannot for the life of me figure out what is going on with continuous variables. I've spent hours searching google without much luck, and I would appreciate if someone could shed some light on this for me. Thanks!


title: Analysis of intervals between events

text: I have been investigating the possibility of using the interval between uncommon events to test for changes in the frequency of such events over time.   As an example, say that the event is breaking a record in some sporting competition. This might occur at most a few times a year, and the data segmentation problem (whether an event falls within a certain interval of observation) causes the usual events-per-year analysis, usually a GLM with a Poisson link, to be strongly affected by even a single event without a very long series.   I first thought of sampling different observation intervals, and then decided to discard the intervals entirely and see if the intervals between events were associated with their serial position. This seemed to be a more powerful test in the data I had. Apparently this was an accepted procedure some 50-60 years ago, but I have not been able to find much in the recent literature. I'm pretty sure that this is related to analyzing the frequency of extreme weather events and the like, but I am not familiar with that field.   Anyone out there know whether this type of test has been superceded by something?


title: Is Jensen's inequality applicable for two variables?

text: This is the Jensen's inequality I saw in my textbook: $$E{ f(X) } \geq f( E(X) ),$$ where $f$ is a convex function. Is this also applicable for two random variables--independent or otherwise-- like this: $$E{ f(X,Y) } \geq f( E(X,Y) )?$$


title: About 2 unit root tests and null hypothesis

text: I have been looking at unit root testing. Specifically 2 tests:   1. **The ADF test.** The ADF (augmented Dickey Fuller) test has the null hypothesis that "the time series has a unit root" (meaning that the time series is not weakly stationary)    2. **The KPSS test.** The KPSS is a (Lagrange Multiplier type) test which has the null hypothesis that "the time series is weakly stationary" (so no unit root here). What I am looking for when applying those tests is to say whether a process is stationary, so I can do things with it. Somebody told me that because of how the null hypothesis was specified, on average one test would give more stationary processes than the other. Could anybody explain to me why ? This is what I think: because the ADF test has as $H_0$ "the process has a unit root", what I can say is "With a given confidence, say 99%, I cannot reject the fact that this process has a unit root"; because of this I can be pretty sure that a process is not stationary, but this is not the same that saying that a process is stationary, so the ADF test should "keep" more processes as stationary than the KPSS test ...


title: Significance testing for a group of samples

text: For instance, I have a bar chart - in which I have 4 samples/bars. for these 4 bars, I want to check statistical significance for all combinations.   1. How can I do it? Can unpaired t-test be used - to compare 2 samples at a time? or this is not correct. What is then the appropriate way in which all the samples/ bars can be compared against each other at the same time.   2. Is there any easy-quick way to nicely represent it in a chart? thanks


title: How to test for simultaneous equality of choosen coefficients in logit or probit model?

text: How to test for simultaneous equality of choosen coefficients in logit or probit model ? What is the standard approach and what is the state of art approach ?


title: Fit of a normal distribution to a one-dimensional dataset in R

text: I've got a set of (continuous) values from a measurement, where each object should be either positive or negative, and I know that the values of the "negative" objects should be approximately normally distributed. I've been using a `k-means`-based algorithm in `R` to cluster the data and thus to classify the objects. Here is a typical example of the clustering result, shown as a density plot: ![enter image description here](http://i.stack.imgur.com/aB0Ob.png) As you can see, although the algorithm does distinguish two populations, it includes a part of the positive objects into the negative cluster (the blue one). You can also see that there are some extremely negative outliers, causing the blue curve to be skewed to the right. These outliers are artefacts, resulting from the method the values are measured and calculated, and I'd like to find a cut-off to exclude them in future. Now, I'd like to fit a normal distribution curve the the negative population and to define a cut-off for those artefacts (outliers) based on this fit as well as to estimate some characteristics of the distribution. I was trying to do this by extracting only those objects that fall into the blue cluster and then fitting a normal distribution on them, but obviously this doesn't work out, mainly due to the bump on the right of the blue population. So I need another classification algorithm that would a priori assume the existence of two Gauss distributions and fit them (in the best case taking into account the possibility of the existence of outliers). However, so far I haven't been able to find an algorithm that would be reliable and fast enough (I've got around 20,000-50,000 objects), although I'm sure there must be a simple way in `R` to do so. So in summary, I'm looking for a classiciation algorithm implemented in `R` to fit a normal distribution to the "negative" objects of my measurement, ideally taking into account that outliers do exist. Does anyone has a suggestion how to accomplish this?


title: Arima model for non-negative data

text: I have been reading a tutorial for an introduction to time series. It contains a dataset, with an $Arima(2,0,0)$ forecast along with a 80% and 95% prediction interval. It looks like this: ![enter image description here](http://i.stack.imgur.com/cZMTu.png) This model was chosen due to the behaviour of the partial autocorrelation function ("it is zero after lag 2"). However, as the tutorial later explains, there is a clear problem with this modelling approach. The prediction interval predicts some negative values, which would be invalid for this data. (The data is assumed to be $\geq 0$). The forecast errors also have undesirable properties. **My question: How would you address this issue?** Some of my thoughts:   * transforming the data $x$ using $\log(x+1)$. Is this a respected method? Any pros / cons I should know about?   * any predicted values $<0$, naively set to zero. It feels, to me, like that there should be a better choice than this.


title: Fitting Inspection Time Data to Curve

text: I have a set of Inspection Time data that consists of two variables: SOA (stimulus onset asynchrony) and Accuracy. This data refers to how accurately you can perform on a simple stimuli discrimination task at a given duration of exposure. The accuracy data is binary (0,1) and the SOA data is in increments of .01 (10 milliseconds).               SOA = [.01, .02, .02, .03, .02, .01.....]     Accuracy = [1, 0, 1, 1, 0, 1.....]      What I need to do is estimate the SOA at 75% accuracy. I've been told that this is done by fitting the data to a curve (y-axis accuracy, x-axis SOA) and then taking the SOA value at 75%. Unfortunately I have mixed information about what kind of curve to use: sigmoid, gauss, logistic, and cumulative normal curve. So far I've managed to fit a binary logistic curve to the data using R, but there seems to be no method of estimating SOA at 75% accuracy. Could someone please suggest a method for accomplishing my task? Any help would be greatly appreciated.


title: Problem related to OLS estimators

text: The problem is: > Suppose we fit a model Y = XA βA + ε. However, the true model is Y = XA βA + > XB βB + ε. A is kA x 1 and B is kB x 1. Show that the OLS estimates of βA > will still be equal if XTAXB = 0. My attempt at a solution: Under OLS, the estimate of β is (XTX)-1XTY. So I was thinking for the first model, we just have X = XA, but in the second model, we have X = XA + XB. So I was trying to expand (XTX)-1XTY but subbing in the previous values for X, and then seeing if anything stood out, but I wasn't having any luck. Any help or tips would be appreciated


title: Do all these estimates of kurtosis and skewness have the same (asymptotic) distribution under normal sample distribution?

text: I have seen five types of estimates of kurtosis and skewness:   * three from http://stats.stackexchange.com/a/84057/1005   * one from page 9 of Analysis of Financial Time Series by Ruey S. Tsay   * one from http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm For an iid normal sample, do all those types of estimates of kurtosis and skewness have the same distribution for each sample size, and have the same asymptotic distribution? In tests for normality using skewness and kurtosis (such as D'Agostino's K-squared test, and Jarque–Bera test),   * does it matter which type of estimates for skewness and kurtosis are used in building the test statistics?    * Will the null distributions of the test statistics based on those types of estimates for skewness and kurtosis be the same? Thanks!


title: Contour plot for a likelihood function

text: For the likelihood associated with the mixture model `1/4*Normal(mu1,1)+3/4*Normal(mu2,1)` I would like to have a figure as like as following one for a simulated sample of 500 observation from `mu1=0` and `mu2=2.5`. My `R` code is:               da=rbind(rnorm(10^2),2.5+rnorm(4*10^2))       like=function(mu1,mu2){       sum(log((.25*dnorm(da-mu1)+.75*dnorm(da-mu2))))}       mmu1=seq(-2,5,length.out=length(da))       mmu2=seq(-2,5,length.out=length(da))       post=outer(mmu1,mmu2,like)         image(mmu1,mmu2,post,xlab=expression(mu1),ylab=expression(mu2))      contour(mmu1,mmu2,post,add=T)       which contains the error:               "Error in dnorm(da-mu1):      dims [product 800] do not match the length of object"      I've changed the length of `mmu1` & `mmu2`, but nothing seemed to be happening. Could someone help me? ![enter image description here](http://i.stack.imgur.com/dKyiH.png)


title: What does it mean that a statistic $T(X)$ is sufficient for a parameter?

text: I am having a hard time understanding what a sufficient statistic actually helps us do. It says that Given $X_1, X_2, ..., X_n$ from some distribution, a statistic $T(X)$ is sufficient for a parameter $\theta$ if $P(X_1, X_2, ..., X_n|T(X), \theta) = P(X_1, X_2, ..., X_n|T(X))$. Meaning, if we know $T(X)$, then we cannot gain any more information about the parameter $\theta$ by considering other functions of the data $X_1, X_2, ..., X_n$. I have two questions: 1) It seems to me that the purpose of $T(X)$ is to make it so that we can calculate the pdf of a distribution more easily. If calculating the pdf yields a **_probability measure_** , then why is it said that we cannot "gain any more information about the parameter $θ$"? **In other words, why are we focused on $T(X)$ telling us something about $\theta$ when the pdf spits out a probability measure, which isn't $\theta$?** 2) When it says: "we cannot gain any more information about the parameter θ by **_considering other functions of the data $X_1,X_2,...,X_n$._** ", what other functions are they talking about? Is this akin to saying that if I randomly draw $n$ samples and find $T(X)$, then any other set of $n$ samples I draw give $T(X)$ also? Thanks!!


title: How to model zero-inflated continuous response using categorical predictors - preferably resulting in multiplicative parameters

text: I'm having trouble finding a suitable model for predicting the AVG value (revenue in cents) of a single click on a product on a large e-commerce site. (assuming a click leading directly to a purchase has value equal to its revenue) Using only categorical (perhaps ordered) predictors, how do I best model a continuous response variable -- most likely Revenue Per Click (RPC) -- which is:   * non-negative   * heavily zero-inflated      * depends on data granularity availability. The % data mass at zero is: ~45% with aggregation by unique predictor combinations (I have this), but in reality it would be ~97% for click-level granularity (don't have this)   * has a very large right skew      * aggregated by unique predictor combinations, the non-zero values average nearly 400 with SD a few times larger than this! (violates Poisson??) For various reasons I'd really like to end up with a set of multiplicative parameters so that for any combination of categorical predictors, I can multiply these parameters (or a transform of them) to estimate a mean response. Further details:   * Predictors will eventually include some or all of these categoricals:      * "product type" (typ, can be as few as 4 types)     * "clicker country" (cc, appx 30 countries/regions)     * "product score" (val, 10+ ordered values) FYI: I'm actually estimating the revenue per click for thousands of products. There are about as many total purchases in the data as products, though many products have no purchases (usually having very few clicks) & none will have purchases for anywhere near EVERY categorical predictor combination (ex: 4 x 30 x 10). I suppose that makes this a bit sparse. Any ideas how to handle this? My initial urge was to take the log of RPC to get multiplier parameters:               log(RPC) ~ typ+cc+val, weights=Clicks      But this fails with so many RPC==0. I next considered adding a small constant to ensure the log is better behaved:               log(RPC+k) ~type+cc+val, weights=Clicks      But whether adding tiny or large constant values, the parameters predicted wildly different/inaccurate values after exponentiating and subtracting the k. :( This led me to my most recent idea:   * Can I just use Poisson (loglinear) regression & pretend the continuous RPC value is a count?   * If so, does Poisson regression in R allow for weights (by the Clicks column) ? If not, the model will be heavily biased. Thanks for your insight!


title: Fisher information and the "surface area of the typical set"

text: On page 331 of "Elements of Information Theory" (1991), author says that while entropy is related to the volume of the typical set, Fisher information is related to the surface area of the typical set, but I can't find anything more on this...can anyone explain this connection?


title: I just installed the latest version of R. What packages should I obtain?

text: **Duplicate thread:** What R packages do you find most useful in your daily work? Are there any R packages that are just plain good to have, regardless of the type of work you are doing? If so, what are these packages? If not, what packages do you find the most useful?


title: Getting residuals from (kernel) canonical correlation analysis

text: I'm playing around with kernel Canonical Correlation Analysis, as implemented in the R package kernlab. Is there a simple way to extract the residuals after fitting? (Is this a well-defined quantity?) Here's some sample code:               library(kernlab)     x = as.matrix(rnorm(100))     y = as.matrix(.5 * x + rnorm(100,0,.5))     fit = kcca(x,y)      fit contains "kcor" (length = 10 since the default for kcca is 10 canonical components), "xcoef" (100 by 10) and "ycoef" (100 by 10).


title: How many people bought wine?

text: Rephrased a problem trying to solve for work in terms of people buying wine, also included progress made so far. **Set-up:** Customers enter a winery with the option of buying bottles of wine. Those who chose not to purchase can leave through one exit while those who wish to purchase go to the cash register. After exiting the register they are asked a survey "How many bottles of wine did you just purchase?". All have to answer. They can give any number between zero and their true value, but not more. If they answer 'zero' nothing is recorded. Any other answer is added to the survey total along with incrementing the number of positive respondents. Knowing:   * $T$ : Actual number of bottles sold (known, from inventory)    * $R$ : Total bottles recorded sold from survey    * $A$ : Total respondents that gave non-zero answer    * Find $X$, total number of actual customers.   * Also given $U$ : Theoretical maximum for answer.  _**U** can be explained that the maximum number of customers can be either the total number of bottles of wine on the shelf (one bottle per customer), or the capacity of the store in terms of customers._ **Small illustrative example:** ![aa](http://i.imgur.com/YDJ79OM.png)   * $T = 12$ bottles actually sold (all squares)    * $R = 7$ bottles recorded being sold (black squares)   * $A = 3$ people gave non-zero answers (brackets with at least one black square)   * $X = 4$ actual number of customers (all brackets; unknown) **Progress so far:** At first tried simple proportion $\hat{X}=\frac{A T}{R} $ but found it over- estimated the true population. This led to exploring what properties the estimate should have. Found at least six:   1. The minimum $X$ can be is $A$. _(At least $A$ people bought wine.)_   2. The maximum $X$ can be is the minimum of $U$ and $T-R+A$. Define this to be $V$. _(Reason: If we know that $A$ people bought $R$ bottles then the remaining $T-R$ bottles could in theory be bought by $T-R$ individual people. Take the minimum as $U$ could be lower due to store capacity)_   3. $X$ can only be $A$ when $R=T$ _(All bottles were bought, nothing to estimate)_   4. $X$ can only be the maximum when $R=A$ _(Every respondent said they bought one bottle)_   5. If they actually purchased a bottle, they purchased at least one.   6. If they are recorded as purchasing a bottle, they purchased at least one. [1] and [2] imply it could be a weighted average between $A$ and $V$ for some value $w$ $$\hat{X}=Aw + V(1-w)$$ [3] and [4] imply the weight could be linear $$w=\frac{R-A}{T-A}$$ when ($R=A$ then $\hat{X}=V$) and ($R=T$ then $\hat{X}=A$). [5] and [6] imply the weight should be modified. Each customer who bought wine, bought at least one; and each customer who gave a non-zero answer in the survey also bought at least one.Out of $R$ recorded only $R-A$ are 'free', the first $A$ were needed to assign 1-1 to each customer. Same argument suggests the denominator should be $T-X$. With $X$ unknown substitute with $\hat{X}$ **Summary:** Define   * $V=Min[U,T-R+A]$    * $w=\frac{R-A}{T-\hat{X}}$   * Solve for $\hat{X}$ in the equation: $\hat{X}=Aw + V(1-w)$ The new estimate gave significantly better answers than simple proportion estimate. **Questions**   1. Is this the best estimate that can be made with given information? (Am I missing an assumption or property which would make it better?)   2. Is there a way to estimate the variance?


title: Identifying outliers for non linear regression

text: I am doing research on the field of functional response of mites. I would like to do a regression to estimate the parameters (attack rate and handling time) of the Rogers type II function. I have a dataset of measurements. **How can I can best determine outliers?** For my regression I use the following script in R (a non linear regression): (the dateset is a simple 2 column text file called `data.txt` file with `N0` values (number of initial prey) and `FR` values (number of eaten prey during 24 hours):               library("nlstools")     dat <- read.delim("C:/data.txt")         #Rogers type II model     a <- c(0,50)     b <- c(0,40)     plot(FR~N0,main="Rogers II normaal",xlim=a,ylim=b,xlab="N0",ylab="FR")     rogers.predII <- function(N0,a,h,T) {N0 - lambertW(a*h*N0*exp(-a*(T-h*N0)))/(a*h)}     params1 <- list(attackR3_N=0.04,Th3_N=1.46)     RogersII_N <-  nls(FR~rogers.predII(N0,attackR3_N,Th3_N,T=24),start=params1,data=dat,control=list(maxiter=    10000))     hatRIIN <- predict(RogersII_N)     lines(spline(N0,hatRIIN))     summary(RogersII_N)$parameters      For plotting the calssic residuals graphs I use following script:               res <- nlsResiduals (RogersII_N)     plot (res, type = 0)     hist (res$resi1,main="histogram residuals")         qqnorm (res$resi1,main="QQ residuals")     hist (res$resi2,main="histogram normalised residuals")         qqnorm (res$resi2,main="QQ normalised residuals")     par(mfrow=c(1,1))     boxplot (res$resi1,main="boxplot residuals")         boxplot (res$resi2,main="boxplot normalised residuals")      ### Questions   * **How can I best determine which data points are outliers?**   * **Are there tests I can use in R which are objective and show me which data points are outliers?**


title: machine learning to predict equations / parameters of equations

text: I'm not sure if machine learning is the best way to do this, but I'm interested in seeing if this problem is feasible. Normally you use machine learning for classification. ie.given the size of a tumor and the color(darkness) of the tumor, determine whether its lethal or benign. However, say you want to figure out how the tumor has been growing/how it will end up in the future, given that you know its current size and current color. Can you use machine learning to predict an equation or curve that describes its size, for the duration of the tumor's life? This is how my data looks:               Subject   SIZE  COLOR   GROWTH     1          5      10    3.4/(1+t^-1.2)     2          8      62    3.8/(1+t^-1.5)     3          15     77    6.5/(1+t^-3)     4          40     94    8.6/(1+t^-15)     …          …      … …


title: negative binomial GLM with time elapsed as a predictor

text: Suppose you have a response variable, number of transactions a customer makes in the 24 hour period following acquisition, which follows a negative binomial distribution. Let's call the response _trans_. Now it's easy to model the number of transactions the customer makes in this immediate 24 hour period, or _trans_ , using Library(MASS) in r and running a negative binomial GLM. I can get a nice fit by thoughtful selection of predictors. Now I want to model the number of transactions a customer makes in _any_ 24 hour period as a function of the number of days elapsed since the customer was acquired (call this predictor _t_ ). Empirically, as time elapsed increases, the 24 hour transaction rate decreases, although the rate of decrease is a decreasing function of time. Plugging number of days elapsed, or _t_ , into the aforementioned negative binomial GLM runs in to all sorts of problems. I can get a good fit by adding a quadratic term _t^2_ but I cannot project beyond the number of days in the data set without getting illogical results; clearly adding a quadratic term is just curve fitting without capturing the underlying relationship of 24 hour transaction rates vs time elapsed since acquisition. Is there any better way to approach this problem? I will add more details if this is unclear.


title: K-S test correcting, or is a $\chi^2$ acceptable?

text: I have a data set (21 binned values) which I have fitted to a Gaussian with IDL ($\mu=3.825$, $\sigma =0.0377$). I have tried to find the $\chi^2$, getting $21.14$ for $20$ d.f.. Here my understanding is a little shakey; I think this means I can be 60% confident that the distribution is normal (please correct me if I'm wrong, every information source seems to say something different at this point!). I went on to try and use the K-S test, with peculiar results; when attempted manually, I found $D=0.09958$ (without ordering min->max), when I ordered values I got $D=0.12$. When given to IDL to calculate it gave $D=0.143$, with a 0.97 probability that it fits the normal. These D-values seem to point to very high probability, and I suspect I've muddled something. Have I too few data bins? Should I need to go so far as the K-S if $\chi^2$ gives me a fair result? I'm supplying the data in case it makes more sense than me. Observed Frequency, Expected frequency (Gauss.)               0,  0.090119938     0,  0.202768898     0,  0.42521436         1,  0.831075789        0,  1.51390793        2,  2.570303594        2,  4.067199812        0,  5.998362667        13, 8.245102466        9,  10.5629572        13, 12.61249796        12, 14.03598269        12, 14.5583        10, 14.07358146         14, 12.68015971       8,  10.64807109        8,  8.333804074         4,  6.079134555         2,  4.133009011         0,  2.618888894         0,  1.54665668


title: Training a decision tree against unbalanced data

text: I'm new to data mining and I'm trying to train a decision tree against a data set which is highly unbalanced. However, I'm having problems with poor predictive accuracy. The data consists of students studying courses, and the class variable is the course status which has two values - Withdrawn or Current. Age Ethnicity Gender Course ... Course Status In the data set there are many more instances which are Current than Withdrawn. Withdrawn instances only accounting for 2% of the total instances. I want to be able to build a model which can predict the probability that a person will withdraw in the future. However when testing the model against the training data, the accuracy of the model is terrible. I've had similar issues with decision trees where the data is dominated by one or two classes. What approach can I use to solve this problem and build a more accurate classifier? Thanks in advance for your help!


title: Repeated measures mixed model using lmer in R

text: I’m hoping to get some guidance in specifying a mixed model using the lme4 package in R. The study is quite straightforward. It’s a repeated measures design with pre/post measurements on the variable of interest where half of the participants received the treatment (a physical therapy program) and the other half did not. The research was carried out in groups at 8 different hospitals. Because the research was carried out in groups, I’m treating the data as multilevel with Time at level 1, Individual at level 2, and Group at level 3. The structure of the data can be seen below:               'data.frame':   1856 obs. of  6 variables:      $ ID           : num  1 1 2 2 3 3 4 4 5 5 ...      $ Group        : Factor w/ 8 levels "Richmond",..: 1 1 1 1 1 1 1 1 1 1 ...      $ Age          : num  12 12 14 14 13 13 18 18 16 16 ...      $ Treatment    : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...      $ Time         : Factor w/ 2 levels "Pre","Post": 1 2 1 2 1 2 1 2 1 2 ...      $ DV           : num  4 7 4 3 4 8 8 9 5 5 ...      I’m trying to answer the following questions:   1. Was there an improvement in the DV over time?   2. Did participants in the treatment program improve more than those in the control group? I’m not expecting any, but I’d also like to know:   3. If there are differences in the DV attributable to age.   4. If there are differences in the DV attributable to the location (group). To that end I have specified the following model treating Time, Treatment and Age as fixed effects and Group as a random effect.               lmer(DV ~ Time + Treatment + Time*Treatment + Age + (1 | ID) + (1 | Group),           data = TLHB.data, REML = FALSE)      However, I’m not sure that this model accounts for the fact that participants were nested in groups so am looking for guidance that my approach so far seems reasonable given my research questions and how to properly specify that the participants within each group were not necessarily independent.


title: Dependence of distribution standard deviation with subtraction of its mode from mean of normal distribution

text: The following question might have relations to this question: Given $\mu$ and $\sigma$ (mean and standard deviation) of normal distribution, find the set of all distributions with their $\mu_1$ and $\sigma_1$ and $Mo_1$ (mode) that $\mu-Mo_1$ is any increasing function of $\sigma_1$. If not a general result, at least providing a clue, references or examples of such distributions are very much appreciated. My question might be wrong or I could put it in a better way.


title: Is it necessary to sample a raster for spatial regression?

text: I'm looking to model land cover change using a variety of environmental predictors (e.g. elevation, rainfall, etc.) stored as raster layers. In most similar studies I've found in the literature the authors sample the raster layers at some subset of the cells before performing the regression. I'm trying to understand if this sampling is done purely for efficiency, i.e. because it would be prohibitively slow to run the analysis over the whole grid? Or, is there a statistical reason that it is better to use a subset of the data, e.g. to reduce spatial autocorrelation?


title: What is the relationship between the GINI score and the log-likelihood ratio

text: I am studying classification and regression trees, and one of the measures for the split location is the GINI score. Now I am used to determining best split location when the log of the likelihood ratio of the same data between two distributions is zero, meaning the likelihood of membership is equally likely. My intuition says that there must be a connection of some sort, that GINI has to have a good foundation in a mathematical theory of information (Shannon) but I don't understand GINI well enough to derive the relationship myself. **Questions:**   * What is the "first principles" derivation of GINI impurity score as a measure for splitting?   * How does the GINI score relate to log of likelihood ratio or other information-theoretic fundamentals (Shannon Entropy, pdf, and cross entropy are part of those)? **References:**   * http://mathworld.wolfram.com/GiniCoefficient.html   * How is the Weighted Gini Criterion defined?   * Mathematics behind classification and regression trees   * http://www.cs.put.poznan.pl/jstefanowski/sed/DM-5-newtrees.pdf   (added)   * http://www.ibe.med.uni-muenchen.de/organisation/mitarbeiter/020_professuren/boulesteix/pdf/gini.pdf   * https://www.youtube.com/watch?v=UMtBWQ2m04g   * http://www.ius-migration.ch/files/content/sites/imi/files/shared/documents/papers/Gini_index_fulltext.pdf   * http://stackoverflow.com/questions/4936788/decision-tree-learning-and-impurity **Shannon's entropy** is described as: $$ H \left(x \right) = \Sigma_{i} P\left(x_{i} \right)\log_{b} P\left(x_{i} \right) $$ Extending this to the multivariate case we get: $$ H \left(X,Y \right)= \Sigma_{x}\Sigma_{y} P\left(x,y \right)\log_{b} P\left(x,y \right) $$ Conditional Entropy is defined as follows: \begin{align} H \left(X|Y \right) &= \Sigma_{y} p\left(x,y \right)\log_{b} \frac {p\left(x \right)} {p\left(x,y \right)} \newline &\text{or,} \newline H \left(X|Y \right) &= H \left(X,Y \right) - H \left(Y \right) \end{align} The log of the ratio of likelihoods is used for abrupt change detection and is derived using these. (I don't have derivation in front of me.) **GINI Impurity:**   * The general form of GINI impurity is $ I = \sum_{i=1}^m f_{i} \cdot \left( 1-f_{i}\right) $ **Thoughts:**   * Splitting is done on a measure of impurity. High "purity" is likely the same as low entropy. The approach is likely related to entropy minimization.   * It is likely that the assumed basis distribution is uniform, or possibly with hand-waving, Gaussian. They are likely making a mixture of distributions.   * I wonder if the Shewhart chart derivation can apply here?   * The GINI Impurity looks like the integral of the probability density function for a binomial distribution with 2 trials, and one success. $P(x=k)= \begin{pmatrix} 2\\ 1\end{pmatrix} p \left( 1-p \right) $  (additional)   * The form is also consistent with a Beta-binomial distribution which is a conjugate prior for a Hypergeometric distribution. Hypergeometric tests are often used to determine which samples are over or under represented in a sample. There is also a relationship to Fisher's exact test, whatever that is (note to self, go learn more about this).


title: Blocked Weighted Bootstrap

text: Bootstrap is a well known resampling method. But I want to know what is blocked weighted bootstrap sampling? Why we need this?


title: P value for interaction term in mixed effects models using lme4

text: I'm analysing some behavioural data using `lme4` in `R`, mostly following Bodo Winter's excellent tutorials, but I don't understand if I'm handling interactions properly. Worse, no-one else involved in this research uses mixed models, so I'm a bit adrift when it comes to making sure things are right. Rather than just post a cry for help, I thought I should make my best effort at interpreting the problem, and then beg your collective corrections. A few other asides are:   * While writing, I've found this question, showing that `nlme` more directly give p values for interaction terms, but I think it's still valid to ask with relation to `lme4`.   * I realise this is an overly long question, and maybe not completely appropriate for the StackExchange format. I've excluded some of the R output to make things a little shorter. Sorry.   * `Livius'` answer to this question provided links to a lot of additional reading, which I'll be trying to get through in the next few days, so I'll comment with any progress that brings. * * * In my data, I have a dependent variable `dv`, a `condition` manipulation (0 = control, 1 = experimental condition, which should result in a higher `dv`), and also a prerequisite, labelled `appropriate`: trials coded `1` for this should show the effect, but trials coded `0` might not, because a crucial factor is missing. I have also included two random intercepts, for `subject`, and for `target`, reflecting correlated `dv` values within each subject, and within each of the 14 problems solved (each participant solved both a control and an experimental version of each problem).               library(lme4)     data = read.csv("data.csv")          null_model = lmer(dv ~ (1 | subject) + (1 | target), data = data)     mainfx_model = lmer(dv ~ condition + appropriate + (1 | subject) + (1 | target), data = data)     interaction_model = lmer(dv ~ condition + appropriate + condition * appropriate + (1 | subject) +          (1 | target), data = data)     summary(interaction_model)      Output:               ## Linear mixed model fit by REML ['lmerMod']     ## ...excluded for brevity....     ## Random effects:     ##  Groups   Name        Variance Std.Dev.     ##  subject  (Intercept) 0.006594 0.0812       ##  target   (Intercept) 0.000557 0.0236       ##  Residual             0.210172 0.4584       ## Number of obs: 690, groups: subject, 38; target, 14     ##      ## Fixed effects:     ##                                Estimate Std. Error t value     ## (Intercept)                    0.2518     0.0501    5.03     ## conditioncontrol               0.0579     0.0588    0.98     ## appropriate                   -0.0358     0.0595   -0.60     ## conditioncontrol:appropriate  -0.1553     0.0740   -2.10     ##      ## Correlation of Fixed Effects:     ## ...excluded for brevity.      ANOVA then shows `interaction_model` to be a significantly better fit than `mainfx_model`, from which I conclude that there's a significant interaction present (p = .035).               anova(mainfx_model, interaction_model)      Output:               ## ...excluded for brevity....     ##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)       ## mainfx_model       6 913 940   -450      901                               ## interaction_model  7 910 942   -448      896  4.44      1      0.035 *     ## ---     ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1      From there, I isolate a subset of the data for which the `appropriate` requirement is met (i.e. `appropriate = 1`), and for it fit a null model, and a model including `condition` as an effect, compare the two models using ANOVA again, and lo, find that `condition` is a significant predictor.               good_data = data[data$appropriate == 1, ]     good_null_model = lmer(dv ~ (1 | subject) + (1 | target), data = good_data)     good_mainfx_model = lmer(dv ~ condition + (1 | subject) + (1 | target), data = good_data)          anova(good_null_model, good_mainfx_model)      Output               ## Data: good_data     ## models:     ## good_null_model: dv ~ (1 | subject) + (1 | target)     ## good_mainfx_model: dv ~ condition + (1 | subject) + (1 | target)     ##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)       ## good_null_model    4 491 507   -241      483                               ## good_mainfx_model  5 487 507   -238      477  5.55      1      0.018 *     ## ---     ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


title: Help in regression interpretation, including interaction terms

text: I m doing a quantitative analysis for my dissertation and I m having issues with the interpretation of my results (I have little knowledge of SPSS). I m doing a panel data analysis with multiple regression. What, I want to understand is when I should keep or remove a control variable? My R2 and adjusted R2 are optimal when I keep my set of variables. But I have a variable (political stability) that, when I add my explanatory variable (FDI) becomes insignificant and another variable (inflation rate) which logically should be negative becomes positive. I have trouble understanding it: what does it mean? And also how to interpret interaction terms? For instance I m studying the effect of FDI on economic growth: both openness of the economy and FDI are significant and positive on FDI, but the association of FDI and openness comes with a negative coefficient. How can I explain this?


title: Good, useful and characteristic experiments for (optimal) statistical design of experiments

text: There are more phenomena to which experimental design may be applied than there are alternative valid design strategies. This should be true, though there are many ways to properly design an experiment. What are the best "problems" that truly demonstrate the value and nuance for the different types of optimal design of experiments? (A, D, E, C, V, phi, ....) Can you provide books, links, articles, references, or at least good empirically driven opinions?


title: How to test if counts differ across temporal categories?

text: I have a homework question in which data is taken from three months of hospital admission data. It is then split into six 4-hour groups (0000–400, 0400–0800, etc.). The data are given with just the hour group and one column containing a number of admits. We are to draw a conclusion about how nurses should be scheduled based on this data. Is there any statistical or graphical method that would result in a significant conclusion, or is this a trick question, and the answer is that there isn't enough data to make a suggestion with any certainty? **_Update:_** the values for the 4 hour window are for the entire 3 month period. So if the value that corresponded to 00:00-04:00 was 80 that would be 80 patients admitted between 00:00 and 04:00 over the last 3 months. If I did a single sample t-test using the values of column 2 and any values that were outside the CI, could I say that more nurses should be allocated to that hour group?


title: Factor analysis with repeated measures

text: Multilevel factor analysis seems to be the technical term for factor analysis with repeated measures, judging from this abstract. To be precise, following Wikipedia's factor analysis notation, the model I want to build is $$x_i =l_{i1} F_1 + \cdots + l_{ik} F_k + z_i + \varepsilon_i$$ where $x_i$ is the $i$th observed variable (already centered and scaled, say), an $n\times 1$ vector. The thing that makes this model different from ordinary factor analysis is the presence of the $n\times 1$ vector $z_i$ on the right- hand side; this is a vector of fixed or random effects that correspond to the repeated measures. Specifically, $z_{i(p)} = z_{i(q)}$ whenever the $p$th and $q$th records come from the same individual. Multiple queries similar to this one exist (here and here). This question is only slightly more general while hopefully also more expository: (A) Where can I find a publicly available and detailed description of multilevel factor analysis? (B) What software exists to do multilevel factor analysis in a pretty straightforward way? Solutions involving R, SAS, Python, or Latent GOLD are of particular interest.


title: Autocorrelation summation identity

text: Let $X_t$ be a weakly stationary process with mean $\mu$ and autocovariance function $\gamma$. How do I show that $$n^{-2}\sum_{i=1}^n \sum_{j=1}^n cov(X_i, X_j)$$ equals $$ n^{-2} \sum_{i-j=-n}^n (n- |i-j|)\gamma(i-j)$$ Thanks


title: Is rstan or my grid approximation incorrect: deciding between conflicting quantile estimates in Bayesian inference

text: I have a model to achieve Bayesian estimates the population size $N$ and probability of detection $\theta$ in a binomial distribution solely based on the observed number of observed objects $y$: $$ p(N,\theta|y)\propto \frac{ \text{Bin}(y|N,\theta)}{N} $$ for $ \left\{N|N\in\mathbb{Z}\land N\ge \max(y)\right\}\times(0,1) $. For simplicity, we assume that $N$ is fixed at the same, unknown value for each $y_i$. In this example, $y=53,57,66,67,73$. This model, when estimated in `rstan`, diverges from the results obtained from a grid approximation of the posterior. I'm trying to pin down why. (Interested readers might find that this question is a follow-on to my answer here.) # `rstan` Approximation For reference, this is the rstan code.               raftery.model   <- "         data{             int     I;             int     y[I];         }         parameters{             real<lower=max(y)>  N;             simplex[2]      theta;         }         transformed parameters{         }         model{             vector[I]   Pr_y;                  for(i in 1:I){                 Pr_y[i] <-  binomial_coefficient_log(N, y[i])                             +multiply_log(y[i],         theta[1])                             +multiply_log((N-y[i]),     theta[2]);             }             increment_log_prob(sum(Pr_y));             increment_log_prob(-log(N));                     }     "     raft.data           <- list(y=c(53,57,66,67,72), I=5)     system.time(fit.test    <- stan(model_code=raftery.model, data=raft.data,iter=10))     system.time(fit     <- stan(fit=fit.test, data=raft.data,iter=10000,chains=5))      Note that I cast `theta` as a 2-simplex. This is just for simplicity. The quantity of interest is `theta[1]`; obviously `theta[2]` is superfluous information. Additionally, $N$ is a real value (`rstan` only accepts real-valued parameters because it is a gradient method), so I wrote a real-valued binomial distribution. ## Rstan results                           mean se_mean       sd   2.5%    25%    50%    75%   97.5% n_eff Rhat     N        1078.75  256.72 15159.79  94.44 148.28 230.61 461.63 4575.49  3487    1     theta[1]    0.29    0.00     0.19   0.01   0.14   0.27   0.42    0.67  2519    1     theta[2]    0.71    0.00     0.19   0.33   0.58   0.73   0.86    0.99  2519    1     lp__      -19.88    0.02     1.11 -22.89 -20.31 -19.54 -19.09  -18.82  3339    1      # Grid Approximation The grid approximation was produced as below. Memory constraints prevent me making a finer grid on my laptop.               theta   <- seq(0+1e-10,1-1e-10, len=1e3)     N       <- round(seq(72, 5000, len=1e3)); N[2]-N[1]     grid    <- expand.grid(N,theta)     y   <- c(53,57,66,67,72)     raftery.prob    <- function(x, z=y){         N       <- x[1]         theta   <- x[2]         exp(sum(dbinom(z, size=N, prob=theta, log=T)))/N     }          post        <- matrix(apply(grid, 1, raftery.prob), nrow=length(N), ncol=length(theta),byrow=F)         post.norm   <- post/sum(post)      I used the grid approximation to produce this display of the posterior density. We can see that the posterior is banana-shaped; this kind of posterior can be problematic for euclidian metric HMC. (The severity of the banana shape is actually suppressed here since $N$ is on the log scale.) If you think about the banana shape for a minute, you'll realize that it must lie on the line $\bar{y}=\theta N$. (Additionally, the grid approximation displayed in this graph is not normalized for reasons of clarity - else the banana is a little too narrow to clearly make out.) ![posterior over a grid](http://i.stack.imgur.com/kiR50.png) ## Grid approximation results               do.call(cbind, lapply(c(0.025, .25, .5, .75, .975), function(quantile){         approx(y=N, x=cumsum(rowSums(post.norm))/sum(post.norm), xout=quantile)     }))       [,1]     [,2]     [,3]     [,4]     [,5]         x 0.025    0.25     0.5      0.75     0.975        y 92.55068 144.7091 226.7845 443.6359 2475.398      # Discussion The 97.5% quantile for $N$ is much larger in my `rstan` model than it is for the grid approximation, but its quantiles are similar to the grid approximation otherwise. I interpret this as indicating that the two methods are generally in agreement. I do not know how to interpret the discrepancy in the 97.5% quantile, though. I've developed several possible explanations for what might be accounting for the divergence between the grid approximation and the results from `rstan` HMC-NUTS sampling, but I'm uncertain how to understand if one, both or neither explanation is correct.   1. Rstan is wrong and the grid is correct. The banana-shaped density is problematic for `rstan`, especially as $N$ drifts off towards $+\infty$, so these tail quantities are not trustworthy. We can see from the plot of the posterior over the grid that the tail is very sharp at larger values $N$.   2. Rstan is correct and the grid is wrong. The grid makes two approximations which may undermine the results. First, the grid is only a finite set of points over a subspace the posterior, so it is a rough approximation. Second, because it's a finite subspace, we're falsely declaring there to be 0 posterior probability over values $N$ larger than our largest grid value for $N$. Likewise, `rstan` is better at getting into the tails of the grid, so its tail quanitles are correct.  * * * I needed more space to clarify a point from Juho's answer. If I understand correctly, we can integrate $\theta$ out of the posterior to obtain the beta- binomial distribution: $$ p(y|N,\alpha,\beta)={N\choose y} \frac{\text{Beta}(y+\alpha, N-y+\beta)}{\text{Beta}(\alpha,\beta)} $$ In our case, $\alpha=1$ and $\beta=1$ because we have a uniform prior on $\theta$. I believe that the posterior should then be $p(N|y)\propto N^{-1}\prod_{i=1}^K p(y_i|N,\alpha=1,\beta=1)$ where $K=\#(y)$ because $p(N)=N^{-1}$. But this appears to wildly diverge from Juho's answer. Where have I gone wrong?


title: MRQAP in multiple samples (Social Network Analysis)

text: I have a set of samples, each depicting a (full) social network of coworkers. The questionnaire, which was used to gather the data, was exactly the same for all the samples; just the respondents and the number of respondents is different. For hypothesis-testing, I use multiple regression (MRQAP). For now, I did this for each individual sample on it's own. However, I would also be interested in the "pooled" results. Is it possible to run an MRQAP over a number of samples simultaneously (where the algorithm of course knows that individuals of different samples have no chance of connecting)? Or, alternatively, is it possible to run the analyses for each sample (as I have done) and then pool the results in the end (similar to what is done during a multiple imputation procedure)?


title: How can I replicate R forecasts for seasonal ARIMA?

text: I have obtained the following estimations and forecasts in R for a seasonal ARIMA(1, 0, 1)(1, 0, 1)[7] > model1 Series: PO ARIMA(1,0,1)(1,0,1)[7] with zero mean Coefficients:                     ar1      ma1    sar1     sma1           0.9895  -0.8241  0.9974  -0.9551     s.e.  0.0053   0.0223  0.0018   0.0136       sigma^2 estimated as 0.07273:  log likelihood=-109.35        AIC=228.7   AICc=228.76   BIC=252.96      > forecast(model1, h=2) Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 t+1 **1.404053** 1.0584363 1.749670 0.8754777 1.932629 t+2 1.266133 0.9158214 1.616444 0.7303778 1.801888 However, considering the back errors and observed values from below, I cannot seem to replicate the forecast for t+1 which is 1.404053, using the formula:               x[t+1]=0.9895*x[t]+0.9974*x[t-6]-0.9895*0.9974*x[t-7]-(-0.8241)*e[t]-(-0.9551)*e[t-6]+(-0.8241)*(-0.9551)*e[t-7]      Instead of 1.404053 I get 1.34755 Point Errors Fitted Observed ... ... t-8 0.091543793 1.439935124 1.5314789 t-7 -0.146540485 1.40181299 1.2552725 t-6 -0.031449518 1.235569501 1.20412 t-5 -0.008707829 1.263980334 1.2552725 t-4 -0.009736316 1.472134314 1.462398 t-3 0.079273123 1.477029378 1.5563025 t-1 0.064970255 1.440179723 1.50515 t 0.135555386 1.444228211 1.5797836 Can you please help me understand what I am doing wrong? Thank you!


title: Can I use mixed ANOVA with 2 groups and 2 time points

text: I have 2 groups (biomaker high vs. low) and I measured quality of life (QoL) score (continuous variable) on 2 time-points. I want to run interaction statistics (group*time) to assess the effect of interaction on QoL score. My question is, can I use Mixed ANOVA on SPSS for this purpose? I have done it and it looks fine, but a colleague of mine mentioned that I can not use mixed ANOVA when it is 2 X 2 factors, explaining that Mixed ANOVA needs more than 2 time-points, e.g. 3 or 4 time-points. Any advice please?


title: Meaning of R's coxph output

text: I am getting below output from `coxph` function in R.               Concordance= 0.581  (se = 0.024 )     Rsquare= 0.06   (max possible= 0.995 )     Likelihood ratio test= 11.52  on 3 df,   p=0.009216      Can somebody explain what do `11.52` and `0.009216` mean in this output? I assume `11.52` is the prediction of overall survival and `0.009216` is the significance of that prediction, but I am not quite sure. Thanks!


title: Autoregressive distributed lag models ADL(p,q) determining amount of lags

text: I would like to know how I can determine the appropriate amount of lags in Matlab or another statistical package. I'm getting confused with VAR models and ARMAX models all the time and I'm a little stuck to be honest. Introduction on ADL models at wikipedia: http://en.wikipedia.org/wiki/Distributed_lag Example of my issue: I have one measure based on prices and one measure based on news. I expect a linear relationship between both but I expect the news based measure to be lagged to the price effect therefore I want to try to do a linear regression with the lagged news measures to the non-lagged price measure to check how the coherence (R-squared) changes. I expect the lag to be 3-5 days so that is what I want to check (also some autocorrelation could be present). I would like to do this without manually adding regressors (general- to-specific/specific-to-general) but with some model selection criteria for the number of lags if this is possible for the ADL(p,q). Still remains unclear to me how I should handle this.


title: R caret - how to combine principal components and non-principal-component predictors?

text: Consider the following example. Suppose I want to model `x` in `DF2` using everything else. While it makes sense to do `pca` pre-processing on variable 1 to 10, I want to keep `z` the factor variable as is and simply center and scale `y` and `v`. My first question is whether this approach makes sense at all.               DF = data.frame(x=rep(c("a","b","c"),each=3),                  y=c(1,3,6), v=1:9,                  z=rep(c('m','n','l'),3))     DF          i <- 1:10     s <- sapply(i, function(x) rep(x*5, 9))     s          DF2 <- cbind(DF, s)     DF2     str(DF2)      If it does make sense, how can I do it efficiently? I think one approach is to use `preProcess` on variable `1` to `10` and recreate the data frame using the resulting principal components and `y`, `v`, and `z`. But that doesn't seem very efficient. At the same time, It is often advised such as in this post that pre-processing is done during resampling in the `train` function. But I haven't figured out how to specify the different pre-processing methods for different predictors. Thanks in advance for any insight.


title: Detecting a significant change in the difference between two rates from T1 to T2

text: I asked this previously but I failed to activate my account so I was not able to respond to comments and my question was never answered. Here is attempt 2: I have two rates of compliance for two different populations:                Group      N         Compliant     Rate (Compliant/N)        A       696          155             22.27%        B       1058         442             41.78%      Assuming these responses conform to a Poisson distribution and the N-s do not change between time 1 and time 2, I want to know how to determine the amount of change in the difference score that is needed for that change to be considered statistically significant. So if the difference at the current time (T1) is 19.51%, what is the magnitude of the change needed to reach statistical significance at T2?


title: Is adding a function of time to a ARMA model still a linear time series model?

text: I wonder if the following model is a linear or nonlinear time series model: $$ X_t = d(t) + \sum_{i=1}^p a_i X_{t-i} + \sum_{i=0}^q b_i e_{t-i} $$ when $d(t) = ct$, $d(t) = e^t$, $d(t)=sin(t)$, and $d(t) = t^3$? According to the definitons in Tsay's Financial Time Series, I guess they are linear because $E(X_t|F_t), F_t=\{X_{t-1}, \dots, X_0\}$ is linear in $X_{t-1}, \dots, X_0$, and $Var(X_t |F_t)$ doesn't depend on $t$. But $X_t$ may be not linear in $t$, so I am not sure. p.s. There are two definitions from the book (I just used the second one), which I guess are equivalent? Definition 1: ![enter image description here](http://i.stack.imgur.com/OxWmR.png) ![enter image description here](http://i.stack.imgur.com/QIpCX.png) Definition 2: ![enter image description here](http://i.stack.imgur.com/klPVk.png)


title: Suitable data transformation in multiple GLS regression of raw variables against PCs

text: Im trying to make a multiple regresion in which the dependent variable and one of the factors were originally positive variables but were centered and scaled for easire interpretation of effect sizes. The other factors are principal components derived from a PCA procedure based on correlation and thus are also centered and scaled. The model has no interactions 7 predictors and 84 observations. The residuals of the test are not normal, and thus id like to transform something to achieve normality in the residuals. I wonder what would be the best approach for this since i cannot use sqrt or og because some of the data are negative. I consulted this post: Transforming variables for multiple regression in R But it seems I cant apply the power transformation to a GLS model made with the nmle package. Rank transformation of both factors and response variable would be acceptable in this case? I experimentally tried ranking everything and got similar results, except for one of the factors that appeared to be significant after ranking. It thus seems i did not lose power by ranking. Thanks in advance. Agus


title: Detecting heteroscedasticity - Can I use Breusch-Pagan Test on binary logistic regression?

text: I'm currently testing a (binary) logistic regression model, which seems to have at least some issues with multicollinearity. Now I don't really trust the data anymore and would like to also test it on heteroscedasticity. I found some information on Breusch-Pagan Test on the internet, but I could not find an answer to the question if this test also applys on Maximum-Likelihood- Methods, as it is usually mentioned in the context of OLS. So, can I apply the Breusch-Pagan Test on my model? Related to this question: How could I plot for heteroscedasticity-detection? I found this thread, but due to the binary nature of my dependent variable, the plot does not really work and unfortunately I'm pretty novice on plotting with R. Thanks in advance!


title: Predicted values for fixed effect quantile regression

text: I'm currently working with the method proposed by Koenker (2004) and Lamarche(2010) on fixed effects for quantile regression, for this I'm using the RQPD code in R. I would like to get the predicted values for each firm after running the RQPD code. The problem is that, after the running the code is only possible to retrieve the residuals. I know that the observed values of Y minus the residuals will give me the predicted values; but I need to run some simulations and for that I need to know how to get the predicted values "by hand" based on the β0 and β1 coefficients. I tried to get the predicted values as in the OLS: β0 + β1 (x) + the individual effect of the corresponding firm, using the following code:               require(rqpd) # it works only for version 2.15.3     set.seed(10)     m <- 3        n <- 10        s <- as.factor(rep(1:n,rep(m,n)))        x <- exp(rnorm(n*m))                X <- cbind(1,x)     u <- x*rnorm(m*n) + (1-x)*rf(m*n,3,3)     a <- rep(rnorm(n),rep(m,n))           y <- a + u     fit <- rqpd(y ~ x | s, panel(lambda = 1))     summary(fit)     fitcoe<-dcoef(fit)          res<-residuals(fit) # list of 100 residuals, I assumed the following:     res25<-res[1:30] # residuals for the 0.25 quantile     res50<-res[31:60] # residuals for the 0.50 quantile     res75<-res[61:90] # residuals for the 0.75 quantile     FE<-res[91:100] # individual effects for each firm          fitted1<-y-res25 # fitted values for the 0.25 quantile     fit1.1<-fitcoe[1,]+ fitcoe[2,]*x[1]+FE[1] # example for 1st observation      Unfortunately this is not working, for example for the first observation: when I substract the residuals from the observed values of Y, I get a predicted value of -0.8690 and when I calculate the predicted values by hand (fit1.1), I get a value of -1.2748. Do I miss something on my way to calculate the predicted values?


title: How to interpret basis function that yields vector in machine learning algorithm?

text: I'm struggling to understand what $\phi(x_{N+1})$ is in this excerpt of an algorithm (namely Linear Bayesian Regression embedded in other algorithm): $c_i = \gamma_i / \sum^L_{j} \gamma_j$ $V_i^{N+1} = ((V_i^N)^{-1}+\beta \phi(x_{N+1})^Tc_i \phi(x_{N+1}))^{-1}$ $\theta_i^{N+1} =V_i ^ {N+1} ((V_i^N)^{-1}\theta_i^N+\beta \phi(x_{N+1})^Tc_i r)$ The authors say: > "...where $\phi(x_i)$ denotes the feature vector of $x_i$ (...) The choice > of $\phi(x)$ depends on the task. In our experiments however, a Gaussian > basis function where the center is given by the augmented state $x_i$ proved > to be a good choice" I know what a Gaussian basis function is, but how does it yield a vector? I think that users with statistics and machine learning background will probably understand right away what $\phi(x)$ is here, but let me know if you need more context.


title: Clarification in information geometry

text: This question is concerned with the paper Differential Geometry of Curved Exponential Families-Curvatures and Information Loss by Amari. The text goes as follows. Let $S^n=\{p_{\theta}\}$ be an $n$-dimensional manifold of probability distributions with a coordinate system $\theta=(\theta_1,\dots,\theta_n)$, where $p_{\theta}(x)>0$ is assumed... We may regard every point $\theta$ of $S^n$ as carrying a function $\log p_{\theta}(x)$ of $x$... Let $T_{\theta}$ be the tangent space of $S^n$ at $\theta$, which is, roughly speaking, identified with a linearized version of a small neighborhood of $\theta$ in $S^n$. Let $e_i(\theta), i=1,\dots,n$ be the natural basis of $T_{\theta}$ associated with the coordinated system... Since each point $\theta$ of $S^n$ carries a function $\log p_{\theta}(x)$ of $x$, it is natural to regard $e_i(\theta)$ at $\theta$ as representing the function $$e_i(\theta)=\frac{\partial}{\partial\theta_i}\log p_{\theta}(x).$$ I don't understand the last statement. This appears in section 2 of the above mentioned paper. How's the basis of the tangent space are given by the above equation? It would be helpful if someone in this community familiar with this kind of material can help me understand this. Thanks. * * * **Update 1:** Although I agree that (from @aginensky) if $\frac{\partial}{\partial\theta_i}p_{\theta}$ are linearly independent then $\frac{\partial}{\partial\theta_i}\log p_{\theta}$ are also linearly independent, how these are members of the tangent space in the first place is not very clear. So how can $\frac{\partial}{\partial\theta_i}\log p_{\theta}$ be considered as basis for the tangent space. Any help is appreciated. **Update 2:** @aginensky: In his book Amari says the following: Let us consider the case where $S^n=\mathcal{P}(\mathcal{X})$, the set of all (strictly) positive probability measures on $\mathcal{X}=\{x_0,\dots,x_n\}$, where we regard $\mathcal{P}(\mathcal{X})$ as a subset of $\mathbb{R}^{\mathcal{X}}=\{X\big|X:\mathcal{X}\to \mathbb{R}\}$. In fact, $\mathcal{P}(\mathcal{X})$ is an open subset of the affine space $\{X\big |\sum_x X(x)=1\}$. Then the tangent space $T_p(S^n)$ of $S^n$ at every point can naturally be identified with the linear subspace $\mathcal{A}_0=\{X\big |\sum_x X(x)=0\}$. For the natural basis $\frac{\partial}{\partial\theta_i}$ of a coordiante system $\theta=(\theta_1,\dots,\theta_n)$, we have $(\frac{\partial}{\partial\theta_i})_{\theta}=\frac{\partial}{\partial\theta_i}p_{\theta}$. Next, let us take another embedding $p\mapsto \log p$, and identify $S^n$ with the subset $\log S^n:=\{\log p\big |p\in S^n\}$ of $\mathbb{R}^{\mathcal{X}}$. A tangent vector $X\in T_p(S^n)$ is then represented by the result of operating $X$ to $p\mapsto \log p$, which we denote by $X^{(e)}$. In particular we have $(\frac{\partial}{\partial\theta_i})_{\theta}^{(e)}=\frac{\partial}{\partial\theta_i}\log p_{\theta}$. It is obvious that $X^{(e)}=X(x)/p(x)$ and that $$T_p^{(e)}(S^n)=\{X^{(e)}\big |X\in T_p(S^n)\}=\{A\in \mathbb{R}^{\mathcal{X}}\big |\sum_x A(x)p(x)=0\}.$$ My question: If both $\frac{\partial}{\partial\theta_i}$ and $(\frac{\partial}{\partial\theta_i})^{(e)}$ are basis for the tangent space then would this not contradict the fact that $T_p$ and $T_p^{(e)}$ are distinct and $\frac{\partial}{\partial\theta_i}^{(e)}\in T_p^{(e)}$? I guess there seems to be an association between ($S^n,T_p$) and $(\log S^n,T_p^{(e)})$. If you can clarify this, it would be of great help. You may give it as an answer.


title: When conducting a t-test why would one prefer to assume (or test for) equal variances rather than always use a Welch approximation of the df?

text: It seems like when the assumption of homogeneity of variance is met that the results from a Welch adjusted t-test and a standard t-test are approximately the same. Why not simply always use the Welch adjusted t?


title: Monte Carlo Based A-Priori Power Calculation for Logistic Regression

text: I have been wanting to get a Monte Carlo based power calculation working for Logistic Regression cases for a little bit. I have put together a workflow with some assisstnace and I wanted to ask the community if this solution is on the right track and help to highlight any possible issues. I will list the code and then I will pose a few specific questions after the code. The code should be executable with no further inputs.               library(parallel)          # Define a sample experiment ----------------------------------------------          if(TRUE) {       X1 = c('L1', 'L2', 'L1', 'L1', 'L2', 'L2', 'L2', 'L1', 'L2', 'L2', 'L1',               'L1', 'L1', 'L2', 'L2', 'L2', 'L2', 'L1', 'L1', 'L1')       X2 = c(0, 50, 0, 100, 100, 50, 100, 50, 0, 0, 0, 50, 0, 0, 0, 100, 100,               100, 100, 100)       data = data.frame(X1, X2)       X = model.matrix(~ (X1 + X2)^2 + I(X2^2),data=data)       }               # Specify some helper functions -------------------------------------------          logit=function(p){       return(log(p/(1-p)))     }          ilogit=function(x){       return(exp(x)/(1+exp(x)))     }          workerfunction = function(c, X, pbar, n, i) {       # Simulate a series of responses       y = rbinom(n, 1, pbar)            # Execute the fit with this simulated data.       # I am using the glm.fit function since I have already generated       # the design matrix       # Question: Should this have the link function listed?       tfit = glm.fit(X, y, family = binomial())            # Extract the p-value for the model term in question       return(summary.glm(tfit)$coefficients[i, 4])     }          # Perform the power calculations ------------------------------------------          cat('Performing Power Calculations')          N = 1000     p0 = 0.7     ES = 0.2     alpha = 0.2          ES_Logit = logit(p0 + ES) - logit(p0)     NTerms = ncol(X)     n = nrow(X)          useParallel = FALSE          #tic()     if(useParallel) {       NWorkers = 8       cl = makeCluster(NWorkers, type = "PSOCK")     }          power = as.numeric(rep(NA, NTerms))     for(i in 2:NTerms) {            # Vector to store the information regarding the p-value of the model term under study       termpvals = as.numeric(rep(NA, N))            # Fit a straight line to the data so that we can determine what the p values should        # be at the different levels of the factor of interest         m = ES_Logit / (max(X[,i]) - min(X[,i]))       b = logit(p0) - (m * min(X[,i]))            # Fix all of the model term coefficients to zero with the exception of the       # model term of interest.  betas[1] should be the y-intercept and betas[x]       # should be the slope that we determined before       betas =  c(b, rep(0, NTerms-1))       betas[i] = m            # Determine the correct probabilities for each of the different levels of       # the model term of interest       pbar = ilogit(X %*% betas)            # Loop over each of the different pseudo-experiments       # The call to the workerfunction is strange here.  The first argument to the worker is supposed       # to be the element of the array that we are looking at, but we dont care.  The workerfunction       # function ignores the first argument.  The remaining arguments to the worker function are        # specified in the Xapply function are the FUN argument (X, pbar, n, i)       if(useParallel) {         termpvals = parLapply(cl, 1:N, FUN = workerfunction, X, pbar, n, i)       } else {           termpvals = lapply(1:N, FUN = workerfunction, X, pbar, n, i)       }            # Remove any NA values from the termpval vector       termpvals = termpvals[!is.na(termpvals)]            # Calculate the power as the fraction of elements in the termpval that is less than alpha       power[i] = sum(termpvals < alpha) / length(termpvals)     }          if(useParallel) {       stopCluster(cl)     }     #toc()          # Display the power calculation results -----------------------------------          cat('Power Calculation Results:')     for(i in 1:ncol(X)) {       cat(sprintf('  %30s: %.2f\n', colnames(X)[i], power[i]))     }      While the code above functions and gives me results, I want to ensure that they are good defensible results. The specific questions besides "Is this correct" are:   1. In the workerfunction function, should the call to the glm.fit include the binomial family argument. I believe that it should since we are fitting binomial data but I wanted to double check.   2. Does this solution effectively handle the difference between categorical and continuous quantities?   3. In the event that I have a categorical factor with more than three levels, how do I determine power for that case? I lean on comparing a model with the dummy terms included and without the dummy terms included.


title: Combining multiple probability distributions

text: Let say I have a 2 samples from the same population that is unknown               Sample 1: n=100 mean=50 range= 0-100          Sample 2: n=40 mean= 70 range=50-140      How would I combine these two data sets to create a sample that better reflects the population for calculating the probability of events. Would I use something like a Gaussian mixture model? I am guessing GMM from this video. But I am not sure this is what I would use


title: Calculating Probability

text: Hi looking for advice on how to calculate the probability of the following. If someone has a 4% chance of catching a disease (so 4 people out of 100) Then if they caught the disease they have a 15% chance of it not showing up with tests, so 15 out of 100 who have caught the disease won't be diagnosed with the diseased. What is the overall probability of this happening combined? Regards, Mike


title: What algorithm should I use to cluster a huge binary dataset into few categories?

text: I have a large (650K rows * 62 columns) matrix of binary data (0-1 entries only). The matrix is mostly sparse: about 8% is filled. I would like to cluster it into 5 groups - say named from 1 to 5. I have tried hierarchical clustering and it was not able to handle the size. I have also used hamming distance based k-means clustering algorithm, considering the 650K bit vectors of length 62. I did not get proper results with any of these. Please help.


title: Convert a categorical variable to a numerical variable prior to regression

text: I am doing a project to estimate students' final graduation GPAs based on several variables. I have students' first year GPAs, high school GPAs, their race, where they come from, and their ACT score, and so on. I have two questions:   1. How to convert race into numbers, I know I can just assign white to be 1, Black to be 2, Asian to be 3, but it may cause some problem that make my result not significant, so how do I convert the race into numbers to make my model more accurate?   2. How do I find which factor make the most contribution to estimate students final GPA, so I can put more weight on it?


title: Repeated measures analysis... or not?

text: Not sure if dealing with repeated measures analysis...? Measuring X at 3 seasons (summer, winter, spring/fall) during a year – might be high summer/low winter, high winter/low summer, consistently high, consistently low, etc. Y (health outcome) measured at same 3 times. Time between measures not fully equal between participants. Y may or may not change within the short time frame. Main question is whether seasonal variation in X is associated with Y. ie. would someone with consistent-X in all seasons have healthier Y vs. someone with more variability in X? Might just be multiple linear/logistic regression… or is it repeated measures, or ANCOVA…? And then what if in some cases X appears to increase (or decrease) in a season but it's not due to season, ie. rather a permanent change over time... which might lead to a change in Y? But real interest is the association between seasonal variation in X and health outcome Y.


title: Forecasting daily data with trend, yearly, day of the week, and moving holiday effects

text: I'm expanding a question I posed earlier because I think it was lacking detail. I'm attempting to forecast daily demand for a restaurant that sells take away food, primarily to office workers on their lunch breaks. They are located in the downtown core of a major city. They are only open on workdays - no holidays, no weekends. I'm familiar with models that take into account seasonality and trend - Holt-Winters triple exponential smoothing, for example. I'm also familiar with models that take into account complex seasonality and trend - the TBATS package for R, for example. My problem is that I've identified 8 components that determine sales on a given day:   1. The yearly seasonal component. Sales are lower in the summer, for example, when many office workers are on vacation.   2. The weekly component. Sales very obviously peak on Thursdays (in the absence of other effects - see below)   3. The _Friday-Long-Weekend-Effect_. The weekly sales pattern changes if the coming Friday is a holiday. Wednesday will typically have higher sales, for example.   4. The _Post-Friday-Long-Weekend-Effect_. The weekly sales pattern changes if the week before was shortened due to the Friday being a holiday.   5. The _Monday-Long-Weekend-Effect_. The weekly sales pattern changes in week $t$ if the Monday in week $t+1$ is a holiday. For example, sales are much lower on Fridays preceding Monday-Long-Weekends. Presumably people are leaving the office early and skipping lunch.   6. The _Post-Monday-Long-Weekend-Effect_. The weekly sales pattern changes if the week is shortened due to the Monday being a holiday.   7. The trend component.    8. The noise component. If holidays fell on the same date every year, then the "long-weekend-effects" would be captured in the yearly seasonal component. However, they don't. My first thought was to include dummy variables. For example, let $X_{M+1}$ be the "Monday-long-weekend-effect" component, and $\beta_{M+1}$ be the associated coefficient, for a given day. Then for the Friday preceding a Monday-Long-Weekend, $X_{M+1}=1$, and for a Friday not preceding a Monday- Long-Weekend, $X=0$. I'm only using three years of data, so it would be easy for me to change the $X_{M+1}$ values to 0 or 1 by hand for each year. However, I don't know how to include such dummy variables in models like those that I've mentioned. Any input as to a model that can take into account the components I've mentioned would be greatly appreciated. It seems like I need to capture moving-holiday-effects, day-of-the-week effects, seasonal patterns, and trend, all in one. **Question: Is there a model I can use that can be implemented in R and take into account the components I've listed?** _My background: I'm a forth year mathematics and economics student. I've also taken statistics classes, and I'm using R to perform my analysis. This is for a final report for a forth year data analysis class._


title: Transformations of Variables

text: Suppose that you are have a response variable $Y$ and explanatory variables $X_1$, $X_2$ and $X_3$. If we want to use a quadratic transformation for $X_1$, would we still include $X_1$? In other words, would we have: $$E[Y|X] = \beta_0+\beta_{1}X_{1}^{2} + \beta_{2}X_{2} + \beta_{3}X_{3}$$ or $$E[Y|X] = \beta_0+\beta_{1}X_{1} + \beta_{2}X_{1}^{2} + \beta_{3}X_{2}+\beta_{4}X_{3}$$ If instead we did a logarithmic transformation, then would it just be: $$E[Y|X] = \beta_0+\beta_{1} \log(X_1) + \beta_{2}X_{2} + \beta_{3}X_{3}$$


title: Estimate just the constant coefficient in logistic regression

text: How do I calculate the constant coefficient in logistic regression manually, i.e without having to use a calculator? My model is $g(Y) = X \beta + \alpha$ is it possible to calculate just the constant parameter $\alpha$ without performing the full regression fit?


title: Generating rules to obtain a given categorical distribution... is it possible?

text: I'm working on a problem, I was wondering if there are any methods available to do the following. I have a data set with information on people (continuous and categorical data). I have 3 categories to assign to these people {A, B, C}, and I want to try and develop some rules of how to assign each person a category. However the tricky part is, the 'target' I'm looking for is a distribution. So for example I would like to be able to pass in {.2, .05, .75} as the distribution of {A, B, C} and then have an algorithm try out different rules and splits of the data (similar to a decision tree) until it achieves the given distribution. Are there any machine learning type methods suited to this problem? Is there a name for this type of problem? I wasn't finding anything with my googling... Thanks edit: Just for a little more information, the way I'm currently going to approach it is random assignment, give everyone a random number, if it's below .2 then give Category A and so on. Not ideal, but the fact that it's random is a major issue for the problem. Although I would much prefer to have a set of rules that I could use instead.