title: When is R squared negative?

text: My understanding is that R squared cannot be negative as it is the square of R. However I ran a simple linear regression in SPSS with a single independent variable and a dependent variable. My SPSS output give me a negative value for R-squared. If I was to calculate this by hand from R then R squared would be positive. What has SPSS done to calculate this as negative?               R=-.395     R squared =-.156     B (un-standardized)=-1261.611      Code I've used:               DATASET ACTIVATE DataSet1.      REGRESSION /MISSING LISTWISE /STATISTICS COEFF OUTS R ANOVA                 /CRITERIA=PIN(.05) POUT(.10) /NOORIGIN                 /DEPENDENT valueP /METHOD=ENTER ageP      I get a negative value. Can anyone explain what this means? Thanks. ![Negative RSquared](http://i.stack.imgur.com/fCCDi.png) ![enter image description here](http://i.stack.imgur.com/h1cBS.png)


title: How to decide bootstrap number of runs?

text: I am using bootstrap for my simulation. The number of the population is flexible for each case, and the sample size is decided by a certain percentage. For example, I have a 10,000 population, and I decide to use 10% for each iteration of bootstrap, so the sample size is 1,000. In practice, I found it is hard to decide how many times to run the bootstrap is enough. With less simulation, the results appear insufficiant, while with a large number of simulation they are purely redundant. May I know if there is a method that can help me to decide the number of iterations to run?


title: Effect of missing data and outliers on least square estimation

text: Why is it that "missing data" and "outliers" can affect the performance of least square estimation?


title: Number of trials required from a binomial distribution to get the desired odds

text: This question may have been asked before, but I couldn't find it. So, here goes. From about 3000 data points that can be characterized as "wins" or "losses" (binomial), it turns out that there are 52.8% dumb luck wins. This is my dependent variable. I also have some additional data that may help in predicting the above wins that could be considered an independent variable. The question I'm trying to answer is: If my independent variable can be used to predict 55% wins, how many trials are required (giving me 55% wins) for me to be 99% sure that this wasn't dumb luck? The following R code is purposely hacky so I can see everything that is happening.               #Run through a set of trial sizes     numtri <- seq(2720, 2840, 1)          #For a 52.8% probability of dumb luck wins, and the current trial size,     #calculate the number of wins at the 99% limit     numwin <- qbinom(p=0.99, size=numtri, prob=0.528)          #Divide the number of wins at the 99% limit by the trial size to     #get the percent wins.     perwin <- numwin/numtri          #Plot the percent wins versus the trial size, looking for the     #"predicted" 55% value. See the plot below.     plot(numtri, perwin, type="b", main="Looking for 0.55")     grid()          #Draw a red line at the 55% level to show which trial sizes are correct     abline(h=0.55, lwd=2, col="red")          #More than one answer? Why?........Integer issues     head(numtri)     head(numwin)     head(perwin)      From the graph, the answer is: 2740 <= numtri <= 2820 As you can guess, I'm also looking for the required number of trials for 56% wins, 57% wins, 58% wins, etc. So, I'll be automating the process. Back to my question. Does anyone see a problem with the code, and if not, does anyone have a way to cleanly sneak up on the right and left edges of the "answer"? ![enter image description here](http://i.stack.imgur.com/2CSHV.jpg) Edit (02/13/2011) ================================================== Per whuber's answer below, I now realize that to compare my alternate 55% wins (which came from the data) with my 52.8% "dumb luck" null (which also came from the data), I have to deal with the fuzz from both measured values. In other words, to be "99% sure" (whatever that means), the 1% tail of BOTH proportions needs to be compared. For me to get comfortable with whuber's formula for N, I had to reproduce his algebra. In addition, because I know that I'll use this framework for more complicated problems, I bootstrapped the calculations. Below are some of the results. The upper left graph was produced during the bootstrap process where N (the number of trials) was at 9000. As you can see, the 1% tail for the null extends further into the alternate than the 1% for the alternate. The conclusion? 9000 trials is not enough. The lower left graph is also during the bootstrap process where N was at 13,000. In this case, the 1% tail for the null falls into an area of the alternate that is less than the required 1% value. The conclusion? 13,000 trials is more than enough (too many) for the 1% level. The upper right graph is where N=11109 (whuber's calculated value) and the 1% tail of the null extends the right amount into the alternate, aligning the 1% tails. The conclusion? 11109 trials are required for the 1% significance level (I am now comfortable with whuber's formula). The lower right graph is where I used whuber's formula for N and varied the alternate q (at both the 1% and 5% significance levels) so I could have a reference of what might be an "acceptable zone" for alternates. In addition, I overlaid some actual alternates versus their associated number of data points (the blue dots). One point falls well below the 5% level, and even though it might provide 72% "wins", there simply aren't enough data points to differentiate it from 52.8% "dumb luck". So, I'll probably drop that point and the two others below 5% (and the independent variables that were used to select those points). ![enter image description here](http://i.stack.imgur.com/sEsZx.jpg)


title: Why we always put log() before the joint pdf when we use MLE(Maximum likelihood Estimation)?

text: Maybe this question is simple, but I really need some help. When we use the Maximum Likelihood Estimation(MLE) to estimate the parameters, why we always put the log() before the joint density? To use the sum in place of product? But why? The wikipedia said it will be convenient. Why? Thank you.


title: Why would SVD be 'unstable' if you don't standardize your data first?

text: I'm reading an article about Direct Linear Transformation which processes data using _SVD_ , and the data set is standardized so that it has _zero_ mean and _unit_ standard deviation (n.b., some people call this transformation 'normalized'). In that article, it is said that for **numerical reasons** , the data set has to be normalized to make sure the DLT algorithm stable. I wonder what are those numerical reasons? Does it mean if I have 2 data sets of the same objects, if no standardization is done, the DLT algorithm will give out 2 numerically different results? **UPDATE** In computer vision, Direct Linear Transform (DLT for short) is a method of homography estimation, it solves the overdetermined linear system via SVD $$Ah=b$$to find a solution $h$ under constraint $\|h\|=1$. Actually it finds the least square solution which minimize $\|Ah - b\|$. I understand the basic idea of this algorithm, but it is recommended to **standardize the data set** before applying DLT on it, and here is a intro about how to do the standardization. It is lectured that data standardization is important to DLT, without standardization the results from DLT are not stable, and the reason is that DLT is quite dependent on the coordinates of the points in the data set. I wonder why? Just because DLT involves solving the linear system using SVD and $A$ might be singular?


title: Improved Monte-Carlo method vs. hit-and-miss method

text: I do not understand _Which is more accurate, the hit-and-miss method or the improved Monte-Carlo method?_ Here it is written that that the hit-and-miss has a higher variance but they showed `(variance of improved Monte-Carlo) - (variance of hit-and-miss method) > 0` $$\sigma^2_M-\sigma^2_H>0$$ Doesn't it imply improved Monte-Carlo method has a higher variance? Can you please explain it experimentally to me with an easy example (such as for a small integration)?


title: Degrees of freedom for a weighted average

text: And yet another for the long list of degrees-of-freedom questions! Given an i.i.d. sample $x_1,..., x_n$ from an arbitrary real-valued distribution with expectation $\mu$, the sample mean can be written as $$\bar x = \sum_i \frac{1}{n} x_i.$$ One can say that the degrees of freedom for $\bar x$ is $n-1$, as there are $n$ data points and one parameter being estimated. Now suppose that some twisted statistician decides that the weights $1/n$ are not good enough, and replaces $\bar x$ with $$\bar x_w = \sum_i w_i x_i,$$ where $w = (w_1, ..., w_n)$ is a pre-determined (i.e., not related to the data) vector of nonnegative weights with $\sum_i w_i = 1$. What is the degrees of freedom for $\bar x_w$? References to versions of this problem in peer-reviewed work would be of particular interest.


title: Neural Nets: computing weight deltas is a bottleneck

text: After profiling my Neural Nets' code I've realized that the method, which computes the weight changes for each arc in the network (`-rate*gradient + momentum*previous_delta - decay*rate*weight`), already given the gradient, is the bottleneck (55% inclusive samples). Is there any trick to compute these values in a efficient manner?


title: Negative coefficient in model stacking

text: I've some data I'm modeling with a couple of different approaches. Random Forests on one side, a logistic regression on the other side. I wanted to create an ensemble model using the predicted probabilities from the RF and regression models. Using cross validation I'm running a third model (a logistic regression) that uses the predicted probabilities from the first two models as inputs. In the examples I've seen on Kaggle the coefficients tend to be positive, and are interpreted as weights of the ensemble model. However one of my two coefficients ends up being negative. What's the intuition to explain that? And what is it telling me from a practical perspective? My hunch is that it is telling me that the model with the negative coefficient is actually detrimental to the ensemble model, but it is just a guess. Thanks


title: Analysis of a time series with a fixed and random factor in R

text: I have a dataset which I am not sure how to analyse. The dataset came from the following experiment: I grew plants (2 different types) and measured their height at different time point (each plant was followed individually). I had in total, 3 boxes in which I grew the plants and, in each box, 3 plants of each type. I took the measurement at 4 different time points. So, if I am not wrong, the plant type is a fixed factor and the box a random one. Here is how it is structured:               time, box, type, height     1, 1, 1, 1.2     1, 1, 1, 1.3     1, 1, 1, 1.1     1, 1, 2, 1.4      1, 1, 2, 1.5     1, 1, 2, 1.6     ...     2, 1, 1, 1.2     2, 1, 1, 1.3     ...     1, 2, 1, 1.2     1, 2, 1, 1.3          I would like to know what is the correct way to check if there is a **difference between the different plant types** using R So far, what I have done is this:               lme1 <- lme(height ~ type, random= ~ 1|box, data=mydata)     anova(lme1)      but I do not know how to include the time variable in the analysis.. Here is the plot of the height evolution with time, for the different plant type. Each line is a plant. ![enter image description here](http://i.stack.imgur.com/5R27I.png)


title: Dealing with dropping non-responders in clinical trials

text: I'm wondering what are some good approaches for analyzing trials where non- responders are dropped out at some set point. This is not my field, so I may be unaware of standard things! Let's say there are 100 patients, divided 50/50 into treatment and control groups (given by an indicator variable Trt). We are measuring a normally- distributed outcome Y. The trial goes for 10 weeks (with a baseline Y taken at week 0), then the patients whose scores have not improved by (say) at least 20% (considered non-responders) are removed from the study (to seek other treatments) and are no longer followed. The trial then continues for another 10 weeks. For simplicity, assume that the only missing data is the outcomes for the non-responders after week 10. This is plugged into a longitudinal model with random slopes and intercepts and Trt as a level-two predictor for slope (since the baseline measure shouldn't depend on group assignment). I've seen the following approaches:   1. Ignore the fact that the missing data is not missing at random.   2. Carry forward the final measurement for the non-responders.   3. Do multiple imputation to complete the data set for the non-responders. My understanding is that no one likes 2, and 3 is okay up to standard caveats about imputation. The problem with 1 seems subtler, and I'm not sure how to properly perform inference in that situation. My questions are: is there an acceptable way to do Option 1? Are there other approaches? For example, would it be feasible to include nonresponsiveness in the model (maybe using BUGS)? Any references would be greatly appreciated!


title: When simulating a bivariate normal distribution, why is $\rho$ chosen instead of estimated from the data?

text: In a video lecture, MrProf shows the 3d-plot of a bivariate normal distribution $\mu_{x_1} = \mu_{x_2} = \sigma_1 = \sigma_2 = 1$ and chooses $\rho = 0.5$ . If stick to Mathworld, $\rho$ simply is the Pearson correlation that is $\dfrac{\mbox{cov}(x_1, x_2) }{[ \mbox{var}(x_1) \mbox{var}(x_2)]^\frac{1}{2}}$. I don't understand why one would have control over $\rho$ as in the example given by MrProf. I understood that $\rho$ is an output from the formula as I plug in $x_1$ and $x_2$, not an input. Can anyone throw any light?


title: Why does a finite, irreducible and aperiodic Markov chain with a doubly-stochastic matrix P have a uniform limiting distribution?

text: The theorem is "If a transition matrix for an irreducible Markov chain with a ﬁnite state space S is doubly stochastic, its (unique) invariant measure is uniform over S." If a Markov Chain has a doubly-stochastic transition matrix, I read that its limiting probabilities make up the uniform distribution, but I do not quite understand why. I have been trying to come up with, and locate, an understandable proof for this. But the proofs I find all gloss over details I don't understand, like proposition 15.5 here (why does it work to just use the [1,...1] vectors?) Could someone point me to (or write) a more simple/detailed proof? (Though not part of anything I will hand in at school, it is part of a course I'm taking so I guess I'll tag it with homework in either case.)


title: Cointegration between two stationary processes

text: I have two stationary time series. I would like to check for cointegration between them. Does this make sense, and can I just use Engle-Granger Test (two step) for Cointegration for this?


title: How do I estimate and interpret a forecast made by Minitab?

text: From 43 years of monthly stream flow data, I want to get 11 years ARIMA future forecast using Minitab 15 statistical Software. The difficulty is the procedure to be followed, interpretation of model outputs (Final Estimates of Parameters and Modified Box-Pierce (Ljung-Box) Chi-Square statistics) and other best fit model selection criteria like AIC, BIC etc. Any one who can help me is appreciated. The raw data:               2.76     2.11     1.7     1.25     1.25     2.14     19.27     42.97     21.54     6.95     3.68     2.87     2.18     1.83     1.57     1.38     1.08     2.17     24.52     33.48     24.86     9.86     7.55     6.05     2.61     1.71     1.56     1.01     1.21     2.13     18.08     31.41     22.5     8.18     4.3     2.97     2.2     1.54     1.46     1.35     1.8     2.84     18.33     34.43     20.46     6.25     7.72     4.29     2.43     1.49     1.13     1.2     1.43     2.94     42.58     51.32     30.71     13     6.12     3.95     2.68     2.06     2.01     1.71     1.26     2.37     11.52     26.97     14.76     10.5     7.43     4.33     2.16     1.62     1.55     1.16     1.33     2.16     16.43     25.17     17.43     4.81     3.39     2.29     1.73     1.38     1.46     1.15     1.06     3.11     24.04     35.75     23.49     9.82     4.44     3.13     2.54     1.77     1.7     1.35     1.42     4.65     32.49     27.43     19.49     5.02     4.86     5.02     2.35     1.76     1.8     1.36     1.39     3.04     31.74     51.66     27.06     6.78     4.06     3.21     2.03     1.71     1.66     1.45     2.4     5.38     22.52     38.41     21.69     7.58     3.69     2.64     2.02     1.55     1.36     1.15     1.22     3.5     18.3     41.48     21.6     6.71     5.08     2.85     2.19     1.41     1.28     1.07     1.11     2.73     10.74     21.15     14.17     4.18     2.87     1.75     1.34     0.98     0.89     0.79     1.12     2.47     13.34     46.49     17.1     7.07     3.34     2.32     1.84     1.3     1.21     0.96     2.23     3.18     25.22     35.36     27.7     8.21     3.94     2.7     2.17     2.11     1.61     1.23     1.26     3.9     49.69     96.27     67.47     17.04     6.28     6.01     3.48     3.36     3.06     2.5     3.52     14.98     28.34     44.81     24.12     7.49     5.12     3.13     2.29     1.55     1.47     1.09     1     2.88     19.7     24.31     15.24     14.28     4.27     2.53     2.42     1.63     1.43     1.46     1.63     8.79     33.17     32.96     26.45     8.96     4.56     3.43     2.49     1.7     1.48     1.17     1.53     3.3     14.07     29.24     19.26     7.79     3.14     2.29     1.65     1.27     1.16     1.15     1.16     7.56     34.35     34.35     12.84     6.6     2.8     2.06     1.73     1.39     1.43     0.96     1.81     4.73     42.59     35.78     26.99     7.51     2.64     1.97     1.79     1.52     1.73     1.42     1.67     6.25     26.93     32.81     21.82     7.99     2.6     1.99     1.54     1.11     1.12     0.85     0.89     3.56     5.86     39.91     14.07     6.66     2.89     1.91     1.36     0.89     0.89     0.78     0.9     10.72     22.42     21.67     17.46     4.62     2.55     2.16     1.37     0.93     0.84     0.75     1.01     3.84     15.71     32.72     16.24     5.93     3.36     2.1     1.49     0.94     0.85     0.65     0.55     3.57     14.13     15.95     15.58     7.45     2.85     2.19     1.57     1.21     1.71     0.89     2.33     9.44     15.77     19.92     12.85     6.6     3.56     2.41     1.99     1.3     1.17     0.7     1.02     6.34     40.25     35.71     26.34     15.6     5.03     3.17     2.41     1.55     1.47     1.19     1.41     7.93     27.54     32.99     16.63     6.4     3.09     1.87     1.43     1.46     1.35     1.06     1.13     2.08     18.02     56.17     21.99     5.82     2.95     2.09     1.62     1.11     1.02     1.28     1.35     4.47     31.81     55.78     26.49     6.13     4.73     3.49     2.55     1.68     1.55     1.49     1.46     3.17     16.67     40.25     20.07     14.23     7.77     4.15     3.6     2.33     1.8     1.74     2.24     7.87     32.87     25.45     24.42     11.17     5.4     3.32     2.32     1.62     1.48     1.22     1.46     6.09     17.88     23.93     21.77     5.89     4.06     2.85     1.9     1.33     1.32     1.17     1.65     5.15     17.28     29.62     21.1     5.25     3.04     2.29     1.71     1.09     1.21     1.55     3.31     14.33     34.6     42.93     24.74     9.36     4.67     3.3     2.47     1.7     1.63     1.27     2.76     8.08     20.12     28.48     14.6     11.18     6.93     3.42     2.33     1.5     1.33     1.09     2.52     7.16     23.86     34.65     25.07     19.58     6.41     4.27     3.06     1.95     1.62     1.35     2.24     10.19     38.02     33.4     20.03     19.02     5.96     4.27     2.61     1.72     1.52     1.88     1.85     5.4     18.6     53.81     19.31     23.97     8.04     3.92     2.6     1.92     1.79     1.4     1.99     11.96     37.47     54.03     17.56     7.41     4.1     3.06     2.38     1.67     1.64     1.29     1.21     5.55     19.57     32.44     15.63     6.24     3.45     2.11


title: When someone says residual deviance/df should ~ 1 for a Poisson model, how approximate is approximate?

text: I've often seen the advice for checking whether or not a Poisson model fit is over-dispersed involving dividing the residual deviance by the degrees of freedom. The resulting ratio should "approximately 1". The question is what range are we talking about for "approximate" - what's a ratio that should set off alarms to go consider alternative model forms?


title: Using z-standardization to account for covariate

text: I would like to know whether z-standardization is an appropriate way to account for a covariate. Please consider the following dummy example (I am not interested in the interpretation of the result but only in the statistical plausibility): In a sample of 100 subjects relationship between IQ and head size is investigated (one measurement per subject). However IQ has been estimated with IQ Test 'A' in some subjects and with IQ test 'B' in the other subjects. Also there is a significant difference in IQ measured by 'A' as compared to IQ measured by 'B'. Therefore 'zIQ' is obtained by z-standardization of all IQ values obtained by test 'A' as well as z-standardization of all IQ values obtained by 'B'. Is it appropriate to use zIQ to investigate its relationship with head size independently of which test ('A' or 'B') has been used?


title: If $cor(X,\epsilon) \approx 0$ in linear regression, can we conclude $X$ is exogenous?

text: Suppose that we run the simple linear regression $Y = \alpha + \beta X + \epsilon$. I want to test whether the independent variable $X$ is exogenous. If the correlation between the independent variable $X$ and the residual of linear regression $\epsilon$ is almost zero, i.e. $cor(X, \epsilon) \approx 0$, can I then conclude that this simple test suggests that the independent variable $X$ is exogenous?


title: Which t-test to use for a two-group pre- post-test design?

text: I have data on 26 participants (13 from computing and remaining 13 from non computing) who have participated in my research. Each participant is treated with a lab module (Hands on Robotics Session). Now each participant will be evaluated using a rubric on scale of 1 to 4. This experiment has both pre and post test. For my research i want to evaluate the following questions: ### Research question 1   * Null Hypothesis: students do not learn about computational thinking (programming basics and algorithmic thinking) with the help of robotics.   * Alternate Hypothesis: Students learn about computational thinking (programming basics and algorithmic thinking) with the help of robotics. To evaluate the above question, the categories i will be considering are Plan, Implementation and Knowledge gained on a scale of Excellent, Good, Fair and Poor. I think i should use DEPENDENT T-TEST FOR PAIRED SAMPLES. Am i correct? ### Research question 2:   * Null Hypothesis: Participants background (computing and non computing) has no effect in learning about algorithmic thinking with help of robotics   * Alternate Hypothesis:Participants background (computing and non computing) has an effect in learning about algorithmic thinking with help of robotics  I will also evaluate question 2 on a scale of Excellent, Good, Fair and Poor, but with respect to background. ### My Statistical Question **Which T test should I use for each of my research questions?**


title: Funnel plot with asymmetrical point sizes

text: Lets say I want to take the usual funnel plot for meta-analysis, and add another dimension to it, visually changing the points used for each study by a covariate. While it might be easier to change the marker or color for categorical variables, for a continuous variable, this gets a little harder. Lets say we want to see if there's obvious asymmetry not just overall, but by say, length of follow-up. Does anyone know a way in either R or Stata to either change the size of the points plotted - so you essentially have a funnel-bubble plot, or to color them on a continuous gradient? Clearly, this can be done by just modifying a scatterplot, but I'd like to save the steps in manipulating the access to get it looking like a conventional funnel plot.


title: Suitable model for predicting range-bound integer attribute?

text: I need to predict an integer variable in the scale (0-50). I am wondering how it should be modeled:   1. Should zero be predicted as a 0/non-0 categorical classification separately, given that observations are mostly zero-inflated?   2. If I use Linear Regression for the range 1-50, how can I make the prediction variable range-bound (1-50)? I am using scikit-learn but cannot see anything that allows for such a thing.


title: Bayes factor calculation

text: I am a newbie to Bayesian stats...   I came across an article for calculating Bayes factor by rounder etal.. In this article how am I supposed to put t value..n scale r on effect size parameters. If you could help me with a specific example it would be a great help.. I have 2 groups in my sample               group A sample size 110, mean 112.76, SD 30     group B sample size 89, mean 100.4, SD 28      What are `t value` ( is it by routine t test ?) and `scale effect` for `r size` in this case..and how can I calculate them Further how do we interpret the output..accompany article is a bit inaccessible to me because of my non mathematical ground...if anyone could explain..how can I interpret these ratios..?


title: Predictive model for network data

text: Assume a network as a set of data, which are defined by their coordination $(x,y,z)$ and a weight on its edge. Now this data can be used as an input data to predict a single value. In my case, coordination of nodes are mostly fixed, and just by changing subset of node I would like to predict the output. For example, in protein structure, coordination of each amino acid is given (nodes) and atomic bond can be calculated based on physical property and distance of each amino acid,and we would like to see by changing couple of amino acids (either by displacement or replacing by some other amino acid) predict the pH or solubility of the protein which are single values. Also, in gene-gene interaction network, one can do more or less same thing by changing some genes, predict the drug performance and side effect. I appreciate any keyword or literature for such a method.


title: Is significance more important than absolute value?

text: I have three groups and I'm comparing two graphs in each group using the t-test: the absolute difference between two graphs in the first group is 34 with a p=0.024. (*) The absolute difference in the second group is 28 with a p=0.04 (*) and in the third group the absolute difference is 4 but the p value is=0.010 (**) How do I explain this data? The absolute value in the first and second group is higher; however, the significance is less than the comparison in group three where the absolute value difference between two groups is only 4. So, should I say there was more change in group three (**) than in group 1 and 2 although the absolute difference (4) is less in this group? Thank you.


title: Refugee from SPSS having issues with fit.contrast in R

text: I'm trying to wean myself off of SPSS and switch to R for psych statistics and I'm mostly getting there, but I've run into a real roadblock with `fit.contrast` from the gmodels package. As far as I can tell, even the provided example doesn't work properly--either that, or I grossly misunderstand contrasts in R (always a possibility). First let's look at the example code that comes with the documentation for `gmodels::fit.contrast()`:               set.seed(03215)     Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)     Time <- factor(sample(1:3, 1000, replace=TRUE))     y <- rnorm(1000)     data <- data.frame(y, Genotype, Time)          model <- aov( y ~ Genotype + Time + Genotype:Time, data=data )     model.tables(model, "means")          fit.contrast( model, "Genotype", rbind("KO vs WT"=c(-1,1) ), conf=0.95 , df=T)      For me, this produces the following output:                                  Estimate Std. Error   t value Pr(>|t|)   lower CI  upper CI     GenotypeKO vs WT 0.01683876  0.1095764 0.1536714   0.8779 -0.1981888 0.2318664      Now consider the later example code:               model <- aov( y ~ Genotype + Time + Genotype:Time, data=data,                   contrasts=list(Genotype=make.contrasts(cm.G),                                  Time=make.contrasts(cm.T) )     )          summary(model, split=list( Genotype=list( "KO vs WT"=1 ),                                Time = list( "1 vs 2" = 1,                                             "2 vs 3" = 2 ) ) )      This produces this output for me:                                                 Df Sum Sq Mean Sq F value Pr(>F)     Genotype                           1    1.2  1.1687   1.121  0.290       Genotype: KO vs WT               1    1.2  1.1687   1.121  0.290     Time                               2    0.7  0.3677   0.353  0.703       Time: 1 vs 2                     1    0.2  0.1784   0.171  0.679       Time: 2 vs 3                     1    0.6  0.5571   0.534  0.465     Genotype:Time                      2    1.2  0.5760   0.552  0.576       Genotype:Time: KO vs WT.1 vs 2   1    0.3  0.3265   0.313  0.576       Genotype:Time: KO vs WT.2 vs 3   1    0.8  0.8256   0.792  0.374     Residuals                        994 1036.6  1.0429        I'm completely flummoxed here. The p-values for the two tests on Genotype should match, but they don't! The p-values are not the same and the F-value isn't the square of the t-value. Also, the value computed for the contrast estimate (0.0168) is NOT the difference in the two means... why not? Any help in understanding what fit.contrast is doing here would be greatly appreciated!


title: Standard error and p-values of exponential smoothing weights

text: Is there any justfification for producing a standard error of a single exponentially weighted coefficient? If yes, how can we interpret the p-value? * * * **Background** I use SAS ETS to estimate a single exponential smoothing model. When i select the optimize option for the estimation of the smoothing weight I get the estimated value (optimal), a standard error, a t-value and a p-value. How is the standard error calculated (e.g. what is the formula for it) and what is the menaing of the p-value. In the users guide it says that: > The standard errors associated with the smoothing weights are calculated > from the Hessian matrix of the sum of squared, one-step-ahead prediction > errors with respect to the smoothing weights used in the optimization > process. As far as I know the optimal weight is calcualted with a non-linear optimization procedure (heuristcs search with objective function to be minimzed the MSE or somthing similar). Based on this the statistical concepts are not applicable in this case (e.g. if we repeat the calculations in another sample then in the 95% of the cases etc.) So what is the meaning of the p-value and the standard error calcualtion? Do we have to check whether the smoothing coefficient is significant or not or we just disregard the standard error and the p-values? Is there any use of the standard error, the t-value and the p-value?


title: Should Kevin Durant be benched?

text: Kevin Durant is a world-class basketball player. He averages roughly 2 fouls for every 40 minutes of game play. Per NBA rules, you must stop playing after picking up 6 fouls. Supposing Kevin Durant has 4 fouls with 20 minutes remaining in a game: a) if he tries to play all 20 remaining minutes, what are the odds he acquires two more fouls and has to leave the game? b) if he were to sit on the bench for 5 minutes (leaving 15 remaining), then what would the odds be? c) if not 5, is there an ideal number for B in generating a greater likelihood of attaining more minutes than A?


title: Multiple linear regression, multi-period data and p-values

text: I would like to better understand the limitations of the MLR in regards to multiperiod data. Lets say I am looking at financial data for 10 companies, and each of those 10 companies has 5 periods worth of data, can I treat all 100 datapoints the same, or do I need to break it down into individual years? Additionally, I am noticing that removing a factor will cause another factors p-value to sporadically rise by a large number. My impression of the way a MLR works is that each factor should be independent of the each other, and while the coefficient is supposed to be divided by p number of factors, its relating stats like p-value and t-stat should remain constant if the data for that factor doesn't change Thanks for the clarity!


title: What are the assumptions of ordinal mixed effects logistic regression?

text: Specifically, what are the assumptions of ordinal mixed effect logistic regression performed with the "ordinal" package in R? I just got knobbled by a reviewer because these weren't clearly stated in a paper. The same reviewer asked me how I could "account for repeated measures". How would you explain the ordinal(ity) assumptions and the ME/ hierarchical model in a relatively parsimonious way?


title: Sparse Additive Generative Models (SAGE) v/s LDA (Latent Dirichlet Allocation)

text: Intuitively, how are Sparse Additive Generate Models (SAGE) by Eisenstein, different from Multinomial Dirichlet distributions. I understood that in this model distributions are added in logarithmic space. However, It is not clear why it can prevail over existing topic models. Does log space has any special characteristic which makes it more suitable, especially for sparse data.


title: Statistical description of how two distributions are different

text: I am have a collection of events recorded in two separate conditions, and I have measured the duration of each one of them. I would like to know two things: 1) are the distributions of durations different between the two conditions, and 2) if so, how are they different? For example, maybe the biggest difference between the two conditions is that there are far more 10 ms long events in condition #1 than in condition #2. I feel confident in using a Kolmogorov-Smirnov test to determine whether the distributions are significantly different, but am stuck on how to describe what the differences actually are using a statistical test or collection of test statistics. I could histogram the data and run some sort of ANOVA, but that seems a little course and inelegant. Or maybe it is fair to plot the distributions on top of each other and simply read off the difference at areas of interest. Any advice on how to tackle this problem? Below is some example data. Here a Kolmogorov-Smirnov test comparing these two distributions returns a very small p-value (p < 0.000001), so I know that the two distributions come from different populations. Based on how I _think_ the KS test works, this also tells me that where the two distributions are widest apart in the horizontal plane (around 0.9 on the y-axis) the two distributions are significantly different (there are more events with durations around 60 ms in group 1 than in group 2). But what about everywhere else in the distribution? Is there a way for me to know whether the difference in the amount of events that are 10 ms long is also statistically significant? ![enter image description here](http://i.stack.imgur.com/6XIxB.png)


title: How are these probabilities being calculated

text: This is an excerpt from the BBC website: > A recent study into speed dating habits concluded that if men and women go > to an evening and have 22 separate dates, men are keen to see about five > women again, while women would only choose to see two again, on average. > > That means that for every offer a woman makes, she has roughly a 50-50 > chance that the man will want to see her again too. > > But for every offer a man makes, he only has a one in five chance that the > desire to meet again is reciprocated. I have read the story a few times, but I don't see how the probabilities (of a successful date) for men and women are being calculated. Anyone care to shed some light on this?


title: How to replicate large well powered mixed effects model with a smaller sample?

text: Some edits made... I have a dataset which other researchers have used mixed effects modelling with to come up with a nice set of associations. I also have a much smaller dataset which is the same variables but from a different country. The first dataset is plenty powerful enough (350 individuals from 30 locations) but the second dataset is not very powerful at all (70 individuals from 15 locations). Is there some way to take what I know about the first data set and use it to boost the power of the second dataset? I've pooled all the data and analysed it and replicated the results, but I fear that is because the first dataset is just swamping the second one. So far I've come up with: 1) Analyse first dataset, then pool, analyse again, and compare AIC to see if it gets worse (seems a bit crude). 2) Analyse both datasets using a dummy variable for first/ second dataset (they're from different countries). This seems like a good idea but I'm not convinced it would eliminate the whole problem of the first swamping the second. 3) It sounds a bit like Bayesian statistics, but I don't know the first thing about Bayesian statistics, so if it is this then I'd better hit the library. Any and all help as ever gratefully received! Clarifications: 1) dependent variable same for both countries 2) oh dear, not sure I understand this question. All the variables are the same, so I have a list of their data, with a country code, and then the same variables again, with a different country code. All observations are individuals clustered in wards clustered in countries. Does this answer? 3) Really I'm just trying to replicate the significant linear associations which they found (or, I guess, highlight where they are different). So e.g. they found that size of ward improved outcomes, we would like to replicate that, along with some other associations.


title: SUM of AICc vs SUM of LL and then AICc

text: I have a technical question concerning calculating AICc for two possible models. For the data set I am working with there are 12 subjects and 10 phases of the experiment. Two different models, a Gamma distribution and an Inverse Gaussian distribution were fit to the data of each individual subject for each phase. Using a function I've found useful in exploratory unimodal distribution fitting (downloaded it off of Matlab's File Exchange) I got AICc values for each subject for each phase. Thus, within each phase of the experiment there are 12 AICc values for each model, one per subject. My question is as follows: Given the AICc values are already calculated for each individual subject, is it okay to sum the AICc values across subjects for each phase to yield to AICc values, one for the gamma distribution and one for the Inverse Gaussian and then calculate the difference or weights based off those two sums? The alternative is to get the Log-likelihood for each animal, sum across animals, and then calculate AICc for each model within each phase. However, it seems to me that this is equivalent to my first solution... Additionally, to determine the best model for the entire dataset (i.e. across phases) does the summing seem appropriate? Does anyone have any insight into this situation (getting the Log-likelihoods is doable, I'd need to edit the function, but the AICc are the default output of the function).


title: Assessing the representativeness of population sampling

text: I am looking for some suggestions about assessing the representativeness of a particular dataset I am analyzing. In this dataset I am looking at the relationship between two variables (e.g., X and Y) in a population that is split into five distinct blocks. The main problem is that the data is based upon reports from the public, so some blocks have much more data than others. The goal is to assess whether the relationship between X and Y differs between the blocks, but also to determine how reliable such estimates are given that we do not have a truly random sample of the overall population. Any suggestions appreciated. Thanks


title: Measuring correlation of trained neural networks

text: I am training an artificial neural network (backpropagation, feed-forward) with non-normal distributed data. Beside the root mean squared error, literature does often suggest the Pearson correlation coefficient for evaluating the quality of the trained net. But, is the Pearson correlation coefficient reasonable, if the training data is not normally distributed? Would it not be more reasonable, to use a rank-based correlation measure, e.g. Spearman rho?


title: Logistic regression with directional data as IV

text: I am looking for good references on using directional data (measure of direction in degrees) as an independent variable in regression; ideally, it would also be useful for hierarchical nonlinear models (the data are nested). I am also interested in directional data more generally. I have found a text by Mardia, which I am going to get, but wondered if there were good articles. I am more interested in practical articles about how to deal with this type of data than in theorems and proofs, or formal statements of distributions and such. Thanks UPDATE I have got the Mardia text, which is quite comprehensive. After some more reading, I may be back with more questions.


title: What to do following poor fit statistics for a confirmatory factor analysis?

text: ### Context I have got some problems with my doctoral dissertation. My thesis is Investigating Secondary Primary School Teachers' Organizational citizenship behaviours through their perceptions about organizational culture and their organizational trust levels. I have a sample of 871 teachers. I have three instruments but they were developed by other researchers and they were used in some other studies. I have been trying to analyse my data using structural equation modelling. However while trying to do confirmatory factor analysis, only one instrument (organizational citizenship) was ok. The other two intruments did not give the first order confirmatory analysis. The RMSEA values were around 0.100. Chi- square was too high, and chi-square divided by degrees of freedom was also too high. ### Question   * What should I do with my instruments and CFAs?    * Would item parcelling or removing items from these instruments be a good idea?    * Or should I skip structural equation modelling and go on with regression?


title: Overall type I error when repeatedly testing accumulating data

text: I have a question about **group sequential methods**. According to Wikipedia: > In a randomized trial with two treatment groups, classical group sequential > testing is used in the following manner: If n subjects in each group are > available, an interim analysis is conducted on the 2n subjects. The > statistical analysis is performed to compare the two groups, and if the > alternative hypothesis is accepted, the trial is terminated. Otherwise, the > trial continues for another 2n subjects, with n subjects per group. The > statistical analysis is performed again on the 4n subjects. If the > alternative is accepted, then the trial is terminated. Otherwise, it > continues with periodic evaluations until N sets of 2n subjects are > available. At this point, the last statistical test is conducted, and the > trial is discontinued But by repeatedly testing accumulating data in this fashion, the type I error level is inflated... If the samples were independent of one another, the overall type I error, $\alpha^{\star}$, would be > $\alpha^{\star} = 1 - (1 - \alpha)^k$ where $\alpha$ is the level of each test, and $k$ is the number of interim looks. But the samples are not independent since they overlap. Assuming interim analyses are performed at equal information increments, it can be found that (slide 6) ![enter image description here](http://i.stack.imgur.com/iFT4L.png) Can you explain me how this table is obtained?


title: Testing significance of overlap in R

text: _I have got an answer for it from Spacedman. But I am not entirely satisfied with the answer as it does not give me any sort of value (p or z value). So I am re-framing my question and posting it again. No offences to Mr.Spacedman._ I have a dictionary of say 61000 elements and out of this dictionary, I have two sets. Set A contains 23000 elements and Set B contains 15000 elements and an overlap of Set A and Set B gives 10000 elements. How can I estimate a p-value or z value to show that this overlap is significant and is not occuring by chance or vice versa. What I have been suggested till now includes MonteCarlo simulation methods. Is it possible to have an analytical method. Thank you in advance.


title: Help with Johansen procedure to check the cointegration

text: I'm using **ca.jo** function in URCA package to use Johansen procedure to check cointegration. I need help to understand the result of that test, this is my example:               > jo <- ca.jo(z, type='eigen', ecdet='none', spec='longrun')     > summary(jo)          ######################      # Johansen-Procedure #      ######################           Test type: maximal eigenvalue statistic (lambda max) , with linear trend           Eigenvalues (lambda):     [1] 0.335639191 0.001256000          Values of teststatistic and critical values of test:                     test 10pct  5pct  1pct     r <= 1 |   1.26  6.50  8.18 11.65     r = 0  | 408.52 12.91 14.90 19.19          Eigenvectors, normalised to first column:     (These are the cointegration relations)                    x.l2       y.l2     x.l2 1.0000000 1.00000000     y.l2 0.3328167 0.01557135          Weights W:     (This is the loading matrix)                     x.l2         y.l2     x.d -0.01751432 -0.002265876     y.d -2.90081864  0.007172311      where **z** is:               x <- diffinv(rnorm(1000))     y <- 2.0-3.0*x+rnorm(x,sd=5)     z <- ts(cbind(x,y))      Ok, now I have to check the critical values of the test. I see that we can reject r = 0 (ttest is above them), then for the second(r = 1) we CAN'T reject so the result should be that THERE IS one cointegrated vector but I don't understand Why R=1 what it does mean? I have a matrix with 2 column why I need to check r = 1 if we can reject r = 0 that mean there is one contegrated vector? Thank you in advance!


title: Can I account for an unknown starting period while estimating reliability?

text: I have been going through this article on how to estimate the reliability of a device using Weibull Analysis. The article clearly states that in order to do this, I need to have information in the following form:               #   Age_A   Age_B     1    26      16     2    31      24     3    35      30     4    38      36     5    41      42     6    44      48     7    47      55     8    51      63     9    55      76     10   70      126     Avg. 44      52      From my understanding, this age must have been measured from some **starting period** that is used as the experiment starting time. Now, if I am not really sure of the experiment starting point and have observations of failures across time, is there any way I can still carry out the Weibull Analysis to draw some conclusions? To put it clearly, let us say that in the middle of nowhere I decide to see which type is more reliable while both types have been operational for an unknown period of time. My idea was to fix a time period `t` and then start calculating the age relative to this `t`. Of course, the serious drawback is for devices that have been operational for a long time that fail right after `t` in which case the age is miscalculated at a small number. Any suggestions?


title: How do I interpret this Weibull plot?

text: I was exploring Weibull analysis for understanding reliability of two specific specimen. I used the R package, "weibulltoolkit" and "survival" to get the plot in question. The dataset is big so I am not posting it.               X = read.table("./test", header=T)     d <- data.frame(ob=X, state=1)     s = Surv(d$ob, d$state)     plot.wb(s)      From my understanding, both plots indicate bad fits! This could be an artifact of having too many data-points in the system or could be due to seasonal- effects or other unknown effects. Looking at the plot, what am I supposed to understand? Specifically,   * Plot 1 (using 4000 points) seems to indicate the presence of two different behaviors (0.05 - 5 and 5-1e+05). Does this mean that I am better off splitting the dataset into two and analyzing them separately? If so, what is the implication of this?   * Same seems to be going for Plot 2 (using 400 points) i.e. I would be better of splitting this dataset as well. What is a good way of comparing these two datasets? Any suggestions? **EDIT:** Upon Googling, I came across something called Multiple failure modes though I am not sure how to do this in R. Any suggestions would be greatly appreciated ### PLOT 1 ![enter image description here](http://i.stack.imgur.com/bEoiU.png) ### PLOT 2 ![enter image description here](http://i.stack.imgur.com/1umaL.png)


title: Uncertainty on parameter of unknown function

text: I have a collection of data points, which are counts as a function of length. Each data point is the result of many trials and has an error bar. The function `Counts(length)` is unknown. ![enter image description here](http://i.stack.imgur.com/wDP42.png) In the end the measurement that I need to quote is the length $X_m$ at which counts are a maximum (so in this case $X_m=0.9m$). Since the data are noisy, there is an uncertainty associated with $X_m$ (i.e., $X_m$ looks like it could be 0.6m-1m). If I knew the function that described these data then I could get an error on the fit parameters, but I have no information on what the function should be, so how to extract an uncertainty? I suppose I could create a hundred data sets like this, find the $X_m$ for each one and find the uncertainty that way, but it seems like that method ignores the information given by the error bars on all the points. Is there a formalism for figuring out the uncertainty on a parameter of an unknown function?


title: P-values for quantile regression with replicate weight via the 'survey' withReplicates function

text: I have been able to successfully apply BRR estimated weights to `rq()` using the `survey` and `quantreg` packages in R. The model comes out fine, I am able to obtain the $\theta$ and SEs, but I also want the $p$-values. Unlike with `rq()` where one is able get the $t$ and $p$-value etc, this approach does not provided the $p$-values but only the $\theta$ and SE. R script:               mod<-withReplicates(repwtx,quote(coef                        (rq(logA~ dose+ chtcm + cage +raceth + sex,                            tau=0.9, weights=.weights,alpha=0.05))))      and here's the output I get               >mod                      theta     SE     (Intercept)  1.3314065 0.1809     dose         0.3255046 0.0367     chtcm        0.0322069 0.0100     cage        -0.0516363 0.0154     raceth      -0.0082423 0.0196     sex          0.2314573 0.2887


title: Combining between and within-subjects designs?

text: I have a two by two design. Ideally I would run the whole thing within- subjects (all participants get all conditions) but the procedure is too long. I could run one factor between-subjects (each participants gets one condition) and one factor within-subjects(each participant gets both conditions), but then I have more power to see the effect for one factor vs. the other. I was thinking about running it such that one participant gets factor A within and another participant gets factor B within. That way I have to recruit fewer participants and I can see within subject and between subject effects for each factor. But I'm not sure how this would affect the analyses. Advice would be appreciated.


title: Comparing means of two simple time series

text: I'm trying to look for difference in timing (ie. earlier/later) in a variable measured at regular intervals between two groups. This seems like a simple experimental design, and working in R, I'm able to visualize the data in a way that makes sense to me, but somehow I'm getting confused when it comes to testing for significance. The data consist of weekly measurements of number of flowers for each individual, within and outside of the greenhouse. To take a small example:               expand.grid(week=(1:6),treatment=c("greenhouse","outside"),individual=1:2)->df     c(0,3,10,2,0,0,0,0,0,2,18,0,0,1,19,0,0,0,0,0,1,2,15,1)->flowers     data.frame(cbind(df,flowers))->df      Visually,               qplot(week,flowers,data=df,facets=treatment~.)      If my interest is simply to determine whether there's a significant difference in the time of flowering between the treatments; should I be doing a repeated measures ANOVA and looking at the interaction? Simplifying (?) the problem even further, what if I remove the quantity of flowers, and just consider how many individuals are flowering? So the summarized data would be               ddply(df, .(treatment,week), function(d) length(d[d$flowers>0,"flowers"]))->indiv      Which looks like this:                qplot(week,V1,data=indiv,facets=treatment~.)      Here, my first thought was that I can just think of these as two distributions, and compare with a t-test; however, only individuals and not individualsxweek are independent, so perhaps this should also be a repeated measures ANOVA? Or do I need to venture into the world of more complex time- series math? Thanks for your help! **Edit** As an update, I'm now also considering failure-time / survival analysis as a possible appropriate method.


title: Remove for systematic error (ANOVA)

text: Having a systematic 'error' (due to a 3 leveled external/experimental factor) influencing the variable I want to predict with my model. How do I best find a sound model that represents the underlying relation of interest, not fitting these 'errors'?


title: Likelihood Ratio of two-sample Uniform Distribution

text: Consider two uniform distributions: $$f \left( x, \theta_i \right) =\begin{cases} \frac{1}{2\theta_i} \quad -\theta_i<x<\theta_i, -\infty<\theta_i<\infty \\\ 0 \quad \text{elsewhere} \end{cases} $$ $i=1,2$ We want to test $H_0:\theta_1=\theta_2 $ against $H_1:\theta_1 \neq \theta_2 $ The order statistics of two independent samples from the respective distributions are $X_1< X_2<\ldots <X_{n_1}$ and $Y_1<Y_2< \ldots< Y_{n_2} $   1. Derive the Likelihood Ratio ,$\Lambda$   2. Find the distribution of $-2\log \Lambda$ when $H_0$ is true. * * * We assume the sample sizes are equal and put $n=n_1=n_2$ I can deal with part $1$. The mles for this distribution is $\hat{\theta_x} = \max \\{-X_1,X_n \\}$ and $\hat{\theta_y}=\max\\{-Y_1,Y_n \\}$ because we wish to make $\theta$ as small as possible given that the likelihood is a decreasing function of it. Thus, the Likelihood Ratio becomes: $$\frac{\hat{\theta_x}^{n} \hat{\theta_y}^{n} }{\left[ \max \\{ \hat{\theta_x}, \hat{\theta_y} \\} \right]^{2n} } $$ Examining the two different cases here, we can see that it equals: $$ \left[ \frac{\min \\{ \hat{\theta_x },\hat{\theta_y } \\}} {\max \\{ \hat{\theta_x}, \hat{\theta_y }\\}} \right]^{n}$$ Now, to finish part $2$ we have to derive the joint distribution of the numeraror and the denominator. According to my book: if we let $U=\min \\{ \hat{\theta_x}, \hat{\theta_y} \\}$ and $V= \max \\{\hat{\theta_x}, \hat{\theta_y} \\}$, the joint pdf is $$g(u,v)=2n^2 u^{n-1} v^{n-1} / \theta^{2n}\quad 0<u<v<\theta $$ where $n=n_1=n_2$ This is precisely the part I do not understand. Where does this distribution come from? My initial thoughts were that this is a joint pdf of maximum and minimum order statistics but this does not seem to be it. I would appreciate some help comprehending that part. I can finish the exercise afterwards. it's only that intermediate step that I find puzzling. The answer by the way is $\chi^{2} (2) $ Thank you in advance. **EDIT** I know that the joint pdf expression I have thrown at you is baffling but that's all I am given. If there is another way to show that the distribution of $-2 \log \Lambda$ is $\chi^{2} $, please let me know.


title: Is ARIMA(1,0,0)+xreg for level shift the same as linear regression model with level shift adjustment and lag1 term?

text: I have a time series with a level shift. Thus, when treating it with an ARIMA model, I use `arima(1,0,0)+xreg`. The `xreg` is a dummy variable for the level shift. And then I use linear regression: $$ y(t)=β_0 + β_1y(t-1) + β_3{\rm levelshift} + μ(t), $$ where ${\rm levelshift}$ is the same as `xreg` in ARIMA. I think these two models are exactly the same, however, `eviews` gives the different results even if I use least square methods for both estimations. Could anyone explain that? Is the different results due to `xreg` or `y(t-1)`?


title: Simple asymmetric ARCH vs asymmetric ARCH

text: I am working with GARCH models. As I try to familiarize myself with Stata's options I have run into a number of questions I am having problems answering. What is the difference between SAARCH and AARCH? If EGARCH and AARCH are both tools to deal with asymmetry, why prefer one vs the other?


title: Independence of residuals in a computer-based experiment/simulation?

text: I conducted a computer-based assessment of different methods of fitting a particular type of model used in the palaeo sciences. I had a large-ish training set and so I randomly (stratified random sampling) set aside a test set. I fitted $m$ different methods to the training set samples and using the $m$ resulting models I predicted the response for the test set samples and computed a RMSEP over the samples in the test set. This is a single **run**. I then repeated this process a large number of times, each time I chose a different training set by randomly sampling a new test set. Having done this I want to investigate if any of the $m$ methods has better or worse RMSEP performance. I also would like to do multiple comparisons of the pair-wise methods. My approach has been to fit a linear mixed effects (LME) model, with a single random effect for **Run**. I used `lmer()` from the **lme4** package to fit my model and functions from the **multcomp** package to perform the multiple comparisons. My model was essentially               lmer(RMSEP ~ method + (1 | Run), data = FOO)      where `method` is a factor indicating which method was used to generate the model predictions for the test set and `Run` is an indicator for each particular **Run** of my "experiment". My question is in regard to the residuals of the LME. Given the single random effect for **Run** I am assuming that the RMSEP values for that run are correlated to some degree but are uncorrelated between runs, on the basis of the induced correlation the random effect affords. Is this assumption of independence _between_ runs valid? If not is there a way to account for this in the LME model or should I be looking to employ another type of statical analysis to answer my question?


title: Perpendicular offsets in a weighted least squares regression

text: Perpendicular offset least square fitting has a lot of advantages compared to the native least square fitting scheme. The following figure illustrates the difference between there, and for a more detailed comparison of these two methods, we refer to here. ![enter image description here](http://i.stack.imgur.com/Aue8i.gif) Perpendicular offset least square fitting, however, is not robust to outliers( points that are not supposed to be used for model estimation). Therefore, I am now considering to use a weighted perpendicular offset least square regression method. The method has two steps:   1. Calculate the weighting factor for each points that are going to be used for line estimation;   2. Perform perpendicular offset in a weighted least square regression scheme.  For the time being, my biggest problem comes from step 2. Suppose the weighting factors are given, how can I get the formula to estimate the parameters of the line? Many thanks! **EDIT:** Based on the kind suggestion of @MvG I have implemented the algorithm in MATLAB:               function  line =  estimate_line_ver_weighted(pt_x, pt_y,w);     % pt_x  x coordinate     % pt_y  y coordinate     % w     weighting factor               pt_x = pt_x(:);     pt_y = pt_y(:);     w    = w(:);               % step 1: calculate n     n = sum(w(:));          % step 2: calculate weighted coordinates      y_square = pt_y(:).*pt_y(:);     x_square = pt_x(:).*pt_x(:);     x_square_weighted = x_square.*w;       y_square_weighted = y_square.*w;       x_weighted        = pt_x.*w;     y_weighted        = pt_y.*w;          % step 3: calculate the formula     B_upleft = sum(y_square_weighted)-sum(y_weighted).^2/n;     B_upright = sum(x_square_weighted)-sum(x_weighted).^2/n;     B_down = sum(x_weighted(:))*sum(y_weighted(:))/n-sum(x_weighted.*pt_y);     B = 0.5*(B_upleft-B_upright)/B_down;          % step 4: calculate b     if B<0         b       = -B+sqrt(B.^2+1);     else         b       = -B-sqrt(B.^2+1);     end          % Step 5: calculate a     a = (sum(y_weighted)-b*sum(x_weighted))/n;          % Step 6: the model is y = a + bx, and now we transform the model to      % a*x + b*y + c = 0;     c_ = a;     a_ = b;     b_ = -1;          line = [a_ b_ c_];      The result is as good as we can expect, which is illustrated in the following script:               %% Procedure 1: given the data     pt_x = [   692   692   693   692   693   693   750];     pt_y = [ 919         971        1022        1074        1126        1230        1289];          % Procedure 2: draw the point       close all; figure; plot(pt_x,pt_y,'b*');          % Procedure 3: estimate the line based on the weighted vertical offset     % least square method.      weighting = ones(length(pt_x),1);      weighting(end) = 0.01;  % we give the last point a low weighting because obvously it is an outlier      myline =    estimate_line_ver_weighted(pt_x,pt_y,weighting);       a = myline(1); b = myline(2); c= myline(3);           % Procedure 4: draw the line      x_range = [min(pt_x):0.1:max(pt_x)];      y_range = [min(pt_y):0.1:max(pt_y)];      if length(x_range)>length(y_range)             x_range_corrspond = -(a*x_range+c)/b;             hold on; plot(x_range,x_range_corrspond,'r');      else             y_range_correspond = -(b*y_range+c)/a;             hold on; plot(y_range_correspond,y_range,'r');      end      The following figure corresponds to the above script: ![enter image description here](http://i.stack.imgur.com/IIz2k.png).


title: ANOVA with uncertain values - attempt to redefine node

text: I'm trying to do one-way ANOVA on uncertain (latent) input value `YReal[i]` (each input value is not exactly known, but it is given as normal distribution with mean `Y[i]` and standard deviation `SE[i]`). But, I got stuck. I tried the following model:               model {          # likelihood     for (i in 1:N) {         YReal[i] ~ dnorm(alpha[group[i]], eps_tau) # one-way ANOVA expression         YReal[i] ~ dnorm(Y[i], 1/(SE[i]*SE[i])) # uncertainty of the input value     }          # priors          eps_tau ~ dgamma(0.01, 0.01)     eps_sigma <- 1 / sqrt(eps_tau)          for (j in 1:no_group) {         alpha[j] ~ dnorm(0, 0.0001)     }          }      but I got message from jags > Attempt to redefine node YReal[1] I tried to rewrite the model as               YReal[i] <- alpha[group[i]] + eps[i]     eps[i] ~ dnorm(0, eps_tau)     YReal[i] ~ dnorm(Y[i], 1/(SE[i]*SE[i]))      or               eps[i] <- YReal[i] - alpha[group[i]]     eps[i] ~ dnorm(0, eps_tau)     YReal[i] ~ dnorm(Y[i], 1/(SE[i]*SE[i]))      But I got the same error. Is there any bayesian tool or trick how to get around that and achieve the goal? Thanks. P.S.: I tried this:               YReal[i] ~ dnorm(alpha[group[i]], eps_tau)     Y[i] ~ dnorm(YReal[i], 1/(SE[i]*SE[i]))      but I'm not even remotely sure it is equivalent to the original model (and that it does what I need)


title: Bonferroni correction: control vs. groups?

text: I'm trying to understand how to set up a Bonferroni correction on several different groups and compare it to the control group. The groups and observations are as follows (with Group 0 being the control): Group $\hspace{5mm}$ Observations   0$\hspace{15mm}$ 97.76$\hspace{5mm}$ 102.56$\hspace{5mm}$ 96.08$\hspace{5mm}$ 125.12   1$\hspace{15mm}$ 91.28$\hspace{5mm}$ 129.20$\hspace{5mm}$ 90.80$\hspace{5mm}$ 72.32   2$\hspace{15mm}$ 67.28$\hspace{5mm}$ 85.76$\hspace{5mm}$ 95.60$\hspace{5mm}$ 73.28   3$\hspace{15mm}$ 80.24 $\hspace{5mm}$ 64.88 $\hspace{5mm}$64.88 $\hspace{5mm}$ 78.56   4$\hspace{15mm}$ 96.08$\hspace{5mm}$ 98.24 $\hspace{5mm}$77.84 $\hspace{5mm}$ 95.36   5$\hspace{15mm}$ 57.68 $\hspace{5mm}$ 89.84$\hspace{5mm}$ 98.48 $\hspace{5mm}$ 92.72   6$\hspace{15mm}$ 68.72 $\hspace{5mm}$85.28 $\hspace{5mm}$68.72 $\hspace{5mm}$ 74.24 I know typically the Bonferroni correction works by multiplying probabilities by the total number of probabilities and then comparing the new value to the critical value. However, I am not familiar with how to perform a Bonferroni with several groups (and subsequent observations within each group).


title: Cluster detection in a distribution

text: > **Possible Duplicate:**   >  Determine different clusters of 1d data from database I am trying to detect clusters, i.e. segments on a linear scale, with visibly higher concentration of distribution elements. Consider the following simplified distribution               VAR01=3     VAR02=17     VAR03=18     VAR04=21     VAR05=23     VAR06=44     VAR07=61     VAR08=89     VAR09=90     VAR10=92     VAR11=94     VAR12=117      If you put these values on a diagram, you could see that there is a visible concentration (cluster) between 17-23 and 89-94 with VAR01, VAR06, VAR07 and VAR12 being lone outliers. What are the best techniques in programmatically recognizing such clusters in a distrubution?


title: What to do if the number of parameters and observations are almost equal?

text: I have a data set of 10 observations on 9 bivariate random variables (one is response & others are covariates). I want to fit a linear regression on the whole data set. I want to know if there is any special technique to do so. I can divide the data into two equal groups. Please suggest something for that situation also. Thanks in advance for any kind of help/suggestion.


title: Not usual piecewise linear regression

text: I know that with fixed breakpoints (for the x variable), piecewise linear regression is easy and linear. However, for a variety of reasons, I would like to fix the breakpoints for the y variables. One reason is that there are natural candidates for the breaks for the y variable, but I have no clue about the x. Because there is a lot of noise, I do not want to just fit the reversed model.   I can write down the expression to be minimized, but it is nonlinear least squares, and the gradients look messy, involving the d/dr(1/r), etc. I want to avoid this. Does anyone know of a good solution to this problem? Maybe I have a missed a way of rewriting it to be linear?


title: Format for distance matrix for running hclust in R

text: I am using R for a simple hierarchical clustering method for finding protein sequence similarities. I already have a distance matrix computed. I am using the `hclust` method to run it (thus I dont want to use `dist` method). Could anyone tell me what format does my distance matrix need to be for `hclust` to run? I used the `read.xls` method to read my Excel file. I need only the lower triangular portion of the matrix, and this is what I am sending to `hclust` to be clustered.               d= as.matrix(mat)     hc =hclust(d, method="single")      However this is the error I get:               Error in if (n < 2) stop("must have n >= 2 objects to cluster") :      argument is of length zero


title: Representing the anova's interaction in R

text: I read in the book The basics of S and S-plus that "A:B" functions as "Interaction between A and B".   However, in other book I read "There is never an interaction effect for two factors when one is nested within the other".   But aov(y~a/b) would have the same results as aov(y~a+a:b). So it seems to me a:b could mean interaction, but also could mean the effect of a factor nested within another. In other words, A:B" definition as "Interaction between A and B" is incomplete.   Am I right? How to explain all this?


title: What is the distribution of this nearly-Hotelling statistic?

text: Let $X$ be an $n \times l$ matrix, and $F$ an $n \times p$ matrix, with the rows of $X$ and $F$ drawn _i.i.d._ from multivariate Gaussians. (The independence applies to rows: the $X$ and $F$ may be correlated.) Let $\Sigma_X, \Sigma_F, \Sigma_{XF}$ be the covariances of the two Gaussians and the covariance of the two Gaussians with each other. (It may be simpler to think of a single multivariate Gaussian which I have partitioned into two parts.) Given $p$-variate vector $v$, I would like to compute the distribution of $$v^{\top}(F^{\top}F)^{-1}F^{\top}X\left(X^{\top}X - X^{\top}F(F^{\top}F)^{-1}F^{\top}X\right)^{-1}X^{\top}F(F^{\top}F)^{-1}v$$ I am usually looking at the case where one of the columns of $F$ is all ones (which we can model as a Gaussian with a very small variance, if need be). In the case where $p = 1$ and $F$ is all ones, then we recover Hotelling's $T^2$ (up to scaling of $v$, which I ignore here). In fact I arrived at this monster by asking the question: given $l$-variate vector $w$, letting $y=Xw$, then regressing $y$ against $F$ to get $\hat{\beta}$, and computing the t-statistics associated with $v$ ( _i.e._ $(\hat{\beta}^{\top}v)/(\hat{\sigma}\sqrt{v^{\top}(F^{\top}F)^{-1}v})$), what is the maximal value of the t-statistic over all $w$?


title: Testing coefficient equality across two Tobit models

text: I have a tobit model the looks like `Y = cons + B*(X) + error`. This model is estimated twice on two different samples. I want to be able to test if, say, B is different across the two models (the sample A model vs the sample B model). I have done this before for logit models using two approaches that yielded the same results. My issue is that these two approaches are not matching when using the tobit model. Approach 1: I stack the two samples. If SA is a dummy that is 1 when sample A and SB is a dummy that is one when Sample B then I estimate: `Y = SA + SB + B1*(SA*X) + B2*(SB*X) + error` then I test `B1==B2` and done. Approach 2: Here I use something called seemingly unrelated models. I am using stata, and the procedure involves the following: estimate `Y = cons + B*X + error` for sample A estimates store ... estimate `Y = cons + B*X + error` for sample B estimates store ... combine the estimates using:               suest sample_A_model sample_B_model      then test the coefficients. Approaches 1 and 2 show different results. Note that if I set up the regressions so that no values are censored, the results are identical. The problem arises when there are censored values. Any thoughts? Thanks in advance!


title: Elliptic regression, basic conceptual question

text: I'm considering circular regression and elliptic regression on a computational and conceptual basis. If we fit an ellipse to our data then we deal with the principal components as reference for the main and for the secondary axis; we set the lengthes of the axes according to the sqrt of the eigenvalues. The idea behind this is, that in the first principal component we have the maximal variance and the second the minimal variance. Accordingly we construct the confidence-ellipse this way. Now I was thinking: assume we have a huge number of normal distributed data, slightly correlated such that we have a nice elliptic shape. If we rotate to the principal components we get the two values for the length of the axes by the variance in each direction. But if we would think of, say, 40 angular segments of the ellipse, consecutively ordered around the origin, then shouldn't the shape of the confidence-curve follow the variance in each direction/each of the 40 segments? What would such a curve in the limit look like? Suddenly, as I think about this I' getting doubt it would be an ellipse...


title: Should I use Pearson Correlation or Linear Regression to show the linearity of two sets of numbers？

text: I have two sets of acceleration data in terms of time. In one set the acceleration increases more linearily than another set. How can I show the linearity difference of the two set of data best? Should I use Pearson Correlation or Linear Regression? If I get r = 0.8 and 0.98 by Pearson Correlation, how I can interpretate the result? Do they have very big difference on the linearity based on 0.8 and 0.98? Thanks.


title: Is it discriminant analysis?

text: I have samplings of one-dimensional data of two classes: $A$ and $B$. I have to predict the posterior probability of class $A$. $$\tilde{P}(A|x) \approx \frac{N_A\tilde{f_A}(x)}{N_A\tilde{f_A}(x) + N_B\tilde{f_B}(x)}.$$ $f_A$, $f_B$ are the estimation of densities of classes $A$ and $B$. I estimate the densities as normal distributions, whose $\mu$ and $\sigma$ are average values and standard deviations of all the points of each class ($A$ and $B$ accordingly). Will it be correct to say, that I'm using `discriminant analysis` in this case?


title: Estimating variance of center-censored Normal samples

text: I have normally-distributed processes from which I get small samples ( _n_ typically 10-30) that I want to use to estimate variance. But frequently the samples are so close together that we can't measure individual points near the center. I have this vague understanding that we should be able to construct an efficient estimator using ordered samples: E.g., if I know the sample contains 20 points, and that 10 are clustered near the center too tightly to measure individually, but I have discrete measurements of 5 on either tail, is there a standard/formulaic approach for estimating the process variance that makes optimal use of such samples? (Note that I don't think I can just weight the center average. For example, it's possible for 7 samples to cluster tightly while another three are asymmetrically skewed to one side but close enough we can't tell that without more tedious single sampling.) If the answer is complicated any tips on what I should be researching would be appreciated. E.g., is this an order-statistic problem? Is there likely to be a formulaic answer, or is this a computational problem? **Updated detail:** The application is analysis of shooting targets. A single underlying sample is the point of impact ( _x,y_ ) of a single shot on the target. The underlying process has a symmetric bivariate normal distribution but there is no correlation between axes, so we are able to treat the { _x_ } and { _y_ } samples as independent draws from the same normal distribution. (We could also say the underlying process is Rayleigh-distributed, but we can't measure the sample Rayleigh variates because we can't be certain of the coordinates of the "true" center of the process, which for small _n_ can be significantly distant from the sample center ($\bar{x}$, $\bar{y}$).) We are given a target and the number of shots fired into it. The problem is that for _n_ >>3 precise guns will typically shoot a "ragged hole" surrounded by distinct shots. We can observe the _x_ \- and _y_ -width of the hole, but we don't know _where_ in the hole the non-distinct shots impacted. Here are some examples of more problematic targets: ![\[Sample target with n=10\]](http://i.stack.imgur.com/RTk1K.jpg) Sample target with n=100 (Granted, in an ideal world we would change/switch targets after each shot and then aggregate the samples for analysis. There are a number of reasons that's often impractical, though it is done when possible.) **Further Notes** following WHuber's clarifications in the comments: Shots produce target holes that are of uniform and known diameter. When a shot is outside any "ragged group" we know the projectile radius and so we can measure the precise center $x_i$. In each "ragged group" we can discern some number of peripheral "balls" and again mark the precise center of those outer shots based on the known projectile radius. It is the _remaining_ "center-censored" shots that we only know impacted _somewhere_ in the interior of a "ragged group" (which is typically -- and if necessary let us assume -- one per target). To facilitate solution I believe it will be easiest to reduce this to a set of one-dimensional samples from the normal, with a central interval of width _w_ > _d_ , where _d_ is the projectile diameter, containing _c_ < _n_ "censored" samples.


title: How to report the standard deviation of a sample of averages?

text: If you have a 100x100 grid of values and you want to report the standard deviation of an area that takes up a 20x20 portion of the grid, this is not problem. But what is the correct way to calculate the standard deviation of the 20x20 area if all the values in the 100x100 grid are themselves averages and therefore also have an associated (common) standard deviation? I haven't dealt with such a spatial statistics problem for a while and I think I have forgotten some terminology...


title: Difference in meaning of these terms: Dataset vs Corpus

text: I have a question pertaining to definition of terms that I couldn't find answer for. What's the difference between Dataset and Corpus? I've seen them being used almost interchangeably. My understanding is that Corpus (meaning collection) is broader and Dataset is more specific (in terms of size, features, etc). Please let me know what you think. I apologize in advance if this isn't the right forum for this question. The reason I have a doubt is because my adviser and I are having a confusion on what is most suitable for our paper titled: "Techniques to extract XX from multimedia corpus/dataset"


title: Warranty claims count prediction

text: I have weekly data of number of warranty claims from week of manufacturing. I want to use this historical data and predict the number of claims that may arise in next few weeks. I have read about count data models like Poisson regression and NB regression but they expect the counts in the range of 0~30 (is this correct?). In my case the claims can be of the range 0 to several thousands. What would be an ideal way to solve this problem?


title: Fit a regression line by using `MATLAB`

text: I have the following data                 Individual     Heart rate     Weight    Hours of exercise per week            1           72            134             3.2            2           81            201             3.5            3           60            156             7.1            4           82            148             2.4            5           75            170             1.2      Now i have to fit a regression line by using `MATLAB`. But which one is independent variable ? Hours of exercise per week? And is that Heart rate & Weight response variables ? Then how can i fit `one` regression line ? If i have two response or two independent variables, won't i get `two` regression line ?


title: Positively constrained parameters in LAD LASSO

text: I am looking for an implementation of LAD LASSO method with constraints on coefficients, such as non-negative beta coefficients for example. $\min_{\beta} \left| \mathbb{y} - \mathbb{X} \boldsymbol{\beta}\right|_1 + \lambda \left|\boldsymbol{\beta}\right|_1 $ Are there any packages available with such properties? Also interested are there any methods which combines elastic net types of penatly with LAD? Thanks


title: kernel for a (semi-) metric space

text: Let's say I have a metric space $(\mathcal{X}, d)$. Is there any kernel function that I can use with SVM? If we change the RBF kernel a little bit, we have $k(x,y) = e^{-d(x,y)^2}$. Is this a valid kernel? Additionally, what if $d$ is semi-metric?


title: Is there a statistical application that requires strong consistency?

text: I was wondering if someone knows or if there exists an application in statistics in which strong consistency of an estimator is required instead of weak consistency. That is, strong consistency is essential for the application and the application would not work with weak consistency.


title: How to calculate absolute fit indices (RMSEA, GFI...) from relative ones (AIC, BIC...)?

text: I have conducted an IRT analysis with Conquest in order to compare two models (1-dimensional vs. 8-dimensional) applied to a given data set (41 items of a questionnaire, N=195). Comparing the deviance scores (as well as resulting AIC/BIC scores), I know that the 8-dimensional model is superior in a relative sense. How do I know if this model describes the data well enough in an absolute sense? My idea would be to calculate fit indices like GFI, NNFI etc. but I don't have a clue how to get them from my Conquest output.


title: Is there any problem with a Bayesian setting that gives positive prior probability on malformed events?

text: I have a model class $M$. Each $\theta \in M$ defines a probability distribution $p(x,y | \theta)$ such that $x \in A$ and $y \in B$. I am interested only in probability distributions of $\theta$ such that their support is contained in $A_0 \subseteq A$ and $B_0 \subseteq B$ -- other distributions that are not of this form are "malformed" (they define a distribution over too permissive space). Now, I have a prior $\pi(\theta)$. This prior can give positive probability to "malformed" probability distributions. But now I assume that I observe $x_1,\ldots,x_n$, where $x_i \in A_0$. I also assume that $p(y | x,\theta)$ for _any_ $\theta$ has a support in $B_0$, if $x \in A_0$. **If I do Bayesian inference with this prior $\pi$ and $x_1,\ldots,x_n$ to compute the posterior $p(\theta, y_1,\ldots,y_n | x_1,\ldots,x_n)$ -- is there any practical or theoretical problem here?** Supposedly, my prior can give positive probability to "bad" $\theta$, but it seems to me that since $x_1,\ldots,x_n$ are always from $A_0$ in practice, there aren't too many problems with the posterior either, because $y_1,\ldots,y_n$ are always going to be from $B_0$.


title: how to read y axis in kernel density graph

text: I need to understand how to read kernel density graphs. How do you come up with the values in y-axis?![enter image description here](http://i.stack.imgur.com/8CRLT.jpg)


title: sampling question

text: I have a question about inclusion criteria to obtain a representative sample of people with "sustained" or uncontrolled high BP. A colleague proposed that I look at health history data retrospectively (i.e. three years back), and include only members who have a BP of 140 over 90 consistently for all three years. I had asked what about those who only had a measurement of 140 over 90 in two or one years? And if they weren't consecutive? Although the rationale sounds logical, I would like a more research-based method. Thank you in advance.


title: Distribution function of an indicated random variable

text: _A random variable $X$ has a distribution function $F$. Let $Y\triangleq XI_{(a,b)}(X)$ with $1 <a<b$. Find $G$, the c. d. f. of $Y$._ $$\begin{align} G(y) &\triangleq \mathbb P(Y\le y) =\\\ &= \mathbb P(X\le y)\,I_{(a,b)}(X) + \mathbb P(0\le y)\,I_{(a,b)^\complement}(X) =\\\ &= I_{(a,b)}(X)\,F(y) + I_{(a,b)^\complement}(X)\,I_{[0,\infty)}(y)\quad. \end{align}$$ I'm stuck in the indicators of $X$. It occurred to me to write $I_{(a,b)}(X)$ as $\mathbb P(a<X<b)$, but that's just $E(I_{(a,b)}(X))$. The book gives the answer $$\begin{align} G(y) = [1 - &F(b) + \mathbb P(X=b) - F(a)]\,I_{[0,a)}(y) +\\\ &+\; [1 + F(y) - F(b) + \mathbb P(X=b) - 2F(a)]\,I_{[a,b)}(y) + I_{[b,\infty)}(y)\quad. \end{align}$$ How do I arrive at it?


title: Hints on this computer vision / machine learning problem

text: I've been working for a while on a pet problem. The task is to identify and segment out the dark lines and possibly the wiggly ones too. I'm not looking for anyone to solve this problem for me...I'm just at the point where I feel I'm losing my footing and could benefit from some ideas! I've so far focussed on neural networks with varying degrees of not-super- success. The nns haven't given much more success than some threshholding. I have tried edge detection but the results look messy. A shrink, laplace, expand showed some promise. Hough transforms...it seems they should work but not yet(for me, anyway). thanks.![example of source image from which to extract the dark lines](http://i.stack.imgur.com/xJdF9.png)


title: hypothesis testing in hierarchical trees

text: Suppose I have an hierarchical tree organising people according to their music preferences. Each person then has 3 extra properties A, B and C (e.g. parent's supportiveness on musical education). I wanted to select people from A and test across several tree distances if their friends tend to be more A or C then the friends from group B. The goal of the toy example is to find distant influences in friend's space (tree) between music preferences and parent's supportiveness on musical education not just of person in question but also of the friends. I was wondering if a Pearson's chi-square test would be appropriate in this scenario. As I see it, there are two cohorts: Cohort1: people from group A. Cohort2: people from group B. Then, for a given tree distance (dist=2,3,4 ... ), test wether cohort1 friends are A or not A more than cohort2. Possible problems: The different cohorts are composed of different people (A!=B) but at a given tree distance cohort1 and cohort2 will have friend's in common (e.g., at distance 4 across tree branches, a friend from group A of a person34 from group A, cohort1, can also be a friend of person38 from group B- cohort2). Does this mean independency between rows of my contingency table does not hold? Thanks!


title: 6-sided die roll average value

text: > A die is rolled 100 times. Find the probability that the mean value would be > between 3 and 3.5. I've found the expected value which is 3.5 without idea what to do next.


title: Predicted probabilities from a multinomial regression model using zelig and R

text: I have a multinomial model estimated with the zelig package in R. Whenever I try to use the setx() command, I get an error message saying there is more than one mode. So in stead of using Zelig, I thought I would do it the hard way. I used the instructions here, but I am not sure if I trust these instructions as they give me some weird results.   * Could anyone explain how to get predicted probabilities from a multinomial logit in R? * * * PS: here are the coefficients from a simple model (there are 6 choice catagories):                structure(c(3.68021133487111, -0.903496528862169, -1.56339830041814,      -1.13238307296064, -1.67706243532044, -0.177585202845615, 0.0611115470557421,      -0.0458373863009504, 0.0881133593132653, -0.0686190052488972,      0.0163917121907627, 0.0165232098847022, 0.0373815294869855, -0.0353209839724262,      -0.00698911507852077), .Names = c("(Intercept):1", "(Intercept):2",      "(Intercept):3", "(Intercept):4", "(Intercept):5", "Deviance:1",      "Deviance:2", "Deviance:3", "Deviance:4", "Deviance:5", "Votes:1",      "Votes:2", "Votes:3", "Votes:4", "Votes:5"))      and here are some values for the independent variables:               structure(c(0.71847390030784, 1.01838748408701, 1.01838748408701,      1.20499277373001, 0.71847390030784, 1.20499277373001, 0.56393315893118,      1.20499277373001, 0.71847390030784, 0.56393315893118, 2, 5, 4,      7, 10, 27, 5, 9, 17, 16), .Dim = c(10L, 2L))


title: logistic regression with slack

text: I have a Bayesian model that tries to predict a binary variable that I am modelling as logistic regression. The training data have lot of wrong labels*. Therefore, I think I might need to introduce some sort of slack (similar to slack variable that SVM). How can I do this? If I model logistic regression ignoring errors in training data, is logistic regression robust? *The logistic regression training-data is sampled from another Bayesian model therefore there is large variance in my Gibbs sampling.


title: Is there a graphical representation of bias-variance tradeoff in linear regression?

text: I am suffering from a blackout. I was presented the following picture to showcase the bias-variance tradeoff in the context of linear regression: ![Polynomial model for data, simple and complex case](http://i.stack.imgur.com/8RlJk.png) I can see that none of the two models is a good fit - the "simple" is not appreciating the complexity of the X-Y relation and the "complex" is just overfitting, basically learning the training data by heart. However I completely fail to see the bias and the variance in these two pictures. Could someone show this to me? PS: The answer to Intuitive explanation of the bias-variance tradeoff? did not really help me, I would be glad if someone could provide a different approach based on the above picture.


title: test for comparing means of two time series

text: I wonder if there is a test for comparing means of two different time series (both measured at the same person or two different persons)? Would be very glad if someone can help me!


title: SAS Code for Survival Analysis

text: In PROC PHREG, how do you set a continuous variable at a certain value as the reference level? For example, suppose **x = 3.5, 3.6, 4.3, 5.4** and we want **x = 6** to be the reference level. How would you do this in SAS? Would this be correct:                     data new;              set old;             x1 = (x=3.5);             x2 = (x= 3.6);             x3 = (x = 4.3);             x4 = (x= 5.4);                 proc phreg data = new;               class x ( ref = '6');              model time*censor(0) = x1 x2 x3 x4 /r1;             run;       **Edit.** Actually, could one just use the point estimate of the coefficient of $x$? In other words suppose the estimated coefficient for $x$ is $0.512$. Then the hazard ratio between $x = 3.5$ and $x = 6$ would be $\exp(0.512(3.5-6))$, the hazards ratio between $x = 3.6$ and $x = 6$ would be $\exp(0.512(3.6-6))$ etc...? So there is no need to even form these groups? Just do the following:                   proc phreg data = new;            model time*censor(0) = x /r1;         run;


title: Reporting non-causal relationships

text: What is the appropriate way of reporting significant regression coefficients of a multiple regression when all variables have been obtained at the same measurement occasion? Specifically, do I imply a causal relationship when I report that A significantly predicts B or A has a significant effect on B?


title: We're trying to analyze which monsters drop what in an MMO

text: After killing a monster there's a percent chance certain items drop that are unique to the monster, but sometimes the monsters are in groups and we only get the drops after killing all of them. Monster A, B, C sometimes appear like [A,B] A with B but never [A] alone, but sometimes as [C,A]. What is it called solving this type of problem, for example isolating what monster A drops with a confidence interval? I think it's a linear regression or a multivariate linear regression but I'm not sure what mathematical tools I would use to solve this. Just stats or probability alone? I've toyed with R but my stats is not that good anymore. I'd like to know if this is within the realm of a problem a beginner can be able to solve with maybe a few weeks of study.


title: Partialling or regressing out a categorical variable?

text: Occasionally I see in literature that a categorical variable such as sex is “partialled” or “regressed” out in (fixed-effects or mixed-effects) regression analysis. I'm troubled with the following practical issues involved in such a statement: (1) Usually the coding method is not mentioned in the paper. Such a variable has to be coded with quantitative values, and I feel the sensible way should be effect coding (e.g., male = 1, female = -1) so that partialling can be achieved with other effects interpreted at the grand mean of both sex groups. A different coding may render a different (and unwanted) interpretation. For example, dummy coding (e.g., male = 0, female = 1) would leave other effects associated with males, not the grand mean. Even centering this dummy-coded variable might not work well for their partialling purpose if there is unequal number of subjects across the two groups. Am I correct? (2) If the effect of such a categorical variable is included in the model, examining its effects first seems necessary and should be discussed in the context because of its consequence on the interpretation of other effects. What troubles me is that sometimes the authors don't even mention the significance of sex effect, let alone any model building process. If the sex effect exists, a natural follow-up question is whether any interactions exist between sex and other variables in the model? If no sex effect and no interactions exist, sex should be removed from the model. (3) If sex is considered of no interest to those authors, what is the point of including it in the model in the first place without checking its effects? Does the inclusion of such a categorical variable (and costing one degree of freedom on the fixed effect of sex) gain anything for their partialling purpose when sex effect exists (my limited experience says essentially no)?


title: Comparing ROC from Logistic Regression with Miss-Classification rate from Neural Network

text: Can I compare the ROC output from a Logistic Regression model with the Miss- Classification rate output from the Neural Network for the same dataset. Are these apples to apples comparison ? ~..


title: RTM and study design issue

text: I was asked to estimate "regression to the mean" (RTM) error for the purpose of measuring program efficacy for reducing spending among the highest spending individuals. Given the constraints of the data (administrative data set), I have identified the following issues:   1. selection bias: patients choice determines whether a program is given; and,   2. unable to separate program from RTM effect: lack of a "true" comparison group Given these constraints, I foresee many challenges to this project. I have, however, read this article, and I want to know if I could conduct something similar to meet the requirements of this ask. Specifically, I am wondering if the following procedure will work:   1. Estimate the RTM for all data, and adjust the outcome variable (spending in $)   2. Identify the "highest" spenders, which was proposed as top 10% spenders   3. Assign "treatment to those who were in the program and assign "control" to those who were not in the program (control)   4. Assess differences while adjusting for "similar" characteristics to try and control for some of the selection bias, while accepting some residual bias (MANOVA?)  The most important step right now would be to estimate a value for RTM; however, I would like to also know what are some options I could take in order to solve the longer-term problem. Thanks in advance.


title: How to use a point estimator on an interval?

text: I got part a and part b just fine, but I am confused on how to do part (c), (d), and (e). I don't really understand what the question is asking, nor do I understand how to do that sort of interval estimation. Any help is greatly needed and appreciated! Thanks! Consider the following sample of observations on coating thickness for low- viscosity paint (“Achieving a Target Value for a Manufacturing Process: A Case Study,” J. of Quality Technology, 1992: 22–26): 0.83 0.88 0.88 1.04 1.09 1.12 1.29 1.31 1.48 1.49 1.59 1.62 1.65 1.71 1.76 1.83 Assume that the distribution of coating thickness is normal (a normal probability plot strongly supports this assumption). a.) Calculate a point estimate of the mean value of coating thickness, and state which estimator you used. Answer = 1.3481234 using x-bar b.) Calculate a point estimate of the median of the coating thickness distribution, and state which estimator you used. Answer = 1.395 using ~x c.) Calculate a point estimate of the value that separates the largest 10% of all values in the thickness distribution from the remaining 90%, and state which estimator you used. [Hint: Express what you are trying to estimate in terms of μ and σ.] d.) Estimate P(X < 1.5), i.e., the proportion of all thickness values less than 1.5. [Hint: If you knew the values of μ and σ, you could calculate this probability. These values are not available, but they can be estimated.] e.) What is the estimated standard error of the estimator that you used in part (b)?